Antoine Deleforge

Yann Traonmilin

PHASE UNMIXING : MULTICHANNEL SOURCE SEPARATION WITH MAGNITUDE CONSTRAINTS

Keywords: Informed Source Separation, Phase Retrieval, Semidefinite Programming

We consider the problem of estimating the phases of K mixed complex signals from a multichannel observation, when the mixing matrix and signal magnitudes are known. This problem can be cast as a non-convex quadratically constrained quadratic program which is known to be NP-hard in general. We propose three approaches to tackle it: a heuristic method, an alternate minimization method, and a convex relaxation into a semi-definite program. The last two approaches are showed to outperform the oracle multichannel Wiener filter in under-determined informed source separation tasks, using simulated and speech signals. The convex relaxation approach yields best results, including the potential for exact source separation in under-determined settings.

INTRODUCTION

Let M sensors record K complex signals through linear instantaneous mixing. The noisy observation y ∈ C M is expressed as:

y = As 0 + n, (1)
where s 0 ∈ C K is the source vector, n ∈ C M is the noise vector and A ∈ C M ×K is the mixing matrix. This model is very common in signal processing and occurs, for instance, when looking at a single time-frequency bin of the discrete short-time Fourier domain. In that case, each entry a m,k of A may be viewed as a frequency-dependent gain and phase offset from source k to sensor m. The classical problem of multichannel source separation consists in estimating the source signals s 0 given one or several observations y. Even without noise, this problem is fundamentally ill-posed. This is always true in the blind case where A is unknown, but also when A is perfectly known and full-rank, as long as M < K (under-determined setting). Due to these ambiguities, source separation under various prior knowledge on A, s 0 or n is a long-standing and still active research topic. Often, specific structures are imposed on A based on physical [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF][START_REF] Kounades-Bastian | A variational em algorithm for the separation of moving sound sources[END_REF] or learned [START_REF] Deleforge | Variational EM for binaural sound-source separation and localization[END_REF] models of signal propagation. It is also quite common to add statistical assumptions on source and noise signals. For instance, signals may be assumed pairwise statistically independent such as in independent component analysis (ICA [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]), or non-Gaussian and non-white such as in TRINI-CON [START_REF] Buchner | Trinicon: A versatile framework for multichannel blind signal processing[END_REF]. If signals are assumed wide-sense stationary, their respective observed (image) covariance can be estimated and used to compute the well-known multichannel Wiener filter (e.g. [START_REF] Cohen | Speech processing in modern communication: challenges and perspectives[END_REF]), which is then the optimal linear filter in the least squared error sense. Although very powerful, the stationarity assumption is unrealistic for many signals of interest such as speech or music in audio. For this reason, Wiener filtering is often used in combination with time-frequency varying source variances and a Gaussian [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF][START_REF] Vincent | Underdetermined instantaneous audio source separation via local gaussian modeling[END_REF] or alpha-stable [START_REF] Liutkus | Generalized wiener filtering with fractional power spectrograms[END_REF] signal model. To avoid over-parameterization, variances are either assumed to be known such as in informed source separation [START_REF] Ozerov | Informed source separation: source coding meets source separation[END_REF][START_REF] Rohlfing | Nmf-based informed source separation[END_REF], or to be provided by a low-dimensional model such as nonnegative matrix factorization (NMF) [START_REF] Kounades-Bastian | A variational em algorithm for the separation of moving sound sources[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] or more recently deep neural networks [START_REF] Aditya | Multichannel music separation with deep neural networks[END_REF]. A common property of all Wiener-filter-based approaches is that they rely on a good model of source variances, while the source phases are left unconstrained and estimated from observations. The oracle Wiener-filter corresponds to the case where instantaneous source variances and mixing matrices are known.

