A Signal Subspace Rotation Method for Localization of Multiple Wideband Sound Sources

TL;DR

This paper introduces a novel wideband sound source localization method using signal subspace rotation, which improves accuracy and reduces computational complexity without needing prior knowledge.

Contribution

The paper presents a new wideband DOA estimation technique based on signal subspace rotation, extending narrowband algorithms efficiently to wideband signals.

Findings

Effective in localizing multiple wideband sound sources

Reduces computational complexity compared to traditional methods

Demonstrates robustness and accuracy in experiments

Abstract

In this paper, the problem of extending narrowband multichannel sound source localization algorithms to the wideband case is addressed. The DOA estimation of narrowband algorithms is based on the estimate of inter-channel phase differences (IPD) between microphones of the sound sources. A new method for wideband sound source DOA estimation based on signal subspace rotation is present. The proposed algorithm normalizes the narrowband signal statistics by rotating the estimated signal subspace to the wideband counterpart in the eigenvector domain. Then the wideband DOA estimate can be obtained by estimating the normalized IPD from these wideband signal statistics. In addition to requiring less computational complexity compared to repeating the narrowband algorithms for all relevant frequencies of wideband signals, the proposed method also does not require any additional prior knowledge. The experimental results demonstrate the efficacy and the robustness of the proposed method.

Tables

Table 1. Table 1 : Single white noise source

	SNR=10dB		SNR=0dB
Algorithm	MAE	SDE	MAE	SDE
hist-ESPRIT	${2.10}^{\circ }$	${2.47}^{\circ }$	${3.44}^{\circ }$	${2.47}^{\circ }$
CSS	${3.81}^{\circ }$	${5.58}^{\circ }$	${5.52}^{\circ }$	${7.34}^{\circ }$
Proposed	${1.41}^{\circ }$	${1.62}^{\circ }$	${1.68}^{\circ }$	${1.75}^{\circ }$

Table 2. Table 2 : Two simultaneous white noise sources

	SNR=10dB		SNR=0dB
Algorithm	MAE	SDE	MAE	SDE
hist-ESPRIT	${2.57}^{\circ }$	${5.64}^{\circ }$	${3.28}^{\circ }$	${7.47}^{\circ }$
CSS	${6.8}^{\circ }$	${9.67}^{\circ }$	${6.5}^{\circ }$	${12.22}^{\circ }$
Proposed	${2.82}^{\circ }$	${3.55}^{\circ }$	${2.75}^{\circ }$	${5.62}^{\circ }$

Table 3. Table 3 : Two simultaneous talkers, SNR ∈ [ − 5 , 10 ] absent 5 10 \in[-5,10] dB

Algorithm	MAE	SDE
hist-ESPRIT	${4.28}^{\circ }$	${6.47}^{\circ }$
CSS	${6.85}^{\circ }$	${15.37}^{\circ }$
Proposed	${5.75}^{\circ }$	${10.62}^{\circ }$

Equations

$$X_p(f_i)=\sum _{q=1}^Q{A}_{p,q}(f_i)S_q(f_i)+N_p(f_i),$$

$$\mathbf{R}(f_i)=\mathbf{A}(f_i){\mathbf{R}}_s(f_i){\mathbf{A}}^{\text{H}}(f_i)+{\mathbf{R}}_n(f_i),$$

$$\begin{array}{r}\begin{array}{c}A(f_i{)}_{p,q}=e^{-j2\pi f_i\Delta t_{p,q}},\\ \Delta t_{p,q}\stackrel{\text{def}}{=}(p-1)\Delta d{c}^{-1}\sin {\theta }_q.\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{c}\mathbf{A}(f_1{)}^{\circ f_1^{-1}}=\mathbf{A}(f_2{)}^{\circ f_2^{-1}},\end{array}\end{array}$$

$$\begin{array}{r}\begin{array}{c}{f}_a=\lfloor \frac{c}{2\Delta d\sin \theta }\rfloor ,\end{array}\end{array}$$

