PSD estimation FFT vs Welch

246 views (last 30 days)
Chris
Chris on 27 Mar 2012
Commented: CduP on 9 Apr 2020
I am attempting to compute PSD estimate for a signal, but am having trouble determining how to scale the output from my FFT and pwelch estimates. An example of my code is as follows:
%%Create Ensemble of Sample Signals
%Data Parameters
fs=400; %sample rate (Hz)
DelT=1/fs; %time between samples
ns=5000; %sample length
ms=10; %number of ensembles
%Time Vector
t_vector=(0:ns-1)*DelT;
%Signal (Sinusoids with zero mean random noise)
x=zeros(ms,ns);
for ims=1:ms
x(ims,:)=25*cos(50*2*pi*t_vector+rand*2*pi)+...
20*cos(40*2*pi*t_vector+rand*2*pi)+...
50*cos(14*2*pi*t_vector+rand*2*pi)+...
30*cos(20*2*pi*t_vector+rand*2*pi)+...
25*cos(10*2*pi*t_vector+rand*2*pi)+...
20*cos(100*2*pi*t_vector+rand*2*pi)+...
50*cos(90*2*pi*t_vector+rand*2*pi)+...
30*cos(80*2*pi*t_vector+rand*2*pi)+...
25*cos(70*2*pi*t_vector+rand*2*pi);
% +1*randn(1,length(t_vector));
end
plot(t_vector(1:200),x(:,1:200))
%%Power Spectral Density (Welch's Method)
for ims=1:ms
XWELCH(ims,:)=pwelch(x(ims,:),[],[],fs);
end
XbarWELCH=mean(XWELCH,1);
figure
plot(XbarWELCH)
title('PSD pwelch Method')
%%Power Spectral Density (Fourrier)
%frequency vector
f_vector=(0:ns/2-1)*fs/ns;
for ims=1:ms
XFFT(ims,:)=abs(fft(x(ims,:)*hamming(length(ns)))).^2;
end
XbarFFT=mean(XFFT,1);
% XbarFFToneside=XbarFFT(1:floor(length(XbarFFT)/2)).^2*2*DelT/ns;
XbarFFToneside=XbarFFT(1:floor(length(XbarFFT)/2))*2*DelT/ns;
figure
plot(f_vector,abs(XbarFFToneside))
title('PSD fourier Method')
If anyone can explain to me how my pwelch and FFT PSD estimates should be scaled in order to give correct and approximately matching results, it would be greatly appreciated. Also, I want to better understand the physical units of the PSD for an input signal in terms of acceleration.

Answers (1)

Paul Pacheco
Paul Pacheco on 30 Mar 2012
Here's an example showing how to calculate the PSD with FFT. The results are identical to the results of using pwelch.m.
Fs = 1024;
t = (0:1/Fs:1-1/Fs).';
x = sin(2*pi*t*200);
Nx = length(x);
% Window data
w = hanning(Nx);
xw = x.*w;
% Calculate power
nfft = Nx;
X = fft(xw,nfft);
mx = abs(X).^2;
% Normalize by window power. Multiply by 2 (except DC & Nyquist)
% to calculate one-sided spectrum. Divide by Fs to calculate
% spectral density.
mx = mx/(w'*w);
NumUniquePts = nfft/2+1;
mx = mx(1:NumUniquePts);
mx(2:end-1) = mx(2:end-1)*2;
Pxx1 = mx/Fs;
Fx1 = (0:NumUniquePts-1)*Fs/nfft;
[Pxx2,Fx2] = pwelch(x,w,0,nfft,Fs);
plot(Fx1,10*log10(Pxx1),Fx2,10*log10(Pxx2),'r:');
legend('PSD via FFT','PSD via pwelch')
You can find more details about scaling the FFT in following link:
HTH
-Paul
  2 Comments
CduP
CduP on 9 Apr 2020
Hi! I'm just wondering why you haven't divided X by the length of the signal?

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!