rootmusic.m

<统计智能信号处理>一书的matlab源程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: rootmusic.m
%
% AUTHOR: Steve Kogon
% 
% DATE: January 18, 1999
%
% DESCRIPTION: This file forms an estimate of the frequency spectrum using 
%              root-MUSIC algorithm (Barabell 1983). 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%-----------------------------------------------------------
% 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] = rootmusic(x,P,M)
% x = signal
% P = number of complex exponentials
% M = time-window length
% Nfft = number of FFT frequencies

% Generate data matrix
N = length(x) - M + 1;
X = zeros(N,M);
for n = 1:M
   X(:,n) = x((1:N)+ (M-n));
end
R = (1/N)*X'*X;

% 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
A = zeros(2*M-1,1);
for n = 1:(M-P)
   A = A + conv(Q(:,M-(n-1)),conj(Q(M:-1:1,M-(n-1)))); 
end
r_A = roots(A);
[i_min] = find(abs(r_A) < 1);
r_A_min = r_A(i_min);
freq = angle(r_A_min)/(2*pi);
[r_A_order,index] = sort(abs((abs(r_A_min) -1)));
fest = freq(index(1:P));         % frequency estimates