$$\mathbf{R}(f_i)\mathbf{U}(f_i)=\mathbf{U}(f_i)\mathbf{\Lambda }(f_i),$$

$$\mathbf{U}(f_i)=\left[{\mathbf{U}}_s(f_i)|{\mathbf{U}}_n(f_i)\right].$$

$${\mathbf{U}}^{{}^{\prime }}(f_i)=\mathbf{U}(f_i{)}^{\circ f_i^{-1}}.$$

$${\mathbf{U}}_{ss}^{{}^{\prime }}=\sum _i\beta (f_i){\mathbf{U}}_s^{{}^{\prime }}(f_i),$$

$$\beta (f_i)=\text{trace}\{{\mathbf{\Lambda }}_s(f_i)\},$$

$${\mathbf{U}}_{ss2}^{{}^{\prime }}={\mathbf{U}}_{ss1}^{{}^{\prime }}\mathbf{\Phi },$$

$$\mathbf{\Phi }=\text{diag}\{e^{-j2\pi \Delta t_1},...,e^{-j2\pi \Delta t_P}\}$$

$${\mathbf{E}}_s={\mathbf{U}}_{ss}^{{}^{\prime }}\mathbf{T},$$

$$\begin{array}{r}\begin{array}{c}{\mathbf{E}}_{s1}={\mathbf{U}}_{ss1}^{{}^{\prime }}\mathbf{T},\\ {\mathbf{E}}_{s2}={\mathbf{U}}_{ss2}^{{}^{\prime }}\mathbf{T}={\mathbf{U}}_{ss1}^{{}^{\prime }}\mathbf{\Phi T},\\ {\mathbf{E}}_{s2}={\mathbf{E}}_{s1}\mathbf{\Psi },\end{array}\end{array}$$

$${\mathbf{R}}^{{}^{\prime }}(f_i)={\mathbf{U}}^{{}^{\prime }}(f_i)\mathbf{\Lambda }(f_i){\mathbf{U}}^{{}^{\prime }}(f_i{)}^{-1}.$$

$${\mathbf{R}}^{{}^{\prime \prime }}=\sum _i\beta (f_i){\mathbf{R}}^{{}^{\prime }}(f_i).$$

$${\mathbf{R}}_1^{{}^{\prime \prime }}={\mathbf{U}}_s^{{}^{\prime }{f}_1{f}_0^{-1}}(f_0)\mathbf{\Lambda }(f_0){\mathbf{U}}_s^{{}^{\prime }{f}_1{f}_0^{-1}}(f_0{)}^{-1}$$

$$\begin{array}{r}\begin{array}{c}{\mathbf{R}}_{i+1}^{{}^{\prime \prime \prime }}=\beta (f_i){\mathbf{U}}_s^{{}^{\prime }}(f_i)\mathbf{\Lambda }(f_i){\mathbf{U}}_s^{{}^{\prime }}(f_i{)}^{-1}+{\mathbf{R}}_i^{{}^{\prime \prime }}\\ {\mathbf{R}}_{i+1}^{{}^{\prime \prime \prime }}\stackrel{\text{def}}{=}{\mathbf{U}}_s^{{}^{\prime \prime \prime }}{\mathbf{\Lambda }}^{{}^{\prime \prime \prime }}{\mathbf{U}}_s^{{}^{\prime \prime \prime }-1}\\ {\mathbf{R}}_{i+1}^{{}^{\prime \prime }}={\mathbf{U}}_s^{{}^{\prime \prime \prime }{f}_{i+2}{f}_{i+1}^{-1}}{\mathbf{\Lambda }}^{{}^{\prime \prime \prime }}{\mathbf{U}}_s^{{}^{\prime \prime \prime }{f}_{i+2}{f}_{i+1}^{-1}}\end{array}\end{array}$$

$$f_{a_0}\stackrel{\text{def}}{=}\lfloor \underset{\theta }{\mathrm{argmin}}{f}_a\rfloor =\lfloor \frac{c}{2\Delta d}\rfloor .$$

