Estimate Power Spectral Density (PSD)

Hello! I want to estimate the PSD of a signal from the CWT coefficients. How can I do it? Can you provide me with some algorithm to do this? I have Matlab R2016a installed. Thank you so much!

[Manikanta Aditya]

Hey,

Check this example code, I feel this can help.

% Assuming 'cwtCoeffs' is your matrix of CWT coefficients and 'Fs' is your sampling frequency

% 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, which is the PSD
psd = scalogram / sum(scalogram);

% Plot the PSD
frequencies = Fs * (0:(length(psd)/2))/length(psd);
plot(frequencies, 10*log10(psd(1:length(frequencies))));
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB/Hz)');

[Gisela Clemente]

Thanks for your answer. But when I compare your CWT-PSD with the clasical FFT-PSD, they aren`t similar. It is the code:

% Constants and Signal Definition

Fs = 1000; % Sampling frequency in Hz

duracion = 1; % Duration in seconds

N = Fs * duracion; % Total number of samples

delta = zeros(1, N); % Initialize delta function

delta(round(N/2)) = sqrt(N); % Impulse in the center

scales=1:16;

cwtCoeffs = cwt(delta, 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, which is the PSD

psd = scalogram / sum(scalogram);

% Plot the PSD

frequencies = Fs * (0:(length(psd)/2))/length(psd);

plot(frequencies, 10*log10((psd(1:length(frequencies)))));

xlabel('Frequency (Hz)');

ylabel('Power Spectral Density (dB/Hz)');

% FFT-based PSD

X = fft(delta);

Pxx1_delta = 2*abs(X(1:N/2)).^2/(N*Fs); % PSD computation

frecuencias = linspace(0, Fs/2, N/2); % Frequency vector for FFT

hold on

plot(frecuencias, 10*log10(Pxx1_delta))