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?