Estimate PSD from CWT

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

Answer 1

[Drishti]

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.

https://www.mathworks.com/help/releases/R2024a/signal/ref/pwelch.html

https://www.mathworks.com/help/releases/R2024a/wavelet/ref/wavedec.html

I hope this helps!