- pwelch: https://www.mathworks.com/help/signal/ref/pwelch.html
- fft: https://www.mathworks.com/help/matlab/ref/fft.html
periodgram vs pwelch not showing similar results
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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))))
0 comentarios
Respuestas (1)
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:
Hope this helps.
0 comentarios
Ver también
Categorías
Más información sobre Spectral Estimation en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!