- 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
5 views (last 30 days)
Show older comments
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 Comments
Answers (1)
Binaya
on 8 Oct 2024
Edited: Binaya
on 8 Oct 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 Comments
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!