# Power Spectrum estimate using pwelch function

7 views (last 30 days)
Matthew Casiano on 9 Feb 2021
I am looking into using the pwelch function to estimate the power spectrum. As I understood, the power spectrum is simply the PSD multiplied by the frequency bandwidth (sampling rate divided by the block size). So I decided to check and compare the pwelch estimate with the PSD*BW approach. I found that for a rectangular window (no window), the methods match exactly. But if I use a window function, then the methods do not match. The ratio that they are off isn’t the window correction factor, but maybe it is related: for example with a hann window the estimated power spectra differ by a factor of 1.5x and for a hamming window it is a factor of 1.36x. Has anyone encountered anything like this or can it be explained? Attached is an example overlaying the two methods for the rectangular window and hamming window.
%% Inputs
fs=10000; % sampling rate (sps)
t1=1; % specified start time (s)
t2=9; % specified end time (s)
BS=1024; % block size (pts)
nfft=BS; % DFT block size (pts)
OL=0.5; % percent overlap (%/100)
%% Example Data
TimeData=0:1/fs:10; % time vector (s)
Data=rand(1,length(TimeData)); % data vector (EU)
tStart=TimeData(1); % data start time
%% PSD/PS Calculation
BW=fs/BS; % bandwidth (Hz)
df=fs/nfft; % frequency resolution (Hz)
ptsOL=floor(OL*BS); % number of points for overlap (pts)
PSDStartIndex=round((t1-tStart)*fs)+1; % find index closest to specified PSD start time
PSDEndIndex=round((t2-tStart)*fs)+1; % find index closest to specified PSD end time
[pxx_rect1,f_rect1]=pwelch(Data(PSDStartIndex:PSDEndIndex),ones(1,BS),ptsOL,nfft,fs);
[pxx_rect2,f_rect2]=pwelch(Data(PSDStartIndex:PSDEndIndex),ones(1,BS),ptsOL,nfft,fs,'power');
[pxx_hamming1,f_hamming1]=pwelch(Data(PSDStartIndex:PSDEndIndex),hamming(BS),ptsOL,nfft,fs);
[pxx_hamming2,f_hamming2]=pwelch(Data(PSDStartIndex:PSDEndIndex),hamming(BS),ptsOL,nfft,fs,'power');
%% Plotting
% Power Spectrum Plot - Rectangular Window
figure(1)
semilogy(f_rect1,pxx_rect1*BW)
hold on
semilogy(f_rect2,pxx_rect2)
hold off
grid on
set(gca,'xlim',[0 fs/2])
set(gca,'ylim',[10^-4 10^-3])
xlabel(gca, 'Frequency (Hz)')
ylabel(gca, 'Amplitude (EU-rms)^2')
legend({'PSD*BW','Power'}, 'Location','NorthEast')
title(gca,{['\fontsize{14}' 'Power Spectrum - Rectangular Window' '\fontsize{10}']},'fontweight','bold')
% Power Spectrum Plot - Hamming Window
figure(2)
semilogy(f_hamming1,pxx_hamming1*BW)
hold on
semilogy(f_hamming2,pxx_hamming2)
hold off
grid on
set(gca,'xlim',[0 fs/2])
set(gca,'ylim',[10^-4 10^-3])
xlabel(gca, 'Frequency (Hz)')
ylabel(gca, 'Amplitude (EU-rms)^2')
legend({'PSD*BW','Power'}, 'Location','NorthEast')
title(gca,{['\fontsize{14}' 'Power Spectrum - Hanning Window' '\fontsize{10}']},'fontweight','bold')

R2020b

### Community Treasure Hunt

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

Start Hunting!