PWELCH vs PSD
95 ビュー (過去 30 日間)
古いコメントを表示
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 件のコメント
ramon mata
2021 年 5 月 27 日
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?
採用された回答
Honglei Chen
2011 年 10 月 19 日
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 件のコメント
PIYUSH MOHANTY
2019 年 4 月 1 日
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
その他の回答 (5 件)
Paul Pacheco
2012 年 3 月 30 日
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
2011 年 10 月 19 日
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 件のコメント
PIYUSH MOHANTY
2019 年 4 月 1 日
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
Chad
2013 年 9 月 4 日
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
0 件のコメント
Honglei Chen
2012 年 3 月 29 日
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.
0 件のコメント
参考
カテゴリ
Help Center および File Exchange で Parametric Spectral Estimation についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!