Magnitude and phase of the signal using cwt

15 visualizzazioni (ultimi 30 giorni)
Giovanni Pirro
Giovanni Pirro il 27 Mag 2024
Risposto: Vinay il 2 Set 2024
Hi,
I have a signal in time that I want to analize in the frequency domain using wavelet transform (to balance time and frequency resolution). When I check the magnitude scalogram that is automatically generated using the function cwt, the value of the magnitude that i read is somehow corresponding to the value I am expecting from the time signal. When I try to extract this values computing abs and angle of the coefficients I get from cwt and compare them with expected magnitude and phase, I see there is some scaling factor. I tried to play with the sampling frequency, and it seems to be dependent also on this, but I could not figure out looking into the documentation how to get rid of this scaling factor.
coefs = cwt(signal, scales, wavelet);
mag=abs(coefs);
ph=angle(coefs);
I am taking the magnitude and phase at the freqeuncy where I have the peak, of course, using max(). That's the one I want to extract.
  1 Commento
Abhimenyu
Abhimenyu il 12 Giu 2024
Hello,
It will really help if you could share the full code.
Thanks

Accedi per commentare.

Risposte (1)

Vinay
Vinay il 2 Set 2024
Hii Giovanni,
The discrepancy arises because the Continuous Wavelet Transform (CWT) inherently includes a scaling factor related to the wavelet's energy and the sampling frequency. This scaling affects the magnitude of the coefficients, making them not directly comparable to the time-domain amplitude.
We have to normalize the coefficients based on the wavelet's properties and the sampling rate using the below steps.
  • Calculate the center frequency of the Morse wavelet using the parameters gamma and beta.
  • Perform the CWT on the signal, obtaining wavelet coefficients and corresponding frequencies.
  • Normalize the wavelet coefficients by multiplying them with the square root of the ratio of frequency to sampling frequency
  • This normalization ensures that the energy distribution across scales is consistent, making the wavelet analysis less dependent on the sampling frequency.
Kindly refer to the following documentations
% sinusoidal signal
Fs = 2000; % Sampling frequency
t = 0:1/Fs:1-1/Fs; % Time vector from 0 to 1 second
f0 = 50; % Frequency of the sinusoid
amplitude = 1;
signal = amplitude * sin(2 * pi * f0 * t);
% Define Morse wavelet parameters
gamma = 3; % Gamma parameter
beta = 50; % Beta parameter
% center frequency of the Morse wavelet
centerFreq = (beta^(1/gamma)) / (2 * pi);
% CWT using the Morse wavelet
[cfs, freqs] = cwt(signal, 'morse', Fs, 'WaveletParameters', [gamma, beta]);
% scaling factor
scales = centerFreq ./ (freqs / Fs);
%
% Normalization of coefficients
normalizedCfs = cfs .* sqrt(freqs(:) / Fs);
% normalized magnitude scalogram
figure;
imagesc(t, freqs, abs(normalizedCfs));
axis xy;
xlabel('Time (s)');
ylabel('Frequency (Hz)');
title('Normalized CWT Magnitude Scalogram');
colorbar;
% frequency with the peak normalized magnitude
[~, idx] = max(max(abs(normalizedCfs), [], 2));
peakFreq = freqs(idx);
% Extract original magnitude and phase at the peak frequency
exactMag = abs(cfs(idx, :));
exactPhase = angle(cfs(idx, :));
% Extract normalized magnitude and phase at the peak frequency
normalizedMag = abs(normalizedCfs(idx, :));
normalizedPhase = angle(normalizedCfs(idx, :));
% Find the maximum magnitude and corresponding phase for exact and normalized values
[maxExactMag, maxExactIdx] = max(exactMag);
maxExactPhase = exactPhase(maxExactIdx);
[maxNormalizedMag, maxNormalizedIdx] = max(normalizedMag);
maxNormalizedPhase = normalizedPhase(maxNormalizedIdx);
% Display the results
fprintf('Peak Frequency: %.2f Hz\n', peakFreq);
fprintf('Maximum Exact Magnitude: %.2f\n', maxExactMag);
fprintf('Phase at Maximum Exact Magnitude: %.2f radians\n', maxExactPhase);
fprintf('Maximum Normalized Magnitude: %.2f\n', maxNormalizedMag);
fprintf('Phase at Maximum Normalized Magnitude: %.2f radians\n', maxNormalizedPhase);

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by