PSD estimation FFT vs Welch

234 views (last 30 days)
Chris
Chris on 27 Mar 2012
Commented: Raphael Krier on 24 Jun 2022
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.
  1 Comment
Raphael Krier
Raphael Krier on 24 Jun 2022
Hi there, Welch and FFT are very different by nature. The FFT is used to get the spectral estimate over the netire signal but it is sensitive to non stationarity. So if you want to have a better estimate for signal with non stationary components, use Welch. Welch spectra breaks down the signal in segment and use a hanning function. The number of segments in Welch depends on the overlap you use and the length of your segment. From this you can define your degrees of freedom (DOF). The welch spectra has a lower frequency resolution fompare to a simple fft as you are averaging a series of fft obtained on each segment. The DOF is conventionaly used to tell how statistically robust is your welch sepctra (e.g. DOF>16 is considered as a good start in ocean wave statistics).
The following script breaks down the welch method. so you can see how it is rlated to the FFT. This is based on welch (1967).
I hope this helps
samprate=fs; %sampling frequency
nfft=length(Waveburst); % enter here any time series quantity (here we have 2048 sampes rows in one burst column)
% get the corresponding time vector (timestamp of the record)
dt=1/fs;%Time or space increment per sample
t = (0:nfft-1)/fs; %Time of 1 burst in decimal minutes
% segmentation in averaging windows
xwin=512; %length of averaging segment
noverlap=xwin*0.5; % number of samples corresponding to a 50% overlap (0.5)
K=(nfft-noverlap)/(xwin-noverlap); % Number of segment that data were divided in
% obtain modified periodogram for each segment
HanningW=hanning(xwin); %generate the anning indow for the given segment length
centreseg=noverlap; % initial centre of the segment
for i=1:K % get the fft on each modified periodogram. This will ony work for 50% overlap and if the signal length is a multibple of the window length
Ak(:,i)=(1/xwin).*fft((Waveburst(centreseg-noverlap+1:centreseg+noverlap,1).*HanningW),xwin); % periodogram are obtained by multiplying the hanning window to each data segment then we realised and fft to obtaine modified periodogram
centreseg=centreseg+noverlap; %shift the centre of the segment at the end of the previous segment;
end
U=(1/xwin).*sum((HanningW).^2); % we need to correct the spectra estimate as we applied a window
Ik=(xwin/U).*abs(Ak).^2; % spectra estimate of each periodogram
%average the periodograms to obtain thepriodogram of the entire signal.
%So we obtain the magnitude of each frequency component in the signal.
%However when k>N/2 the sequences reapeat (2 sided spectra)
P=(1/K).*sum(Ik,2)';
% apply Nyquist limit to avoid repetion and obtained 1 sided spectra
Nyquist = xwin/2+1;
P=P(1:Nyquist); % truncate the spectral estimate within the boundary of the Nyquist
P(2:end-1) = P(2:end-1)*2; % now to obtain the amplitude of each frequency, we need to multiply the magnitude of the 1 sided spectra by 2 as we shiched from 2 sided to one sided spectra
Fx1 = (0:Nyquist-1)*samprate/xwin; % these are the corresponding frequencies for thhe PSD in P
Pxx1 = P/samprate; % the is the power spectra density for each frequency
% compare to welch
figure
plot(Fx1,Pxx1)
hold on
[psd freq]=pwelch(Waveburst(:,1),xwin,noverlap,nfft,samprate); %use hamming window by default
plot(freq,psd,'r')

Sign in to comment.

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
jack star
jack star on 14 May 2016
Hello. Is this program has no overlap or %50 overlap?
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!