Multichannel Linear Prediction for Blind Reverberant Audio Source Separation

TL;DR

This paper introduces a novel use of multichannel linear prediction as a pre-processing step to improve blind source separation in reverberant environments, simplifying the problem and enhancing separation performance.

Contribution

It proposes a new application of MCLP as a pre-processing step, reducing reverberant separation to a non-reverberant problem, validated with real microphone array recordings.

Findings

Pre-processing with MCLP improves source separation in reverberant conditions.

Theoretical reduction of reverberant separation to non-reverberant case.

Experimental validation with real recordings confirms effectiveness.

Abstract

A class of methods based on multichannel linear prediction (MCLP) can achieve effective blind dereverberation of a source, when the source is observed with a microphone array. We propose an inventive use of MCLP as a pre-processing step for blind source separation with a microphone array. We show theoretically that, under certain assumptions, such pre-processing reduces the original blind reverberant source separation problem to a non-reverberant one, which in turn can be effectively tackled using existing methods. We demonstrate our claims using real recordings obtained with an eight-microphone circular array in reverberant environments.

Tables

Table 1. TABLE I: PESQ and SIR scores for Exp.  3 . (LP : Multichannel Linear Prediction, G : Geometric Source Separation, P : Post filter)

PESQ

Equations

$$y(t)\approx x(t)\ast \tilde{h}(t),$$

$$h(t)=\tilde{h}\bigl{(}t+D/c\bigr{)}.$$

$$y_i(t)\approx x(t)\ast h_i(t-{\tau }_i)=x(t-{\tau }_i)\ast h_i(t).$$

$$y_i(t)\approx \sum _{k=1}^K{x}_k(t)\ast h_{i,k}(t-{\tau }_{i,k})=\sum _{k=1}^K{x}_k(t-{\tau }_{i,k})\ast h_{i,k}(t)$$

$$y(n)=x(n)\ast h(n).$$

$$Y(s,\omega )\approx X(s,\omega )H(\omega ),$$

$$Y(s,\omega )\approx X(s,\omega )\stackrel{s}{\ast }\tilde{H}(s,\omega ),$$

$$X(s,\omega )\stackrel{s}{\ast }\tilde{H}(s,\omega )=\sum _{l}X(s-l,\omega )\tilde{H}(l,\omega ).$$

$$Y_{i}(s,\omega)\approx\sum_{k=1}^{K}\,X_{k}(s,\omega)\,\overset{s}{\ast}\,\Bigl{(}\tilde{H}_{i,k}(s,\omega)\,\exp(-j\omega\,\tau_{i,k})\Bigr{)}.$$

$$y_i=x\ast h_i,\text{ for }i=1,2,\dots ,N,$$

$$\left[\begin{array}{c}{\mathcal{Y}}_1(z)\\ ⋮\\ {\mathcal{Y}}_N(z)\end{array}\right]=\left[\begin{array}{c}{\mathcal{H}}_1(z)\\ ⋮\\ {\mathcal{H}}_N(z)\end{array}\right]\mathcal{X}(z).$$

$$\left[\begin{array}{ccc}{\mathcal{G}}_1(z)& \dots & {\mathcal{G}}_N(z)\end{array}\right]\left[\begin{array}{c}{\mathcal{H}}_1(z)\\ ⋮\\ {\mathcal{H}}_N(z)\end{array}\right]=1.$$

$$\begin{array}{rl}{\mathcal{G}}_r(z)& =1-z^{-d}{\mathcal{U}}_r(z)\text{ for a reference index }r,\\ {\mathcal{G}}_i(z)& =-z^{-d}{\mathcal{U}}_i(z)\text{ for }i=1,\dots ,N,i\ne r,\end{array}$$

$$\underset{u_1,\dots ,u_N}{\min }P\left(y_r(s)-\sum _{i=1}^N\sum _{l=0}^L{y}_i(s-d-l)u_i(l)\right),$$

$$\left[\begin{array}{ccc}{\mathcal{G}}_1(z)& \dots & {\mathcal{G}}_N(z)\end{array}\right]\left[\begin{array}{c}{\mathcal{H}}_1(z)\\ ⋮\\ {\mathcal{H}}_N(z)\end{array}\right]=c,$$

$$h_i(s)={\tilde{H}}_{i,1}(s,\omega )\exp (-j\omega {\tau }_{i,1}).$$

$${\hat{y}}_r(s)=h_r(0)\exp (-j\omega {\tau }_{r,1})x(s).$$

$$h_i(0)=h_{i^{\prime }}(0)\text{ for any distinct pair }(i,i^{\prime }),$$

$${\hat{y}}_r(s)=c^{\prime }\exp (-j\omega {\tau }_{r,1})x(s),$$

$${\tilde{H}}_i(0,\omega )={\tilde{H}}_{i^{\prime }}(0,\omega ),$$

$$\sum _i{g}_i(n)\ast h_i(n)=c_r\delta (n)$$

$$y_i=\sum _{k=1}^K{x}_k\ast h_{i,k},\text{ for }i=1,2,\dots ,N.$$

$$\left[\begin{array}{c}{\mathcal{Y}}_1(z)\\ ⋮\\ {\mathcal{Y}}_N(z)\end{array}\right]=\underset{\mathbf{H}}{\underbrace{\left[\begin{array}{ccc}{\mathcal{H}}_{1,1}(z)& \dots & {\mathcal{H}}_{1,K}(z)\\ & \ddots & \\ {\mathcal{H}}_{N,1}(z)& \dots & {\mathcal{H}}_{N,K}(z)\end{array}\right]}}\left[\begin{array}{c}{\mathcal{X}}_1(z)\\ \dots \\ {\mathcal{X}}_K(z)\end{array}\right]$$

$$\left[\begin{array}{ccc}{\mathcal{G}}_{1,1}(z)& \dots & {\mathcal{G}}_{1,N}(z)\\ & \ddots & \\ {\mathcal{G}}_{K,1}(z)& \dots & {\mathcal{G}}_{K,N}(z)\end{array}\right]\left[\begin{array}{ccc}{\mathcal{H}}_{1,1}(z)& \dots & {\mathcal{H}}_{1,K}(z)\\ & \ddots & \\ {\mathcal{H}}_{N,1}(z)& \dots & {\mathcal{H}}_{N,K}(z)\end{array}\right]=I,$$

$$\left[\begin{array}{ccc}{\mathcal{G}}_1(z)& \dots & {\mathcal{G}}_N(z)\end{array}\right]\left[\begin{array}{ccc}{\mathcal{H}}_{1,1}(z)& \dots & {\mathcal{H}}_{1,K}(z)\\ & \ddots & \\ {\mathcal{H}}_{N,1}(z)& \dots & {\mathcal{H}}_{N,K}(z)\end{array}\right]=c^{(r)},$$

