%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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
% 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; % 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;
% Find local maxima (values of Rbar that are larger than their neighbors)
z1 = Rbar(2:(Nfft-1)) - Rbar(1:(Nfft-2));
z2 = Rbar(2:(Nfft-1)) - Rbar(3:Nfft);
peak_index = find((z1 > 0) & (z2 > 0)) + 1;
if Rbar(1) > Rbar(2)
peak_index = [1
peak_index];
end
if Rbar(Nfft) > Rbar(Nfft-1)
peak_index = [peak_index
Nfft];
end
Npeaks = length(peak_index);
f_peaks = f(peak_index);
% Determine the P largest peaks (taken from local maxima)
[dummy,fest_index] = sort(Rbar(peak_index));
fest = f_peaks(fest_index(Npeaks:-1:(Npeaks-P+1))).';
