MUSIC - Maple Help

# Online Help

###### All Products    Maple    MapleSim

SignalProcessing

 MUSIC
 performs the Multiple Signal Classifier (MUSIC) Method on a signal

 Calling Sequence MUSIC( signal, options )

Parameters

 signal - 1-D rtable or list of data, which must contain two or more elements.

Options

 • compiled: (optional) Either true or false, determines if internal auxiliary procedures will be compiled. The default is false.
 • dimension: (optional) Positive integer for the dimension of the signal subspace, which must be less than the size of the signal. The default is NULL.
 • plotoptions: (optional) List containing any additional options to be passed to the plotting command when creating the pseudo-power spectrum. The default is [].
 • samplerate: (optional) Positive numeric value for the sampling rate. The default is 1.0.
 • size: (optional) The number of points to be used to construct the pseudo-power spectrum, which cannot be smaller than the size of the signal. The default is 4*numelems(signal).
 • threshold: (optional) Positive number no larger than 1.0 for the minimum value for a normalized eigenvalue to be considered part of the signal subspace. The default is NULL.
 • window: (optional) Either a list, name, or string, specifies the windowing command to be applied to the signal. The default is "none" (for no windowing to be applied). If a list is passed, the first element provides the name of the windowing command, and any remaining terms are passed as options to the command.
 • windownormalization (optional) Either true or false, indicates if the windowing function is to be normalized. The default is true.
 • output: (optional) The type of output. The supported options are:
 – eigenvalues: Returns a Vector, of float[8] datatype and length $m=\mathrm{numelems}\left(\mathrm{signal}\right)$, containing the normalized (with respect to the maximum) eigenvalues of the auto-correlation Matrix.
 – eigenvectors: Returns a Matrix, of complex[8] datatype and size $m$ by $m$, containing the eigenvectors of the auto-correlation Matrix.
 – frequencies: Returns a Vector, of float[8] datatype and length $n=\mathrm{size}$, containing the frequencies.
 – peaks: Returns a Matrix of float[8] datatype containing the detected signal frequencies, in the first column, and their heights in the pseudo-power spectrum, in the second column. The rows are sorted by decreasing heights.
 – plot: Returns a scaled plot which displays the pseudo-power spectrum $\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%.}\left(10,\mathrm{log10}\left(P\right)\right)\right]\right)$ versus the frequencies $F$.
 – power: Returns a Vector of float[8] datatype containing the pseudo-power spectrum.
 – record: Returns a record with the previous options.
 – list of any of the above options: Returns an expression sequence with the corresponding outputs, in the same order. The default is [frequencies,power].

Description

 • The Multiple Signal Classifier (MUSIC) Method is an algorithm used to estimate frequencies in a noisy signal. To this end, the data is separated into a signal subspace and a noise subspace, where the signal subspace is orthogonal to the uncorrelated noise subspace. Now, for any frequency $\mathrm{\omega }$ that is present in the signal, the sum of the projections of the steering Vector $S$ of size $m$, where $m$ is the size of the signal and

${S}_{i}={ⅇ}^{I\left(\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%/}\left(2\mathrm{Pi},r\right)\right]\right)\right)\left(i-1\right)\mathrm{\omega }}$

 for each $i$, where $r$ is the sample rate, will be close to zero. Thus, we can construct a pseudo-power spectrum at uniformly-spaced discrete frequencies, with denominator that is small near a signal frequency, and thus the locations of the peaks tell us the frequencies.
 • To estimate the frequencies:
 1 Start with the 1-D signal $X$ of size $m$.
 2 Compute the mean $\mathrm{\mu }$ of $X$.
 3 Compute the auto-correlation Vector $A$ of $X-\mathrm{\mu }$, which is of size $m$.
 4 Take the auto-correlation Matrix $C$ of $A$ to be the square Matrix of size $m$ with elements given by:

${C}_{i,j}=\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%.}\left(0.5,{A}_{\mathrm{modp}\left(j-i,m\right)+1}+\mathrm{conjugate}\left({A}_{\mathrm{modp}\left(i-j,m\right)+1}\right)\right)\right]\right)$

 Here, $C$ is Hermitian, positive definite, Toeplitz, and circulant and, consequently, the eigenvalues are real-valued and positive.
 5 Determine the sorted (in decreasing size) eigenvalues Vector $U$ and eigenvectors Matrix $V$, with column $j$ of $V$ being the eigenvector associated with eigenvalue ${U}_{j}$. The first $d$ eigenvalues and eigenvectors correspond to the signal subspace, where $d$ is specified apriori, and the remaining $m-d$ eigenvalues and eigenvectors correspond to the noise subspace.
 6 Define the frequencies Vector $F$ of size $n$, where $n$ is the length of the pseudo-power spectrum, element-wise as