$$c_k^{(r)}={\mathcal{H}}_{r,k}(\infty ),\text{ for }k=1,2,\dots ,K.$$

$$h_{i,k}(s)={\tilde{H}}_{i,k}(s,\omega )\exp (-j\omega {\tau }_{i,k}),$$

$$\sum _{i=1}^N{g}_i(s)\ast h_{i,k}(s)=c_k^{(r)}\delta (s)$$

$$\sum _{i=1}^N{g}_i(s)\ast y_i(s)=c^{\prime }\sum _{k=1}^K{x}_k(s)\exp (-j\omega {\tau }_{r,k}),$$

$${\hat{y}}_r(s)=y_r(s)-\sum _{i=1}^N\sum _{l=0}^{L-1}{y}_i(s-d-l){\hat{u}}_i^{(r)}(l).$$

Full text

I Introduction

In anechoic environments, audio source separation with a microphone array can be achieved by beamforming and post-filtering [7]. On the other hand, in reverberant environments, reverberation causes significant degradation in the quality and intelligibility of the beamformer output, especially if the room impulse responses (RIR) are not taken into account or are not known precisely. This is a serious limitation because in practice, RIRs are usually not known with precision and they can vary wildly with respect to the positions of the microphone/source pair.

A related problem addressing reverberation without explicit RIR knowledge is called blind dereverberation [17]. A successful class of methods for blind dereverberation fall under the framework of multi-channel linear prediction (MCLP) [16, 11, 10]. These methods require a microphone array and employ linear prediction to extract information about the properties of reverberation. In this paper, we propose to employ MCLP, which is primarily used for single source dereverberation, as a pre-processing step in source separation. We show theoretically and demonstrate numerically that such a scheme not only dereverberates the sources but also helps in the subsequent separation.

In order to make this statement more precise, let us briefly discuss an MCLP formulation. Consider an array of microphones, recording a source in a reverberant room. Recent MCLP formulations [16, 11, 10] select one of the microphones as the reference and use delayed versions of the whole set of observations to linearly predict the reference signal. The residual of prediction, which may also be interpreted as the ‘innovations’ component, forms the dereverberated estimate of the source. In our scenario, we have multiple sources recorded with a microphone array and we would like to separate and dereverberate the sources. For this, suppose we apply MCLP for each possible choice of the reference microphone. If we have $N$  microphones in the array, this results in $N$ hopefully dereverberant but mixed observations. We say ‘hopefully dereverberant’ because the use of MCLP has previously been considered and justified solely for the case of a single source [11]. Further, even if the MCLP outputs thus obtained are dereverberated, in order to achieve separation, we have to know the relation between the recordings. In this paper, we show that, under a mild assumption on the RIRs, MCLP actually suppresses reflections other than the direct source and preserves the phase of the direct source. Therefore, the MCLP preprocessing step essentially reduces the reverberant separation problem to an anechoic one. This greatly simplifies the separation step because we can now use the known array geometry to estimate the direction of arrival and perform separation using this estimate. We argue that this effect is not dependent on the following DOA estimation and separation algorithms. In that sense, the proposed MCLP stage complements currently available source separation algorithms, improving their performance. Consequently, we do not have any reservations about which method to use following the MCLP stage. Nevertheless, for the sake of completeness, we include a brief description of the geometric source separation (GSS) algorithm [18] followed by a post-filter.

Related Work

Literature on dereverberation/source separation with a microphone array is vast. Here, we focus on blind 111We use the term ‘blind’ to refer to methods that do not require to know the RIRs. However, the blind methods we consider may use the known array geometry, which is independent of the characteristics of the environment. methods for dereverberation and source separation for arrays with a known geometry.

For source separation, provided the directions of arrival (DOA) of the sources are known, and reverberation is limited, beamformers like the minimum variance distortionless response (MVDR) beamformer, along with a post-filter [2, 7, 23] can be effective. In MVDR, one aims to minimize the energy of the reconstructed signals, while preserving the signals from the DOAs. However, reverberation causes a source to contribute to different directions, reducing the effectiveness of such approaches. Geometric source separation (GSS) [18] addresses this problem by also enforcing statistical independence of the reconstructed sources. This is achieved by selecting the beamformer weights so as to minimize a cost function that consists of the sum of a term imposing geometric constraints based solely on the array geometry and a term that promotes the statistical independence of the sources (motivated by blind source separation techniques [9]). In practice, GSS outputs are still reverberant for moderate to highly reverberant rooms. GSS has been complemented with a post-filter in [28, 27] to devise a high performance source-separation method that can operate in real time. An alternative approach, based on the assumption that the activity of the different sources do not overlap in the time-frequency domain [31], has been proposed in [24]. However, the model in [24] assumes that the room impulse responses are relatively short, compared to the analysis window size used in the STFT.

Despite its simplicity, blind dereverberation is a challenging problem due mainly to the characteristics of room impulse responses. For a general overview and discussion of relatively recent methods, we refer to the edited book [17]. Here, we focus on a successful class of methods that can be collected under the heading ‘multichannel linear prediction’ (MCLP), as they are more directly related to the content of our work. These methods usually fit in a schema as follows. Given multiple observations $y_{i}$, for $i=1,2,\ldots,M$, one of the observations is chosen as the reference. Then, this reference observation is linearly predicted by the delayed version of all of the observations. Once the linear prediction coefficients are determined, this information has been used differently in prior work. For instance, in [12], the authors use the coefficients to estimate the power of late reverberations, which are subsequently reduced by spectral subtraction. In more recent work [16, 11, 10], the residual of linear prediction is taken as the dereverbed audio signal. Prior work also differentiates based on the criteria used for determining the weights used in linear estimation. In [12], the authors propose to select the linear estimation weights so as to minimize the energy of the residual. In [16], the source signal is modeled as a time-varying Gaussian process and the weights are selected so as to maximize the likelihood of such a model. More recently, [11, 10] employ a sparsity- based model in the time-frequency domain and select the weights so that the residual minimizes a sparsity promoting cost function.

An interesting approach that employs MCLP to achieve suppression of late reverberation and source separation is presented in [26]. The authors use MCLP to suppress late reverberation, using a statistical framework similar to that in [16]. In contrast to [16], a time-varying scenario with multiple sources is considered and the MCLP weights are determined along with the statistics of the late reverberant part of the sources in an iterative manner.

Contribution