Full text

**Index Terms— ** Sound source localization, wideband, signal subspace rotation

1 Introduction

Sound source localization is an important component in many multichannel signal processing systems aiming, e.g., at source tracking, signal separation, enhancement and noise suppression [1, 2] . Traditional sound source localization algorithms such as GCC-PHAT [3] estimates the time delay of arrival (TDOA) between a pair of microphones to localize a single source. SRP-PHAT [4] as an extension can further localize multiple sources simultaneously. Another group of algorithms that resolves the simultaneous multiple source localization problem is based on high-resolution subspace techniques, such as MUSIC [5] and ESPRIT [6] . For highly reverberant senarios, the Independent Component Analysis (ICA) based algorithms are proposed  [7, 8].

ESPRIT [6] exploits the algebraic properties of the spatial covariance matrix. This method features good performance for narrowband signals. However, it is not directly applicable to wideband signals. Hence, it process the wideband signals by estimating the individual narrowband results in the Short Time Fourier Transform (STFT) domain. Post-processing schemes such as the histogram method only exploit the narrowband estimation results while the mutual information between different frequency bins is not considered. Especially for known type of the sound sources, such as speech, the process based on clustering STFT bins such as [9, 10] can further improve the accuracy.

Many studies focus on extensions of narrowband localization algorithms to wideband signals. The methods of [11, 12] transform signals to cylindrical/spherical harmonics for circular/spherical microphone arrays. Then the rotational invariances are normalized through frequencies and the further estimation via ESPRIT is based on wideband signals. A Coherent Signal Subspace (CSS) method that derives the so-called focusing matrices based on steering vectors was introduced in [13]. The focusing matrices aim to adapt each narrowband signal covariance matrix to generalize to the wideband case. Theoretical analysis and experiments show that the focusing matrices lead to a wideband MUSIC algorithm [13], and in [14, 15, 16], the focusing matrices extend ESPRIT to wideband signals. As pointed out in [17, 18], the focusing matrices method adapts the covariance matrices via linear transformations that are dependent on frequency and DOAs. Therefore, these extensions require a priori knowledge of DOA estimates and on the array manifold for initializing the transformation matrices. To avoid this requirement, the algorithm in [17] estimates an AR model for the transmission channel of the wideband signals as a priori knowledge for wideband localization. The performance of this algorithm crucially depend on how well the AR model is estimated. However, the model that they proposed and the conventional models are not suitable for non-stationary scenarios.

A novel method that adapts the narrowband ESPRIT to wideband signals is introduced in this work. The key idea is the rotation of the eigenvectors that span the signal subspace by the corresponding frequency, and to reconstruct the covariance matrices for each frequency bin based on the rotated eigenvectors. In comparison with the focusing matrices approach, the proposed method utilize the wideband signal second-order statistics while does not require a priori knowledge of the DOA which is typically obtained by repeating narrowband MUSIC [16]. Experiments with simulated and real recordings show the advantages of the new method regarding performance and computational complexity.

2 Signal model

We assume a microphone array which fits the requirements for using ESPRIT, i.e., to localize $Q$ sources. We use $P$ microphones ($P>Q$) that form a linear array with a uniform spacing $\Delta d$ and with mutually uncorrelated zero-mean sensor noise. The environment is assumed with free-field and far-field conditions. The source-microphone model in the STFT domain is defined as

[TABLE]

where ${X}_{p}(f_{i})$ and ${N}_{p}(f_{i})$ denote the observation and the noise at the $p$-th microphone in the STFT domain at frequency bin $f_{i}$, respectively, ${S}_{q}(f_{i})$ denotes the $q$-th source signal, and $A_{p,q}(f_{i})$ denotes the frequency response of the propagation path which are given by the rows of ${A(f_{i})}$ for $q$-th source arriving at the $p$-th microphone.