In this paper, we introduce a slight shift of this paradigm by replacing the prior knowledge on instantaneous source variances by a prior knowledge on instantaneous source magnitudes. Note that both quantities are related, the former being the maximum-likelihood estimate of the latter for non-stationary Gaussian signals. We refer to this problem as phase unmixing and focus here on the oracle case where the source magnitudes and the mixing matrix A are exactly known, while only source phases need to be estimated. Applications are hence informed source separation or situations where good magnitude and mixing models are available. We show that multichannel phase unmixing can be expressed as a non-convex quadratically constrained quadratic minimization problem, which is known to be NP-hard in general and hard to solve in practice. We propose three different approaches to tackle it: a heuristic approach, an alternated minimization approach, and a convex relaxation of the problem into a semi-definite program (SDP). Some preliminary theoretical insights and a detailed experimental study on simulated data are presented to compare the proposed methods to the oracle multichannel Wiener filter. A task of informed multichannel speech source separation is also performed. The proposed convex scheme yields particularly encouraging results, including stability to noise and the potential for exact source separation in under-determined settings. Related work. The considered problem of phase unmixing is related but not to be confused with the problem of phase retrieval, which has triggered considerable research interest over the past 30 years [START_REF] James R Fienup | Phase retrieval algorithms: a comparison[END_REF][START_REF] Yonina C Eldar | recent advances in phase retrieval[END_REF] and has recently regained momentum thanks to novel methodologies [START_REF] Emmanuel | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF][START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF]. In phase retrieval, only the magnitudes of y are observed while s 0 is completely unknown. This problem occurs in applications such as adaptive optics [START_REF] Robert | Phase retrieval and diversity in adaptive optics[END_REF] or X-ray crystallography [START_REF] Harrison | Phase problem in crystallography[END_REF], where the phases of the Fourier transform are intrinsically lost during measurement. Phase retrieval is also sometimes used to find phase estimates which are consistent [START_REF] Le | Consistent wiener filtering for audio source separation[END_REF] with magnitudes estimates of a single-channel signal of interest. This consistency may be imposed by the properties of the shorttime Fourier transform [START_REF] Le | Consistent wiener filtering for audio source separation[END_REF][START_REF] Jaganathan | Stft phase retrieval: Uniqueness guarantees and recovery algorithms[END_REF] or sinusoidal source models [START_REF] Magron | Phase reconstruction of spectrograms with linear unwrapping: application to audio signal restoration[END_REF]. In contrast, the proposed framework solely relies on complex multichannel observations to perform recovery and does not require any structure on phases.

PHASE UNMIXING

Under model (1), we consider the problem of estimating the phases of s 0 ∈ C K given strictly positive magnitudes b = |s 0 |, the multichannel observation y ∈ C M and the mixing matrix A ∈ C M ×K , which we assume full rank. A natural approach is to minimize the Euclidean norm of the residual:

ŝ = argmin s As -y 2 2 s.t. |s k | 2 = b 2 k , k = 1 . . . K. (ΦLS)
We refer to this problem as phase least-squares (ΦLS), because without constraints, it becomes a standard least-squares problem. Least squares has infinitely many solutions in the under-determined case (K > M ) and a unique closedform solution ŝLS = A † y otherwise, where {•} † denotes the Moore-Penrose pseudo-inverse. On the other hand, (ΦLS) is an instance of quadratically constrained quadratic program (QCQP). These problems are non-convex, and solving them or even finding whether they have a solution is NP-hard in general [START_REF] Audet | A branch and cut algorithm for nonconvex quadratically constrained quadratic programming[END_REF]. While branch-and-bound methods exist to solve non-convex QCQPs [START_REF] Audet | A branch and cut algorithm for nonconvex quadratically constrained quadratic programming[END_REF][START_REF] Nikolaos | Baron: A general purpose global optimization software package[END_REF], they are extremely slow in practice 1 . A generally exact and efficient solution to (ΦLS) is thus most likely out-of-reach, but we propose in the following three practical approaches to tackle it.

Normalized multichannel Wiener filter. The multichannel Wiener filter (MWF) is one of the most widely used methods in signal processing [START_REF] Cohen | Speech processing in modern communication: challenges and perspectives[END_REF]. One interpretation of MWF is that it is the maximum a posteriori estimator of s 0 given y, assuming that source and noise signals are zero-mean complex circular-symmetric Gaussian (e.g. [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]). For an i.i.d. noise with variance σ 2 n and independent sources with vari-

ances b 2 1 , . . . , b 2 K , the MWF estimate ŝMWF is 2 ŝMWF = argmin s 1 σ 2 n As -y 2 2 + K k=1 |s k | 2 /b 2 k (2) = (σ 2 n Diag{b} -2 + A H A) -1 A H y (3)
where {•} H denotes Hermitian transposition. While (ΦLS) assumes that source magnitudes are known, MWF assumes that source variances are known, which is the same quantity of prior information. Note that adding the magnitude constraint

|s k | 2 = b 2 k to (2)

recovers (ΦLS).

A simple heuristic approach to (ΦLS) is thus to normalize ŝMWF by changing its magnitudes to b while keeping its phases unchanged. We refer to this as normalized multichannel Wiener filtering (NMWF). Alternated minimization. A second approach to solve (ΦLS) is by alternated minimization w.r.t. each coordinate s i until convergence, i.e., by coordinate descent. The Lagrangian of (ΦLS) writes:

L(s, λ) = As -y 2 2 + K k=1 λ k (|s k | 2 -b 2 k ). (4
) Finding the zeros of the derivatives of (4) w.r.t. to the real and imaginary parts of s i and λ i yields

s i = b i y -A :,i c s i c , a i /| y -A :,i c s i c , a i |, (5)
where a i is the i-th column of A, s i c is s deprived of its i-th element and A :,i c is A deprived of its i-th column. Note that ( 5) is almost surely well-defined (see 3 for detailed derivations). Given an initial guess s (0) ∈ C K , repeatedly applying (5) for i = 1 . . . K until convergence yields the phase unmixing by alternated minimization (PhUnAlt) method, described in Alg. 1. PhUnAlt converges because the nonnegative residual error r (p) decreases at each iteration. Since (ΦLS) is not convex, it does not generally converge to a global minimum but to a local minimum which depends on the initial guess.

Lifting scheme. We first note the following identity:

As -y 2 2 = A -y s 1 2 2 = Ax 2 2 = x H A H Ax = trace A H Axx H = trace (CX), (6) where
A = [A, -y], x = [s, 1] ∈ C K+1 , X = xx H ∈ C (K+1) 2 and C = [A, -y] H [A, -y] ∈ C (K+1) 2
. We can now consider the following convex relaxation of (ΦLS):

X = argmin X trace (CX) s.t. diag{X} = b, X 0 (PhUnLift)
where b = [b .2 , 1] and PhUnLift stands for phase unmixing by lifting. Note that the rank-1 constraint on X has been removed. Using [START_REF] Cohen | Speech processing in modern communication: challenges and perspectives[END_REF] and observing that the diagonal of X contains the squared magnitudes of s, it is easy to see that if X = xx H is a rank-1 solution of (PhUnLift), then ŝ = x1:K /x K+1 = X 1:K,K+1 is a global solution of (ΦLS). Hence, the NP-hard quadratic problem (ΦLS) has been relaxed to a simple convex, linear semi-definite program (SDP). SDPs have been extensively studied over the past decades, and many methods are available to solve them efficiently [START_REF] Wen | Alternating direction augmented lagrangian methods for semidefinite programming[END_REF][START_REF] Wen | Block coordinate descent methods for semidefinite programming[END_REF]. Following [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF], we propose to use the particularly inexpensive block-coordinate descent (BCD) method of [START_REF] Wen | Block coordinate descent methods for semidefinite programming[END_REF]. Algorithm 2 shows the method adapted to (PhUnLift). Each iteration only requires two matrixvector multiplications, which makes iterations of PhUnAlt and PhUnLift of comparable computational complexity.

Algorithm 1 PhUnAlt

Input:

y ∈ C M , A ∈ C M ×K , b ∈ R K + , s (0) ∈ C K . Output: Source estimate ŝ with |ŝ| = b.
1: p := 0; r (0) := +∞; 2: repeat 3: p := p + 1; r (p) := y -As (p) 2 2 ; // Residual error 8: until (r (p-1) -r (p) )/r (p) < 10 -3 9: return s (p) The problem is that solutions of (PhUnLift) are not necessarily rank-1. To understand why this relaxation may still be a viable approach to phase unmixing, let us first draw a connection to related works. Without loss of generality, b can be set to 1 by changing C to Diag{ b}CDiag{ b}. Then, (PhUnLift) has the same form as the SDP PhaseCut, recently introduced in [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF] for phase retrieval. The name PhaseCut was chosen because the real rank-1 counterpart of the problem is known to be equivalent to the classical NPhard MaxCut problem in graph theory. An SDP relaxation of MaxCut was proposed 20 years ago [START_REF] Michel | Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming[END_REF], and since then many extensions have been developed. For the phase retrieval application, it was showed through a connection to the method PhaseLift [START_REF] Emmanuel | Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming[END_REF] that PhaseCut does yield a rank-1 solution with high-probability when some stringent conditions on C are verified [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF]. In the presence of noise, solutions are no longer rank-1 but some stability results are available. Unfortunately, none of these results can be directly transposed to the phase unmixing problem. Indeed, the PhaseCut/PhaseLift equivalence only occurs when C has a specific form involving an orthogonal projection matrix [START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF], which is not the case for phase unmixing. We provide here a first stability theorem for PhUnLift in the determined case only (K ≤ M ). A proof of this theorem is available in the supplementary material 3 : Theorem 1 Let y = As 0 + n, b = |s 0 |, A be full-rank and K ≤ M . Let ŝ be the output of Algorithm 2. We have:

for i = 1 → K do 4: s (p) i := b i y -A :,i c s (p) i c , a i /| y -A :,i c s (p) i c , a i |;
ŝ -s 0 2 ≤ 2 √ 2 σ min (A) n 2 (7)
where σ min (A) is the smallest singular value of A.

In other words, PhUnLift recovers the true source vector up to an error proportional to the noise level. In particular, a rank-1 solution and exact recovery is obtained in the noiseless case. Note that for K ≤ M , bounds similar to [START_REF] Vincent | Underdetermined instantaneous audio source separation via local gaussian modeling[END_REF] can also be obtained for the least-squares, MWF and NMWF estimates.

For the more interesting under-determined case, a different approach is needed because then σ min (A) = 0. A theoretical extension of theorem 1 to K > M and additional properties on A is likely intricate to obtain, although numerical results of section 3 do suggest that this is possible. [START_REF] Wen | Block coordinate descent methods for semidefinite programming[END_REF] (in fact, ν = 0 worked well in practice). Output: Source estimate ŝ with |ŝ| = b.

Algorithm 2 PhUnLift (Block-coordinate descent) Input: y ∈ C M , A ∈ C M ×K , b ∈ R K + , ν > 0, typically small
1: C = [A, -y] H [A, -y]; 2: p := 0; r (0) := +∞; X (0) := I K+1 ; 3: repeat 4:

for i = 1 → K do 5: z := X (p) i c ,i c C i c ,i ; γ := z H C i c ,i ; X (p+1) i c ,i := 6: X (p+1)H i,i c := -bi-ν γ z for γ > 0, 0 otherwise; 7:
end for 8:

p := p + 1; r (p) := trace (CX (p) ); // Residual error 9: until (r (p-1) -r (p) )/r (p) < 10 -3 10: return ŝ with the same phases as X n used to generate the noise in all experiments (interestingly, this is not needed by PhUnAlt or PhUnLift). Three initializations are considered for PhUnAlt: random phases (PhUnAlt), the output of NMWF (NMWF+) or the output of PhUnLift (PhUnLift+). Moreover, a brute-force approach (PhUnAlt*5) is considered, which picks the PhUnAlt estimate with smallest residual out of 5 randomly initialized runs. Phase unmixing problems are generated by randomly picking all the elements of A, s 0 and n from i.i.d. zero-mean complex Gaussian distributions of respecive standard deviations σ A , σ s and σ n . In each experiment, σ A and σ s are randomly uniformly picked in [0, 2] while σ n is adjusted to the desired signal-to-noise-ratio SNR = As 0 2 2 /(M σ 2 n ). For each considered combination of (M, K, SNR), all methods are ran on 1000 random tests. as a function of M when K = M (determined case), under low-noise conditions (SNR = 60dB). As predicted by Theorem 1, near-exact recovery is possible with PhUnLift, and MWF and NMWF yield similarly low errors. PhUnLift+ does not improve over PhUnLift suggesting that the global minimum of (ΦLS) is already reached, while NMWF+ leads to the same solution as PhUnLift. PhUnAlt and PhUnAlt*5 perform relatively less well due to local-minimum convergence. Fig. 1(b)-(c) shows the same experiment in under-determined settings. PhUnLift seems to yield solutions sufficiently close to the global optimum to achieve again near-exact reconstruction with PhUnLift+ for M sufficiently large. This does not seem to be the case with other methods. This is further illustrated in Fig. 1(d values of (M, K). Here, exact means a relative error lower than 10 -8 . 100% exact recovery seems possible with PhUn-Lift+ in a number of under-determined cases where PhUnAlt only achieves around 80%. Fig. 1(f) illustrates that the error of PhUnLift is proportional the noise error when K ≤ M , as predicted by Theorem 1. This is also true for MWF, NMWF and NMWF+, but not for PhUnAlt due to local-minima. In the under-determined setting showed in Fig. 1(g), stability to noise is less obvious. PhUnLift and PhUnLift+ perform best, closely followed by PhUnAlt*5. In general, the fair results obtained with PhUnAlt*5 suggests that the number of localminima is often not too high, making multiple initialization of PhUnAlt a feasible approach. From a computational point-ofview, the non-iterative methods MWF and NMWF are much faster but also perform less well. The computational times of other methods depend on the number of iterations needed for convergence. PhUnAlt generally converges in a few hundred iterations. The same is observed for PhUnLift, except in under-determined cases with high SNRs (≥ 30dB), where tens of thousands of iterations are often needed. This calls for using an SDP solver with faster convergence rate.