Our proposal in this paper is of a complementary nature and is likely to improve the performance of existing separation algorithms. We present an inventive way to employ MCLP as a pre-processing step for blind source separation using a microphone array. We show theoretically that, under certain assumptions, such a pre-processing step converts the reverberant source separation problem into a non-reverberant one, which in turn can be solved by existing methods such as [24, 18, 27] outlined above. Such an approach is inherently different than previous uses of MCLP which aim to dereverb a single source as in [12, 16, 11, 10]. Our proposal also differs from that of [26] because [26] requires the source statistics to cancel the late reverberation and thus resorts to a joint estimation of the linear prediction weights and source statistics. This approach in [26] requires an iterative procedure and is handled with an EM algorithm. In contrast, we propose to employ the MCLP step once, in order to reduce the reverberant problem to an approximately non-reverberant one.

In addition to theoretical justification, we demonstrate the validity of the claims using real recordings. For the sake of completeness, we include a description of GSS along with a simple post-filter, to be applied following the MCLP step, in order to achieve separation. However, after MCLP, other methods that assume short RIRs, or operate under a free propagation model may also be used, as noted above.

Notation

Sequences and continuous functions of time are denoted with a small letter as in $h(n)$ or $h(t)$, respectively. The DFT and STFT of sequences are denoted with a capital letter as $H(\omega)$ and $H(s,\omega)$ respectively – which one is referred to will be clear from the context. The $z$-transform of sequences are denoted with calligraphic letters as in $\mathcal{H}(z)$. We refer to $\mathcal{H}(z)$ as ‘causal’ if it is a polynomial of $z^{-1}$. For a causal $\mathcal{H}(z)$, we denote the constant $h(0)$ also as $\mathcal{H}(\infty)$.

Throughout the paper, $K$ denotes the number of sources and $N$  denotes the number of microphones in the array.

II Description of the Model

Throughout the paper, we are interested in far-field processing. Thus, the distance of the source to the microphones is assumed to be significantly greater than the array size. Our proposed formulation/algorithm is in the STFT domain. However, for simplicity of exposition, we first introduce a model in the continuous-time domain and then present an approximation of this model in the STFT domain.

II-A A Time-Domain Observation Model

We introduce the time-domain model in three stages, where complexity increases with each stage. We start with the case of a single source, single microphone. In the second stage, we consider the single source, multiple microphone case. Finally, we discuss the multiple source, multiple microphone case, which is the main scenario of interest in this paper.

II-A1 Single Source, Single Microphone

For a given microphone-source pair, let $x(t)$ denote the signal produced by the source. The observation $y(t)$ can be approximated as,

[TABLE]

where $\tilde{h}(t)$ is the room impulse response (RIR) for the given source-microphone positions. Now suppose that the distance from the source to the microphone is $D$ meters and the speed of sound in the environment is $c$ m/sec. Then, the acoustic signal emitted by the source reaches the microphone after a delay of $D/c$ seconds. We declare the distance compensated room impulse response (DCIR) to be

[TABLE]

Notice that for practical purposes, we can work with the DCIR instead of the RIR for a single microphone-source setup. For multiple microphones or sources, we need to take into account relative distances, as discussed next.

II-A2 Single Source, Multiple Microphones

Consider a setup consisting of a microphone array and a source as shown in Fig. 1. Suppose we set the center of the array as the origin for a coordinate system. For this origin, let the vector $p_{i}$ denote the position of the the $i^{\text{th}}$ microphone. Also, let $u_{x}$ denote the unit vector in the direction of the source. Finally, let $d_{i}=\langle p_{i},u_{x}\rangle$. Note that $d_{i}$ denotes the coordinate of the projection of the $i^{\text{th}}$ microphone’s position to the subspace spanned by $u_{x}$. The induced time-delay is $\tau_{i}=-d_{i}/c$. Notice that $\tau_{i}$ may be positive or negative, i.e., may represent an advance or delay. Finally, let the DCIR from the source to the $i^{\text{th}}$ microphone be denoted as $h_{i}(t)$. Then, we model the observation at the $i^{\text{th}}$ microphone, namely $y_{i}$, as,

[TABLE]

The far field assumption amounts to the following : If the distance of the source to the center of the array is $D$, then the distance of the source to the $i^{\text{th}}$ microphone is approximately $D-d_{i}$. Under the far-field assumption, the described model preserves the ability of the array to differentiate different directions. We next add another layer to this model by considering multiple sources.

II-A3 Multiple Sources, Multiple Microphones

Suppose now that there are $K$ sources and $N$ microphones. Let the $k^{\text{th}}$ source signal be denoted by $x_{k}(t)$. Also, let the unit vector in the direction of the $k^{\text{th}}$ source be denoted by $u_{k}$. For the $i^{\text{th}}$ microphone, positioned at $p_{i}$, we set $\tau_{i,k}=\langle p_{i},u_{k}\rangle/c$. Finally, let $h_{i,k}(t)$ denote the DCIR between the $k^{\text{th}}$ source and the $i^{\text{th}}$ microphone. Then, we model the observation at the $i^{\text{th}}$ microphone as,

[TABLE]

Notice that, the absolute delay from the $k^{\text{th}}$ source is not included in the model. However, this is not a careless omission, because this delay is in fact not known and is not estimated in practice. Therefore, the described model is suitable for practical far-field processing scenarios.

II-B An Approximate Model in the STFT Domain

We work in the STFT domain for reasons to be briefly clarified in the sequel. We are specifically interested in the representation of an LTI system in the STFT domain. Consider a convolution equation as,

[TABLE]

Let $Y(s,\omega)$, $X(s,\omega)$ denote the STFTs of $y$ and $x$ respectively, obtained using a window $w(n)$. We differentiate between two cases, following the discussion in [13] (see Lemma 1 in [13]).

1. (I)

If the support of $h$ is contained in the support of $w$, then

[TABLE]

where $H(\omega)$ denotes the DFT of $h$. 2. (II)

If the support of $h$ is significantly larger than the support of $w$, then

[TABLE]

for some $\tilde{H}(s,\omega)$, where $\overset{s}{\ast}$ denotes convolution with respect to the time variable $s$, defined as

[TABLE]

We remark that $\tilde{H}(s,\omega)$ is not necessarily the STFT of $h$, but is chosen so as to make the approximation as close as possible – see [20] for a discussion.

We remark that, as a special case of (I), if $\tau$ is small enough, then the STFT of $x(n-\tau)$ can be approximated as $X(s,\omega)\,\exp(-j\omega\,\tau)$. Therefore, the model introduced in (4) can be approximated in the STFT domain as

[TABLE]