${F}_{k}=\frac{\left(k-1\right)r}{n}$

 where $r$ is the sample rate.
 7 The pseudo-power spectrum is the Vector of size $n$ with elements given by the following:

${P}_{k}=\frac{1}{{\sum }_{j=d+1}^{m}{\left|{\sum }_{i=1}^{m}\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%.}\left(\mathrm{exp}\left(I\mathrm{%/}\left(2\mathrm{Pi},r\right)\left(i-1\right){F}_{k}\right),{V}_{i,j}\right)\right]\right)\right|}^{2}}$

 The pseudo-power spectrum is then normalized with respect to the maximum value.
 8 The peaks of the pseudo-power spectrum indicate, ideally, the frequencies present in the signal.
 • If the number of frequencies to be determined is known, it is recommended that the value of $d=\mathrm{dimension}$ be chosen to be double this number.
 • The unsorted eigenvalues of a circulant Matrix C (which must be square) can be computed as the unnormalized Fast Fourier Transform (that is, using the FFT command with normalization=none) of the first (or any) row of the Matrix. The normalized eigenvectors, in fact, depend only on the size of C and not the specific elements. Specifically, the Matrix V of normalized eigenvectors, with the columns being the eigenvectors associated with the eigenvalues computed using the FFT, are given by the following formula:

${V}_{i,j}=\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%.}\left(\frac{1}{\mathrm{sqrt}\left(m\right)},\mathrm{exp}\left(I\mathrm{%/}\left(2\mathrm{Pi},m\right)i\left(j-1\right)\right)\right)\right]\right)$

 where $m$ is the row and column dimension of C.
 • The aforementioned properties of circulant matrices can be used to speed up the computation of the pseudo-power spectrum. Let $U$ be the Vector of unsorted eigenvalues computed as the unnormalized FFT of the first row of the auto-correlation Matrix $C$, and let $V$ be the Matrix of eigenvectors, with elements defined above. Define the $m$ by $n$ Matrix $S$ by:

${S}_{i,k}={ⅇ}^{I\left(\mathrm{Typesetting}:-\mathrm{_Hold}\left(\left[\mathrm{%/}\left(2\mathrm{Pi},r\right)\right]\right)\right)\left(i-1\right){F}_{k}}$

 Now, the inner sum of the denominator of ${P}_{k}$ is the Inverse Fast Fourier Transform (IFFT) of column $k$ of $S$. If $T$ is the Matrix formed by first taking the IFFT of each column of $S$, and then setting to zero the rows that correspond to the $d$ largest eigenvalues of $U$, we can write the pseudo-power spectrum as the following:

${P}_{k}=\frac{1}{{\sum }_{j=1}^{m}{\left|{T}_{j,k}\right|}^{2}}$

 • The value of size determines the frequency resolution of the pseudo-power spectrum. Moreover, when the signal is real-valued, the maximum search frequency is samplerate/2, otherwise the maximum search frequency is samplerate.
 • If signal is of type AudioTools:-Audio, then the sample rate stored in the attributes of signal will override the sample rate passed in the samplerate parameter.
 • Internally, the dimension d of the signal subspace is determined as follows:
 – If dimension=NULL and threshold=NULL, then d=m-1, where m is the size of the signal.
 – If dimension<>NULL, then d=dimension, even when threshold<>NULL.
 – If dimension=NULL and threshold<>NULL, then d is the number of eigenvalues no smaller than threshold.
 • In the event that one the points of the pseudo-power spectrum is HFloat(infinity), that is an exact signal frequency was detected at one of the discrete frequencies, then the value is replaced with double the value of the highest finite value present.
 • The pseudo-power spectrum is not intended to match the heights of the actual power spectrum. Rather, what is important is the locations of the most prominent peaks.
 • When MUSIC is called in the session for the first time with compiled=true, the internal auxiliary routines which otherwise would be called with evalhf will instead be compiled. If MUSIC has already been called in the session with compiled=true, it will continue to use the compiled versions of the internal commands, even if called with compiled=false.
 • Use of the compiled=true option is recommended when the sizes of the signal and pseudo-power spectrum are large, or the MUSIC command is to be called numerous times. In these scenarios, the cost of the overhead for compilation is likely worth it.
 • Maple will attempt to coerce the provided signal to a 1-D Vector of either float[8] or complex[8] datatype, and an error will be thrown if this is not possible. For this reason, it is most efficient for the passed input to use this datatype.
 • The input signal cannot have an indexing function, and must use rectangular storage.
 • The MUSIC command is not thread safe.
 • As the underlying implementation of the SignalProcessing package is a module, it is also possible to use the form SignalProcessing:-MUSIC to access the command from the package. For more information, see Module Members.

Examples

 > $\mathrm{with}\left(\mathrm{SignalProcessing}\right):$
 > $\mathrm{with}\left(\mathrm{Statistics}\right):$

