--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/music.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,262 @@
+function [music_data,msg] = music(x,p,varargin)
+%MUSIC  Implements the heart of the MUSIC algorithm of line spectra estimation.
+%   MUSIC is called by both PMUSIC and ROOTMUSIC.
+%   
+%   Inputs:
+%     
+%     x - vector or matrix. If vector it is a signal, if matrix it may be either a data
+%         matrix such that x'*x=R, or a correlation matrix R.
+%     p - scalar or two element vector. If scalar, it indicates the dimension of the 
+%         signal subspace. If vector, p(2) is a threshold used to determine the 
+%         aforementioned dimension.
+%     nfft - (optional) to be used only with PMUSIC. A scalar indicating the number of
+%            points used in the evaluation of the pseudospectrum.
+%     Fs - (optional) a scalar specifying the sampling frequency. If omitted, we work
+%          in rad/sample; if empty it defaults to 1 Hz.
+%     nw - (optional) a scalar or vector indicating either the order of the correlation
+%          matrix or (when a vector) a window whose length is the order of the matrix
+%          and whose values are used to window each column of the data matrix.
+%     noverlap - (optional) a integer indicating the number of samples to overlap from
+%                column to column.
+%     strings - Optional input strings are: 'corr', 'EV' and range ('half' or 'whole').
+%
+%   Outputs:
+%
+%     msg - a possible error message.
+%
+%     music_data - a structure with the following fields:
+%          
+%          noise_eigenvects - a matrix whose columns are the noise subspace eigenvectors.
+%          signal_eigenvects - a matrix whose columns are the signal subspace eigenvectors.
+%          eigenvals - the eigenvalues of the correlation matrix.
+%          p_eff - the effective dimension of the signal subspace.
+%          nfft - number of points used to evaluate the pseudospectrum (only used in PMUSIC).
+%          Fs - sampling freq.
+%          range - string indicating whether 'half' or the 'whole' pseudospectrum should be
+%                  computed. (Only used in PMUSIC.)
+%          EVFlag - flag, 0 = MUSIC method; 1 = EigenVector method.
+
+%   Author(s): R. Losada
+%   Copyright 1988-2006 The MathWorks, Inc.
+%   $Revision: 1.1 $  $Date: 2010/02/18 11:16:20 $ 
+
+%   References:
+%     [1] Petre Stoica and Randolph Moses, Introduction To Spectral
+%         Analysis, Prentice-Hall, 1997, pg. 15
+%     [2] S. J. Orfanidis, Optimum Signal Processing. An Introduction. 
+%         2nd Ed., Macmillan, 1988.
+
+
+xIsReal = isreal(x);
+msg = '';
+music_data = [];
+
+if isempty(p),
+   msg = 'The signal subspace dimension cannot be empty.';
+   return
+end
+
+[opts,msg] = music_options(xIsReal,p,varargin{:});
+if ~isempty(msg),
+   return
+end
+
+% Compute the eigenvalues and eigenvectors of the correlation matrix
+[eigenvals,eigenvects] = computeeig(x,opts.CorrFlag,opts.CorrMatrOrd,opts.nw,opts.noverlap,opts.window,opts.EVFlag);
+
+% Determine the effective dimension of the signal subspace
+p_eff = determine_signal_space(p,eigenvals);
+
+% Separate the signal and noise eigenvectors
+signal_eigenvects = eigenvects(:,1:p_eff);
+noise_eigenvects = eigenvects(:,p_eff+1:end);
+
+% Generate the output structure
+music_data.noise_eigenvects = noise_eigenvects;
+music_data.signal_eigenvects = signal_eigenvects;
+music_data.eigenvals = eigenvals;
+music_data.p_eff = p_eff;
+music_data.nfft = opts.nfft;
+music_data.Fs = opts.Fs;
+music_data.EVFlag = opts.EVFlag;
+music_data.range = opts.range;
+
+
+%--------------------------------------------------------------------------------------
+function [options,msg] = music_options(xIsReal,p,varargin)
+%MUSIC_OPTIONS   Parse the optional inputs to the MUSIC function.
+%   MUSIC_OPTIONS returns a structure, OPTIONS, with the following fields:
+%
+%   options.nfft         - number of freq. points at which the psd is estimated
+%   options.Fs           - sampling freq. if any
+%   options.range        - 'onesided' or 'twosided' pseudospectrum (they correspond to
+%                          'half' and 'whole' respectively, but are returned as is by
+%                          psdoptions.m
+%   options.nw           - number of columns in the data matrix 
+%   options.noverlap     - number of samples to overlap
+%   options.window       - a vector with window coefficients
+%   options.CorrFlag     - a flag indicating whether the input is a correlation matrix
+%   options.EVFlag       - flag, 0 = MUSIC method ; 1 = EigenVector method
+%   options.CorrMatrOrd  - order of the correlation matrix to be used in computations
+
+
+
+% Assign Defaults
+msg = '';
+options.nw = [];
+options.noverlap = [];
+options.window = [];
+options.nfft = 256;
+options.Fs = [];
+options.CorrFlag = 0;
+options.EVFlag = 0;
+% Determine if frequency vector specified
+freqVecSpec = false;
+if (length(varargin) > 0 && isnumeric(varargin{1}) && length(varargin{1}) > 1)
+    freqVecSpec = true;
+end    
+
+if xIsReal && ~freqVecSpec,
+   options.range = 'onesided';
+else
+   options.range = 'twosided';
+end
+
+[options,msg] = psdoptions(xIsReal,options,varargin{:});
+
+if length(options.nfft) > 1,
+    if strcmpi(options.range,'onesided')
+        warning(generatemsgid('InconsistentRangeOption'),...
+            'Ignoring the ''onesided'' option. When a frequency vector is specified, a ''twosided'' PSD is computed');
+    end
+    options.range = 'twosided';
+end
+
+% psdoptions doesn't handle this field, assign it separetely
+options.CorrMatrOrd = 2*p(1);
+
+%-----------------------------------------------------------------------------------------
+function [eigenvals,eigenvects] = computeeig(x,CorrFlag,CorrMatrOrd,nw,noverlap,window,EVFlag)
+%COMPUTEEIG  Compute eigenvalues and eigenvectors of correlation matrix.
+%
+%   Inputs:
+%      
+%     x        - input vector or matrix
+%     CorrFlag - (flag) indicates whether x is a correlation matrix
+%     nw       - (integer) length of the rows of the data matrix 
+%                (only used if x is vector)
+%     noverlap - (integer) overlap between the rows of the data matrix
+%                (used in conjunction with nw) 
+%     window   - (vector) window to be applied to each column of data
+%                 matrix  (not used if x is a correlation matrix)
+%     EVFlag   -  True if eigenvector method, false if MUSIC.
+%   
+%
+%   Outputs:
+%    
+%     eigenvals
+%     eigenvects
+%
+%     If x is a matrix,  
+%          If CorrFlag = 1, input x is a correlation matrix, we compute the
+%          eigendecomposition and order the eigenvalues and eigenvectors.
+%
+%     If x is a vector,
+%          a data matrix is formed by calling corrmtx unless a custom nw
+%          and noverlap are specified. In that case, we use buffer to form
+%          the data matrix. 
+%
+%     If window is not empty, each row of the data matrix will be
+%     multiplied by the window.
+
+% Determine if the input is a matrix
+xIsMatrix = ~any(size(x)==1);
+
+if xIsMatrix && CorrFlag,
+    % Input is Correlation matrix
+
+    % Compute the eigenvectors and eigenvalues
+    [E,D] = eig((x+x')/2); % Ensure hermitian
+    [eigenvals,indx] = sort(diag(D),'descend');
+    eigenvects = E(:,indx);
+
+else
+    if xIsMatrix
+        % Input is already a data matrix
+        [Mx,Nx] = size(x); % Determine size of data matrix
+        if EVFlag && (Nx > Mx),
+            errmsg = 'The number of columns in the data matrix cannot exceed the number of rows.';
+            error(generatemsgid('invalidDataMatrix'),errmsg);
+        end
+
+    else
+        % x is a vector
+        x = x(:); % Make it a column
+        if isempty(nw),
+            x = corrmtx(x,CorrMatrOrd-1,'cov');
+        else
+            if EVFlag && nw > (ceil((length(x)-nw)/(nw-noverlap))+1),
+                errmsg = sprintf(['The segment length and overlap specified result in\n',...
+                    'a data matrix with more columns than rows.']);
+                error(generatemsgid('invalidDataMatrix'),errmsg);
+            end
+            Lx = length(x);
+            x = buffer(x,nw,noverlap,'nodelay');
+            if Lx <= nw,
+                error(generatemsgid('invalidSegmentLength'),'The segment length, NW, must be smaller than the signal length.');
+            end
+            x = x'./sqrt(Lx-nw); % Scale appropriately such that X'*X is a scaled estimate of R
+        end
+    end
+    if ~isempty(window),
+        % Apply window to each row of data matrix
+        if length(window) ~= size(x,2),
+            error(generatemsgid('InvalidDimensions'),'Window length must equal the number of columns in the data matrix.');
+        end
+        window = repmat(window(:).',size(x,1),1);
+        x = x.*window;
+    end
+
+    % Compute the eigenvectors and eigenvalues via the SVD
+    [U,S,eigenvects] = svd(x,0);
+    eigenvals = diag(S).^2; % We need to square the singular values here
+end
+
+
+%--------------------------------------------------------------------------------------------
+function p_eff = determine_signal_space(p,eigenvals)
+%DETERMINE_SIGNAL_SPACE   Determines the effective dimension of the signal subspace.
+%   
+%   Inputs:
+%
+%     p         - (scalar or vector) signal subspace dimension 
+%                 (but may contain a desired threshold).
+%     eigenvals - (vector) contains the eigenvalues (sorted in decreasing order)
+%                 of the correlation matrix
+%
+%   Outputs:
+%
+%     p_eff - The effective dimension of the signal subspace. If a threshold
+%             is given as p(2), the signal subspace will be equal to the number
+%             of eigenvalues, NEIG, greater than the threshold times the smallest
+%             eigenvalue. However, the dimension of the signal subspace is at most
+%             p(1), so that if NEIG is greater than p(1), p_eff will be equal to
+%             p(1). If the threshold criteria results in an empty signal subspace,
+%             once again we make p_eff = p(1).
+
+
+% Use the signal space dimension or the threshold to separate the noise subspace eigenvectors
+if length(p) == 2,
+   % The threshold will be the input threshold times the smallest eigenvalue
+   thresh = p(2)*eigenvals(end); 
+   indx = find(eigenvals > thresh);
+   if ~isempty(indx)
+      p_eff = min( p(1), length(indx) );
+   else
+      p_eff = p(1);
+   end
+else
+   p_eff = p;
+end
+
+% [EOF] - music.m