function [w_i,powers,w_mse,p_mse] = rootmusic(x,p,varargin)%ROOTMUSIC Computes the frequencies and powers of sinusoids via the% Root MUSIC algorithm.% W = ROOTMUSIC(X,P) returns the vector of frequencies W of the complex% sinusoids contained in signal vector X. W is in units of rad/sample.% P is the number of complex sinusoids in X. If X is a data matrix,% each row is interpreted as a separate sensor measurement or trial.% In this case, X must have a number of columns larger than P. You can% use the function CORRMTX to generate data matrices to be used here. %% W = ROOTMUSIC(R,P,'corr') returns the vector of frequencies W, for a% signal whose correlation matrix estimate is given by the positive% definite matrix R. Exact conjugate-symmetry of R is ensured by forming% (R+R')/2 inside the function. The number of rows or columns of R must% be greater than P.%% If P is a two element vector, P(2) is used as a cutoff for signal and% noise subspace separation. All eigenvalues greater than P(2) times% the smallest eigenvalue are designated as signal eigenvalues. In % this case, the signal subspace dimension is at most P(1).%% F = ROOTMUSIC(...,Fs) uses the sampling frequency Fs in the computation% and returns the vector of frequencies, F, in Hz.%% [W,POW] = ROOTMUSIC(...) returns in addition a vector POW containing the% estimates of the powers of the sinusoids in X.%% EXAMPLES:% s1 = RandStream.create('mrg32k3a');% n=0:99; % s=exp(i*pi/2*n)+2*exp(i*pi/4*n)+exp(i*pi/3*n)+randn(s1,1,100); % X=corrmtx(s,12,'mod'); % Estimate the correlation matrix using% % the modified covariance method.% [W,P] = rootmusic(X,3); % % See also ROOTEIG, PMUSIC, PEIG, PMTM, PBURG, PWELCH, CORRMTX, SPECTRUM.% Reference: Stoica, P. and R. Moses, INTRODUCTION TO SPECTRAL ANALYSIS,% Prentice-Hall, 1997.% Author(s): R. Losada% Copyright 1988-2008 The MathWorks, Inc.% $Revision: 1.1 $ $Date: 2010/02/18 11:16:00 $%%%%%%%%%%%%%%%%%%%%%%%%%% Added function to compute approx. MSE for the case of a unique sinusoid% % REFERENCES: Rao, B. Performance Analysis of Root-Music% IEEE Trans. Acoust. Speech and Sig. Proc. 37, 1989%% VERSION: $Id: rootmusic.m,v 1.1 2010/02/18 11:16:00 miquel Exp $% % M Nofrarias 12/02/2010%error(nargchk(2,5,nargin,'struct'));xIsReal = isreal(x);% Check for an even number of complex sinusoids if data is realif xIsReal && rem(p,2), error(generatemsgid('InvalidDimensions'),'Real signals require an even number p of complex sinusoids.');endnfft = []; % Root Music doesn't use nfft, but the parser needs itvarargin = {nfft,varargin{:}};[md,msg] = utils.math.music(x,p,varargin{:});if ~isempty(msg), error(generatemsgid('SigErr'),msg); end% Find the Complex Sinusoid Frequenciesw_i = compute_freqs(md.noise_eigenvects,md.p_eff,md.EVFlag,md.eigenvals);% Estimate the noise variance as the average of the noise subspace eigenvaluessigma_w = sum(md.eigenvals(md.p_eff+1:end))./size(md.noise_eigenvects,2);% Estimate the power of the sinusoids[powers] = compute_power(md.signal_eigenvects,md.eigenvals,w_i,md.p_eff,sigma_w,xIsReal);% Compute MSE[w_mse,p_mse] = compute_mse(sigma_w,powers,length(x));% Convert the estimated frequencies to Hz if Fs was specifiedif ~isempty(md.Fs), w_i = w_i*md.Fs./(2*pi); w_mse = w_mse*(md.Fs./(2*pi))^2;end%---------------------------------------------------------------------------------------------function w_i = compute_freqs(noise_eigenvects,p_eff,EVFlag,eigenvals)%Compute the frequencies via the roots of the polynomial formed with the noise eigenvectors%% Inputs:%% noise_eigenvects - a matrix whose columns are the noise subspace eigenvectors% p_eff - signal subspace dimension% EVFlag - a flag indicating of the eigenvector methos should be used% eigenvals - a vector with all the correlation matrix eigenvalues. % However, we use only the noise eigenvalues as weights % in the eigenvector method.%% Outputs:%% w_i - frequencies of the complex sinusoids% compute weightsif EVFlag, % Eigenvector method, use eigenvalues as weights weights = eigenvals(end-size(noise_eigenvects,2)+1:end); % Use the noise subspace eigenvalueselse weights = ones(1,size(noise_eigenvects,2));end% Form a polynomial D, consisting of a sum of polynomials given by the product of% the noise subspace eigenvectors and the reversed and conjugated version.D = 0;for i = 1:length(weights), D = D + conv(noise_eigenvects(:,i),conj(flipud(noise_eigenvects(:,i))))./weights(i);endroots_D = roots(D);% Because D is formed from the product of a polynomial and its conjugated and reversed version,% every root of D inside the unit circle, will have a "reflected" version outside the unit circle.% We choose to use the ones inside the unit circle, because the distance from them to the unit% circle will be smaller than the corresponding distance for the "reflected" root.roots_D1 = roots_D(abs(roots_D) < 1);% Sort the roots from closest to furthest from the unit circle[not_used,indx] = sort(abs(abs(roots_D1)-1)); %#oksorted_roots = roots_D1(indx);% Use the first p_eff roots to determine the frequenciesw_i = angle(sorted_roots(1:p_eff));%-----------------------------------------------------------------------------------------------function [powers] = compute_power(signal_eigenvects,eigenvals,w_i,p_eff,sigma_w,xIsReal)%COMPUTE_POWER Solves the system of linear eqs. to calculate the power of the sinusoids.%% Inputs:%% signal_eigenvects - the matrix whose columns are the signal subspace eigenvectors% eigenvals - a vector containing all eigenvalues of the correlation matrix% w_i - a vector of frequency estimates of the sinusoids% p_eff - the dimension of the signal subspace% sigma_w - the estimate of the variance of the white noise% xIsReal - a flag indicating wether we have real or complex sinusoids%% Outputs:%% powers - a vector that contains the power of each sinusoid%This is just the solution of a linear system of eqs, Ax=b% For real sinusoids, the system of eqs. has half the number of unknownsif xIsReal, w_i = reshape(w_i,2,length(w_i)./2); w_i = w_i(1,:); % Use only the positive freqs. w_i = w_i(:); p_eff = p_eff./2;end% Form the A matrixif length(w_i) == 1, % FREQZ does not compute the gain at a single frequency, handle this separately A = polyval(signal_eigenvects(:,1),exp(1i*w_i));else for n = 1:p_eff, A(:,n) = freqz(signal_eigenvects(:,n),1,w_i); endendA = abs(A.').^2;% Form the b vectorb = eigenvals(1:p_eff) - sigma_w;% The powers are simply the solution to the set of eqs.powers = A\b;%--------------------------------------------------------------------------function [w_mse,p_mse] = compute_mse(sigma_w,powers,N)% implements eq.30 in Reference L = 1; % one element array p_mse = 12 * (sigma_w/(powers*N*L^2)); % first term of eq.30 in paper is to pass from frequency to DOA % this sigma_w^2 could be wrong w_mse = 12/(2*L)* (sigma_w^2/(powers*N*L^2)); % [EOF] rootmusic.m