The power spectral density matrix for all Q source signal components at the $i$-th frequency bin $f_{i}$ is denoted as $\mathbf{R}_{s}(f_{i})$ and the power spectral density matrix for the observed noise is denoted as $\mathbf{R}_{n}(f_{i})$. The narrowband component covariance matrix $\mathbf{R}(f_{i})$ of the microphone signals can be expressed as

[TABLE]

where $*^{\text{H}}$ denotes the Hermitian transpose, $\mathbf{A}(f_{i})$ denotes a $P\times Q$ matrix that captures the steering vectors at the frequency bin $f_{i}$. The element of $\mathbf{A}(f_{i})$ at the $p-{\text{th}}$ row and $q-{\text{th}}$ column, $A(f_{i})_{p,q}$, can be expressed as

[TABLE]

Therefore, the steering vectors are frequency-dependent and the relationship between different frequency bins $f_{1},f_{2}$ can be expressed as,

[TABLE]

where $\{*\}^{\circ f^{-1}}$ describes element-wise exponentiation by $f^{-1}$ and $f_{1},f_{2}$ denote any two frequencies below the spatial aliasing frequency $f_{a}$,

[TABLE]

where $\lfloor*\rfloor$ denotes a function that returns the next lower frequency bin.

For our proposed evolution of ESPRIT to wideband source localization, the first step is to estimate the narrowband steering vector. Then the IPDs of each source are obtained from the estimation result. Therefore, the IPDs are depending on frequency and ESPRIT has to be repeated for each narrowband to localize wideband sources, as it is also the case for other narrowband localization algorithms  [5, 19]. A solution to estimate a unique IPD for each source through the frequency bins is to rotate the estimated steering vector based on (4) as it is shown in the following.

3 Proposed Approach

The proposed approach considers a novel way for adapting ESPRIT to a single or multiple wideband sources. In case of a single source scenario, for each narrowband component, the proposed approach rotates the eigenvectors that span the signal subspaces by normalizing the IPD. In case of multiple sources, it reconstructs covariance matrices from the rotated eigenvectors.

To obtain an estimate of the signal subspace, the least-squares (LS) criterion is conventionally employed to find $Q$ vectors that describe the signal subspace and $P-Q$ noise vectors to represent the noise subspace in an LS sense (see, e.g., [5, 6]). For $Q$ independent target sources, the $P\times P$ matrix $\mathbf{A}(f_{i})\mathbf{R}_{s}(f_{i})\mathbf{A}^{\text{H}}(f_{i})$ is at least of rank $Q$ and positive semidefinite by construction. Therefore, by taking the eigenvalue decomposition of $\mathbf{R}(f_{i})$, the eigenvectors corresponding to the largest $Q$ eigenvalues are assumed to be optimum to span the signal subspace at the frequency bin $f_{i}$ (see, e.g., [6]).

The eigenvalue decomposition of $\mathbf{R}(f_{i})$ can be expressed as

[TABLE]

where $\mathbf{U}$ denotes the eigenvector matrix and $\boldsymbol{\Lambda}$ is a diagonal matrix that contains the corresponding eigenvalues. The eigenvector spans the signal subspace is denoted as $\mathbf{U}_{s}$, and the eigenvector for the noise subspace is denoted as $\mathbf{U}_{n}$,

[TABLE]

By the relationship defined in (4), to rotate the estimated subspaces such that they become frequency-independent for all frequency subbands, the estimated signal subspace rotation is defined as

[TABLE]

Then the frequency component of the IPDs are assumed to cancelled.

ESPRIT is based on source subspace analysis [6]. For the single-source ESPRIT, the estimated source subspace is described by the vector $\mathbf{U}^{{}^{\prime}}_{s}(f_{i})$. By weighting and summing the rotated eigenvectors that span the source subspace $\mathbf{U}^{{}^{\prime}}_{s}(f_{i})$, the estimated wideband source subspace $\mathbf{U}^{{}^{\prime}}_{ss}$ can be obtained as

