pmusic (Signal Processing Toolbox)

Power spectrum estimate using MUSIC eigenvector method.

Syntax

[Pxx,f] = pmusic(x,p)
[Pxx,f] = pmusic(x,[p thresh])
[Pxx,f] = pmusic(x,[p thresh],nfft,Fs,window,noverlap)
[Pxx,f] = pmusic(x,...,'corr')
[Pxx,f] = pmusic(x,...,'ev')
[Pxx,f,evects,svals] = pmusic(x,...)

Description

pmusic estimates the power spectral density (PSD) of a signal or correlation matrix using Schmidt's eigen-analysis method [1]. The name MUSIC is an acronym for MUltiple SIgnal Classification. The eigenvector method, which uses eigenvalue weighting, is also supported [2]. The calling syntax is similar to that of pwelch, which also performs spectrum estimation. pwelch uses the classical FFT-based approach while pmusic performs eigen-analysis of the signal's correlation matrix.

[Pxx,f] = pmusic(x,p)

and

[Pxx,f] = pmusic(x,[p thresh])

return Pxx, the power spectrum estimate, and f, a vector of frequencies at which the PSD is estimated. x is the input signal, where:

A row or column vector represents one observation of the process output (for example, one "signal")
A rectangular (possibly square) array assumes that each column of x is a separate observation of the process output (for example, each column is one output of an array of sensors, as in array processing)
A square matrix, given the trailing argument 'corr', represents a correlation matrix

The second argument is a one- or two-element vector, either p or [p thresh]. If only p is specified, the signal subspace dimension is p. If [p thresh] is specified, thresh is multiplied by

_min, the smallest eigenvalue; eigenvalues below the threshold

_min*thresh are assigned to the noise subspace. In this case, p is the maximum dimension of the signal subspace.

NOTE
pmusic must assign eigenvectors to the noise and signal subspaces, but this is very difficult to do in practice. The two parameters p and thresh are provided for flexibility and control.

[Pxx,f] = pmusic(x,[p thresh],nfft,Fs,window,noverlap)

specifies the FFT length nfft (default is 256) and the sampling frequency for the signal Fs (default is 1). If Fs is specified, the output frequency vector f is scaled by this value. If the input signal is real-valued, the frequency range is 0 to Fs/2; for the complex case, it is 0 to Fs. window is a scalar specifying the rectangular window length, or a vector giving the actual window coefficients. noverlap, used in conjunction with window, is a scalar that gives the number of points by which to overlap successive windows.

[Pxx,f] = pmusic(x,...,'corr')

forces x to be taken as a correlation matrix. In this case, the arguments window and noverlap are ignored.

[Pxx,f] = pmusic(x,...,'ev')

selects the eigenvector variant of the MUSIC estimator. See the "Algorithm" section below for an explanation of how this is different from the MUSIC method.

[Pxx,f,evects,svals] = pmusic(x,...)

returns two additional arguments. evects is a matrix of eigenvectors spanning the noise subspace (one per column). svals is either a vector of singular values (squared) from svd or a vector of eigenvalues of the correlation matrix when the 'corr' option is present.

Remarks

The input x can be a vector or a matrix. x can be interpreted as signal data or as a correlation matrix, in one of three ways:

x is a vector of signal values (row or column). In this case, the dimension of the eigenvectors must be given. This is done either by taking the default value of 2*p or by specifying a window length using window.
x is a rectangular (m-by-n, possibly square) matrix. In this case, each column of x is a separate observation signal that enters into the SVD analysis, n is the number of observations, and the dimension of the eigenvectors is equal to m, the length of a column.
x is a square matrix and the trailing 'corr' is present. x is treated as a correlation matrix. In this case, the matrix must have only real, nonnegative eigenvalues.

The inputs p and thresh can determine the number of noise eigenvectors in one of four ways:

If thresh < 1, or if it is unspecified, the number of eigenvectors spanning the signal subspace will be equal to p. p must be an integer satisfying 0 p < n, where n is the dimension of the eigenvectors. This dimension n is the column length in the data matrix case, the matrix size in the correlation matrix case, or the window length for signal vectors. The value of thresh is unused.
If p n, thresh must be at least 1. thresh is used as the multiplier to determine an absolute threshold for splitting the eigenvalues between the signal and noise subspaces:
If thresh < 1, there will be no noise eigenvectors. This case is not allowed and gives the following error message:
When p < n and thresh 1, p specifies the maximum number of signal eigenvectors. However, the threshold test specified by thresh can also take eigenvectors from the signal subspace and assign them to the noise subspace.

