PWELCH vs PSD

106 views (last 30 days)
Kevin
Kevin on 18 Oct 2011
Commented: ramon mata on 27 May 2021
Hello,
I am trying to update some code that uses the deprecated function PSD to use PWELCH instead; however, I am getting completely different results when I am using what seems to be the same parameters. Does anyone know why? Here are what I believe to be the equivalent function calls:
[pxx,f] = pwelch(data, hanning(1024), [], 1024, 250);
[pxx,f] = psd(data,1024,250,hanning(1024));
Where: data = signal in a vector
hanning(1024) = window vector
1024 = nfft, number of points in the fft
250 = Fs, the sampling rate
[ ] = noverlap, defaults to 50% window overlap
Matlab has removed all help information for the PSD function, and instead says to use its functional equivalent pwelch, so I don't have anyway of looking up what the original documentation says about the function's inputs and outputs. Could the outputs be scaled differently? Is the PSD calculated differently between the two functions?
Any help would be greatly appreciated.
Kevin
  1 Comment
ramon mata
ramon mata on 27 May 2021
Sorry, could you give us a theorical background about this? i've been looking for one, but i can't find a theorical background about this scale problem and what is DC? If Nyquist frecuency is too high why we have several problems at the begins?

Sign in to comment.

Accepted Answer

Honglei Chen
Honglei Chen on 19 Oct 2011
Hi Kevin,
The result of psd is not correctly scaled, therefore you should use pwelch. For spectral density, the result should be scaled by the sampling frequency, which is not performed by psd.
If you look at the two results, the f vector should be the same. If you take the ratio of two pxx, you can that most samples, the differ by a factor of 125, which is basically half of sampling frequency. This is because the resulting spectrum is one-sided. The two end points differs by a factor of 250. This is due to the fact that these two points correspond to DC and Nyquist frequency and should not be doubled even if it's one-sided.
BTW, for your case, there is really no overlap happening because your data and window length are the same.
HTH
  2 Comments
Kevin
Kevin on 27 Mar 2012
Thank you very much Honglei for your response. However, I have a question about your last comment about there being no overlap. You said the data and the window length are the same, but I never explicitly said how long the data was. Did you mean that there is no overlap because the nfft is the same size as the window? Or, did you think that my variable "data" is the same length my window? The length of "data" is actually at least 10,000 samples long. I just want to make sure that I am using the correct parameter values to make sure that my samples have 50% overlap.
Thanks again!
PIYUSH MOHANTY
PIYUSH MOHANTY on 1 Apr 2019
Hi Kevin,
Thanks for sharing the knowledge. Can you please let me know the resource; where I can have good knowledge about the dependednce of PSD or pwelch on sampling frequency. I have carried out some experiments and have acquired good data bout foundations of structures. I need to work on finding out natural frequency of the foundations.
Cheers,
Piyush

Sign in to comment.

More Answers (5)

Paul Pacheco
Paul Pacheco on 30 Mar 2012
Here's a small tweak to Rick's code which takes into consideration the overlap value and the fact that the DC and Nyquist values should not be scaled. Also, I scaled the results of psd.m instead of pwelch.m. Now Pxx and Pyy are identical.
load handel;
N = 1024;
[Pxx,Fx] = psd(y,N,Fs,hanning(N));
Pxx(2:end-1) = Pxx(2:end-1)*2;
Pxx = Pxx/Fs;
[Pyy,Fy] = pwelch(y, hanning(N),0, N, Fs);
plot(Fx,10*log10(Pxx),Fy,10*log10(Pyy));
legend('psd','pwelch');
-Paul

Rick Rosson
Rick Rosson on 19 Oct 2011
Honglei is correct. I created a simple test script to compare the two, taking the scaling factor into account.
Please try the following:
load handel;
N = 1024;
[Pxx,Fx] = psd(y,N,Fs,hanning(N));
[Pyy,Fy] = pwelch(y, hanning(N), [], N, Fs);
Pyy = Pyy*Fs/2;
figure;
plot(Fx,10*log10(Pxx),Fy,10*log10(Pyy));
legend('psd','pwelch');
HTH.
Rick
  1 Comment
PIYUSH MOHANTY
PIYUSH MOHANTY on 1 Apr 2019
Hi Rick,
Thanks for sharing the knowledge. Can you please let me know the resource; where I can have good knowledge about the dependence of PSD or pwelch on sampling frequency. I have carried out some experiments and have acquired good data about foundations of structures. I need to work on finding out natural frequency of the foundations.
Cheers,
Piyush

Sign in to comment.


