# PSD estimation FFT vs Welch

246 views (last 30 days)
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.

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 CommentsShowHide 1 older comment
CduP on 9 Apr 2020
Hi! I'm just wondering why you haven't divided X by the length of the signal?