[TABLE]

where $\beta(f_{i})$ denotes a frequency-dependent weighting factor. The eigenvalues of the sources are relevant to the signal power at a certain frequency and can be assumed to reflect the reliability of the signal subspaces estimates. Therefore, the weighting function is chosen as

[TABLE]

where $\text{trace}\{*\}$ denotes a function which returns the trace of the matrix and $\boldsymbol{\Lambda}_{s}(f_{i})$ denotes the eigenvalue matrix of the signal subspace.

Following ESPRIT [6], the submatrices satisfy the invariance relation

[TABLE]

where $\mathbf{U}^{{}^{\prime}}_{ss1}$ and $\mathbf{U}^{{}^{\prime}}_{ss2}$ denote the vectors contain the first and last $P-1$ elements of $\mathbf{U}^{{}^{\prime}}_{ss}$, respectively. Since the combined vector $\mathbf{U}^{{}^{\prime}}_{ss}$ is assumed to span the wideband signal subspace $\mathbf{E}_{s}$, it holds that

[TABLE]

where $\mathbf{T}$ is a non-singular matrix. Therefore, the subspaces of the two subarrays can be defined as

[TABLE]

where $\mathbf{\Psi}=\mathbf{T}\mathbf{\Phi}\mathbf{T}^{-1}$. $\mathbf{\Psi}$ can then be obtained from (14) by applying a standard least-squares or total least-squares solver. By realizing that the eigenvalues of $\mathbf{\Psi}$ are the diagonal elements of $\mathbf{\Phi}$ the locations of the sources can be estimated.

For the multi-source ESPRIT, each estimated narrowband source subspace is described by the columns of the matrix $\mathbf{U}^{{}^{\prime}}_{s}(f_{i})$. The order of the columns of the matrix is generally not known and may be different for different frequency bins. Therefore, the estimated signal subspaces cannot be combined using (9).

A simpler and more robust solution than the source subspace identification is to reconstruct estimated signal subspaces back to the form of a $P\times P$ covariance matrix like $\mathbf{R}^{{}^{\prime}}(f_{i})$. This can be achieved by the inverse process of eigenvalue decomposition (6) using the frequency component cancelled matrix $\mathbf{U}^{{}^{\prime}}(f_{i})$,

[TABLE]

With the same weighting factor definition in (10), the reconstructed covariance matrices are calculated by

[TABLE]

Therefore, the eigenvectors that span the estimated signal subspaces are remixed and accumulated through frequencies in $\mathbf{R}^{{}^{\prime\prime}}$. By using the conventional narrowband ESPRIT, wideband signal subspace can be separated back and wideband DOAs can then be estimated.

4 Implementation

In the estimated signal subspace rotation step (8), when $f$ gets larger, a finer quantization is required to limit the effect of quantization errors on the DOA estimation. An iterative accumulation method is proposed to solve this numerical sensitivity problem. In each iteration, the estimated signal subspace from the $i$-th $(i\in\mathbb{N_{+}})$ frequency bin is rotated to adapt the next frequency bin ($i+1$) to reconstruct the covariance matrix $\mathbf{R}^{{}^{\prime\prime}}_{i+1}$. After weighting and summing to $\mathbf{R}^{{}^{\prime\prime}}_{i}$, the new accumulated covariance matrix $\mathbf{R}^{{}^{\prime\prime\prime}}_{i+1}$ is rotated to the next frequency bin ($\mathbf{R}^{{}^{\prime\prime}}_{i+1}$) for the next iteration. The rotation step is then as small as the power of $f_{i+2}f_{i+1}^{-1}$. The iteration process can be expressed as

Initialization:

[TABLE]

In each iteration:

[TABLE]

Because of the spatial aliasing problem, the iteration continues until the lowest aliasing frequency, denoted by $f_{a_{0}}$,

[TABLE]

