MUSIC (algorithm)

Last updated January 22, 2024

MUSIC (MUltiple SIgnal Classification) is an algorithm used for frequency estimation ^[1]^[2]^[3] and radio direction finding.^[4]

History

In many practical signal processing problems, the objective is to estimate from measurements a set of constant parameters upon which the received signals depend. There have been several approaches to such problems including the so-called maximum likelihood (ML) method of Capon (1969) and Burg's maximum entropy (ME) method. Although often successful and widely used, these methods have certain fundamental limitations (especially bias and sensitivity in parameter estimates), largely because they use an incorrect model (e.g., AR rather than special ARMA) of the measurements.

Pisarenko (1973) was one of the first to exploit the structure of the data model, doing so in the context of estimation of parameters of complex sinusoids in additive noise using a covariance approach. Schmidt (1977), while working at Northrop Grumman and independently Bienvenu and Kopp (1979) were the first to correctly exploit the measurement model in the case of sensor arrays of arbitrary form. Schmidt, in particular, accomplished this by first deriving a complete geometric solution in the absence of noise, then cleverly extending the geometric concepts to obtain a reasonable approximate solution in the presence of noise. The resulting algorithm was called MUSIC (MUltiple SIgnal Classification) and has been widely studied.

In a detailed evaluation based on thousands of simulations, the Massachusetts Institute of Technology's Lincoln Laboratory concluded in 1998 that, among currently accepted high-resolution algorithms, MUSIC was the most promising and a leading candidate for further study and actual hardware implementation.^[5] However, although the performance advantages of MUSIC are substantial, they are achieved at a cost in computation (searching over parameter space) and storage (of array calibration data).^[6]

Theory

MUSIC method assumes that a signal vector, $\mathbf {x}$ , consists of $p$ complex exponentials, whose frequencies $\omega$ are unknown, in the presence of Gaussian white noise, $\mathbf {n}$ , as given by the linear model

\mathbf {x} =\mathbf {A} \mathbf {s} +\mathbf {n} .

Here $\mathbf {A} =[\mathbf {a} (\omega _{1}),\cdots ,\mathbf {a} (\omega _{p})]$ is an $M\times p$ Vandermonde matrix of steering vectors $\mathbf {a} (\omega )=[1,e^{j\omega },e^{j2\omega },\ldots ,e^{j(M-1)\omega }]^{T}$ and $\mathbf {s} =[s_{1},\ldots ,s_{p}]^{T}$ is the amplitude vector. A crucial assumption is that number of sources, $p$ , is less than the number of elements in the measurement vector, $M$ , i.e. $p<M$ .

The $M\times M$ autocorrelation matrix of $\mathbf {x}$ is then given by

\mathbf {R} _{x}=\mathbf {A} \mathbf {R} _{s}\mathbf {A} ^{H}+\sigma ^{2}\mathbf {I} ,

where $\sigma ^{2}$ is the noise variance, $\mathbf {I}$ is $M\times M$ identity matrix, and $\mathbf {R} _{s}$ is the $p\times p$ autocorrelation matrix of $\mathbf {s}$ .

The autocorrelation matrix $\mathbf {R} _{x}$ is traditionally estimated using sample correlation matrix

{\widehat {\mathbf {R} }}_{x}={\frac {1}{N}}\mathbf {X} \mathbf {X} ^{H}

where $N>M$ is the number of vector observations and $\mathbf {X} =[\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{N}]$ . Given the estimate of $\mathbf {R} _{x}$ , MUSIC estimates the frequency content of the signal or autocorrelation matrix using an eigenspace method.

Since $\mathbf {R} _{x}$ is a Hermitian matrix, all of its $M$ eigenvectors $\{\mathbf {v} _{1},\mathbf {v} _{2},\ldots ,\mathbf {v} _{M}\}$ are orthogonal to each other. If the eigenvalues of $\mathbf {R} _{x}$ are sorted in decreasing order, the eigenvectors $\{\mathbf {v} _{1},\ldots ,\mathbf {v} _{p}\}$ corresponding to the $p$ largest eigenvalues (i.e. directions of largest variability) span the signal subspace ${\mathcal {U}}_{S}$ . The remaining $M-p$ eigenvectors correspond to eigenvalue equal to $\sigma ^{2}$ and span the noise subspace ${\mathcal {U}}_{N}$ , which is orthogonal to the signal subspace, ${\mathcal {U}}_{S}\perp {\mathcal {U}}_{N}$ .

