SignalProcessing
MUSIC
performs the Multiple Signal Classifier (MUSIC) Method on a signal
Calling Sequence
Parameters
Options
Description
Examples
References
Compatibility
MUSIC( signal, options )
signal
-
1-D rtable or list of data, which must contain two or more elements
dimension: Positive integer for the dimension of the signal subspace, which must be less than the size of the signal. The default is NULL.
plotoptions: List containing any additional options to be passed to the plotting command when creating the pseudo-power spectrum. The default is [].
samplerate: Positive numeric value for the sampling rate. The default is 1.0.
size: 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: 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: 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: Either true or false, indicates if the windowing function is to be normalized. The default is true.
output: The type of output. The supported options are:
eigenvalues: Returns a Vector, of float[8] datatype and length m=numelems⁡signal, 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=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 10·log10⁡P 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].
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 ω 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
Si=ⅇI⁢2⁢πr⁢i−1⁢ω
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:
Start with the 1-D signal X of size m.
Compute the mean μ of X.
Compute the auto-correlation Vector A of X−μ, which is of size m.
Take the auto-correlation Matrix C of A to be the square Matrix of size m with elements given by:
Ci,j=0.5·Amodp⁡j−i,m+1+Amodp⁡i−j,m+1&conjugate0;
Here, C is Hermitian, positive definite, Toeplitz, and circulant and, consequently, the eigenvalues are real-valued and positive.
Determine the sorted (in decreasing size) eigenvalues Vector U and eigenvectors Matrix V, with column j of V being the eigenvector associated with eigenvalue Uj. 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.
Define the frequencies Vector F of size n, where n is the length of the pseudo-power spectrum, element-wise as
Fk=k−1⁢rn
where r is the sample rate.
The pseudo-power spectrum is the Vector of size n with elements given by the following:
Pk=1∑j=d+1m⁡∑i=1m⁡ⅇI⁢2⁢πr⁢i−1⁢Fk·Vi,j2
The pseudo-power spectrum is then normalized with respect to the maximum value.
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=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:
Vi,j=1m·ⅇI⁢2⁢πm⁢i⁢j−1
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:
Si,k=ⅇI⁢2⁢πr⁢i−1⁢Fk
Now, the inner sum of the denominator of Pk 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:
Pk=1∑j=1m⁡Tj,k2
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 of 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.
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.
with⁡SignalProcessing:
with⁡Statistics:
Real-valued signal
Consider the following real-valued signal without noise:
m≔28:
fs≔100.0:
a≔5.0,3.0,−7.0:
f≔10.25,10.40,20.35:
X≔Vector⁡m,seq⁡add⁡ai⁢sin⁡fi⋅2⁢π⁢k−1fs,i=1..3,k=1..m,datatype=float8:
SignalPlot⁡X,samplerate=fs,color=blue,title=Pure Signal Plot,font=Verdana,15,size=800,400
Periodogram⁡X,samplerate=fs,color=blue,title=Pure Signal Periodogram,font=Verdana,15,size=800,400
Now, construct a Vector of noise, and add it to the pure signal:
N≔Vectorcolumn⁡Sample⁡RandomVariable⁡Normal⁡0,0.5,m:
Y≔X+N:
SignalPlot⁡Y,samplerate=fs,color=blue,title=Noisy Signal Plot,font=Verdana,15,size=800,400
Periodogram⁡Y,samplerate=fs,color=blue,title=Noisy Signal Periodogram,font=Verdana,15,size=800,400
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≔210:
R≔MUSIC⁡Y,samplerate=fs,dimension=d,size=n,output=record:
Peaks≔Rpeaks
Rplot
We can also specify a threshold for the eigenvalues, as opposed to the dimension of the signal space:
t≔0.10:
E,P≔MUSIC⁡Y,samplerate=fs,threshold=t,size=n,output=eigenvalues,peaks
Complex-valued signal
Consider a complex-valued signal with noise:
t≔t:
f≔cos⁡5⁢t+3⁢I⁢sin⁡15⁢t
a≔0:
b≔2⁢π:
fs,X,Y≔GenerateSignal⁡f,t=a..b,m,noisedeviation=1.0,noisetype=additive,output=samplerate,puresignal,signal
Now, compare the periodogram of the pure signal with the pseudo-power spectrum of the noisy signal:
PowerSpectrum⁡X,variety=signal,samplerate=fs,output=periodogram,powerscale=dB,periodogramoptions=color=blue,title=Pure Signal Periodogram,font=Verdana,15
d≔4:
n≔4⁢m:
Schmidt, Ralph O. "Multiple Emitter Location and Signal Parameter Estimation." IEEE transactions on antennas and propagation. Vol. 34, No. 3 (1986), pp. 276-280.
Starting with Maple 2023, external C code is used for the auxiliary procedures that were formerly implemented in Maple and could be compiled by passing the compiled option. The compiled option is now deprecated, but it is still accepted for backwards compatibility.
The SignalProcessing[MUSIC] command was introduced in Maple 2021.
For more information on Maple 2021 changes, see Updates in Maple 2021.
The SignalProcessing[MUSIC] command was updated in Maple 2023.
See Also
SignalProcessing[GenerateSignal]
SignalProcessing[Periodogram]
SignalProcessing[PowerSpectrum]
SignalProcessing[SignalPlot]
SignalProcessing[Welch]
Statistics
Download Help Document