Magnitude and phase of the signal using cwt
15 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentario
Respuestas (1)
Vinay
el 2 de Sept. de 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 comentarios
Ver también
Categorías
Más información sobre Continuous Wavelet Transforms en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!