music.m

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: music.m
%
% AUTHOR: Steve Kogon
% 
% DATE: January 18, 1999
%
% DESCRIPTION: This file forms an estimate of the frequency spectrum using 
%              MUSIC algorithm (spectral version), (Schmidt 1986).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%-----------------------------------------------------------
% 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] = music(x,P,M,Nfft)
% x = signal
% P = number of complex exponentials
% M = time-window length
% Nfft = number of FFT frequencies
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
Qbar = zeros(Nfft,1);
for n = 1:(M-P)
   Qbar = Qbar + abs(fftshift(fft(Q(:,M-(n-1)),Nfft))).^2;
end
Rbar = 1./Qbar;
[dummy,findex] = sort(Rbar);                  % sort pseudo-spectrum
fest = f(findex(Nfft:-1:(Nfft-P+1)));         % frequency estimates