Note that for $M=p+1$ , MUSIC is identical to Pisarenko harmonic decomposition. The general idea behind MUSIC method is to use all the eigenvectors that span the noise subspace to improve the performance of the Pisarenko estimator.

Since any signal vector $\mathbf {e}$ that resides in the signal subspace $\mathbf {e} \in {\mathcal {U}}_{S}$ must be orthogonal to the noise subspace, $\mathbf {e} \perp {\mathcal {U}}_{N}$ , it must be that $\mathbf {e} \perp \mathbf {v} _{i}$ for all the eigenvectors $\{\mathbf {v} _{i}\}_{i=p+1}^{M}$ that spans the noise subspace. In order to measure the degree of orthogonality of $\mathbf {e}$ with respect to all the $\mathbf {v} _{i}\in {\mathcal {U}}_{N}$ , the MUSIC algorithm defines a squared norm

d^{2}=\|\mathbf {U} _{N}^{H}\mathbf {e} \|^{2}=\mathbf {e} ^{H}\mathbf {U} _{N}\mathbf {U} _{N}^{H}\mathbf {e} =\sum _{i=p+1}^{M}|\mathbf {e} ^{H}\mathbf {v} _{i}|^{2}

where the matrix $\mathbf {U} _{N}=[\mathbf {v} _{p+1},\ldots ,\mathbf {v} _{M}]$ is the matrix of eigenvectors that span the noise subspace ${\mathcal {U}}_{N}$ . If $\mathbf {e} \in {\mathcal {U}}_{S}$ , then $d^{2}=0$ as implied by the orthogonality condition. Taking a reciprocal of the squared norm expression creates sharp peaks at the signal frequencies. The frequency estimation function for MUSIC (or the pseudo-spectrum) is

{\hat {P}}_{MU}(e^{j\omega })={\frac {1}{\mathbf {e} ^{H}\mathbf {U} _{N}\mathbf {U} _{N}^{H}\mathbf {e} }}={\frac {1}{\sum _{i=p+1}^{M}|\mathbf {e} ^{H}\mathbf {v} _{i}|^{2}}},

where $\mathbf {v} _{i}$ are the noise eigenvectors and

\mathbf {e} ={\begin{bmatrix}1&e^{j\omega }&e^{j2\omega }&\cdots &e^{j(M-1)\omega }\end{bmatrix}}^{T}

is the candidate steering vector. The locations of the $p$ largest peaks of the estimation function give the frequency estimates for the $p$ signal components

{\hat {\omega }}=\arg \max _{\omega }\;{\hat {P}}_{MU}(e^{j\omega }).

MUSIC is a generalization of Pisarenko's method, and it reduces to Pisarenko's method when $M=p+1$ . In Pisarenko's method, only a single eigenvector is used to form the denominator of the frequency estimation function; and the eigenvector is interpreted as a set of autoregressive coefficients, whose zeros can be found analytically or with polynomial root finding algorithms. In contrast, MUSIC assumes that several such functions have been added together, so zeros may not be present. Instead there are local minima, which can be located by computationally searching the estimation function for peaks.

Dimension of signal space

The fundamental observation MUSIC and other subspace decomposition methods are based on is about the rank of the autocorrelation matrix $\mathbf {R} _{x}$ which is related to number of signal sources $p$ as follows.

If the sources are complex, then $M>p$ and the dimension of the signal subspace ${\mathcal {U}}_{S}$ is $p$ . If sources are real, then $M>2p$ and dimension of the signal subspace is $2p$ , i.e. each real sinusoid is generated by two base vectors.

This fundamental result, although often skipped in spectral analysis books, is a reason why the input signal can be distributed into $p$ signal subspace eigenvectors spanning ${\mathcal {U}}_{S}$ ( $2p$ for real valued signals) and noise subspace eigenvectors spanning ${\mathcal {U}}_{N}$ . It is based on signal embedding theory ^[2]^[7] and can also be explained by the topological theory of manifolds.^[4]

Comparison to other methods

MUSIC outperforms simple methods such as picking peaks of DFT spectra in the presence of noise, when the number of components is known in advance, because it exploits knowledge of this number to ignore the noise in its final report.

Unlike DFT, it is able to estimate frequencies with accuracy higher than one sample, because its estimation function can be evaluated for any frequency, not just those of DFT bins. This is a form of superresolution.