Such a model is preferred over the time domain model because prior or penalty functions that enforce (structured) sparsity of the audio signals are typically defined in the STFT (or some other time-frequency) domain. If we use the time domain convolutional model (4) in an iterative algorithm, we have to work in the time domain to undo the effects of reverberation, and the STFT domain to take advantage of the signal prior. Switching from the time domain representation to the STFT representation at each iteration increases the computational cost significantly. Working with a convolutional model directly in the STFT domain, such as (9), allows to carry out the iterations in the STFT domain, eliminating the need to apply the forward and inverse STFT operations at each iteration.

III Multi-Channel Linear Prediction for Dereverberation

In this section, we discuss the implications of using multi-channel linear prediction (MCLP) as a preliminary step for the source separation problem. We derive conditions under which MCLP transforms the original blind separation/dereverberation problem to a non-blind one.

We start the discussion with a single source and then consider extensions to multiple sources.

III-A MCLP for a Single Source

Consider a signal $x$ that is observed via different channels as

[TABLE]

where $h_{i}$ denotes the impulse response of the $i^{\text{th}}$ channel. In the $z$-domain, this can be written as,

[TABLE]

We would like to find causal filters $\mathcal{G}_{1}(z)$, …, $\mathcal{G}_{N}(z)$ such that

[TABLE]

MINT [15] ensures that such filters exist provided that $\mathcal{H}_{i}$’s are co-prime, i.e., they do not have a common non-constant factor. Even though it is relieving to know that such filters exist, finding the filters is still a challenging problem when $h_{i}$’s are not known. MCLP addresses this challenge.

In MCLP, the idea is to estimate $x$ by subtracting from $y_{r}$ (for a chosen reference index $r$) delayed weighted combinations of $y_{i}$’s. This is equivalent to using in (12) a set of filters of the form

[TABLE]

for a positive integer $d$ (determining the amount of delay) and causal $\mathcal{U}_{k}$’s. In [11], the selection of the filters is guided by the assumption that the original $x$ is sparse. More specifically, suppose $y_{i}(s)$ denotes the 1D time-series derived from a frequency band of the observations, that is $y_{i}(s)=Y_{i}(s,\omega)$ for some $\omega$. The filters $u_{i}$ are found by solving the following minimization problem.

[TABLE]

where $P(\cdot)$ is a sparsity promoting function. We will have more to say about $P$, but first we would like to stress a point about the choice (13). MINT ensures the existence of dereverberating filters but does not guarantee that they can be of the form (13). Below, we propose modified conditions which ensure the existence of such $G_{i}$’s. Another issue is the role of the reference index $r$. It turns out that repeating the linear prediction for different $r$ values leads to an interesting outcome, summarized in Prop. 2 below. In order to derive Prop. 2, we need an auxiliary result (Prop. 1) and an assumption.

Proposition 1**.**

Suppose that the $N-1$ polynomials $\mathcal{H}_{i}$, $i\in\{1,\ldots,N\}\setminus\{r\}$, are co-prime. Then, for $d=1$, there exist a set of filters of the form (13) such that

[TABLE]

for a constant $c$. Further, for a set of causal filters of the form (13), if (15) holds, then $c=\mathcal{H}_{r}(\infty)$.

Proof.

See Appendix A. ∎

Consider now the application of this result to the model (9). Assume that there is a single source, i.e., $K=1$. For a fixed frequency $\omega$, let

[TABLE]

Under the hypotheses of Prop. 1, if, for $d=1$, MCLP can recover a set of filters satisfying (15), then Prop. 1 implies that,

[TABLE]

Now, if further,

[TABLE]

then,

[TABLE]

for some constant $c^{\prime}$, independent of $r$ (actually, $c^{\prime}=h_{i}(0)$, for any $i$). We remark that the difference of (17) and (19) is that the rhs of (19) depends only on the direction of arrival, whereas this is not the case in (17). In other words, MCLP eliminates reflections from different directions.

The assumption in (18) states that, around the origin, distance compensated impulse responses are similar. This assumption is key to the development in this paper and we state it explicitly for later reference. We will verify this assumption experimentally in Section V.

Assumption 1**.**

For a given source - microphone array pair, let $\tilde{H}_{i}(s,\omega)$ denote the STFT representation of the distance compensated impulse response for the $i^{\text{th}}$ microphone. Then,

[TABLE]

for $i\neq i^{\prime}$ and all $\omega$.

We state the foregoing observations as a proposition.

Proposition 2**.**

For a fixed frequency $\omega$, let $h_{i}(s)=\tilde{H}_{i,1}(s,\omega)\,\exp(-j\omega\,\tau_{i,1})$. Assume also that the collection of polynomials $\mathcal{H}_{i}$, $i\in\{1,\ldots,N\}\setminus\{r\}$ are co-prime. Then, for $d=1$, there exist a set of filters of the form (13) such that

[TABLE]

for a constant $c_{r}$.

Further, for a set of filters of the form (13), if (21) holds, and Assumption 1 is in effect, then $c_{r}=c^{\prime}\,\exp(-j\omega\tau_{r,1})$ for a constant $c^{\prime}$, independent of $r$. ∎

Next, we extend the discussion to the case of multiple sources.

III-B MCLP for Multiple Sources

Consider an observation model with $K$ sources $x_{1}$,…$x_{K}$ and $N>K$ observations defined as,

[TABLE]

In the $z$-domain, this can be written as,

[TABLE]

For this model, MINT also ensures that we can find filters $\mathcal{G}_{k,i}$, for $k=1,\ldots,K$, $i=1,\ldots,N$ such that

[TABLE]

provided that the (Smith) canonical form of $\mathbf{H}$ (see Appendix B or [8], Sec.VI.2, Defn.3) has a non-zero constant (i.e., no polynomial higher than zero order) diagonal. However, unlike the single source case, we cannot use filters of the form (13) to achieve (24). We have the following extension in this case.

Proposition 3**.**

For a fixed reference microphone index $r\in\{1,\ldots,N\}$, let $\bar{\mathbf{H}}_{r}$ be the matrix obtained by removing the $r^{\text{th}}$ row of $\mathbf{H}$. If the greatest common divisor of all the $K$-minors of $\bar{\mathbf{H}}_{r}$ is a non-zero constant, then, for $d=1$, there exist a set of filters of the form (13) such that

[TABLE]

for a constant vector $c^{(r)}$ (i.e., containing no polynomial of order greater than zero).

Further, for a set of filters of the form (13), if (25) holds, then the components of $c$ satisfy,

[TABLE]

Proof.

See Appendix C. ∎

Remark 1**.**

For an $N\times K$ matrix $\mathbf{H}$ with $K\leq N$, it can be shown that, the greatest common divisor of all the $K$-minors of $\mathbf{H}$ is a non-zero constant if and only if the Smith canonical form of $\mathbf{H}$ has a non-zero constant diagonal. This equivalence is shown in Prop. 8 in Appendix B. Therefore, we can replace the hypothesis in Prop. 3 with an alternative, stated in terms of the Smith canonical form of $\bar{\mathbf{H}}_{r}$.