Examples

This example analyzes a signal vector xx, assuming that two real signals are present in the signal subspace. In this case, the dimension of the signal subspace is 4 because each real sinusoid is the sum of two complex exponentials:

nn = 0:199;
xx = cos(0.257*pi*nn) + sin(0.2*pi*nn) + 0.01*randn(size(nn));
[PP,ff] = pmusic(xx,4);

This example analyzes the same signal vector xx with an eigenvalue cutoff of 10% above the minimum. Setting p = Inf forces the signal/noise subspace decision to be based on thresh. Use eigenvectors of dimension 7 and a sampling frequency Fs of 8 kHz:

[PP,ff] = pmusic(xx,[Inf,1.1],[],8000,7); % window length = 7

With the third and fourth outputs, by plotting the zeros of the noise-eigenvector polynomials, it is possible to create a "Root-MUSIC" algorithm, as the following zplane plot illustrates:

[PP,ff,v_noise] = pmusic(xx,4);
for kk = 1:size(v_noise,2)
   rr(:,kk) = roots(v_noise(:,kk));
end
zplane(rr)

Assume that RR is a square correlation matrix (for example, 7-by-7):

RR = toeplitz(cos(0.1*pi*[0:6])) + 0.1*eye(7);
[PP,ff] = pmusic(RR,4,'corr');

Make an observation matrix xx that is rectangular (100-by-7):

xx = reshape(cos(0.257*pi*(0:699)),7,100) + 0.1*randn(7,100);
[PP,ff] = pmusic(xx,4);

Use the same signal, but let pmusic form the 100-by-7 data matrix using its window and overlap inputs. In addition, use a longer FFT:

yy = xx(:);
[PP,ff] = pmusic(yy,4,512,[],7,0);

If we set p = 0, all the eigenvectors are assigned to the noise subspace. 'ev' specifies the eigenvector weighting. This turns out to be equivalent to MVDL (Capon's MLM):

[PP,ff] = pmusic(RR,0,'ev','corr');

Algorithm

The MUSIC estimate is given by the formula

where N is the dimension of the eigenvectors and v_kis the k-th eigenvector of the correlation matrix of the input signal. The integer p is the dimension of the signal subspace, so the eigenvectors v_k used in the sum correspond to the smallest eigenvalues and also span the noise subspace. The vector e(f) consists of complex exponentials, so the inner product

amounts to a Fourier transform. The second form is preferred for computation because the FFT is computed for each v_k and then the squared magnitudes are summed.

In the eigenvector method, the summation is weighted by the eigenvalues

_k of the correlation matrix:

The function relies on the svd matrix decomposition in the signal case, and it uses the eig function for analyzing the correlation matrix. If SVD is used, the correlation matrix is never explicitly computed, but the singular values are the

_k.

Diagnostics

There must be at least one output argument and at least two inputs; otherwise, pmusic stops and gives one of the following error messages:

Must have at least 1 output argument.
Must have at least 2 input arguments.

The first argument must be a full matrix, otherwise pmusic gives the following error message:

Input signal or correlation cannot be sparse.

If the second argument was entered as an empty matrix, or if it has more than two elements, or if it has negative or non-integer elements, pmusic gives one of the following error messages:

P cannot be empty.
Second input must have only 1 or 2 elements.
P must be an integer.
Second input must contain non-negative entries.

If the value of p is too large with respect to the eigenvector dimension, and thresh is less than 1, no eigenvectors can be assigned to the noise subspace and the algorithm fails. In this case, pmusic gives the following error message:

Noise subspace dimension cannot be zero.

If the 'corr' parameter is used, then the first input must be a square correlation matrix. If it is not, pmusic gives the following error message:

Correlation matrix (R) is not square.

The correlation matrix is then checked for validity; if it fails, pmusic gives the following error message:

Correlation matrix (R) has negative or complex eigenvalue.

`lpc`	Linear prediction coefficients.
`pburg`	Power spectrum estimate using the Burg method.
`pcov`	Power spectrum estimate using the covariance method.
`pmcov`	Power spectrum estimate using the modified covariance method.
`pmtm`	Power spectrum estimate using the multitaper method (MTM).
`prony`	Prony's method for time domain IIR filter design.
`pwelch`	Estimate the power spectral density (PSD) of a signal using Welch's method.
`pyulear`	Power spectrum estimate using Yule-Walker AR method.