Its chief disadvantage is that it requires the number of components to be known in advance, so the original method cannot be used in more general cases. Methods exist for estimating the number of source components purely from statistical properties of the autocorrelation matrix. See, e.g. ^[8] In addition, MUSIC assumes coexistent sources to be uncorrelated, which limits its practical applications.

Recent iterative semi-parametric methods offer robust superresolution despite highly correlated sources, e.g., SAMV ^[9]^[10]

Other applications

A modified version of MUSIC, denoted as Time-Reversal MUSIC (TR-MUSIC) has been recently applied to computational time-reversal imaging.^[11]^[12] MUSIC algorithm has also been implemented for fast detection of the DTMF frequencies (Dual-tone multi-frequency signaling) in the form of C library - libmusic^[13] (including for MATLAB implementation).^[14]

Related Research Articles

In mathematics, the discrete Fourier transform (DFT) converts a finite sequence of equally-spaced samples of a function into a same-length sequence of equally-spaced samples of the discrete-time Fourier transform (DTFT), which is a complex-valued function of frequency. The interval at which the DTFT is sampled is the reciprocal of the duration of the input sequence. An inverse DFT (IDFT) is a Fourier series, using the DTFT samples as coefficients of complex sinusoids at the corresponding DTFT frequencies. It has the same sample-values as the original input sequence. The DFT is therefore said to be a frequency domain representation of the original input sequence. If the original sequence spans all the non-zero values of a function, its DTFT is continuous, and the DFT provides discrete samples of one cycle. If the original sequence is one cycle of a periodic function, the DFT provides all the non-zero values of one DTFT cycle.

In computational mathematics, an iterative method is a mathematical procedure that uses an initial value to generate a sequence of improving approximate solutions for a class of problems, in which the n-th approximation is derived from the previous ones.

In engineering, a transfer function of a system, sub-system, or component is a mathematical function that models the system's output for each possible input. It is widely used in electronic engineering tools like circuit simulators and control systems. In simple cases, this function can be represented as a two-dimensional graph of an independent scalar input versus the dependent scalar output. Transfer functions for components are used to design and analyze systems assembled from components, particularly using the block diagram technique, in electronics and control theory.

H_∞methods are used in control theory to synthesize controllers to achieve stabilization with guaranteed performance. To use H_∞ methods, a control designer expresses the control problem as a mathematical optimization problem and then finds the controller that solves this optimization. H_∞ techniques have the advantage over classical control techniques in that H_∞ techniques are readily applicable to problems involving multivariate systems with cross-coupling between channels; disadvantages of H_∞ techniques include the level of mathematical understanding needed to apply them successfully and the need for a reasonably good model of the system to be controlled. It is important to keep in mind that the resulting controller is only optimal with respect to the prescribed cost function and does not necessarily represent the best controller in terms of the usual performance measures used to evaluate controllers such as settling time, energy expended, etc. Also, non-linear constraints such as saturation are generally not well-handled. These methods were introduced into control theory in the late 1970s-early 1980s by George Zames, J. William Helton , and Allen Tannenbaum.

In linear algebra, a circulant matrix is a square matrix in which all rows are composed of the same elements and each row is rotated one element to the right relative to the preceding row. It is a particular kind of Toeplitz matrix.

A sensor array is a group of sensors, usually deployed in a certain geometry pattern, used for collecting and processing electromagnetic or acoustic signals. The advantage of using a sensor array over using a single sensor lies in the fact that an array adds new dimensions to the observation, helping to estimate more parameters and improve the estimation performance. For example an array of radio antenna elements used for beamforming can increase antenna gain in the direction of the signal while decreasing the gain in other directions, i.e., increasing signal-to-noise ratio (SNR) by amplifying the signal coherently. Another example of sensor array application is to estimate the direction of arrival of impinging electromagnetic waves. The related processing method is called array signal processing. A third examples includes chemical sensor arrays, which utilize multiple chemical sensors for fingerprint detection in complex mixtures or sensing environments. Application examples of array signal processing include radar/sonar, wireless communications, seismology, machine condition monitoring, astronomical observations fault diagnosis, etc.

