# MUSIC implementation in Pyroomacoustics library

I am currently using the multiple signal classification (MUSIC) algorithm provided by the Pyroomacoustics library in order to estimate direction of arrival angles and it works pretty well for both narrowband signals and wideband signals.

However, after looking into the source code, I have few questions regarding the implementation. The code snippet below is taken from the Pyroomacoustics library. Here the main issue I have is the tensor denoted by X, which is defined as the STFT of the input signal with shape (M x n_fft/2+1 x num_frames), with M being the number of microphones.

Why do we need this tensor X in the first place? The MUSIC algorithm provided by the original paper (and anywhere else as far as I know) does not use any frequency representation. It only calculates the autocorrelation matrix in time domain and uses the steering vector which is known due to the geometry of the setup then follows it by eigenvalue decomposition.

    def _process(self, X):
"""
Perform MUSIC for given frame in order to estimate steered response
spectrum.
"""
# compute steered response
self.Pssl = np.zeros((self.num_freq, self.grid.n_points))
C_hat = self._compute_correlation_matricesvec(X)
# subspace decomposition
Es, En, ws, wn = self._subspace_decomposition(C_hat[None, ...])
# compute spatial spectrum
identity = np.zeros((self.num_freq, self.M, self.M))
identity[:, list(np.arange(self.M)), list(np.arange(self.M))] = 1
cross = identity - np.matmul(Es, np.moveaxis(np.conjugate(Es), -1, -2))
self.Pssl = self._compute_spatial_spectrumvec(cross)
if self.frequency_normalization:
self._apply_frequency_normalization()
self.grid.set_values(np.squeeze(np.sum(self.Pssl, axis=1) / self.num_freq))

# vectorized version
def _compute_correlation_matricesvec(self, X):
# change X such that time frames, frequency microphones is the result
X = np.transpose(X, axes=[2, 1, 0])
# select frequency bins
X = X[..., list(self.freq_bins), :]
# Compute PSD and average over time frame
C_hat = np.matmul(X[..., None], np.conjugate(X[..., None, :]))
# Average over time-frames
C_hat = np.mean(C_hat, axis=0)
return C_hat

# vectorized version
def _subspace_decomposition(self, R):
# eigenvalue decomposition!
# This method is specialized for Hermitian symmetric matrices,
# which is the case since R is a covariance matrix
w, v = np.linalg.eigh(R)

# This method (numpy.linalg.eigh) returns the eigenvalues (and
# eigenvectors) in ascending order, so there is no need to sort Signal
# comprises the leading eigenvalues Noise takes the rest

Es = v[..., -self.num_src :]
ws = w[..., -self.num_src :]
En = v[..., : -self.num_src]
wn = w[..., : -self.num_src]

return (Es, En, ws, wn)

def _compute_spatial_spectrumvec(self, cross):
mod_vec = np.transpose(
np.array(self.mode_vec[self.freq_bins, :, :]), axes=[2, 0, 1]
)
# timeframe, frequ, no idea
denom = np.matmul(
np.conjugate(mod_vec[..., None, :]), np.matmul(cross, mod_vec[..., None])
)
return 1.0 / abs(denom[..., 0, 0])



After looking into various sources related to MUSIC in Wikipedia, and the original paper, I see no FFT operations or spectrogram calculations. In addition, I have also implemented the MUSIC algorithm for a uniform linear array using the description from the paper and it also works properly.

So my questions are:

• Why the implementation above works, and how does it work without the need of a steering vector which is calculated from the information of microphone geometry?
• Is there a way to generate candidate steering vectors for arbitrary microphone array placements?
• Code is hard to decompose. Can you clarify your question by defining all the variables and writing out in equations what the code is doing? I think you have a good question here, but it's hard to help in its current form. Jul 1, 2022 at 11:21
• Yeah, I will try do it soon. Since I am also unsure what is being done in the source code, it might take some time though.
– Avio
Jul 1, 2022 at 12:55
• Also, what is the subspace_decomposition function doing? I assume it's doing SVD or something, but we don't know what its output variables are. Jul 1, 2022 at 13:13
• Still not entirely clear, but I believe $Es$ contains the signal eigenvectors, $ws$ the signal eigenvalues, $En$ the noise eigenvectors, and $wn$ the noise eigenvalues. The object self seems to have a pre-determined property num_src that decides which components are signals and which are "noise". Jul 1, 2022 at 15:55
• num_src denotes the number of sources, which is in line with MUSIC algorithm since it requires the number of sources active in the frame. I did not put the entire code in the question but the link github.com/LCAV/pyroomacoustics/blob/master/pyroomacoustics/doa/… has the full class. The main issue I have is the theoretical part and the fact that we are taking the autocorrelation of the STFT frames. I have never seen this before. I mean, the FFT of autocorrelation is PSD, but what exactly is the autocorrelation of FFT which is done here?
– Avio
Jul 1, 2022 at 17:40

## Why the Code Works

It's important to note that the MUSIC algorithm works on the microphone array dimension only. So it doesn't care how the data on the receive channels (microphones in this case) has been processed as long as the same thing is done to every channel.

In your case, it appears that MUSIC is being combined with a STFT on each channel so that you can find the direction of arrival (DOA) of signals in a particular frequency band. Thus, the data in each channel has been STFTed, giving a spectrogram on each channel, but it won't matter to MUSIC.

The code takes the data in the frequency bins of interest, and forms an estimate of the channel correlation matrix, $$C$$, which is $$M$$-by-$$M$$. It then performs an eigen-decomposition of $$C$$ in order to use the eigenvectors corresponding to low eigenvalues as a basis for a noise (anti-signal) sub-space.

The way MUSIC works is by estimating a basis for the the noise sub-space (non-signal of interest), and then projecting steering vectors for every direction into that noise sub-space. Directions that have large signal components (rather than noise) will have very small values in the noise basis, so you can find them by taking the reciprocal of the projection, and finding peaks: \begin{align} C &= \text{Channel covariance matrix, estimated by averaging over time bins} \\ E &= \text{Eigenvectors of C, in ascending order} \\ \text{num_src} &= \text{Number of signals} \\ E_s &= \text{last num_src eigenvectors of E} \\ R_{noise} &= \text{Noise covariance matrix} = I - E_s \cdot E_s^H \\ V &= \text{mode_vec = Array steering vectors for all directions of interest} \\ d &= \text{denom} = V^HR_{noise}V \\ \end{align}

Plotting $$1/d$$ would give peaks at the directions corresponding to signals. For another answer that explains MUSIC fairly clearly, see this question.

## How to Calculate Steering Vectors for an Arbitrary Array

Let $$X_{array}$$ be a 3-by-M matrix of 3-D coordinate locations for the M microphones, and $$\hat D_{dir}$$ be a 3-by-N matrix of N unit vectors that point in the directions of interest (e.g., perhaps equally spaced in azimuth relative to the array from -90 to +90 degrees). Then, the N steering vectors are given by: $$V = e^{j\frac{2\pi}{\lambda}\cdot X_{array}^T \hat D_{dir}}$$

where $$V$$ is M-by-N and $$\lambda$$ is the wavelength.

Interesting fact: The above formula reduces to a discrete Fourier transform (DFT) when the array is a uniform linear array (ULA), and the look directions are chosen to be equally spaced in sine-space (e.g., $$sin(\theta) = [-1,\, -1+\Delta,\, -1+2\Delta,... 1-\Delta]$$). See this answer, for example.