To utilize higher frequencies, the replication method [20] can potentially be useful for detecting the aliasing frequency for each source. Similar to multi-source localization, the DOA can be obtained using the matrix $\mathbf{R}^{{}^{\prime\prime}}(f_{a})$.

5 Evaluation

A set of experiments was performed in order to evaluate the performance of the algorithm using real-world recordings. The evaluation includes comparisons to the narrowband ESPRIT with the histogram method (hist-ESPRIT) [21] and to the CSS method [13].

5.1 Experimental setup

The recordings were captured by a uniform linear array (ULA) with five microphones in a low-reverberation lab ($\text{T60}\approx 0.2s$ of size $9m\times 8m\times 3m$). The microphone model was AKG C562CM. The spacing between microphones $\Delta d$ was $0.044$m. The background noise was white noise, and it was played back via 22 surround speakers to emulate diffuse background noise. The set of sources contains white sources (independent from the background noise) and a set of speech recordings selected from the GRID Corpus [22]. There were two sources located at an angle of $45^{\circ}$ and at $-45^{\circ}$ at a distance of $3m$ from the center of the microphone array. In order to evaluate the robustness of the proposed localization method, three test conditions were used: (1) single white source (2) two competing white sources, (3) two competing talkers. The total length of the recordings is $800$s. The sampling rate was $16000$ samples per second.

The proposed method, hist-ESPRIT [21] and CSS [13] were processed block-wise with $50\%$ overlaps. The length of the block was $1024$ samples. The frequency band used for these algorithms was up to $3800$Hz, below the aliasing frequency.

5.2 Results

The scenario of the first experiments consisted of a single white noise as the source under $10$dB and [math]dB background diffuse white noise. The source signal was played in an alternating fashion from two speakers located at $45^{\circ}$ and $-45^{\circ}$. The result is shown in Table 1. The performance of the algorithm is evaluated by the mean absolute error (MAE) and the standard deviation of the error (SDE). It can be seen from the results that the proposed algorithm is more accurate and more stable under this strong background noise scenario.

With the same background noise conditions, the second experiment features two white noise sources being played by both speakers in a competing fashion. The result is shown in Table 2. The MAE of hist-ESPRIT is slightly lower than the proposed algorithm as the error influence of estimation in low frequencies is slightly lower, while the SDE for the proposed algorithm is clearly superior.

The final experiment used the same scenario as the second experiment, but two simultaneously active speech sources were to be localized. The level of the speech signals changed over time, and the estimated SNR was in the range of $[-5,10]$dB. The result is shown in Table 3. The performance of the proposed algorithm is better than CSS algorithm, but worse than the hist-ESPRIT algorithm, especially on SDE. The speech signal is spectrally sparse and the energy is compacted in the fundamental frequency and its harmonics. In contrast, under the white noise background, the intervals between harmonics in the spectrum are noisier (lower SNR). With the inspection of the narrowband processing, the estimated DOA results from the intervals had large errors (up to $90^{\circ}$). The reason behind the reduction in performance of the proposed algorithm is the inferior robustness of the least squares criterion in comparison to the histogram method. In the white sources experiments above, the SNR is constant for the entire frequency range. Therefore, those experiments had better results.

Through the above experiments, compared to the computing time of hist-ESPRIT algorithm, CSS algorithm (with $0.1^{\circ}$ resolution of the spatial spectrum) was $13\%$ faster, and the proposed algorithm was $22.6\%$ faster.

6 Conclusion

A wideband signal subspace DOA estimation approach is presented. The proposed signal subspace rotation method and the narrowband signal covariance matrix reconstruction method are high and outperform the existing conventional approaches in computational complexity. Additionally, the proposed approach avoids the necessity of additional prior knowledge for extending the narrowband ESPRIT to a wideband scheme. Experiments based on real recordings validate the effectiveness and low computational complexity of the proposed method. As part of the future work, the proposed algorithm should be analyzed in the context of solving the spatial aliasing problem, in order to improve the localization accuracy of the spectrally sparse sources in environments with background noise.