MUSIC implementation in Pyroomacoustics library

Question

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?

Answer 1

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.