How to change psd function to Periodogram
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi there,
I have a code written by psd function and would like to update it using Periodogram/Pwelch. I have tried but I could not find what should be f value as the input parameters. I wonder if you could please help me with that.
Thanks.
function varargout=psd2(varargin)
%PSD2 Power Spectral Density and frequency characteristics.
% PSD2 estimates the power spectral density (PSD), mean power frequency (MPF),
% median power frequency (F50),
% peak frequency (PEAK), and limit frequency (F95) that contains up to 95%
% of the PSD using Welch's averaged periodogram method (PSD Matlab function).
% [MPF,PEAK,F50,F95,F,P]=psd2(X,Fs,NFFT,WINDOW,NOVERLAP,DFLAG)
% [MPF,PEAK,F50,F95,F,P]=psd2(X,Fs) uses default values
% Inputs:
% X: data vector
% Fs: sampling frequency
% NFFT,WINDOW,NOVERLAP,DFLAG: see PSD matlab function for help
% default values: NFFT=1024, WINDOW=256, NOVERLAP=WINDOW/2, DFLAG='mean'
% Outputs:
% MPF: mean power frequency
% PEAK: peak frequency (mode)
% F50: median frequency of the power spectral density
% F95: frequency limit that contains up to 95% of the power spectral density
% F: frequency vector
% P: PSD estimated vector
%
% See also PSD
% Marcos Duarte mduarte@usp.br 11oct1998
% one typo and F50 comment modified slightly by Hyun Gu Kang
if nargin == 2
x=varargin{1};
fs=varargin{2};
if length(x)<1000
nfft=length(x); window=round(nfft/4);
else
nfft=1024; window=nfft/4;
end
noverlap=round(window/2); dflag='mean';
end
if nargin >= 3
if ~isempty(varargin{3})
nfft=varargin{3};
else
if length(x)<1000, nfft=length(x); else nfft=1024; end
end
end
if nargin >= 4
if ~isempty(varargin{4}), window=varargin{4}; else window=round(nfft/4); end
end
if nargin >= 5
if ~isempty(varargin{5}), noverlap=varargin{5}; else noverlap=round(window/2); end
end
if nargin >= 6
if ~isempty(varargin{6}), dflag=varargin{6}; else dflag='mean'; end
end
if nargin==1 | nargin > 6
error('Incorrect number of inputs')
return
end
%power spectral density:
[p,f]=psd(x,length(x),fs,length(x));
%mpf:
mpf=trapz(f,f.*p)/trapz(f,p);
%peak:
[m,peak]=max(p);
peak=f(peak);
%50 and 95% of PSD:
area=cumtrapz(f,p);
f50=find(area >= .50*area(end));
if ~isempty(f50)
f50=f(f50(1));
else
f50=0;
end
f95=find(area >= .95*area(end));
if ~isempty(f95)
f95=f(f95(1));
else
f95=0;
end
varargout{1}=mpf;
varargout{2}=peak;
varargout{3}=f50;
varargout{4}=f95;
varargout{5}=f;
varargout{6}=p;
1 Kommentar
Abhishek Kumar
am 4 Dez. 2020
Bearbeitet: Abhishek Kumar
am 4 Dez. 2020
Dear Hamed, I tried to run your script by replacing psd with pwelch, the code works fine could you please elaborate on what issue you are facing to be able to help you better?
Antworten (0)
Siehe auch
Kategorien
Mehr zu Parametric Spectral Estimation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!