periodgram vs pwelch not showing similar results

4 visualizaciones (últimos 30 días)
Andres Antonio
Andres Antonio el 29 de Mzo. de 2019
Editada: Binaya el 8 de Oct. de 2024
I am having a hard time trying to get similar answers for a power spectral density doing it using the fft and the pwelch function. I know that because of the window used in pwelch, the answers will not be the same, however I know that at least both plots should be similar. Although I can see that for both methods the peaks happen at the same frequencies, the amplitudes are different. One method shows peaks way above the other. Can somebody tell me what is wrong with my code?
Thank you
clear
dt = 0.01;
y=0:dt:1000;
X = sin(5*y)+cos(3*y);
%%% pwelch
x = X;
Fs = 1/dt;
Fn = Fs/2;
PSDx = pwelch(x);
Fv = linspace(0, 1, fix(length(PSDx)))*Fn;
Iv1 = 1:length(Fv)/20;
%%%fft
L = length(X);
fftx = fft(X);
P2 = abs(fftx);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
psdx = (1/(L*Fs))*abs(P1).^2;
f = linspace(0, 1, fix(length(psdx)))*Fn;
Iv2 = 1:length(f)/20;
%%%plot both
plot(f(Iv2),P1(Iv2))
hold on
plot(Fv(Iv1),(abs(PSDx(Iv1))))

Respuestas (1)

Binaya
Binaya el 8 de Oct. de 2024
Editada: Binaya el 8 de Oct. de 2024
Hi Andreas
I understand that you are trying to plot the power spectral density using two different approaches. The 'pwelch' function directly gives you the power spectral density while 'fft' gives you the power spectrum. To obtain the power spectral density from the fft output, appropriate normalization must be applied.
The following modified code includes the changes to calculate power spectral density using both the approaches:
clear
dt = 0.01;
y = 0:dt:1000;
X = sin(5*y) + cos(3*y);
%%% pwelch
x = X;
Fs = 1/dt;
Fn = Fs/2;
[PSDx, Fv] = pwelch(x, [], [], [], Fs);
%%% fft
L = length(X);
fftx = fft(X);
P2 = abs(fftx/L).^2; % Normalize by the length of the signal
P1 = P2(1:floor(L/2)+1);
P1(2:end-1) = 2*P1(2:end-1); % Double the amplitude for positive frequencies
psdx = P1 * Fs; % Scale by the sampling frequency
f = linspace(0, Fn, floor(L/2)+1);
%%% plot both
Iv1 = 1:floor(length(Fv)/20);
Iv2 = 1:floor(length(f)/20);
figure;
plot(f(Iv2), 10*log10(psdx(Iv2)), 'b', 'DisplayName', 'FFT-based PSD')
hold on
plot(Fv(Iv1), 10*log10(PSDx(Iv1)), 'r', 'DisplayName', 'pwelch PSD')
legend show
The spectrum is plotted using the logarithmic scale to make the comparison between the plots more prominent due to the large differences in the magnitude.
You can refer to the following documentation links for more details on:
  1. pwelch: https://www.mathworks.com/help/signal/ref/pwelch.html
  2. fft: https://www.mathworks.com/help/matlab/ref/fft.html
Hope this helps.

Categorías

Más información sobre Spectral Estimation en Help Center y File Exchange.

Productos


Versión

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by