PWELCH vs PSD

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

[ramon mata]

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?

Answer 1

[Honglei Chen]

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

[Kevin]

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]

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

Answer 2

[Paul Pacheco]

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

[zohar]

+1

Answer 3

[Rick Rosson]

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

[PIYUSH MOHANTY]

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

Answer 4

[Chad]

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

Answer 5

[Kevin]

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.

Answer 6

[Honglei Chen]

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.