phd.m

来自「matlab谱分析工具箱」· M 代码 · 共 41 行

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: phd.m
%
% AUTHOR: Steve Kogon
% 
% DATE: January 18, 1999
%
% DESCRIPTION: This file forms an estimate of the frequency spectrum using 
%              Pisarenko harmonic decomposition.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------


function [fest,Rbar] = phd(x,P,Nfft)
% x = signal
% P = number of complex exponentials
% Nfft = number of FFT frequencies
M = P+1;              % time-window length
R = est_corr(x,M);    % estimate correlation matrix

% Compute eigendecomposition and order by descending eigenvalues
[Q0,D] = eig(R);
[lambda,index] = sort(abs(diag(D)));
lambda = lambda(M:-1:1);
Q=Q0(:,index(M:-1:1));

% Compute pseudo-spectrum
f = (-Nfft/2:(Nfft/2-1))/Nfft;                % FFT frequencies
Rbar = 1./abs(fftshift(fft(Q(:,M),Nfft))).^2; % pseudo-spectrum
[dummy,findex] = sort(Rbar);                  % sort pseudo-spectrum
fest = f(findex(Nfft:-1:(Nfft-P+1)));         % frequency estimates