Magnitude and phase of the signal using cwt
15 views (last 30 days)
Show older comments
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 Comment
Answers (1)
Vinay
on 2 Sep 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
- ‘cwt’ function: https://www.mathworks.com/help/wavelet/ref/cwt.html
- ‘abs’ function: https://www.mathworks.com/help/matlab/ref/abs.html
% 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);
0 Comments
See Also
Categories
Find more on Continuous Wavelet Transforms in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!