We finally conduct an informed speech separation task us- For each experiment, both gains and delays are uniformly picked at random such that mixing matrices remain full rank. To save computational time, when a source has its magnitude lower than -40dB in a given time-frequency bin, the source is ignored and assigned a random phase by all methods. Mean signal-to-distortion-ratios (SDRs) calculated with [START_REF] Févotte | Bss eval toolbox user guide-revision 2[END_REF] for each considered method are showed in table 1 (Rand means random phases with correct magnitudes).

PhUnLift and PhUnLift+ outperform the other methods in under-determined settings, while MWF, NMWF, NMWF+, PhUnLift and PhUnLift+ performs similarly for K = M .

CONCLUSION

The problem of oracle phase unmixing, i.e., multichannel source separation with known magnitudes and mixing matrix, was introduced and cast as a non-convex quadratic problem. Three approaches were proposed to tackle it, including a lifting scheme which showed best performance in practice. The proposed methods outperformed the oracle multichannel Wiener filter in under-determined settings. More theoretical investigations are needed to understand why the lifting scheme works well in the latter case. Moreover, the genericity of proposed methods calls for a number of extensions to more realistic scenarios, e.g., estimation of A, more flexible magnitude constraints and additional phase structure.

  ,K+1 and |ŝ| = b 3. EXPERIMENTS AND RESULTS We now compare the efficiency of MWF (3), NMWF, PhUnAlt (Alg.1) and PhUnLift (Alg. 2) on the task of estimating the phases of s 0 given an M -channel mixture y = As 0 + n. All these methods are compared in the oracle setting: the true magnitudes b = |s 0 | and the true mixing matrix A are provided. Moreover, Wiener-filter-based methods are given the true variance σ 2

Fig. 1 (

 1 Fig. 1(a) shows the mean relative error ŝs 0 2 2 / s 0 2 2

  ) and Fig. 1(e), which show the probability of exact reconstruction of PhUnLift and PhUnAlt for different

Fig. 1 .

 1 Fig. 1. (a)-(c) : Mean relative error for a fixed SNR of 60dB. (d)-(e) Probability of exact reconstruction with PhUnLift+ and PhUnAlt (noiseless). (f)-(g): Robustness to noise in determined cases and under-determined cases (M = 2 . . . 10).

ing random 1

 1 second utterances from the TIMIT dataset [27]. The clean speech signals are sub-sampled at 16 kHz and transformed to the short-time Fourier (STF) domain using a 64 ms sliding window with 50% overlap, yielding F = 512 positivefrequency bins and T = 33 time bins. They are then mixed using a global gain g(m, k) ∈ [-5dB, +5dB] and a discrete time-domain delay τ (m, k) ∈ [0, 50] in samples from each source k to each microphone m. The corresponding mixing matrices used in the STF domain are defined by A f (m, k) = 10 g(m,k)/20 exp(jτ (m, k)f /F ) were f = 0 . . . F -1 is the frequency index.

Table 1 .

 1 Mean

	M, K →	2, 2	2, 3	2, 4	4, 4	4, 5	4, 6
	Input	0.49	-2.70 -4.66 -4.51 -5.76 -7.27
	Rand	-6.44 -6.28 -4.99 -4.41 -4.46 -5.28
	MWF	59.6	21.7	17.0	58.6	27.9	25.0
	NMWF	59.9	21.6	16.9	59.4	28.6	25.9
	NMWF+	59.9	22.8	19.2	59.7	34.1	31.6
	PhUnAlt	22.7	17.6	15.3	25.8	21.5	22.0
	PhUnAlt*5	43.5	35.4	21.5	47.7	37.3	33.6
	PhUnLift	59.9	37.6	22.3	59.0	57.3	40.9
	PhUnLift+	59.9	39.6	21.2	58.0	59.3	44.8

SDR (dB) for 1-second M -channel mixtures of K speech sources. Means are over the K sources for each mixture.

Solving one instance of (ΦLS) using the Matlab version of BARON[START_REF] Nikolaos | Baron: A general purpose global optimization software package[END_REF] on a regular laptop takes over a minute with M =

2, K = 3.[START_REF] Kounades-Bastian | A variational em algorithm for the separation of moving sound sources[END_REF] In the under-determined case, the expression (3) is equivalently replaced by Diag{b} 2 A H (ADiag{b} 2 A H + σ 2 n I M ) -1 y for numerical stability.

http://people.irisa.fr/Antoine.Deleforge/ supplementary_material_ICASSP17.pdf.