Array processing is a wide area of research in the field of signal processing that extends from the simplest form of 1 dimensional line arrays to 2 and 3 dimensional array geometries. Array structure can be defined as a set of sensors that are spatially separated, e.g. radio antenna and seismic arrays. The sensors used for a specific problem may vary widely, for example microphones, accelerometers and telescopes. However, many similarities exist, the most fundamental of which may be an assumption of wave propagation. Wave propagation means there is a systemic relationship between the signal received on spatially separated sensors. By creating a physical model of the wave propagation, or in machine learning applications a training data set, the relationships between the signals received on spatially separated sensors can be leveraged for many applications.

<span class="mw-page-title-main">Scoring rule</span> Measure for evaluating probabilistic forecasts

In decision theory, a scoring rule provides a summary measure for the evaluation of probabilistic predictions or forecasts. It is applicable to tasks in which predictions assign probabilities to events, i.e. one issues a probability distribution $as prediction. This includes probabilistic classification of a set of mutually exclusive outcomes or classes.$

Pisarenko harmonic decomposition, also referred to as Pisarenko's method, is a method of frequency estimation. This method assumes that a signal, $, consists of complex exponentials in the presence of white noise. Because the number of complex exponentials must be known a priori, it is somewhat limited in its usefulness.$

The theoretical and experimental justification for the Schrödinger equation motivates the discovery of the Schrödinger equation, the equation that describes the dynamics of nonrelativistic particles. The motivation uses photons, which are relativistic particles with dynamics described by Maxwell's equations, as an analogue for all types of particles.

In statistical signal processing, the goal of spectral density estimation (SDE) or simply spectral estimation is to estimate the spectral density of a signal from a sequence of time samples of the signal. Intuitively speaking, the spectral density characterizes the frequency content of the signal. One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

In linear algebra, eigendecomposition is the factorization of a matrix into a canonical form, whereby the matrix is represented in terms of its eigenvalues and eigenvectors. Only diagonalizable matrices can be factorized in this way. When the matrix being factorized is a normal or real symmetric matrix, the decomposition is called "spectral decomposition", derived from the spectral theorem.

In time series analysis, singular spectrum analysis (SSA) is a nonparametric spectral estimation method. It combines elements of classical time series analysis, multivariate statistics, multivariate geometry, dynamical systems and signal processing. Its roots lie in the classical Karhunen (1946)–Loève spectral decomposition of time series and random fields and in the Mañé (1981)–Takens (1981) embedding theorem. SSA can be an aid in the decomposition of time series into a sum of components, each having a meaningful interpretation. The name "singular spectrum analysis" relates to the spectrum of eigenvalues in a singular value decomposition of a covariance matrix, and not directly to a frequency domain decomposition.

Estimation theory, or estimation of signal parameters via rotational invariant techniques (ESPRIT) is a technique to determine parameters of a mixture of sinusoids in a background noise. This technique is first proposed for frequency estimation, however, with the introduction of phased-array systems in everyday technology, it is also used for angle of arrival estimations.

Common spatial pattern (CSP) is a mathematical procedure used in signal processing for separating a multivariate signal into additive subcomponents which have maximum differences in variance between two windows.

In machine learning, the kernel embedding of distributions comprises a class of nonparametric methods in which a probability distribution is represented as an element of a reproducing kernel Hilbert space (RKHS). A generalization of the individual data-point feature mapping done in classical kernel methods, the embedding of distributions into infinite-dimensional feature spaces can preserve all of the statistical features of arbitrary distributions, while allowing one to compare and manipulate distributions using Hilbert space operations such as inner products, distances, projections, linear transformations, and spectral analysis. This learning framework is very general and can be applied to distributions over any space $on which a sensible kernel function may be defined. For example, various kernels have been proposed for learning from data which are: vectors in, discrete classes/categories, strings, graphs/networks, images, time series, manifolds, dynamical systems, and other structured objects. The theory behind kernel embeddings of distributions has been primarily developed by Alex Smola, Le Song, Arthur Gretton, and Bernhard Schölkopf. A review of recent works on kernel embedding of distributions can be found in.$

Lightfieldmicroscopy (LFM) is a scanning-free 3-dimensional (3D) microscopic imaging method based on the theory of light field. This technique allows sub-second (~10 Hz) large volumetric imaging with ~1 μm spatial resolution in the condition of weak scattering and semi-transparence, which has never been achieved by other methods. Just as in traditional light field rendering, there are two steps for LFM imaging: light field capture and processing. In most setups, a microlens array is used to capture the light field. As for processing, it can be based on two kinds of representations of light propagation: the ray optics picture and the wave optics picture. The Stanford University Computer Graphics Laboratory published their first prototype LFM in 2006 and has been working on the cutting edge since then.