Chad
Chad on 4 Sep 2013
I recreated the pwelch algorithm! :) (ignore the self-test features)
if true
% function varargout = pwelch(data_in,Fs,NFFT,Noverlap,window)
%
%
%SPECTRAL ESTIMATION BY WELCH'S METHOD OF AVERAGED PERIODOGRAMS
%
%Pwelch will take a matrix of time and data vectors and compute the
%corresponding Power Spectral Density estimation.
%
%1) Data_in must contain doubles %data points to be analyzed
%2) Sampling Frequency (Fs) must be a 1x1 double
%3) NFFT must be an integer, preferably a power of 2 %Number of FFT points per bin
%4) Noverlap must be between 0 and 1 %Percentage overlap of
% %each segmented data
%
%Input: data_in, Fs, NFFT, Noverlap
%
%Output: varargout{1} = PSD estimate
% varargout{2} = frequency bins
%
%Example:
%
% t = 0:.1:100;
% x = cos(2*pi*t);
% Fs = 1;
% NFFT = 512;
% Noverlap = .5 %50% overlap
%
%
%pxx = pelch(x,Fs,NFFT,Noverlap); %returns the values of the
% %spectral density estimation%
%
%[pxx,f] = pwelch(x,Fs,NFFT,Noverlap); %returns two vectors where f is an
% %index of bin frequencies
%
%semilogy(f,pxx) %plots the data on a logarithmic scale.
%
%
varargout{1} = [];
varargout{2} = [];
if ~exist('data_in','var')
help(mfilename);
return
end
if ischar(data_in)
help(mfilename);
errordlg('Data input must be of type doubles',mfilename);
return
end
full_file_path = which('pwelch');
if isempty(full_file_path)
errordlg(['Unable to find "',full_file_path,'".'],mfilename);
return
end
if nargin < 5 help(mfilename); errordlg('Not enough input parameters. Please see "pwelch" help above.',mfilename); return end
if nargin > 5 help(mfilename); errordlg('Too many input parameters. Please see "pwelch" help above.',mfilename); return end
if isempty(data_in) help(mfilename);
t = 1:.01:100;
data_in = cos(2*pi*10*t);
NFFT = 512;
disp('**************************************************************');
disp('* SELF-TEST *');
disp('**************************************************************');
disp('* *');
disp('***********See resulting PSD estimate in figure***************');
disp('********PSD of cosine function with frequency 10 Hz***********');
disp('**************************************************************');
[pxx, f] = pwelch_chad(data_in,100,NFFT,0);
semilogy(f,pxx);
return
%error('Data Input Not Right Format. Please See "pwelch" help above.');
end
if NFFT > length(data_in) Pxx1 = []; errordlg('FFT points exceeds data points. Automatically set to number of data points.') NFFT = length(data_in); end
T = 1/Fs; %Sample interval (sec/sample)
N = length(data_in); %number of data points
bin = 1/(NFFT*T); %frequency gap between successive data points (Hz)
if strcmpi(window,'hanning')
w = hanning(NFFT-1); %standard hanning window
elseif strcmpi(window,'hamming')
w = hamming(NFFT - 1);
elseif strcmpi(window,'blackman')
w = blackman(NFFT - 1);
elseif strcmpi(window,'rectangular')
w = rectangular(NFFT - 1);
else
help(mfilename);
return
end
s = min(size(data_in));
w = repmat(w,1,s);
index = 1:NFFT; %index of data in first segment
[m,n] = size(data_in);
if n > m
data_in = data_in.'; %column to rows
%data_in = detrend(data_in,'constant')
%data_in(index,:) = detrend(data_in(index,:),'constant');
%else
%data_in = detrend(data_in,'constant')
%data_in(index,:) = detrend(data_in(index,:),'constant');
end
k = (fix((N-NFFT)/(NFFT*(1-Noverlap)))) + 1; %number of windows to be applied
P1 = zeros(NFFT,s);
% SCN = (k*norm(w)).^2; %normalized scaling factor
for j = 1:k
XW1 = w.*data_in(index,:);
index = index + fix((NFFT*(1 - Noverlap)));
P1 = P1 + (abs(fft(XW1,NFFT))).^2; %frequency bins
end
%Compensate for windowing.
window_mean_sq = mean(w.^2);
window_mean_ext = repmat(window_mean_sq,NFFT,1);
P1 = P1 ./ window_mean_ext;
P1 = P1(1:((NFFT/2) + 1),:);
% Average power by number of windows.
P1 = P1 / k;
% Normalize by sampling rate & window length.
P1 = P1 / (Fs * NFFT);
% Compensate for windowing.
% P1 = P1 / w.^2;
% Account for double-sided nature of FFT.
% (But the end points--zero frequency & Nyquist frequency--aren't doubled.)
P1(2:end-1) = 2 .* P1(2:end-1);
varargout{1} = P1;
%
% Pxx1 = (P1((1:(NFFT/2)+1),:))/(bin*SCN); %PSD estimate
% varargout{1} = Pxx1();
q = 0:fix(NFFT/2);
varargout{2} = q*bin;
%n = max(size(Pxx1));
%F = (0:n-1)*1/NFFT/Fs;
end
end

Kevin
Kevin on 29 Mar 2012
Thank you very much Honglei for your response. However, I have a question about your last comment about there being no overlap. You said the data and the window length are the same, but I never explicitly said how long the data was. Did you mean that there is no overlap because the nfft is the same size as the window? Or, did you think that my variable "data" is the same length my window? The length of "data" is actually at least 10,000 samples long. I just want to make sure that I am using the correct parameter values to make sure that my samples have 50% overlap.

Honglei Chen
Honglei Chen on 29 Mar 2012
Hi Kevin,
I must somehow interpreted that your data is only 1024 samples long, so it is only enough for one window, that's why I said there is no overlap. If you data is 10,000 samples long then yes, the 50% overlap will happen by default. Sorry about the confusion.

Community Treasure Hunt

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

Start Hunting!