When we apply this result to our observation model with $K$ sources, we obtain the following proposition.

Proposition 4**.**

For a fixed frequency $\omega$, let

[TABLE]

for $i\in\{1,\ldots,N\}$, $k\in\{1,\ldots,K\}$. Also, let $\mathbf{H}$ be defined as in (23), and $\bar{\mathbf{H}}_{r}$ be defined as in Prop. 3. Finally, let $y_{i}(s)=Y_{i}(s,\omega)$.

If the greatest common divisor of all the $K$-minors of $\bar{\mathbf{H}}_{r}$ is a non-zero constant, then, for $d=1$, there exist a set of filters of the form (13) such that

[TABLE]

for some constants $c^{(r)}_{k}$.

Further, for a set of filters of the form (13), if (28) holds, and Assumption 1 is in effect, then we have

[TABLE]

where $c^{\prime}$ is a constant independent of $r$.

Proof.

See Appendix D. ∎

This proposition suggests that, in a scenario with multiple sources, MCLP, if successful, suppresses reflections and preserves the direct signal even when there exist multiple sources from different directions. Therefore, after the MCLP step, the recordings appear as if they are recorded with the same microphone array, but in a non-reverberant environment. This in turn allows to read off the directions of the sources with a simple procedure like MUSIC [22]. Once the source directions are determined, we can then achieve separation by using known manifold vectors of the microphone array.

IV Implementation Details

The discussion after Prop. 4 sets the stage for the proposed framework. The main building blocks of the framework are shown in Fig. 2. We first input the $N$ observations $y_{1},\ldots,y_{N}$, recorded in a reverberant environment with unknown reverberation, to an MCLP block. This block applies MCLP for every possible choice of the reference microphone. The output of this block is $N$ dereverberated signals $\hat{y}_{1},\ldots,\hat{y}_{N}$. These signals are input to a DOA estimation algorithm (such as MUSIC [22]) and the directions of the $K$ sources are determined. Then these estimated directions are used to separate the sources in $\hat{y}_{i}$ using a non-blind source separation method. The estimated sources are denoted as $\hat{x}_{i}$ in Fig. 2. In the following, we briefly present the algorithm we used for a practical implementation of the components. However, we remark that following MCLP, the problem turns into a non-reverberant source separation problem (with known array geometry) and one can use any method of choice for this task.

IV-A Multi-Channel Linear Prediction

We use the formulation given in Sec. III-B. We work in the STFT domain and process the frequency channels separately. For a fixed value of $\omega$, we set $y_{i}(s)=Y_{i}(s,\omega)$. Then for each value of $r\in\{1,\ldots,N\}$, the dereverbed estimate $\hat{y}_{r}$ is defined as,

[TABLE]

where $\hat{u}^{(r)}_{i}$ denote filters determined by solving

[TABLE]

We use a weighted $\ell_{1}$ norm as the sparsity promoting function $P(\cdot)$. More precisely, for $u\in\mathbb{C}^{n}$, we set

[TABLE]

where $w\in\mathbb{R}_{+}^{n}$ is a given weight vector. We update the weight vector $w$ as the iterations progress, to improve the sparsity promoting property of $P$.

We remark that for MCLP, a similar formulation has been discussed in [11], along with an algorithm based on the alternating direction method of multipliers (ADMM) algorithm [3]. We employ the Douglas-Rachford (DR) algorithm [14, 5] for solving this problem. To that end, we first define the matrices $X^{(i)}$ for $i\in\{1,\ldots,N\}$ such that

[TABLE]

and as $s$  ranges over the domain of the observation $y_{r}$. Notice that $X^{(i)}$ is Toeplitz and has $L$ columns and as many rows as the number of samples in $y_{r}$. By concatenating $X^{(i)}$’s, we construct the data matrix $X$ as,

[TABLE]

Similarly, we construct the variable vector $u$ as

[TABLE]

where $u^{(i)}$ denotes the column vector with $u^{(i)}_{k}=u_{i}(k)$. Also, let $y^{(r)}$ denote the column vector with $y^{(r)}_{k}=y_{r}(k)$. Using these, we can rewrite the problem (31) as,

[TABLE]

We split variables and introduce $z$ in order to make the formulation feasible for employing DR.

For this formulation, the Douglas-Rachford algorithm is of the form outlined in Algorithm 1. We remark that in this algorithm, at each $M$ iterations, the weights $w$ used to compute $\|\cdot\|_{w}$ are updated, so as to improve the sparsity inducing property of the penalty function, following [11]. We refer to [11] for more intricate update strategies.

IV-B DOA Estimation

Because audio is a wideband signal, narrowband DOA estimation techniques like MUSIC [22], or ESPRIT [21] are not directly applicable. However, simple extensions of MUSIC to wideband scenarios exist and we employ one of them. In a nutshell, the idea is to estimate the DOA for each frequency separately and average the estimates. This approach is referred to as the ‘incoherent signal subspace method’ (ISSM) in the recent literature [32] – see [30, 25] for early work. Despite its simplicity, ISSM is noted to perform well in high SNR [32].

Subspace methods for DOA estimation make use of two essential components, namely the covariance matrix of the data and the manifold vector, determined by the array geometry. For a fixed frequency $\omega$, and time instance $s$ (in the STFT domain), we estimate the correlation matrix $R_{Y}(s,\omega)$ as

[TABLE]

where

[TABLE]

In this definition, $\eta\in[0,1]$ is a smoothing (or forgetting) factor justified by the assumption that the sources move slowly.

Let us now define the manifold vector. Recall from Sec. II-A that $p_{i}$ denotes the position of the $i^{\text{th}}$ microphone with respect to the center of the array. Also, let $u_{\theta}$ denote the unit vector in the direction $\theta$. Finally, let $c$ denote the speed of sound, $f_{s}$ denote the sampling frequency and

[TABLE]

denote the relative delay (in samples) of the $i^{\text{th}}$ microphone for a source in the direction $\theta$. The manifold vector is defined as

[TABLE]

Assuming there are $K$ sources, let $C_{K}(s,\omega)$ denote the $(N-K)\times N$ unitary matrix whose columns are the eigenvectors of $R_{Y}(s,\omega)$ corresponding to the smallest $N-K$ eigenvalues. The DOA function is defined as

[TABLE]

For each time instance $s$, the $K$ dominant peaks of $D(s,\theta)$ with respect to $\theta$ are taken to be the directions of the sources.

IV-C Source Separation

