Hi Gisela,
I see that you are trying to compute the Power Spectral Density (PSD) from the Continuous Wavelet Transform (CWT) coefficients and compare it with the PSD obtained from the Fast Fourier Transform (FFT). I believe the issue might be related to normalization and the fact that the CWT might not handle the impulse signal in the same way as the FFT, leading to a different energy distribution.
The FFT is inherently better suited for handling impulse signals, while the CWT provides a time-frequency representation that might not be as straightforward for such signals. I used the scal2frq function to map scales to frequencies accurately.
The following example code might help you understand this better:
clear all
close all
% Constants and Signal Definition
Fs = 1000; % Sampling frequency in Hz
t = 0:1/Fs:1-1/Fs; % Time vector (1 second duration)
% Signal: Sinusoid with noise
f_sin = 50; % Frequency of the sinusoid
signal = sin(2*pi*f_sin*t) + 0.5*randn(size(t)); % Sinusoid with Gaussian noise
% FFT-based PSD
N = length(signal); % Number of samples
X = fft(signal);
Pxx_fft = (1/(Fs*N)) * abs(X).^2; % PSD computation
Pxx_fft = Pxx_fft(1:N/2+1); % Take only the positive frequencies
Pxx_fft(2:end-1) = 2*Pxx_fft(2:end-1); % Adjust for the single-sided spectrum
frequencies_fft = linspace(0, Fs/2, N/2+1); % Frequency vector for FFT
% Plot FFT-based PSD
figure;
plot(frequencies_fft, 10*log10(Pxx_fft));
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('PSD from FFT');
% CWT-based PSD
scales = 1:128; % Use a broader range of scales for better frequency resolution
cwtCoeffs = cwt(signal, scales, 'db6'); % CWT
% Compute the absolute square of the CWT coefficients
cwtCoeffsSquared = abs(cwtCoeffs).^2;
% Compute the time-averaged power, which gives the scalogram
scalogram = mean(cwtCoeffsSquared, 2);
% Normalize the scalogram to get the relative energy distribution as a function of frequency
totalEnergy = sum(scalogram);
psd_cwt = scalogram / totalEnergy;
% Calculate the corresponding frequencies for the scales used in CWT
frequencies_cwt = scal2frq(scales, 'db6', 1/Fs);
% Plot the PSD
figure;
plot(frequencies_cwt, 10*log10(psd_cwt));
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
title('PSD from CWT');
% Compare both PSDs
figure;
plot(frequencies_cwt, 10*log10(psd_cwt), 'LineWidth', 1.5, 'DisplayName', 'CWT-based PSD');
hold on;
plot(frequencies_fft, 10*log10(Pxx_fft), 'LineWidth', 1.5, 'DisplayName', 'FFT-based PSD');
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');
legend('show');
title('Comparison of PSD Estimations');
You can see the peak at the sinusoidal frequency in both the FFT-based and CWT-based PSDs happen at the same frequency.
Kindly have a look at the following documentation links to have more information on:
- Continuous Wavelet Transform: https://www.mathworks.com/help/wavelet/ref/cwt.html
- Fast Fourier Transform: https://www.mathworks.com/help/matlab/ref/fft.html
Hope that helps!
Balavignesh
