Introduction to Statistical and Adaptive Signal Processing – Random Signals – Spectral Estimation

				%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
				%
				% 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))).';