For source separation, we use the geometric source separation (GSS) method [18, 27] complemented with a simple post-filter. We assume that the sources are stationary, so that the estimated directions are constant with respect to time. Extension to moving sources is straightforward, at least in principle, and will not be further discussed.

IV-C1 Geometric Source Separation

Let $\theta_{i}$, for $i=1,\ldots,K$ denote the directions of the sources determined by the method described in Sec. IV-B. For each $\omega$, we define a matrix $A(\omega)$ as,

[TABLE]

We model the MCLP outputs as

[TABLE]

This model is not exact because despite the significant dereverberation by MCLP, $\hat{Y}$ still contains reflections other than the direct signal.

For this model, GSS proposes to estimate $X(s,\omega)$ as $\hat{W}(\omega)\,\hat{Y}(s,\omega)$, where the $K\times N$ weight matrix $\hat{W}(\omega)$ is determined by solving a minimization problem. To describe the problem, let us denote

[TABLE]

for a given weight matrix $W$. Also, let $R_{S}(\omega)$ denote the empirical autocorrelation matrix of $S(s,\omega)$, obtained by averaging $S(s,\omega)\,S^{*}(s,\omega)$ over the time variable $s$. Given this notation, $\hat{W}(\omega)$ is defined as the solution of

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

Thanks to the differentiability of both terms in the cost function, gradient descent can be applied for obtaining a stationary point. For a fixed $\omega$, we remark that the gradient [4] of the complex valued $C(W)$ is given as [18],

[TABLE]

where $R_{\hat{Y}}$ denotes the empirical autocorrelation matrix of $\hat{Y}$, obtained by averaging $\hat{Y}^{*}(s,\omega)\,\hat{Y}(s,\omega)$ over the time variable $s$.

IV-C2 Post-Filter

The GSS algorithm aims to achieve source separation by considering the best linear combination of the inputs. In practice, separation performance can be improved by post-filtering. Specifically, we view the separated signals as contaminated with noise and eliminate this noise using a soft-threshold, where the time-frequency varying threshold value depends on an estimate of the remaining noise variance for each time-frequency sample.

More specifically, if $\hat{S}_{1}$, …, $\hat{S}_{K}$ denotes the $K$ source estimates obtained by GSS, for each value of the triplet $(i,s,\omega)$, we regard $\hat{S}_{i}(s,\omega)$ as an unbiased estimate of $X_{i}(s,\omega)$ with variance $\sigma_{i}^{2}(s,\omega)$. The post-filtered estimate is formed as,

[TABLE]

where $\alpha$ is a parameter.

In order to realize such a post-filter, we need to estimate $\sigma_{i}^{2}(s,\omega)$. We estimate $\sigma_{i}^{2}(s,\omega)$ as,

[TABLE]

where $a_{\theta_{i},n}(\omega)$ denotes the $n^{\text{th}}$ entry of $a_{\theta_{i}}(\omega)$. While this expression may make intuitive sense, further justification is provided in Appendix E. We also remark that this is a generalized version of the post-filter employed in [1].

V Experiments

Our primary claim in this paper is that an MCLP step suppresses reflections other than the direct source, under a certain assumption (Assumption 1). The experiments in this section are geared towards demonstrating the claims and verifying Assumption 1. There are three experiments. In the first, we verify Assumption 1 using real impulse responses obtained in a classroom. In the second experiment, we demonstrate that an MCLP preprocessing step eliminates reflections other than the direct source, by showing how the DOA function changes after the MCLP step. Finally, the third experiment validates our claim about the effect of employing an MCLP preprocessing step in a source separation problem.

All of the experiments are performed on real data. In order to obtain the data, we used a circular microphone array with eight omni-directional microphones. The configuration and the definition of DOA, namely $\theta$, are as shown in Fig. 3. The microphones are distributed uniformly on a circle of radius 5 cm.

In the experiments, we use an STFT with a smooth window of length 64 ms, and hop-size is 16 ms. The window is selected such that the STFT is self-inverting (i.e., a Parseval frame). The sampling frequency for the audio signals is 16 KHz. For the MCLP block, the parameters (see Algorithm 1) are chosen as $d=2$, $L=30$, $M=50$, $\alpha=0.05$, $\lambda=0.5$. We run the algorithm for 2000 iterations.

Some of the audio signals used in the experiments, as well as those from experiments not contained in the paper can be listened to at “http://web.itu.edu.tr/ibayram/MCLP/”.

Experiment 1** (Verification of Assumption 1).**

In order to verify Assumption 1, we used estimated impulse responses from different DOAs. The impulse responses are estimated using the method described in [6], using loudspeakers placed at known DOAs in the far field. Knowing the DOAs, we computed the distance compensated impulse responses and obtained filters $\tilde{H}_{i}(s,\omega)$ that approximate their action in the STFT domain. For each $s$ and $\omega$ value, we computed the variance of $\tilde{H}_{i}(s,\omega)$’s with respect to $i$ (recall $i$ takes eight different values), which we denote as $v(s,\omega)$. We normalize $v(s,\omega)$ by the total power at the location $s$, $\omega$, defined by $p(s,\omega)=\sum_{i}|H_{i}(s,\omega)|^{2}$. Thus for each $s$ and $\omega$ value, we obtained a number $\hat{v}(s,\omega)=v(s,\omega)/p(s,\omega)$. We remark that Assumption 1 is valid if and only if $\hat{v}(0,\omega)=0$ for all $\omega$. In Fig. 4, we show (in decibels) the fields $\hat{v}(s,\omega)$ for different DOAs, namely [math], $\pi/4$ and $\pi/2$ radians. Observe that for all DOAs and almost all $\omega$ shown, $\hat{v}(0,\omega)$ (the first column in each image) is significantly lower than $\hat{v}(s,\omega)$ for $s\neq 0$ (the remaining columns). Therefore, Assumption 1 is seen to be valid for practical purposes.

Experiment 2** (Direction of Arrival).**

In this experiment, we demonstrate the effect of an MCLP front-end on estimating the DOA.

The source signal is produced with a loudspeaker located at $\theta=\pi/2$ radians in a reverberant environment. The spectrogram of the signal observed at the first microphone is shown in Fig. 5a. Notice that the harmonics are seriously smeared. For this signal, MCLP significantly removes the effects of reverberation – see Fig. 5b.

In order to visualize the effects of MCLP on DOA estimation, we show the DOA function $D(k,\theta)$ (defined in (41)) averaged over the time variable $k$ in Fig. 6, for both the reverberant observations and MCLP outputs. We remark that averaging is meaningful here, since the source is stationary. For comparison, we also normalize the averaged DOA functions so that their peak is unity. The correct position of the source is also indicated with a dotted line. We see that the DOA function is clearly improved by the application of MCLP. The variance is clearly reduced and the function is closer to an impulse located at the correct position.