Real-valued signal

 • Consider the following real-valued signal without noise:
 > $m≔{2}^{8}:$
 > $\mathrm{fs}≔100.0:$
 > $a≔\left[5.0,3.0,-7.0\right]:$
 > $f≔\left[10.25,10.40,20.35\right]:$
 > $X≔\mathrm{Vector}\left(m,\left[\mathrm{seq}\left(\mathrm{add}\left(a\left[i\right]\mathrm{sin}\left(\frac{f\left[i\right]\cdot 2\mathrm{\pi }\left(k-1\right)}{\mathrm{fs}}\right),i=1..3\right),k=1..m\right)\right],'\mathrm{datatype}'='\mathrm{float}\left[8\right]'\right):$
 > $\mathrm{SignalPlot}\left(X,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{color}'='\mathrm{blue}','\mathrm{title}'="Pure Signal Plot",'\mathrm{font}'=\left['\mathrm{Verdana}',15\right],'\mathrm{size}'=\left[800,400\right]\right)$
 > $\mathrm{Periodogram}\left(X,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{color}'='\mathrm{blue}','\mathrm{title}'="Pure Signal Periodogram",'\mathrm{font}'=\left['\mathrm{Verdana}',15\right],'\mathrm{size}'=\left[800,400\right]\right)$
 • Now, construct a Vector of noise, and add it to the pure signal:
 > $N≔\mathrm{Vector}\left[\mathrm{column}\right]\left(\mathrm{Sample}\left(\mathrm{RandomVariable}\left(\mathrm{Normal}\left(0,1\right)\right),m\right)\right):$
 > $Y≔X+N:$
 > $\mathrm{SignalPlot}\left(Y,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{color}'='\mathrm{blue}','\mathrm{title}'="Noisy Signal Plot",'\mathrm{font}'=\left['\mathrm{Verdana}',15\right],'\mathrm{size}'=\left[800,400\right]\right)$
 > $\mathrm{Periodogram}\left(Y,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{color}'='\mathrm{blue}','\mathrm{title}'="Noisy Signal Periodogram",'\mathrm{font}'=\left['\mathrm{Verdana}',15\right],'\mathrm{size}'=\left[800,400\right]\right)$
 • The periodogram of the noisy signal makes it difficult to isolate the three frequencies of the original signal. With the MUSIC command, on the other hand, even in the present of noise, the frequencies are correctly detected:
 > $d≔6:$
 > $n≔{2}^{10}:$
 > $R≔\mathrm{MUSIC}\left(Y,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{dimension}'=d,'\mathrm{size}'=n,'\mathrm{output}'='\mathrm{record}'\right):$
 > $\mathrm{Peaks}≔R\left['\mathrm{peaks}'\right]$
 ${\mathrm{Peaks}}{≔}\begin{array}{c}\left[\begin{array}{cc}{10.5468750000000}& {1.}\\ {10.1562500000000}& {0.881269822689514}\\ {20.3125000000000}& {0.314567906139360}\\ {11.1328125000000}& {3.45849990198759}{×}{{10}}^{{-29}}\\ {9.57031250000000}& {3.45836305214661}{×}{{10}}^{{-29}}\\ {19.7753906250000}& {3.40254084261547}{×}{{10}}^{{-29}}\\ {20.8496093750000}& {3.40232715742633}{×}{{10}}^{{-29}}\\ {11.5234375000000}& {3.32818070232832}{×}{{10}}^{{-29}}\\ {9.17968750000000}& {3.32798685701327}{×}{{10}}^{{-29}}\\ {19.3359375000000}& {3.30069539811017}{×}{{10}}^{{-29}}\\ {⋮}& {⋮}\end{array}\right]\\ \hfill {\text{126 × 2 Matrix}}\end{array}$ (1)
 > $R\left['\mathrm{plot}'\right]$
 • We can also specify a threshold for the eigenvalues, as opposed to the dimension of the signal space:
 > $t≔0.10:$
 > $E,P≔\mathrm{MUSIC}\left(Y,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{threshold}'=t,'\mathrm{size}'=n,'\mathrm{output}'=\left['\mathrm{eigenvalues}','\mathrm{peaks}'\right]\right)$
 ${E}{,}{P}{≔}\begin{array}{c}\left[\begin{array}{c}{1.}\\ {1.00000000000000}\\ {0.592904430279541}\\ {0.592904430279541}\\ {0.232892172398253}\\ {0.232892172398253}\\ {0.0399675246048442}\\ {0.0399675246048442}\\ {0.0365223947082229}\\ {0.0365223947082229}\\ {⋮}\end{array}\right]\\ \hfill {\text{256 element Vector[column]}}\end{array}{,}\begin{array}{c}\left[\begin{array}{cc}{10.5468750000000}& {1.}\\ {10.1562500000000}& {0.881269822689514}\\ {20.3125000000000}& {0.314567906139360}\\ {11.1328125000000}& {3.45849990198759}{×}{{10}}^{{-29}}\\ {9.57031250000000}& {3.45836305214661}{×}{{10}}^{{-29}}\\ {19.7753906250000}& {3.40254084261547}{×}{{10}}^{{-29}}\\ {20.8496093750000}& {3.40232715742633}{×}{{10}}^{{-29}}\\ {11.5234375000000}& {3.32818070232832}{×}{{10}}^{{-29}}\\ {9.17968750000000}& {3.32798685701327}{×}{{10}}^{{-29}}\\ {19.3359375000000}& {3.30069539811017}{×}{{10}}^{{-29}}\\ {⋮}& {⋮}\end{array}\right]\\ \hfill {\text{126 × 2 Matrix}}\end{array}$ (2)

Complex-valued signal

 • Consider a complex-valued signal with noise:
 > $m≔{2}^{8}:$
 > $\mathrm{fs}≔\mathrm{evalhf}\left(\frac{m-1}{2\mathrm{\pi }}\right):$
 > $T≔\mathrm{Vector}\left(m,k↦\mathrm{evalhf}\left(\frac{k-1}{\mathrm{fs}}\right),'\mathrm{datatype}'='\mathrm{float}\left[8\right]'\right):$
 > $X≔\mathrm{Vector}\left(m,k↦5\cdot \mathrm{evalhf}\left(\mathrm{cos}\left(5.0\cdot T\left[k\right]\right)+3\cdot I\cdot \mathrm{sin}\left(15.0\cdot T\left[k\right]\right)\right),'\mathrm{datatype}'='\mathrm{complex}\left[8\right]'\right):$
 > $\mathrm{PowerSpectrum}\left(X,'\mathrm{variety}'='\mathrm{signal}','\mathrm{samplerate}'=\mathrm{fs},'\mathrm{output}'='\mathrm{periodogram}','\mathrm{powerscale}'='\mathrm{dB}','\mathrm{periodogramoptions}'=\left['\mathrm{color}'='\mathrm{blue}','\mathrm{title}'="Pure Signal Periodogram",'\mathrm{font}'=\left['\mathrm{Verdana}',15\right]\right]\right)$
 > $\mathrm{Nr}≔\mathrm{Vector}\left[\mathrm{column}\right]\left(\mathrm{Sample}\left(\mathrm{RandomVariable}\left(\mathrm{Normal}\left(0,1\right)\right),m\right)\right):$
 > $\mathrm{Ni}≔\mathrm{Vector}\left[\mathrm{column}\right]\left(\mathrm{Sample}\left(\mathrm{RandomVariable}\left(\mathrm{Normal}\left(0,1\right)\right),m\right)\right):$
 > $Y≔X+\mathrm{Nr}+I\mathrm{Ni}:$
 > $d≔4:$
 > $n≔{2}^{10}:$
 > $R≔\mathrm{MUSIC}\left(Y,'\mathrm{samplerate}'=\mathrm{fs},'\mathrm{dimension}'=d,'\mathrm{size}'=n,'\mathrm{output}'='\mathrm{record}'\right):$
 > $\mathrm{Peaks}≔R\left['\mathrm{peaks}'\right]$
 ${\mathrm{Peaks}}{≔}\begin{array}{c}\left[\begin{array}{cc}{0.792666220477213}& {1.}\\ {2.37799866143164}& {0.191629384082801}\\ {39.7918442679561}& {0.00199162461744076}\\ {38.2065118270017}& {0.00193605728902122}\\ {0.554866354334049}& {5.85893092353225}{×}{{10}}^{{-30}}\\ {40.0296441340993}& {5.85893092353221}{×}{{10}}^{{-30}}\\ {39.5540444018129}& {5.85845833721967}{×}{{10}}^{{-30}}\\ {1.03046608662038}& {5.85845833721954}{×}{{10}}^{{-30}}\\ {2.14019879528848}& {5.85498889041153}{×}{{10}}^{{-30}}\\ {38.4443116931448}& {5.85498889041150}{×}{{10}}^{{-30}}\\ {⋮}& {⋮}\end{array}\right]\\ \hfill {\text{252 × 2 Matrix}}\end{array}$ (3)
 > $R\left['\mathrm{plot}'\right]$

References

 Schmidt, Ralph O. "Multiple Emitter Location and Signal Parameter Estimation." IEEE transactions on antennas and propagation. Vol. 34, No. 3 (1986), pp. 276-280.

Compatibility

 • The SignalProcessing[MUSIC] command was introduced in Maple 2021.
 • For more information on Maple 2021 changes, see Updates in Maple 2021.

 See Also