Estimate PSD from CWT

5 views (last 30 days)
Gisela Clemente
Gisela Clemente on 3 Apr 2024
Answered: Drishti on 10 Oct 2024
Good afternoon! I want to estimate the power spectral density (psd) using the cwt coefficients. I have read that a good estimator is the global spectrum wavelet (gws). But I have the following problem: the sample rate is 1000Hz, so the original psd goes from 0 to 500hz. When I convert the scales of my interest (from 1 to 16), with the db6 wavelet, to frequencies, these are 45 to 700 Hz. How do I compare the traditional psd with gws? Any algorithm to help me? thank you

Answers (1)

Drishti
Drishti on 10 Oct 2024
Hi Gisela,
I understand that you are trying to estimate Power Spectral Density (PSD) and wants to compare the traditional PSD with Global Wavelet Spectrum (GWS) using ‘db6’ wavelet.
To estimate PSD, you can utilize the ‘pwelch’ function that follows Welch’s method.
Refer to the code snippet below for estimating PSD:
% Compute traditional PSD using Welch's method
[pxx, f] = pwelch(signal, [], [], [], fs);
Also, for the generation of GWS, you can utilize the ‘wavedec’ function with ‘db6’ wavelet.
Refer to the implementation below for better understanding:
% Perform DWT using db6
maxLevel = 16; % Maximum level of decomposition
[coeffs, levels] = wavedec(signal, maxLevel, 'db6');
% Calculate the corresponding frequencies for each scale
center_frequency = 0.666 * (fs / 2);
% Calculate frequencies for each scale
frequencies = center_frequency ./ (2.^(0:maxLevel-1));
% Estimate PSD from DWT coefficients (GWS)
gws = zeros(size(frequencies));
for i = 1:maxLevel
% Reconstruct signal from detail coefficients at each level
reconstructed = wrcoef('d', coeffs, levels, 'db6', i);
% Compute power of reconstructed signal
gws(i) = sum(reconstructed.^2);
end
After calculating the values of both PSD and GWS, you can compare them by simply plotting the values.
For more information, you can refer to the MATLAB Documentation of ‘pwelch’ and ‘wavedec’ functions.
I hope this helps!

Tags

Community Treasure Hunt

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

Start Hunting!