Experiment 3** (Source Separation).**

In this experiment, we demonstrate the improvement gained by employing MCLP as a pre-processing step in a blind separation setup. We use the same microphone array as in Exp. 2. We use the source from Exp. 2 as the first source. The second source is produced by a loudspeaker placed at $\theta=1.75\,\pi$ radians (so that the sources are separated by $0.75\,\pi$ radians). The spectrogram of the observation at the first microphone is shown in Fig. 7a. The output of MCLP when the first microphone is taken as the reference is shown in Fig. 7b.

The average (over time) DOA functions obtained from the original recordings and MCLP outputs are shown in Fig. 8. We observe that MCLP has considerably sharpened the functions. Also, without MCLP, the peaks of the DOA functions are offset by some amount from the actual directions, which would lead to a biased estimate of the true DOA.

The source separation results are shown in Fig. 9. Perceptually, without MCLP, especially for the second source, we can clearly hear the first source, despite the post-filter. Further, reverberation is not eliminated. These artifacts are eliminated using MCLP.

In order to quantify the results, we also report PESQ and SIR scores [29]. PESQ gives an objective perceptual distance between an input speech signal and a reference. In order to obtain meaningful distances, we use the multiple channel recordings of source-1 and source-2, prior to mixing, to extract the reference signals. We take the $i^{\text{th}}$ microphone as reference and apply MCLP to the multiple microphone recordings to obtain the references for each source. We then compute the PESQ scores for the methods described in this paper. We finally average over $i$ to obtain a reliable PESQ score. The resulting scores are tabulated in Table I. We remark that the inclusion of an MCLP step (denoted as LP in the table) improves the PESQ scores noticeably.

The SIR scores for separation with different combinations of the discussed methods are also provided in Table I. We observe that including an MCLP preprocessing step also improves significantly the separation quality measured by SIR.

VI Conclusion

In this paper, we advocated the use of MCLP as a pre-processing step for blind reverberant source separation. We provided theoretical justification for such a scheme and demonstrated numerically that it can significantly improve the performance of readily available separation algorithms.

The current algorithm for MCLP is computationally intensive. However, for a practical scenario with slowly moving sources, if we have approximate MCLP weights at a certain time, we think it would be possible to update the weights in a computationally efficient manner by warm-starting the optimization procedure. Therefore, the proposed scheme has potential for real-time applications as well.

Another issue involves the correlation between the weights obtained for different reference microphones. In the current formulation, MCLP is applied independently for each choice of the reference microphone. However, we would expect some correlation between these weights, especially if the microphones are close with respect to the room dimensions, which is usually the case for a microphone array. We plan to explore the implications of such correlation in future work.

Appendix A Proof of Prop. 1

Without loss of generality, suppose $r=1$ in the following.

Since $\mathcal{H}_{2}$, …, $\mathcal{H}_{N}$ are co-prime, we can find causal polynomials $\tilde{\mathcal{U}}_{2}$, …, $\tilde{\mathcal{U}}_{N}$ such that

[TABLE]

Now set, $\mathcal{U}_{1}(z)=1$ and

[TABLE]

Then for $\mathcal{G}$ as in (13), and $d=1$, we have

[TABLE]

For the converse, observe that for causal $\mathcal{U}_{k}$ and $\mathcal{H}_{k}$, the polynomial

[TABLE]

does not have any constant terms, i.e., $C(\infty)=0$. Observe also that $\mathcal{G}_{1}(\infty)=1$. Thus if (15) holds, then

[TABLE]

Appendix B Canonical Form of a Polynomial Matrix

In this appendix, we briefly present, for the sake of completeness, the (Smith) canonical form of a polynomial matrix, along with some of its properties. For a more detailed exposition, we refer to the fine treatment in Sec. VI.1–3 of Gantmacher’s book [8].

We refer to a quantity or matrix as non-zero constant if it is a polynomial of order zero.

Let us start by introducing three elementary row operations on a polynomial matrix $\mathbf{A}(z)$ :

Multiplication of a row with a non-zero constant. 2. 2.

Addition, to the $j^{\text{th}}$ row, the product of the $i^{\text{th}}$ row and a polynomial $p(z)$. 3. 3.

Interchange of two rows.

These operations can be performed by left multiplication with a square matrix whose determinant is a non-zero constant. We refer to matrices that can be used to perform elementary operations as elementary matrices in this appendix. We also remark that similar operations can be defined on the columns instead of rows. Elementary column operations can be performed by right-multiplication with elementary matrices.

We say that two polynomial matrices $\mathbf{A}$ and $\mathbf{A}^{\prime}$ are equivalent if $\mathbf{A}$ can be transformed to $\mathbf{A}^{\prime}$ by elementary row/column operations.

Let us now introduce the (Smith) canonical form.

Definition 1**.**

An $N\times K$ (with $K\leq N$) diagonal polynomial matrix $\mathbf{A}(z)$ of rank $r$ with diagonal $a_{i}(z)$ for $i\in\{1,\ldots,K\}$ is said to be in (Smith) canonical form if

$a_{i}(z)\neq 0$ if $i\leq r$, $a_{i}(z)=0$ if $i>r$, 2. 2.

$a_{i}(0)=1$ for $i\leq r$, 3. 3.

$a_{i}(z)$ divides $a_{i+1}(z)$ for $i<r$.

A fundamental result is that we can transform any polynomial matrix to one in canonical form via elementary row/column operations.

Proposition 5**.**

Each $N\times K$ polynomial matrix is equivalent to a unique matrix in canonical form.

In general, the canonical form is reached by performing both row and column operations. However, we are interested in row operations only, for the development in this paper. We have the following result.

Proposition 6**.**

Suppose an $N\times K$ matrix $A$ has rank $K$ and the diagonal of its canonical form is constant. Then, $A$ can be reduced to its canonical form via row operations.

This result is a corollary of the discussion in [8] but is not stated explicitly. We provide a proof for the sake of completeness. We need a few intermediate results for the proof. All of the non-obvious claims stated before the proof of Prop. 6 can be found in [8].

Given an $N\times K$ polynomial matrix $A(z)$, let $D_{r}(z)$ denote the greatest common divisor of all minors of $A$ of order $r$, normalized such that the leading coefficient is unity. Notice that $D_{r}$ is divisible by $D_{r-1}$. Therefore, for $r>1$, the ratios $i_{r}=D_{r}/D_{r-1}$ are polynomials. We also set $i_{1}=D_{1}$. It can be shown that $D_{r}$’s are invariant under elementary row/column operations. This property transfers to $i_{r}$ as a consequence of its definition.