The streamline upwind Petrov–Galerkin pressure-stabilizing Petrov–Galerkin formulation for incompressible Navier–Stokes equations can be used for finite element computations of high Reynolds number incompressible flow using equal order of finite element space by introducing additional stabilization terms in the Navier–Stokes Galerkin formulation.

Generalized pencil-of-function method (GPOF), also known as matrix pencil method, is a signal processing technique for estimating a signal or extracting information with complex exponentials. Being similar to Prony and original pencil-of-function methods, it is generally preferred to those for its robustness and computational efficiency.

Tau functions are an important ingredient in the modern mathematical theory of integrable systems, and have numerous applications in a variety of other domains. They were originally introduced by Ryogo Hirota in his direct method approach to soliton equations, based on expressing them in an equivalent bilinear form.

References

↑ Hayes, Monson H., Statistical Digital Signal Processing and Modeling, John Wiley & Sons, Inc., 1996. ISBN 0-471-59431-8.
1 2 Gregor, Piotr (2022). Zastosowanie algorytmu MUSIC do wykrywania DTMF [Application of MUSIC algorithm to DTMF detection] (Thesis) (in Polish). Warsaw University of Technology.
↑ Costanzo, Sandra; Buonanno, Giovanni; Solimene, Raffaele (2022). "Super-Resolution Spectral Approach for the Accuracy Enhancement of Biomedical Resonant Microwave Sensors". IEEE Journal of Electromagnetics, RF and Microwaves in Medicine and Biology . 6 (4): 539–545. doi:10.1109/JERM.2022.3210457. ISSN 2469-7249. S2CID 252792474.
1 2 Schmidt, R.O, "Multiple Emitter Location and Signal Parameter Estimation," IEEE Trans. Antennas Propagation, Vol. AP-34 (March 1986), pp. 276–280.
↑ Barabell, A. J. (1998). "Performance Comparison of Superresolution Array Processing Algorithms. Revised" (PDF). Massachusetts Inst of Tech Lexington Lincoln Lab. Archived (PDF) from the original on May 25, 2021.
↑ R. Roy and T. Kailath, "ESPRIT-estimation of signal parameters via rotational invariance techniques," in IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, pp. 984–995, Jul 1989.
↑ Penny, W. D. (2009), Signal Processing Course, University College London, Lecture notes 1999–2000 academic year
↑ Fishler, Eran, and H. Vincent Poor. "Estimation of the number of sources in unbalanced arrays via information theoretic criteria." IEEE Transactions on Signal Processing 53.9 (2005): 3543–3553.
↑ Abeida, Habti; Zhang, Qilin; Li, Jian; Merabtine, Nadjim (2013). "Iterative Sparse Asymptotic Minimum Variance Based Approaches for Array Processing". IEEE Transactions on Signal Processing. Institute of Electrical and Electronics Engineers (IEEE). 61 (4): 933–944. arXiv: 1802.03070 . Bibcode:2013ITSP...61..933A. doi:10.1109/tsp.2012.2231676. ISSN 1053-587X. S2CID 16276001.
↑ Zhang, Qilin; Abeida, Habti; Xue, Ming; Rowe, William; Li, Jian (2012). "Fast implementation of sparse iterative covariance-based estimation for source localization". The Journal of the Acoustical Society of America. 131 (2): 1249–1259. Bibcode:2012ASAJ..131.1249Z. doi:10.1121/1.3672656. PMID 22352499.
↑ Devaney, A.J. (2005-05-01). "Time reversal imaging of obscured targets from multistatic data". IEEE Transactions on Antennas and Propagation. 53 (5): 1600–1610. Bibcode:2005ITAP...53.1600D. doi:10.1109/TAP.2005.846723. ISSN 0018-926X. S2CID 25241225.
↑ Ciuonzo, D.; Romano, G.; Solimene, R. (2015-05-01). "Performance Analysis of Time-Reversal MUSIC". IEEE Transactions on Signal Processing. 63 (10): 2650–2662. Bibcode:2015ITSP...63.2650C. doi:10.1109/TSP.2015.2417507. ISSN 1053-587X. S2CID 5895440.
↑ "libmusic: A powerful C library for spectral analysis". Data and Signal. 2023.
↑ "libmusic_m : MATLAB implementation". Data and Signal. 2023. MathWorks.