Azzera filtri
Azzera filtri

Estimate the psd of a signal

9 visualizzazioni (ultimi 30 giorni)
Gisela Clemente
Gisela Clemente il 5 Apr 2024
In the following code I calculate the PSD of my signal, a delta, using the FFT. Now, I want to estimate this PSD using the CWT coefficients, using the Global Wavelet Spectrum. But the graphics that remain are not comparable with the original PSD. I need the decomposition to be in the first 16 scales, with the db6 wavelet:
Fs = 1000; % Hz
duracion = 1;
N = Fs * duracion;
t = linspace(-duracion/2, duracion/2, N);
delta = zeros(1, N);
delta(round(N/2)) = sqrt(N);
cantidad_muestras = length(delta);
N = length(delta);
X = fft(delta);
Pxx1_delta = 2*abs(X(1:N/2)).^2/(N*Fs); % PSD
frecuencias = linspace(0, Fs, cantidad_muestras);
frecuencias=frecuencias(1:500);
delta_time=1/Fs;
scales=1:16;
freqs=scal2frq(scales,'db6',delta_time); % computes pseudo frequencies from the scales ..
% Wavelet transform and power spectrum
coefs_base=cwt(delta,scales,'db6'); % computes the Wavelet Transform for dataset at scales
pow_base=(abs(coefs_base)).^2; % computes the Wavelet Time-Frequency Power Spectrum for dataset
% bias correction
bias_corr_pow_base=ones(length(scales),length(t));
for k=1:length(scales)
bias_corr_pow_base(k,:)=pow_base(k,:)/(scales(k));
% after the work of C. Torrence and G. Compo
end
% Time averaged wavelet power spectrum
time_avg_cor_pow_base_spec=mean(bias_corr_pow_base,2);
figure;
plot(freqs,time_avg_cor_pow_base_spec)
title('Bias Corrected Time Averaged Power Spectrum Dataset');
hold on
plot(frecuencias, Pxx1_delta)
return

Risposta accettata

Hassaan
Hassaan il 5 Apr 2024
% 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
% 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
% CWT-based PSD
delta_time = 1/Fs;
scales = 1:16; % Define scales
freqs = scal2frq(scales, 'db6', delta_time); % Convert scales to frequencies
coefs_base = cwt(delta, scales, 'db6'); % CWT
pow_base = (abs(coefs_base)).^2; % Power spectrum from CWT coefficients
bias_corr_pow_base = pow_base ./ scales'; % Correcting for scale bias
time_avg_pow_base_spec = mean(bias_corr_pow_base, 2); % Global Wavelet Spectrum
% Interpolate Wavelet Spectrum to match FFT frequencies
interp_wavelet_spectrum = interp1(freqs, time_avg_pow_base_spec, frecuencias, 'linear', 'extrap');
% Plotting
figure;
loglog(frecuencias, interp_wavelet_spectrum, 'LineWidth', 2, 'DisplayName', 'Interpolated Wavelet PSD');
hold on;
loglog(frecuencias, Pxx1_delta, 'LineWidth', 2, 'DisplayName', 'FFT PSD');
xlabel('Frequency (Hz)');
ylabel('Power');
title('Comparison of FFT and Interpolated Wavelet-Based PSD');
legend('show');
grid on;
-----------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.
Professional Interests
  • Technical Services and Consulting
  • Embedded Systems | Firmware Developement | Simulations
  • Electrical and Electronics Engineering
Feel free to contact me.
  1 Commento
Gisela Clemente
Gisela Clemente il 5 Apr 2024
Thanks for the reply! I don't understand why PSD calculated from the CWT does not converge to the traditional PSD. How do you understand these results?

Accedi per commentare.

Più risposte (0)

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by