Suppose now that the canonical form of $A(z)$ consists of a diagonal matrix with elements $a_{r}(z)$ as in Defn. 1. It turns out that $a_{r}=i_{r}$, giving us an alternative definition of $a_{r}$.

We need a final result before proving Prop. 6, which can be shown by a simple elimination argument.

Proposition 7**.**

Using elementary row operations, any polynomial matrix $B(z)$ can be reduced to an upper triangular matrix $U(z)$ such that, at each column of $U(z)$, if the polynomial on the diagonal is non-zero, then the degree of any off-diagonal polynomial is strictly smaller than the degree of the polynomial on the diagonal.

We are now ready for the proof of Prop. 6.

Proof of Prop. 6.

Suppose the canonical form of $A$ is $B$. Observe that the only non-zero minor of $B$ of order $K$ is given by the product of its diagonals, which is unity. Therefore, the invariant polynomial $D_{K}$ of $A$ is a constant.

Suppose we reduce $A$ to an upper triangular form $U$ as in Prop. 7 using elementary row operations only. Observe that the only non-zero minor of $U$ of order $K$ is given by the product of its diagonals. Let us denote this non-zero minor by $U_{K}(z)$. But since this is the only non-zero minor of order $K$, and elementary matrices have constant determinants, it follows that we must have $U_{K}(z)=c\,D_{K}$ for some constant $c$. Thus, $U_{K}(z)$ is a constant. Therefore, $U$ is a constant diagonal matrix with a non-zero diagonal. Each constant on the diagonal can be made unity by an elementary row operation, yielding the canonical form $B$. ∎

In a final proposition, we show the equivalence stated in Remark 1.

Proposition 8**.**

Suppose $A$  is an $N\times K$  polynomial matrix with $K\leq N$. The greatest common divisor of all the $K$-minors of $A$ is a non-zero constant if and only if the Smith canonical form of $A$ has a non-zero constant diagonal.

Proof.

As in the discussion above, let $D_{r}(z)$ denote the greatest common divisor of all minors of $A$ of order $r$, where the leading coefficient is unity, $i_{r}=D_{r}/D_{r-1}$ for $r>1$, with $i_{1}=D_{1}$. Recall that $i_{r}$’s form the diagonal in the Smith canonical form of $A(z)$.

Suppose now that $D_{K}(z)=1$. This implies that $i_{r}=1$ for all $r\leq K$. Therefore the diagonal of the canonical form of $A$  is the identity matrix.

For the converse, suppose the canonical form of $A$ has a constant diagonal. This means that $i_{r}(z)$’s are constant for $r\leq K$. But this implies that $D_{r}(z)=1$ for $r\leq K$. ∎

Appendix C Proof of Prop. 3

Without loss of generality, suppose $r=1$ in the following.

By Remark 1, which is proved in Prop. 8, we first note that the gcd of the $K$-minors of $\bar{\mathbf{H}}_{1}$ is a non-zero constant if and only if the canonical form of $\bar{\mathbf{H}}_{1}$ has a non-zero constant diagonal. Thus, by Prop. 6, we can find a $K\times(N-1)$ polynomial matrix $\mathbf{V}$ such that

[TABLE]

For $k=1,\ldots,K$, let us define the polynomials

[TABLE]

and set

[TABLE]

Also, set $\mathcal{U}_{1}(z)=0$. Then for $\mathcal{G}$ defined as in (13), we have

[TABLE]

For the converse, observe that, for causal $\mathcal{U}_{k}$ and $\mathcal{H}_{i,k}$, the polynomial

[TABLE]

does not have any constant terms, regardless of $k$, i.e., $C_{k}(\infty)=0$ for all $k$. Observe also that $\mathcal{G}_{1}(\infty)=1$. Thus if (25) holds, then,

[TABLE]

Appendix D Proof of Prop. 4

We remark that (28) follows by expressing (25) in the time-domain.

Let us consider (29). If (25) holds, and $\mathcal{G}_{i}$’s satisfy (13), then by (26), we have

[TABLE]

By Assumption 1 and (27), it follows that

[TABLE]

where $c^{\prime}$  is independent of $r$  and $k$. Using (22), and the properties of convolution, we finally obtain

[TABLE]

Appendix E Motivation for the Post-Filter in Sec. IV-C2

In order to justify the variance expression in (48), we will make use of Prop. 9 below. For simplicity, this proposition uses a seemingly independent notation. We will identify the elements in the proposition with our separation scenario right after the proposition. Recall that a complex valued Gaussian random variable $Z$ is said to be circular [19] and have variance $2\,\sigma^{2}$  if its real and imaginary parts are independent and Gaussian random variables with variance $\sigma^{2}$.

Proposition 9**.**

Suppose $\gamma$ is an unknown constant and we observe $Y_{n}=c_{n}\gamma+u_{n}\,Z_{n}$, for $n=1,\ldots,N$, where $|c_{n}|=1$, $u_{n}\in\mathbb{C}$ and $Z_{n}$’s are circularly symmetric, independent, unit variance, zero-mean Gaussian random variables.

The uniformly minimum variance unbiased (UMVU) estimator for $\gamma$ is of the form

[TABLE]

where

[TABLE]

with

[TABLE]

We also have that,

[TABLE]

is an unbiased estimator for the variance of $\hat{\gamma}$.

Proof (Sketch).

It can be shown, by taking into account the Gaussian pdf of $Y_{n}$, that a complete sufficient statistic for $\gamma$ is

[TABLE]

Observing that $\mathbb{E}(U)=\beta^{-1}\,\gamma$, it follows by the Rao-Blackwell theorem that the UMVUE is as given in (68) for $w_{n}$ as in (69) and $\beta$ as in (70). It can also be checked that the variance of $\hat{\gamma}$  is,

[TABLE]

Let us now set \nu=\sum_{n=1}^{N}|w_{n}|\,\Bigl{|}\,\bigl{(}Y_{n}-c_{n}\hat{\gamma}\bigr{)}\Bigr{|}^{2}. Then,

[TABLE]

∎

In order to relate this proposition with the discussion in Section IV-C2, we make the following identifications : ${Y_{n}\to\hat{Y}_{n}(s,\omega)}$; ${c_{n}\to a_{\theta,n}}$; ${w_{n}\to W_{i,n}}$; ${\hat{\gamma}\to\hat{S}_{i}(s,\omega)}$.

For the $i^{\text{th}}$ source, if we regard all of the remaining sources as noise terms with a Gaussian distribution, then the optimal estimate of the $i^{\text{th}}$ source is obtained by a linear combination of the observations. We assume that GSS approximately finds these optimal weights. In that case, the expression in (48) gives an estimate of the variance of these unwanted components (or, ‘noise’) that remain after GSS.