Changchun Wenmeng Xiong
email: wen-meng.xiong@bjut.edu.cn

Maoshen JIA, Senior Member, IEEE Jing Zhou

José Picheral
email: jose.picheral@centralesupelec.fr

Joint DOA Estimation and Dereverberation Based on Multi-channel Linear Prediction Filtering and Azimuth Sparsity

Keywords: DOA estimation, MCLP, reverberation, sparsity, PALM

Source localization in reverberant environments has been a prominent research topic in the past two decades. In this paper, instead of the commonly employed time-frequency (TF) bin based methods which rely on empirically selected threshold values, we leverage the microphone array signal model comprising an early reverberant component and a late reverberant component, to propose a novel method for the source localization problem in reverberant environments. Our proposed criterion involves the joint removal of the late reverberant component using the multi-channel linear prediction (MCLP) filter, while estimating the directions of arrival (DOAs) of the actual sources using the early component signals. By applying the azimuth sparsity constraint, the true DOA can be estimated with high resolution and free from the interference of the early reflections. To solve the proposed criterion, DOAs, source signals, and MCLP filter coefficients are estimated by alternative iterations. Additionally, we present a source localization criterion specifically designed for the single source scenario as a special case of the multiple sources scenario. Finally, a source number estimation method and a postprocessing procedure are discussed for searching the global solutions to our proposed criteria. Evaluations with both simulated and realistic data demonstrate the advantages of our proposed methods over the baseline methods.

I. INTRODUCTION

Multiple sources localization problem is of great importance for numerous applications such as speech enhancement, source separation, and acoustic imaging. The conventional source localization methods, including the steered response power with phase transform (SRP-PHAT)/beamforming [START_REF] Do | Srp-phat methods of locating simultaneous multiple talkers using a frame of microphone array data[END_REF], the multiple signal classification (MUSIC) [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], and the maximum likelihood [START_REF] Johnson | Music, g-music, and maximum-likelihood performance breakdown[END_REF], have been widely investigated and proven efficient in real applications. However, the performance of these methods degrades severely in environments with high noise and reverberation.

To address the challenge of reverberation, many methods have been developed based on the multi-channel linear prediction (MCLP) filter in time domain or in short time Fourier transform (STFT) domain. The STFT-based approach is adequate for wideband signals due to its parallel implementation at each frequency bin [START_REF] Yoshioka | Integrated speech enhancement method using noise suppression and dereverberation[END_REF] [START_REF] Yoshioka | Generalization of multi-channel linear prediction methods for blind mimo impulse response shortening[END_REF]. Early signals can be obtained by removing the delayed microphone signal components of earlier frame indices filtered with the MCLP filter, where the coefficients of filters are to be estimated on-line [START_REF] Braun | Online dereverberation for dynamic scenarios using a kalman filter with an autoregressive model[END_REF] [START_REF]Linear prediction-based online dereverberation and noise reduction using alternating kalman filters[END_REF][8] or off-line [START_REF] Yoshioka | Integrated speech enhancement method using noise suppression and dereverberation[END_REF] [START_REF] Yoshioka | Generalization of multi-channel linear prediction methods for blind mimo impulse response shortening[END_REF]. The early signals can be obtained for one single microphone or all the microphones in the array [START_REF] Yoshioka | Integrated speech enhancement method using noise suppression and dereverberation[END_REF]. In [START_REF] Yoshioka | Generalization of multi-channel linear prediction methods for blind mimo impulse response shortening[END_REF], the generalized weight prediction error (GWPE) method is proposed, modelling the early signals as time-dependent complex Gaussian and solving a cost function based on the correlation between the early signals iteratively. Various efforts have been made to improve the performance of the MCLP filter. Instead of Gaussian distribution assumption, a sharper distribution such as Lapacian [START_REF] Jukić | Speech dereverberation using weighted prediction error with laplacian model of the desired signal[END_REF] or a distribution with a heavier tail such as Student's t-dsitribution [START_REF] Chetupalli | Late reverberation cancellation using bayesian estimation of multi-channel linear predictors and student's t-source prior[END_REF] have been found more appropriate for describing the statical characteristics of the information bearing signals including speech and music signals. The order of the prediction filter is assumed to be known a-priori. However, a too large order leads to over-estimate the reverberation and distort the desired signals, while a too small order brings more residual reverberation. To solve this problem, the MCLP filter coefficients are modelled with a prior Gaussian distribution in [START_REF] Chetupalli | Late reverberation cancellation using bayesian estimation of multi-channel linear predictors and student's t-source prior[END_REF], and a statistical model for late reverberation is proposed to avoid the over-estimation of the undesired late reverberant component in [START_REF] Jukić | Adaptive speech dereverberation using constrained sparse multichannel linear prediction[END_REF]. As demonstrated by [START_REF] Dietzen | Comparative analysis of generalized sidelobe cancellation and multi-channel linear prediction for speech dereverberation and noise reduction[END_REF], the MCLP method and the beamforming methods are more suitable for predominant late-reverberant environments and predominant noise environments, respectively. Combinations of beamforming and MCLP filter have been introduced for dereverberation in noisy environments [START_REF] Cohen | Combined weighted prediction error and minimum variance distortionless response for dereverberation[END_REF] [START_REF] Hashemgeloogerdi | Joint beamforming and reverberation cancellation using a constrained kalman filter with multichannel linear prediction[END_REF] [START_REF] Nakatani | Jointly optimal denoising, dereverberation, and source separation[END_REF]. To achieve high-efficiency dereverberation in noisy environments, the MCLP can either work alone [START_REF] Huang | Kronecker product multichannel linear filtering for adaptive weighted prediction error-based speech dereverberation[END_REF] or combine with a spatial filter [START_REF] Yang | Robust dereverberation with kronecker product based multichannel linear prediction[END_REF] [START_REF] Yang | Dereverberation with differential microphone arrays and the weighted-prediction-error method[END_REF] [START_REF] Yang | A bilinear framework for adaptive speech dereverberation combining beamforming and linear prediction[END_REF] in the Kronecker product array framework, or combine with the generalized sidelobe cancellation filter in the Kalman filter framework [START_REF] Dietzen | Integrated sidelobe cancellation and linear prediction kalman filter for joint multimicrophone speech dereverberation, interfering speech cancellation, and noise reduction[END_REF]. In [START_REF] Wang | Convolutive prediction for monaural speech dereverberation and noisy-reverberant speaker separation[END_REF], a deep neural network is trained to predict the direct sound from noisyreverberant speech, rather than exploiting the linear filter structure. To the best of our knowledge, few work on the MCLP filter has been proposed for source localization in reverberant environments [START_REF] Xue | Direct-path signal cross-correlation estimation for sound source localization in reverberation[END_REF].

The commonly used source localization methods in reverberant environments are based on the W-disjoint orthogonal (WDO) assumption [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] that at most one source is dominant with respect to other sources at one time-frequency (TF) bin. Extensive research has focused on selecting such bins to enhance the robustness of the source localization methods against reverberations. The single source regions (SSRs) and low reverberant single source (LRSS) regions are selected by coherent tests in [START_REF] Mohan | Localization of multiple acoustic sources with small arrays using a coherence test[END_REF][25] [START_REF] Wu | Multi-source direction-ofarrival estimation in a reverberant environment using single acoustic vector sensor[END_REF]. Single source confidence measure is proposed in [START_REF] Rafaely | Speaker localization using direct path dominance test based on sound field directivity[END_REF]. The direct path dominance (DPD) tests are based on the ratio between the largest and second largest eigenvalues [START_REF] Rafaely | Speaker localization in reverberant rooms based on direct path dominance test statistics[END_REF] or the dominant eigenvector [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF] of the covariance matrix of the microphone signals, or the sound filed directivity [START_REF] Rafaely | Speaker localization using direct path dominance test based on sound field directivity[END_REF]. Compared with the planar and spherical microphone arrays, acoustic vector sensors (AVS) for source localization [START_REF] Wu | Multi-source direction-ofarrival estimation in a reverberant environment using single acoustic vector sensor[END_REF][30] require less space and computation, as the bins are selected by comparing the phases of the real and imaginary parts of the microphone signals without eigenvalue decomposition [START_REF] Rafaely | Speaker localization in reverberant rooms based on direct path dominance test statistics[END_REF]. Further enhancement have been achieved through outlier removal [START_REF] Geng | Multi-level time-frequency bins selection for direction of arrival estimation using a single acoustic vector sensor[END_REF], modelling outliers as an extra Gaussian cluster [START_REF] Jia | Multi-source doa estimation in reverberant environments by jointing detection and modeling of time-frequency points[END_REF], generelizing single source dominance TF bin to single source dominance window [START_REF] Ying | Window-dominant signal subspace methods for multiple short-term speech source localization[END_REF], and deconvolution of the global pseudo-spectrum by iteratively source detection [START_REF] Yang | Multiple sound source counting and localization based on tf-wise spatial spectrum clustering[END_REF] for the sources with small separations. All these TF bin selection based methods require empirically chosen threshold values for different experiment configurations. A too large threshold leads to miss detection of the sources, while a too small threshold brings more outliers.

In our paper, in order to avoid the need for empirical threshold selection, we leverage the signal model composed of the early reverberant component and the late reverberant component. Specifically, the algorithm is designed for multiple sources scenario, where the late reverberant component is removed with the MCLP filter and the directions of arrival (DOAs) are estimated jointly with the early component signals.

To account for the assumptions that the sources are unoverlapped, finite in power, and their direction angles occupy only a fairly small part of the whole angular section of interest, sparsity in the l 1 norm and orthonormalization constraints are imposed on the DOA estimation matrix. By applying the azimuth sparsity constraint, the true DOA can be estimated without the interference of the early reflections. Besides, considering the sparsity of the TF spectrum of information bearing signals, the l 1 norm constraints are imposed both on the source signals and the dereverberated signals. To solve the proposed criterion, the MCLP filter, the DOAs, and the source signals are estimated alternately. Notably, proximal alternating linearizations are introduced for estimating the coupling parameters with non-smooth constraints [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], and the orthogonal Procruste problem framework [START_REF] Everson | Orthogonal, but not orthonormal, procrustes problems[END_REF] is employed to deal with the orthonormalization constraint. In the following, we propose a DOA estimation algorithm for the single source scenario, where the orthonormalization constraint on the DOA estimation matrix is replaced by the constraint on the norm of the DOA estimation vector. This algorithm can be considered a special case of the multiple sources localization problem. At last, the source number estimation method and a postprocessing procedure for searching the global solution to the proposed algorithm are discussed.

The paper is organized as follows: in section II, the microphone array signals model in reverberant environments and the MCLP filter are briefly reviewed. To clarify the motivation of our work, the disadvantages of the cascade of the MCLP filter and the conventional MUSIC are also outlined. In section III and IV, the localization algorithms are proposed for multiple and single source scenarios, respectively. In section V, pre and post-processing procedures for better performance are discussed. Evaluations on both simulated and real-world data illustrate the advantage of our proposed algorithm in section VI. Finally, conclusions are given in section VII.

II. PRELIMINARIES

A. Signal model

Considering that Q far-field wideband acoustic sources s tq , q = 1, ..., Q impinge on M microphones in a reverberant room, the time domain signals received by the microphone array at time t can be given as:

x t (t) = Q q=1 h tq * s tq (t) + b t (t), (1)
where h tq (t) = h tq,1 , ..., h t q,M T ∈ C M ×1 denotes the room impulse response vector from source q to the microphone array at time index t, * denotes the convolution operator, b t (t) denotes an additive noise at the microphone array at time t with variance σ 2 b . In reverberant environments, the signals received at the microphone array in the TF domain are approximately formulated as a convolution of the STFT of the room impulse response h q (n, ω) of order K and the STFT of the signals s q (n, ω) along the time frame axis for each frequency bin [START_REF] Talmon | Relative transfer function identification using convolutive transfer function approximation[END_REF], since the length of the room impulse response is generally much larger than that of the time frame of the STFT. The microphone array signals in the TF domain can be given as:

x(n, ω) = Q q=1 K k =0 h q (k , ω)s q (n -k , ω) + b(n, ω) (2)
where n = 1, ..., N and ω denote the time frame and frequency bin indices, respectively. s q (n, ω) and b(n, ω) indicate the STFT coefficients of the q th source and the additive noise at the microphone array, respectively.

The microphone array signals x(n, ω) can be reformulated as the sum of the early and late reverberant components and the noise:

x(n, ω) = x e (n, ω) + x l (n, ω) + b(n, ω), (3)
where

x e (n, ω) = Q q=1 ke k =0 h q (k , ω)s q (n -k , ω
) is defined as the early reverberant component which contains mainly the direct-path and few early reflections of the room impulse response, and x l (n, ω) = Q q=1 K k =ke h q (k , ω)s q (n -k , ω) is defined as the late reverberant component which contains all other reflections of the room impulse response in the TF domain. The value of k e varies from 0 to 4 according to experimental configurations. In our work, we choose k e = 0 in order to minimize the influence of the early reflections on the DOA estimation performance. The details of x e (n, ω) can be given as:

x e (n, ω) = Q q=1   a(θ q , ω) + Q q =1
γ q,q a(θ q,q , ω)   s q (n, ω), (4) where a(θ q , ω) is called the steering vector of a point source from direction θ q such that a(θ q , ω) = 1, e -jω∆t1(θq) , ..., e -jω∆t M -1 (θq) ∈ C M ×1 , when the first microphone is taken as the reference, ∆t m (θ q ) is the time delay from the m th microphone to the reference microphone. The term a(θ q , ω) in (4) denotes the direct-path signals from the q th sources, and a(θ q,q , ω) in (4) denote the q early reflection from the q th source. The q th early reflection of the q th source can be considered as sources correlated with s q but from DOA θ q,q which is different from the original DOA θ q with an attenuation factor γ q,q due to the absorption and time delay during the reflections. In our work, we conduct the experiments with speech signals, in which case the length of the time frame is 20 to 30 ms according to the short time stationary property, the early reverberant component Q is generally no more than 1.

B. The MCLP filter and the GWPE method

The generalized weighted prediction error (GWPE) method [START_REF] Yoshioka | Generalization of multi-channel linear prediction methods for blind mimo impulse response shortening[END_REF] reduces the late reverberation component correlated to the past time frames with a MCLP filter matrix for each current time frame of the microphone array signals in the TF domain. Defining G(ω) = [g 1 , ..., g m , ..., g M ] ∈ C K l M ×M as the MCLP filter matrix for frequency bin ω, where K l gives the prediction order, g m are the MCLP filters, the dereverberated speech components can be estimated as:

ŷ(n, ω) = x(n, ω) -G H (ω)x(n, ω), (5)
where

x(n, ω) = x T (n -∆ -1, ω), ..., x T (n -∆ -K l , ω) T ,
∆ denotes the delay tap index of the beginning frame of the late reverberation component. Assuming that the estimated speech component ŷ(n, ω) follows a complex multivariate Gaussian distribution such that ŷ(n, ω) ∼ N C (0, Λ n ), n = 1, ..., N , where Λ n denotes the time varying covariance matrix of ŷ(n, ω), the log-likelihood function can be given as:

L(Θ ω ) = N n=1 -log (det(Λ n ))-ŷ H (n, ω)Λ -1 n ŷ(n, ω)+const, (6)
where Θ ω = G H (ω), Λ n is the set of unknown parameters to be estimated. ( 6) can be maximized by an alternating twostep scheme: In the first step, ( 6) is maximized with respect to G H (ω) while Λ n is kept fixed:

Ĝ(ω) = arg max G(ω) N n=1 -ŷ H (n, ω)Λ -1 n ŷ(n, ω), (7)
and in the second step, ( 6) is maximized with respect to Λ n while G H (ω) is kept fixed:

Λn = arg max Λn N n=1 -log (det(Λ n )) -ŷH (n, ω)Λ -1 n ŷ(n, ω).
(8) The algorithm can converge in a few iterations.

C. Combined GWPE and conventional source localization method

One approach we can try for the source localization problem in reverberant environments is to combine a conventional source localization problem with the GWPE method. In this approach, the source localization method is applied to the output ŷ(n, ω) of the GWPE method. In our experiments, the sources are located with the famous MUSIC algorithm [START_REF] Stoica | Music, maximum likelihood, and cramer-rao bound[END_REF]. With the number of source Q assumed to be known, the DOAs are estimated by searching the corresponding steering vectors which are the most orthogonal to the noise subspace such that:

θq = arg max θ 1 a H (θ, ω)V 2 2 , (9)
where the matrix V is the basis of the noise subspace of the covariance matrix of the microphone array signals, • 2 is the l 2 norm operator.

Figure 1 shows the MUSIC pseudo-spectrum with the dereverberated speech components ŷ(n, ω) (legend: "GWPE+MUSIC") and with the original microphone array signals x(n, ω) (legend:"MUSIC"), respectively. Room impulse responses were generated using [START_REF] Habets | Room impulse response generator[END_REF] for an 8m×8m×3m room with a linear microphone array composed of M = 8 microphones spaced by 5cm. T 60 is set to 600 ms and the additive noise is not considered in this simulation. Two incoherent speech sources come from the DOA 90 • and 135 • , respectively. The assumed source number for MUSIC is set to Q = 2 in Fig. 1(a) and Q = 3 in Fig. 1(b). It can be observed that the conventional MUSIC identifies one peak around 50 • both in Fig. 1(a) and 1(b) due to reflections. However, by combining the GWPE and MUSIC, this false peak is substantially eliminated and only a small residual peak remains around 50 • , since in (4) the early reverberant component contains also a few early reflections. The source from the DOA 90 • is omitted in the pseudo-spectra in Fig. 1(a) for both methods, while successfully located with a minor bias in Fig. 1(b) where Q is larger than the true source number. The simulation results highlight that the simple combination of the dereverberation method and the DOA estimation method is insufficient for the source localization problem in the reverberant environments.

III. THE PROPOSED ALGORITHM IN MULTIPLE SOURCES

SCENARIO

In section III and IV, the joint dereverberation and DOA estimation algorithm are investigated for both multiple sources and single source scenarios. For the sake of simplicity, all the algorithms are studied in the same sub-frequency bin in the STFT domain, ω is omitted in the parameters in the following of this paper.

A. The cost function

Assuming that the whole angular sector of interest is discretized by L grid point such that φ k = φ 1 + (k -1)δ, k = 1, ..., L, where δ is the angular separation between two points of the grid. A matrix α ∈ C L×Q can be defined such that the non-zero elements in the q th column of α corresponds to the DOAs of the direct signals and early reflections of the q th source. We also define H = [a(φ 1 ), ..., a(φ k ), ..., a(φ L )] as the matrix composed of the steering vectors corresponding to all the φ k . Therefore, according to (4), a grid model for the early reverberant component can be given as:

x e (n) ≈ Hαs(n), (10)
where the elements in the vector s(n) = [s 1 (n), ..., sQ (n)] T are linearly proportional to the actual source signals

s(n) = [s 1 (n), ..., s Q (n)] T ∈ C Q×1 in the TF domain.
The proposed DOA estimation criterion involves an optimization problem with three constraints:

1) The data fitting term: According to (3), we aim to minimize the error between the actual received microphone array signal x(n) and the signal modelled as the sum of the early and late reverberant components. The early components are given by [START_REF] Chetupalli | Late reverberation cancellation using bayesian estimation of multi-channel linear predictors and student's t-source prior[END_REF] and the late components are given by G H x(n) obtained in section II-B. The data fitting term can be given as:

C 1 = N n=1 x(n) -G H x(n) -Hαs(n) 2 2 , (11)
It is worth noting that ( 11) is also an approximation of the additive noise b(n, ω) in (3). Since the power of the noise is independent with the time index, we omit the covariance matrix of the noise, which is different from [START_REF]Linear prediction-based online dereverberation and noise reduction using alternating kalman filters[END_REF].

2) Super-Gaussian constraint for the early reverberant components and source signals: Opposed to the Gaussian distribution assumption, previous works [START_REF] Jukić | Speech dereverberation using weighted prediction error with laplacian model of the desired signal[END_REF][10][40] [START_REF] Wenmeng | Speech enhancement with robust beamforming for spatially overlapped and distributed sources[END_REF] have illustrated that the super-Gaussian distribution model is more appropriate for information bearing signals such as speech signals. The optimization problem for the super-Gaussian distributed signals generally involves the minimization of l p norm of the targets with 0 < p < 2. Here, we introduce the l 1 norm of the dereverberated speech components x(n)-G H x(n) and the Q original sources s(n):

C 2 = λ z1 N n=1 x(n) -G H x(n) 1 + λ s N n=1 s(n) 1 , (12)
where λ z1 and λ s are positive sparsity penalization parameters,

• 1 denotes the l 1 norm operator. By minimizing [START_REF] Dietzen | Comparative analysis of generalized sidelobe cancellation and multi-channel linear prediction for speech dereverberation and noise reduction[END_REF], we can exploit the super Gaussian distribution property as well as the sparsity of their STFT spectra of the signals.

3) Azimuth sparsity constraint: The l 1 norm of α is minimized by exploiting the spatial sparsity property of the DOAs of the sources. Indeed, our simulations in section VI have shown that if the sparsity weighting parameter is large enough, the elements corresponding to the early reflections diminish, leaving only those corresponding to the direct signals. The estimated DOAs are obtained from the indices of the nonzero elements in each column of α. Besides, we insert the constraint α H α = I under the assumption that the sources are not spatially overlapped and have finite power:

C 3 = α 1 , s.t. α H α = I. (13)
It is important to note that the normalization constraint on the diagonal of α H α aims to place the proposed criterion in the framework of the orthogonal Procrustes problem [START_REF] Everson | Orthogonal, but not orthonormal, procrustes problems[END_REF] in order to theoretically guarantee convergence. For satisfying the constraint on α, the ratio between sq (n) and s q (n) varies with the power of the q th source.

Combining C 1 , C 2 and C 3 in ( 11)-( 13), the parameters for the source localization problem in reverberant environments can be estimated by:

Ĝ, α, ŝ(n) = arg min G,α,s(n) N n=1 x(n) -G H x(n) -Hαs(n) 2 2 + λ α α 1 + λ z1 N n=1 x(n) -G H x(n) 1 + λ s N n=1 s(n) 1 , s.t. α H α = I, (14)
where λ α is the sparsity penalization parameter.

B. Solutions for the cost function

Equation ( 14) can be solved via Gauss-Seidel iteration scheme [START_REF] Ortega | Iterative solution of nonlinear equations in several variables[END_REF], the solution converges under the condition that the minimum is uniquely attained in each step [START_REF] Bertsekas | Parallel and distributed computation: Numerical methods[END_REF]. Since the cost function consists of coupling second-order terms and non-smooth constraints, the linearized Gauss-Seidel method, namely the Proximal Alternating Linearized Minimization (PALM) method [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] is utilised to alternately estimate G, α, and s(n) until convergence is reached. All the three parameters are resolved by sub-iterations rather than simple closed-form expressions. For simplicity, only the indices in sub-iterations are kept and those in the top level iteration are omitted. The alternative iterations are as follows:

1) Solve G: In the branch for estimating G, α and s(n) are assumed constant, and equation ( 14) is simplified as:

Ĝ = arg min G N n=1 x(n) -G H x(n) -Hαs(n) 2 2 + λ z1 N n=1 x(n) -G H x(n) 1 . (15)
By introducing an extra variable z 1 (n), (15) can be reformulated as:

Ĝ = arg min G N n=1 x(n) -G H x(n) -Hαs(n) 2 2 + λ z1 N n=1 z 1 (n) 1 , s.t. z 1 (n) = x(n) -G H x(n). ( 16
)
The augmented Lagrangian of ( 16) can be given as:

L(G, z 1 (n), η G (n)) = N n=1 x(n) -G H x(n) -Hαs(n) 2 2 + λ z1 z 1 (n) 1 + Re η H G (n)(x(n) -G H x(n) -z 1 (n)) + 1 2ρ G x(n) -G H x(n) -z 1 (n) 2 2 , (17)
where Re {•} is the real part operator, ρ G is the penalization parameter of the convex term.

The problem [START_REF] Yang | Robust dereverberation with kronecker product based multichannel linear prediction[END_REF] can be solved via some iterative steps, in the (l + 1) th iteration, the closed-form solutions for G and z 1 are given as:

G (l+1) = N n=1 1 + 1 2ρ G x(n)x H (n) -1 • N n=1 1 + 1 2ρ G x(n)x H (n) + 1 2 x(n)η (l)H G (n) -x(n)s H (n)α H H H - 1 2ρ G x(n)z (l)H 1 (n) , (18)
and:

z (l+1) 1 (n) = S λz 1 /µz 1 z (l) - 1 µ z1 ∇ z V (G (l+1) , z (l) 1 , η (l) G , n) , (19)
where • -1 is the inverse operator of a matrix, and S λz 1 /µz 1 (v) is the soft thresholding operator of the vector v such that:

S λz 1 /µz 1 (v) =    v -λ z1 /µ z1 , if v ≥ λ z1 /µ z1 , v + λ z1 /µ z1 , if v ≤ -λ z1 /µ z1 , 0, otherwise, (20)
where the comparisons and subtractions are applied elementwise, and the comparisons are with the absolute values of the elements in v.

The details for obtaining ( 18) and ( 19) can be found in Appendix A.

2) Solve α and s(n):

In the branch for estimating α and s(n), G is assumed constant. The cost function is constructed by introducing ŷ(n) given by ( 5) into ( 14):

α, ŝ(n) =arg min α,s(n) N n=1 ŷ(n) -Hαs(n) 2 2 + λ α α 1 + λ s N n=1 s(n) 1 + ι(α), (21)
where ι(α) is an indicator function such that:

ι(α) = 0, if α H α = I, +∞, otherwise. (22)
The problem ( 21) consists of a coupling function Hαs(n) as well as non-smooth functions α 1 and s(n) 1 , which can also be solved by the linearized Gauss-Seidel iteration scheme. Proximal terms are introduced to assure convexity, as the problem includes a non-convex constraint in ι(α). In the (l + 1) th iteration, the parameters can be estimated as:

s(l+1) (n) =arg min s(n) N n=1 ŷ(n) -Hα (l) s(n) 2 2 + λ s s(n) 1 + d s 2 s(n) -s(l) (n) 2 2 , (23a)
α (l+1) =arg min α N n=1 ŷ(n) -Hαs (l+1) (n) 2 2 + λ α α 1 + ι(α) + d α 2 α -α (l) 2 2 , (23b)
where d s and d α are penalization parameters of convex terms.

Similarly to (40b) in appendix A, defining that

B(s (l) (n), α (l) ) = ŷ(n) -Hα (l) s(l) (n) 2 2
, and the gradient of B(s

(l) (n), α (l) ) in terms of s (l) (n) is ∇ sB(s (l) (n), α (l) ) = -α (l)H H H ŷ + α (l)H H H Hα (l) s(l) (n),
the closed-form solution of (23a) can be given as:

s(l+1) (n) = S λs/ds s(l) - 1 d s ∇ sB(s (l) (n), α (l) ) . ( 24
)
To solve (23b), we introduce an extra variable z 2 ∈ C 180Q×1 to solve the non-smooth functions separately:

α (l+1) = arg min α N n=1 ŷ(n) -Hαs (l+1) (n) 2 2 + λ α z 2 1 + ι(α) + d α 2 α -α (l) 2 2 , s.t. z 2 = vec(α), (25)
where vec(•) is the vectorization operator of a matrix. The Lagrangian of (25) can be given as:

L α (α (l) , z 2 , η α ) = N n=1 ŷ(n) -Hαs (l+1) (n) 2 2 + ι(α) + λ α z 2 1 + d α 2 α -α (l) 2 2 + Re η H α (vec(α) -z 2 ) + 1 2ρ α α -z 2 2 2 , ( 26
)
where ρ α is the penalization parameter of the convex term, η α is a vector of size LQ × 1.

The problem [START_REF] Wu | Multi-source direction-ofarrival estimation in a reverberant environment using single acoustic vector sensor[END_REF] can be solved by a few iteration steps, in the (j + 1) th iteration step, the closed-form solution for α can be derived in the orthogonal Procrustes problem framework as:

α (l,j+1) = U L Q U R , (27)
and z 2 is given as:

z (j+1) 2 = S λα/dz 2 z (j) 2 - 1 λ α ∇ z2 C(z (j) 2 ) , (28)
where U L Q is the matrix composed of the first Q columns of U L , U L and U R are orthonormal matrix defined in [START_REF] Hu | Decoupled multiple speaker direction-of-arrival estimator under reverberant environments[END_REF], C(z

(j) 2 ) = Re η H α (vec(α) -z (j) 2 ) + 1 2ρα α -z (j) 2 2
2 , the gradient of C(z

(j) 2 ) in terms of z (j) 2 (n) is ∇ z2 C(z (j) 2 ) = -1 2 η α + 1 2ρα (z (j) 2 -α), d z2 is the coefficient of the proximal term.
The details for obtaining [START_REF] Rafaely | Speaker localization using direct path dominance test based on sound field directivity[END_REF] can be found in Appendix B. The algorithm for solving ( 14) is summarized in Algorithm 1. In addition, the time and space complexity of each step in one iteration of the proposed algorithm is given in table I, where the big O notation O(•) is the asymptotic upper bound of the complexity of an algorithm. The real-time applicability and the hardware requirements of the proposed algorithm can be evaluated by the time complexity and the space complexity, respectively. A few iterations even without full convergence are sufficient for a rough estimation of the DOA. However, to achieve a precise DOA estimation with high resolution that equals the angular separation between two adjacent grid points in the candidate DOA sector (namely 1 • in our experiments), about 200 iterations of the outer loop and a few iterations of the inner loop are needed, which requires more computational time than the conventional methods of DOA estimation and dereverberation. Thus, there is a tradeoff between the resolution and the computational cost of the DOA estimation.

IV. THE PROPOSED ALGORITHM IN SINGLE SOURCE

SCENARIO

The single source scenario can be considered a special case of the multiple sources scenario, where the matrix α, the vector s(n) and s(n) shrink to a vector of dimension L×1, the scalar s(n) and s(n), respectively. Since there is no overlap problem of multiple sources in this scenario, the constraint on the orthonormality of α can be replaced by a constraint on the norm of α and only the sparsity penalization of α is kept. The DOA estimation criterion can thus be given as:

Ĝ, α, ŝ(n) = arg min G,α,s(n) N n=1 x(n) -G H x(n) -Hαs(n) 2 2 + λ α α 1 + (α) + λ z1 N n=1 x(n) -G H x(n) 1 + λ s N n=1 |s(n)|, (29)
Algorithm 1 Localization algorithm for multiple sources scenario.

1: Initialization: l = j = 0, λ z1 , λ s, λ α , µ z1 , d s, d α , d z2 , ρ G , ρ α ; 2: repeat Calculate G (l+1) with (18) for solving (40a); (n) with [START_REF] Yang | A bilinear framework for adaptive speech dereverberation combining beamforming and linear prediction[END_REF] for solving (40b); Calculate s(l+1) (n) with [START_REF] Mohan | Localization of multiple acoustic sources with small arrays using a coherence test[END_REF] for solving (23a); Calculate α (l,j+1) with ( 27) for solving (42a); where:

(α) = 0, if α H α = 1, +∞, otherwise. (30)
The framework for solving [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF] is based on the integrated PALM method, which is similar to the multiple sources scenario. The main difference lies in estimation of α. Instead of (23b), α can be estimated by:

α (l+1) =arg min α N n=1 ŷ(n) -Hαs (l+1) (n) 2 2 + λ α α 1 + (α) + d α 2 α -α (l) 2 2 . ( 31
)
No extra augmented Lagrangian is needed here for solving [START_REF] Geng | Multi-level time-frequency bins selection for direction of arrival estimation using a single acoustic vector sensor[END_REF]. That is to say, the steps between line 17 and line 24 in Algorithm 1 are replaced by the following two-step scheme for estimating α (l+1) :

α m =S λα/dα α (l) - 1 λ α N n=1 (-H H ŷ(n)s (l+1) (n) +H H Hα s(l+1) (n) 2 ) , (32)
where the soft thresholding is used for the l 1 norm penalization followed by a projection operator for the unitary norm penalization:

α (l+1) = P c (α m ), (33)
x(n) ∈ C K l M ×1 , s(n) ∈ C Q×1 , α ∈ C L×Q , H ∈ C M ×L , ŷ(n) ∈ C M ×1 equation (18) (19) (40c) (24) (27) (28) (42c) time O(K 2 l M 2 N ) + O(K 3 l M 3 ) O(K l M 2 N ) O(K l M 2 N ) O(QLN ) O(QM N ) + O(Q 2 N ) O(1) O(QL) +O(QLM N ) + O(QK l M N ) +O(M LN ) +O(Q 2 L) + O(QM L) space O(K l M ) + O(QL) O(K l M 2 ) O(K l M 2 ) O(QL) O(QL) + O(M L) O(QL) O(QL) +O(M L) +O(M L) +O(Q 2 )
where:

P c (v) = v, if v H v = δ v , v/ v 2 , otherwise, (34)
where δ v denotes the norm constraints on v.

V. PRE-AND POST-PROCESSING

A. Source number estimation

In ( 14) the source number is assumed a-priori to be known for determining the size of α. The K-means algorithm is commonly used for jointly estimating the source number and the DOAs [START_REF] Hu | Decoupled multiple speaker direction-of-arrival estimator under reverberant environments[END_REF] [START_REF] Jain | Data clustering: 50 years beyond k-means[END_REF]. Here, to improve the efficiency, we employ the eigenvalues of the covariance matrix of the microphone array signals to directly estimate the source number. Numerical experiments have shown that although the reverberations may cause perturbations in the eigenvalue calculations, the eigenvalues for the early and late reflections are generally significantly smaller than those for the actual sources [START_REF] Yang | Multiple sound source counting and localization based on tf-wise spatial spectrum clustering[END_REF]. Hence, assuming that the eigenvalues are in descending order such that λ 1 ≥ λ 2 ... ≥ λ M , the estimated source number Q can be obtained using the following criterion:

λ 1 λ Q > ε > λ 1 λ Q+1 , ( 35
)
where ε is empirically selected and set to 6 when the noise is neglected.

B. Post-processing of DOA estimation

A post-processing step can be added at the output of the DOA estimation algorithm in section III and IV. As the proposed criteria involve non-convex optimization due to αs(n) in ( 14) as well as αs(n) in [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF], it is crucial to choose well-suited initial value of the parameters, especially of α, to avoid local solutions. In our work, we initialize α with J different values α 01 , ..., α 0j , ..., α 0 J , each vector α 0j is composed of contiguous 0 and 1 such that:

α 0j = [0, ..., 0, 1, ..., 1, 0, ..., 0] , (36)
where the region of component 1 covers the candidate DOAs of the sources. The borders of the components 1 of the J vectors are uniformly distributed around the peaks of the pseudo-spectrum of the conventional beamforming for an apriori rough DOA estimation. At last, we select α 0j for which the final solution yields the minimum value of the cost functions ( 14) and ( 29) as the good initial value.

VI. EVALUATION

In this section, evaluations are performed on simulated data and realistic data. In all the experiments, our proposed methods are compared with the baseline methods including the conventional beamforming [START_REF] Do | Srp-phat methods of locating simultaneous multiple talkers using a frame of microphone array data[END_REF], conventional MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] in [START_REF] Jukić | Speech dereverberation using weighted prediction error with laplacian model of the desired signal[END_REF] and the DPD based method [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF]. Here, to select the directpath dominant regions, the local covariance matrix Rx (n, ω) is calculated by averaging the microphone array measurements over D n time frames and D ω frequency bins [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF]

such that Rx (n, ω) = 1 DnDω Dn dn=1 Dω dω=1 x(n -d n , ω -d ω )x H (n - d n , ω-d ω )
, D n and D ω are generally from 2 to 10 and chosen empirically. Three tests reported in [START_REF] Madmoni | Direction of arrival estimation for reverberant speech based on enhanced decomposition of the direct sound[END_REF] are used:

A 1 = (n, ω) : σ 1 (n, ω) σ 2 (n, ω) > T th1 , (37)
A 2 =    (n, ω) : max θ |a H (θ, ω)x(n, ω)| 2 a H (n, ω)a(n, ω) > T th2    , (38)
and

A 3 = (n, ω) : max θ 1 a H (θ, ω)P ⊥ u1(n,ω) 2 > T th3 , (39)
where A r (•), r = 1, 2, 3 means the set of the TF bins that pass these tests and are considered valid for the DOA estimation, σ 1 (n, ω) and σ 2 (n, ω) are the largest and the second largest eigenvalues of Rx (n, ω), u 1 (n, ω) is the dominant eigenvector of Rx (n, ω). P ⊥ u1(n,ω) gives the subspace of Rx (n, ω) which is orthogonal to u 1 (n, ω). The MUSIC algorithm is applied to each selected single source dominant TF bin or region for estimating the DOA of the corresponding dominant source, and then all the potential DOAs are plotted in a histogram and the peaks are considered as the estimated DOAs of the actual sources. As the thresholds T th1 ∼ T th3 require empirical selection depending on the experimental configurations, the performance of DOA estimation is evaluated with all the three tests and the best one is chosen as that of the DPD based method.

A. Evaluations with simulated data

For the simulation setup, the RIR is generated using the image method [START_REF] Habets | Room impulse response generator[END_REF] in a room with size of 8 m × 8 m × 2.6 m. The microphone array composed of 8 microphones is centered at (3, 0.5, 1) m and the inter-element space is 0.03 m along the x axis. Three point sources are placed at (3, 5.5, 1) m for s 1 , (6, 3.5, 1)m for s 2 and (0.4, 2.9, 1)m for s 3 , the corresponding DOAs are 90 • , 135 bin of 875 Hz using the STFT with half overlapping Hann windows of 16 ms. Noise is neglected in most cases, except in Fig. 7, where white Gaussian noise is added to the microphones with different signal to noise ratios (SNRs) and 100 Monte Carlo simulations are conducted for each configuration.

In Fig. 2, the dereverberation performance of the proposed algorithm is investigated in single-source and two-source scenarios. The magnitude of the TF signals at 875 Hz are plotted against the frame index ranging from 650 to 1100, since the power of signals dominates in this region. The legend x 1 represents the signals received at the first microphone in (2), and ŷ1 stands for the first element in the dereverberated signal vector x(n) -G H x(n) in [START_REF] Hashemgeloogerdi | Joint beamforming and reverberation cancellation using a constrained kalman filter with multichannel linear prediction[END_REF]. It can be observed that the indices of the sharp peaks of x 1 and ŷ1 generally coincident. However, compared to x 1 , the tails of the sharp peaks of ŷ1 are cut off in both signal-source and two-source cases. The differences between x 1 and ŷ1 illustrate an approximation of the dereverberation performance of the proposed algorithms. TABLE II: DOA estimation results in Fig. 3.

src separation 80 • 60 • 40 • true DOA θ 1 θ 2 θ 1 θ 2 θ 1 θ 2 55 • 135 • 75 • 135 • 95 • 135 • beamforming 68 • 131 • 64 • 130 • 97 • 130 • MUSIC 57 • 123 • 61 • 124 • NAN 125 • DPD 69.5 • 132.3 • 65 • 132 • 80 • 132.1 • proposed 62 • 134 • 75 • 138 • 98 • 137 •
In Fig. 3 peak and fail to separate the two sources. In contrast, the proposed method successfully separates the two sources due to the sparsity constraint imposed on α. Fig. 3(b) demonstrates that the proposed method achieves more accurate peaks for DOA estimation compared to other baseline methods although all the methods succeed in separating the two sources. In Fig. 3(c), the proposed method exhibits the best estimation result for the second source with DOA θ 2 = 135 • . However, for the first source with DOA θ 1 = 55 • , the peak of α deviates more from the true DOA compared to MUSIC. This phenomenon can be attributed to the non-convexity term in the algorithm [START_REF] Hashemgeloogerdi | Joint beamforming and reverberation cancellation using a constrained kalman filter with multichannel linear prediction[END_REF] which leads to local solutions of α. To verify this, we use a criterion with the same framework of ( 14), but with a fixed a-priori value of α that is good, and estimate only G and s(n). Compared to the results of ( 14), the results of this new criterion yields a smaller value of [START_REF] Jukić | Adaptive speech dereverberation using constrained sparse multichannel linear prediction[END_REF] which is one part of the cost function. In addition, as stated in section III, when the sparsity weight is sufficiently large, the proposed method eliminates early reflections and only estimates the DOAs of the actual sources. The pseudo-spectra intuitively illustrate the resolution of the DOA estimation algorithms. The main lobes of the proposed method are thinner than those of the baseline methods, which implies its higher resolution. Table II presents the numerical results of DOA estimation in Fig. 3, where "NAN" means the algorithm fails in detecting the source, namely there is no obvious peak near the true DOA.

From table II it can be seen that the proposed method achieves higher DOA estimation accuracy than the baseline methods. Fig. 4 to Fig. 6 illustrate the mean absolute errors (MAEs) between the estimated DOAs and the true DOAs for one source, two sources, and three sources, respectively. Several observations can be made from these figures. Firstly, it is evident that the proposed method outperforms the three baseline methods with the smallest MAEs, following behind is the DPD method; secondly, the MAE increases with the source number and the value of T 60 . In particular, the sub-figure for s 1 in Fig. 5 shows that the MAE for beamforming and MUSIC exceeds 20 • degree and approaches 40 • when T 60 exceeds 200ms, which indicates that spurious peaks are detected, rendering the DOA estimation results invalid. Similarly, in the sub figures for s 1 and s 2 in Fig. 6 point. The proposed method achieves superior performance in most cases. However, the simulation results also reveal the limitations of the proposed method: firstly, as the MCLP filter cannot eliminate completely the late reverberation, the performance of the proposed algorithm deteriorates as T 60 increases; secondly, in the sub-figure for s 3 in Fig. 6, a local solution is observed where the performance of the proposed method is slightly inferior to the DPD based methods. In Fig. 7 the mean value of the estimated DOAs are presented as a function of SNR, with the true DOAs indicated by dashed pink lines. Generally, the three baseline methods fail to give a DOA estimation when the SNR is -20 dB. In Fig. 7(a) it can be observed that the influence of the noise decreases significantly when the SNR exceeds 10 dB, which can be supported by Fig. 7(b) -Fig. 7(d), where the estimated DOAs from all the methods exhibit stability with a SNR larger than 10 dB. In this analysis, instead of the MAEs, the estimated and the true DOA lines are plotted separately, since the crossing points of these lines in Fig. 7(b) and Fig. 7(c) illustrate that the noises at the microphones and the room reverberations have opposite effects on the performance the DOA estimation methods. This intriguing phenomenon motivates further theoretical investigations in future studies.

B. Evaluations with real-world data

We utilise the single-and multichannel audio recordings database (SMARD) [START_REF] Nielsen | The single-and multichannel audio recordings database (smard)[END_REF] In Fig. 8 the performance of the four methods with real data is illustrated in both one source and two sources cases. In one source case (Fig. 8(a)) beamforming and MUSIC find one good peak for the true source and one spurious peak near 70 • due to the early reflections. For beamforming the spurious peak is dominant rather than the good peak in the pseudo-spectrum. The DPD-based method provides a more accurate DOA estimation yet a spurious peak also appears due to outliers in the TF bins selecting procedure. The proposed method outperforms the three baseline methods for two reasons. On one hand, spurious peaks are avoided in the pseudo-spectra, since the early reflections are eliminated by the sparsity constraint of α. On the other hand, the DOA estimation accuracy is improved, since the late reverberant component is reduced by the MCLP filter matrix G. In two sources case (Fig. 8(b)), beamforming fails to estimate neither of the two sources, MUSIC can give a rough estimate for the source with DOA 126 • but fails to find the source with DOA 90 • , the DPD-based method gives an accurate estimation for the source wih DOA 90 • but fails for the source with DOA 126 • , exhibiting several spurious peaks due to the outliers as in Fig. 8(a). The proposed method shows a significant advantage over the three baseline methods for the estimation of both the two sources. Similarly to Fig. 3, the main lobes in the pseudo-spectra illustrate the resolution of the algorithms, a thinner main lobe implies a higher resolution of the DOA estimation method. Table III present the numerical results of the DOA estimation for the four methods, where "NAN" means the algorithm fails in detecting the source, namely there is no obvious peak near the true DOA. It can be seen that, although the baseline methods may have higher accuracy in some cases, they are prone to fail in detecting the sources or produce false peaks. The proposed method, on the other hand, can consistently estimate the DOA of the sources accurately without spurious peak. All these experiment results illustrate the robustness of the proposed method.

VII. CONCLUSION

This paper presents a novel approach for joint DOA estimation and dereverberation in reverberant environments. The proposed method effectively removes the late reverberant component using the MCLP filter, allowing for accurate DOA estimation using the early component signals in conjunction. The early reflections are neglected and only the DOAs of the actual sources are estimated when the sparsity weight is large enough. The DOAs, the source signals, and the MCLP filter coefficients are estimated by the linearized alternative iterations, namely PALM. In particular, the orthogonal Procrustes problem framework is introduced to satisfy the orthonormalization constraint on the DOA estimation matrix. Evaluations on both simulated an real-world data demonstrate that compared to the baselines methods, our proposed methods tend to be more accurate and more robust against the reverberations and the noises in case of both multiple sources and single source scenarios.

3 Fig. 1 .

 31 Fig. 1. Pseudo-spectrum of MUSIC with original and dereverberated microphone array signals (sources from θ 1 = 90 • and θ 2 = 135 • ).

until%

  Convergence or up to 10 iterations; 11: % line 12 to line 26 for s(n) and α estimation 12: line 15 for s(n) estimation 15:

16 :%

 16 line 17 to line 24 for α estimation in (23b)

untiluntil

  Convergence or up to 10 iterations; 24:α (l+1) = α (l,j) Convergence or up to 10 iterations; 27: until Convergence or up to 10 iterations;

Fig. 2 .

 2 Fig. 2. Magnitude of STFT signals of microphone array.

Fig. 3 .

 3 Fig. 3. Pseudo-spectrum with different source separations, T 60 = 600 ms.

Fig. 4 .Fig. 5 .

 45 Fig. 4. MAE vs. T 60 , θ = 135 • in one source case.

Fig. 6 .

 6 Fig.4to Fig.6illustrate the mean absolute errors (MAEs) between the estimated DOAs and the true DOAs for one source, two sources, and three sources, respectively. Several observations can be made from these figures. Firstly, it is evident that the proposed method outperforms the three baseline methods with the smallest MAEs, following behind is the DPD method; secondly, the MAE increases with the source number and the value of T 60 . In particular, the sub-figure for s 1 in Fig.5shows that the MAE for beamforming and MUSIC exceeds 20 • degree and approaches 40 • when T 60 exceeds 200ms, which indicates that spurious peaks are detected, rendering the DOA estimation results invalid. Similarly, in the sub figures for s 1 and s 2 in Fig.6beamforming and MUSIC fail to provide any estimations when T 60 exceeds 200ms. In contrast, the DPD based method successfully estimate DOAs for all values of T 60 by utilizing the mask of the single source dominant TF

Fig. 7 .

 7 Fig. 7. Mean value of the estimated DOAs vs. SNR, T 60 = 600 ms.

Fig. 8 .

 8 Fig. 8. Pseudo-spectrum obtained by real world data.

TABLE I :

 I Complexity analyses of Algorithm 1.

parameters

  from Aalborg University for the realworld data experiments. Different types of signals are recorded in a 7.34 m × 8.09 m × 2.87 m box-shaped listening room with a reverberation time of approximately 0.15 seconds. Here male speech samples with the OmniPower 4296 loudspeaker and from the DOAs 90 • and 126 • are selected as sources. A

TABLE III :

 III DOA estimation results in Fig.8.

	src number	1	2
	true DOA	θ 126 •	θ 1 90 •	θ 2 126 •
	beamforming	121 •	NAN	NAN
	MUSIC	126 •	NAN	133 •
	DPD	120 •	90 •	117 •
	proposed	123 •	94 •	122 •

This work was supported by the National Natural Science Foundation of China (Grant No.62101013, No.61831019, No.61971015, and No.62001012).

APPENDIX A

In the (l + 1) th iteration for solving [START_REF] Yang | Robust dereverberation with kronecker product based multichannel linear prediction[END_REF], the parameters can be determined by:

For (40a), a closed-form expression of G (l+1) can be obtained by searching a result such that the gradient of L(G (l) , z

G ) with respect to G (l+1) equals to 0, where • is the conjugate operator, which results in [START_REF] Yang | Dereverberation with differential microphone arrays and the weighted-prediction-error method[END_REF].

For solving (40b), noting that V (G (l+1) , z

, introducing the proximal terms and ignoring the constant terms, (40b) can be derived as:

where < •, • is the inner product operator.

With the soft thresholding operator, a closed-from expression for (41) can be given as in [START_REF] Yang | A bilinear framework for adaptive speech dereverberation combining beamforming and linear prediction[END_REF].

APPENDIX B

In the (j + 1) th iteration for solving [START_REF] Wu | Multi-source direction-ofarrival estimation in a reverberant environment using single acoustic vector sensor[END_REF], the parameters can be determined by:

).

Similarly to ( 24) and (40b), using the linearized method, (42a) can be reformulated as:

where

, and the gradient of

, where M {•} is the operator which reshapes the vector η α and z 2 into a matrix of the same size with α. The (43) can be considered the Frobenius-norm minimization problem with orthonormalization constraint which is also called the orthogonal Procrustes problem [START_REF] Dietzen | Square root-based multi-source early psd estimation and recursive retf update in reverberant environments by means of the orthogonal procrustes problem[END_REF][36], and can be solved by means of the singular value decomposition (SVD).

Defining matrix A such that A = α (l,j) -1 dα ∇ α F (α (l,j) ), the SVD of A can be given as:

where U L ∈ C L×L and U R ∈ C Q×Q are orthonormal matrices, Σ ∈ C L×Q consists of a diagonal matrix of singular values of dimension Q × Q, and a zero-value matrix of dimension (L -Q) × Q. Thus, α (l,j+1) can be given as in [START_REF] Rafaely | Speaker localization using direct path dominance test based on sound field directivity[END_REF].


Pierre Joanne
email: pierre.joanne@sorbonne-universite.fr

Yeranuhi Hovhannisyan

Alexandre Simon

Gaëlle Revet

Romain Diot

Gabriel Friob

Denisa Calin

Zhenlin Li

Anthony Béhin

Karim Wahbi

Gérard Tachdjian

Onnik Agbulut
email: onnik.agbulut@sorbonne-universite.fr

Generation of human induced pluripotent stem cell lines from five patients with Myofibrillar myopathy carrying different heterozygous mutations in the DES gene

A B S T R A C T Myofibrillar myopathy (MFM) is a rare genetic disorder characterized by muscular dystrophy that is often associated with cardiac disease. This disease is caused by mutations in several genes, among them DES (encoding desmin) is the most frequently affected. Peripheral blood mononuclear cells from 5 different MFM patients with different DES mutations were reprogrammed into induced pluripotent stem cells (IPSC) using non-integrative vectors. For each patient, one IPSC clone was selected and demonstrated pluripotency hallmarks without genomic abnormalities. SNP profiles were identical to the cells of origin and all the clones have the capacity to differentiate into all three germ layers.

Resource Table:

Unique

Resource utility

These IPSC lines, generated from five different patients affected by MFM associated with dilated cardiomyopathy (DCM) with five different DES mutations, can be used to study cellular and molecular mechanisms of the MFM. They also can be used as a cell model to develop new therapeutic approaches, especially based on drug-screening.

Resource details

Desmin is the specific intermediate filament of muscle cells, within which it forms a three-dimensional network connecting the various cellular components to each other, including organelles and contractile apparatus [START_REF] Brodehl | Molecular insights into cardiomyopathies associated with desmin (DES) mutations[END_REF]. Mutations in the gene encoding this protein is known to cause desmin-related myofibrillar myopathy (MFM1: OMIM:601419), a heterogeneous disease affecting primarily skeletal muscles (weakness) [START_REF] Joanne | Viral-mediated expression of desmin mutants to create mouse models of myofibrillar myopathy[END_REF]. However, approximately 75 % of patients also present signs of cardiac disease (van Spaendonck-Zwarts et al., 2010) that mainly evolve toward dilated cardiomyopathy and heart failure. At the cell level, mutations of desmin cause several disorganizations of the muscle cell that lead to contractile and metabolic dysfunctions [START_REF] Hovhannisyan | Critical contribution of mitochondria in the development of cardiomyopathy linked to desmin mutation[END_REF][START_REF] Winter | Mutant desmin substantially perturbs mitochondrial morphology, function and maintenance in skeletal muscle tissue[END_REF].

Here we described the generation of human induced pluripotent stem cells from the peripheral blood mononuclear cells (PBMC) of 5 patients affected by MFM associated with a cardiac dysfunction. These five patients have different mutations in DES even if they present similar cardiac symptoms (see Table 1.). These mutations (DES D214-E245Del , DES S46Y , DES P419H , DES E439K , DES E245D ) have been previously described [START_REF] Wahbi | High cardiovascular morbidity and mortality in myofibrillar myopathies due to DES gene mutations: a 10-year longitudinal study[END_REF] and are linked to the development of dilated cardiomyopathy. Five different clones (INSRMi012-C from the patient PC173T19 with the mutation DES D214-E245Del , INSRMi013-A from the patient PC179c1 with the mutation DES S46Y , INSRMi019-A from the patient PC129k8 with the mutation DES P419H , INSRMi020-A from the patient PC130k2c with the mutation DES E439K and INSRMi021-A from the patient PC177 3c14 with the mutation DES E245D ) were derived and are described here (Table 1). In details, PBMC of these patients were reprogrammed using non-integrative strategy (Phenocell SAS, France) and clones for each patient were manually isolated and amplified until passage 20 where they were frozen in liquid nitrogen with CTS PSC Cryomedium (Thermo Fisher Scientific, France). After thawing on Matrigel, these clones were passaged at least one time before analysis. Colonies of all clones demonstrate normal morphology in culture with a high nucleus/cytoplasm ratio (Fig. 1A, scale bar = 100 µm) and express pluripotency markers as shown by immunocytochemistry (Fig. 1B, scale bar = 100 µm) and cytometry (Fig. 1C and Supplemental Data A). The absence of residual reprogramming plasmids was confirmed by PCR. Genomic integrity, gender and identity were confirmed for each clone by G-banding karyotyping (Fig. 1E) and by SNP analysis through the comparison of the iPSCs genotype with that of parental PBMCs to confirm their origin (see Table 1 and Supplemental Data C and D). Moreover, each clone was tested for its capacity to differentiate into the three germ layer lineages using the STEMdiff™ Trilineage Differentiation Kit (Stem Cell Technologies, France). As shown in Fig. 1D, each clone expresses endodermic, mesodermic and ectodermic markers after specific differentiation (expresses as fold change compared to undifferentiated condition) further confirming their pluripotency. The presence of the specific mutations (Table 1) was assessed by sequencing the genomic DNA of each clone and reveal that the expected mutation was correctly observed (Supplemental Data B). Finally, absence of mycoplasma was confirmed using the Mycoplasma check service of Eurofins Genomics.

Materials and methods

Generation and culture of iPSCs

Patient PBMC were reprogrammed by a non-integrative method using a commercial episome-based kit (Cat. #A15960, ThermoFischer Scientific) or sendai-based kit (#A16517, ThermoFischer Scientific). After manually isolated, clones are cultured and amplified under 37 C, 5 % CO 2 , 20 % O 2 on Matrigel® (hESC-qualified, Corning) coated dishes with daily replacement of mTeSR TM 1 (Stem Cell Technologies) and passaged in clumps (1:30 ratio) at approximately 80 % confluency using ReLeSR (Stem Cell Technologies).

Immunocytochemistry

Cells cultured on coverslip coated with Matrigel were fixed with 4 % paraformaldehyde for 5 min and permeabilization were performed with cold methanol (-20 • C) for 5 min. After incubation with 5 % bovine serum albumin for one hour cells were incubated with primary antibodies (see Table 2) overnight at 4 • C. After 45 min incubation with secondary antibodies (see Table 2), nuclei were counterstained with Hoescht 33,342 (Life Technologies). Then, slides were mounted with homemade medium containing Mowiol. Images were taken on an inverted epifluorescence microscope (DMi8, Leica).

Cytometry

Cells were detached with accutase (10 min of incubation at 37 • C). Then, cells were fixed with a 1 % formaldehyde solution for 20 min and permeabilized using 0.1 % TritonX-100. Then, cells were labelled fluorescent-coupled antibodies (see Table 2) for 15 min, centrifuged (300 g, 5 min), washed in 0.5 % BSA and resuspended in the same buffer before analysis on a MACSQuant10 cytometer (Miltenyi Biotec). Analysis was performed using the FlowJo TM v.10.6.1.

Karyotyping

At least 2 passages after thawing cells were karyotyped by Cell Guidance systems (UK) and at least 10 metaphases were counted.

Differentiation assay

STEMdiff Trilineage Differentiation Kit (Cat. # 05230, StemCell Technologies) was used according to the manufacturer's instructions, at a density of 400,000 cells/well for the ectoderm and endoderm conditions and 100,000 cells/well for the mesoderm condition. Differentiated cells were then collected for RNA isolation using TRIzol. RevertAid First Strand cDNA Synthesis Kit was used to generate cDNA using random hexamer primers and expression of lineages markers (Table 2) was performed in SYBR Green Master Mix (Roche) on a LightCycler480 (Roche).

Genomic DNA extraction

Genomic DNA was extracted using the Genelute Mammalian Genomic DNA kit (Sigma) and was amplified using specific primers (Table 2). After migration into an 0.1 % agarose gel, bands were purified using the NucleoSpin Gel and PCR clean-up kit (Macherey Nagel) following the manufacturer's instructions.

Sequencing

Genomic DNA was sequenced by Eurofins Genomics using the primers indicated in Table 2.

SNP analysis

Genomic DNA (at least 500 ng) was analysed for SNP by Integragen (Evry, France) using Illumina-Infinium Core24 technology. Results were analyzed using Genome Studio by comparing the SNP profile of the clone to the profile of the corresponding genomic DNA extracted from PBMC of the patient to verify cell identity and genomic integrity.

Mycoplasma checking

Supernatant medium (500 µL) from 80 to 90 % confluency iPSC were boiled at 95 • C for 10 min, centrifuge at 10,000 g to pellet cellular debris and analyzed by eurofins genomics. Reports were deposited as archive with the manuscript.

Credit authorship contribution statement

Fig

  Fig. 1.

Pierre Joanne :

 Joanne Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Validation, Visualization, Writingoriginal draft, Writingreview & editing.

Table 1

 1 Characterization and validation.

	Classification	Test	Result	Data
	Morphology	Photography Bright	Normal for all	Fig. 1A
		field	INSRMi012-C:
			passage 24
			INSRMi013-A:
			passage 21
			INSRMi019-A:
			passage 20
			INSRMi020-A:
			passage 20
			INSRMi021-A:
			passage 21
	Phenotype	Qualitative analysis by	Normal staining	Fig. 1B
		immunocytochemistry	of Nanog
			(green) and Tra-
			1-80 (red)
			markers
			INSRMi012-C:
			passage 24
			INSRMi013-A:
			passage 21
			INSRMi019-A:
			passage 20
			INSRMi020-A:
			passage 20
			INSRMi021-A:
			passage 21
		Quantitative analysis	Assess more	Fig. 1C and
		by flow cytometry	than 90 %	Supplemental
			positive cells for	Data A
			Nanog, Oct3/4
			and Sox.
			staining of
			Nanog and Tra-
			1-80 markers
			INSRMi012-C:
			passage 24
			INSRMi013-A:
			passage 21
			INSRMi019-A:
			passage 20
			INSRMi020-A:
			passage 20
			INSRMi021-A:
			passage 21
	Genotype	Karyotype (G-	Normal	Fig. 1E
		banding) and	karyotype (46,
		resolution	XY or 46,XX)
			Resolution
			450-500
			INSRMi012-C:
			passage 24
			INSRMi013-A:
			passage 21
			INSRMi019-A:
			passage 23
			INSRMi020-A:
			passage 20
			INSRMi021-A:
			passage 23
	Identity	Microsatellite PCR	Not performed
		(mPCR) OR
		STR analysis	SNP Analysis	Supplemental
			INSRMi012-C:	Data C and D
			passage 23	Original data
			INSRMi013-A:	submitted in
			passage 24	archive with
			INSRMi019-A:	journal
			passage 19
			INSRMi020-A:
			passage 19
			INSRMi021-A:
			passage 20

Table 1 (

 1 continued )

	Classification	Test	Result	Data
	Mutation	Sequencing	Heterozygous	Supplemental
	analysis (IF		mutations	Data B
	APPLICABLE)	Southern Blot OR WGS	Not performed
	Microbiology	Mycoplasma	Mycoplasma	Submitted as
	and virology		testing by RT-	supplementary
			PCR	files
	Differentiation	Directed	Expression of	Fig. 1D
	potential	differentiation	genes from cells
			induced for
			differentiation
			toward
			endoderm,
			mesoderm or
			ectoderm
			germlayers:
			FOXA2, EOMES,
			Brachyury,
			MIXL1, LHX2
			and MAP2
			INSRMi012-C:
			passage 24
			INSRMi013-A:
			passage 21
			INSRMi019-A:
			passage 20
			INSRMi020-A:
			passage 20
			INSRMi021-A:
			passage 21
	Donor	HIV 1 + 2 Hepatitis B,	Not performed
	screening	Hepatitis C
	(OPTIONAL)
	Genotype	Blood group	Not performed
	additional	genotyping
	info	HLA tissue typing	Not performed
	(OPTIONAL)

Table 2

 2 Reagents details.

		Antibodies used for immunocytochemistry/flow-cytometry
		Antibody	Dilution	Company Cat	RRID
				#
	Pluripotency	Rabbit (IgG1)	1:200	Cell Signaling	RRID:
	Markers	anti Nanog		#4903	AB_10559205
		clone D73G4
	Pluripotency	Mouse	1:100	Merck	RRID:
	Markers	monoclonal		#MAB4381	AB_177638
		(IgM) anti Tra-
		1-81
	Pluripotency	PE Mouse anti-	1:10	BD Biosciences	RRID:
	Markers	Human Nanog		Human	AB_2722505
		(clone N31-		Pluripotent
		355)		Stem Cell
				Transcription
				Factor Analysis
				Kit
				#51-9006414
	Pluripotency	PerCP-	1:10	BD Biosciences	RRID:
	Markers	Cy™5.5		Human	AB_2722505
		Mouse anti-		Pluripotent
		Oct3/4 (clone		Stem Cell
		40/oct-3)		Transcription
				Factor Analysis
				Kit
				#51-9006267
	Pluripotency	Alexa Fluor®	1:10	BD Biosciences	RRID:
	Markers	647 Mouse		Human	AB_2722505
		anti-Sox2		Pluripotent
				Stem Cell
				Transcription
				Factor Analysis
				Kit
				#51-9006407
	Isotype control	PE Mouse	1:10	BD Biosciences	RRID:
	antibodies	IgG1, κ		Human	AB_2722505
		Isotype		Pluripotent
		Control		Stem Cell
				Transcription
				Factor Analysis
				Kit
				#51-9006406
	Isotype control	PerCP-	1:10	BD Biosciences	RRID:
	antibodies	Cy™5.5		Human	AB_2722505
		Mouse IgG1, κ		Pluripotent
		Isotype		Stem Cell
		Control		Transcription
				Factor Analysis
				Kit
				#51-9006272
	Isotype control	Alexa Fluor®	1:10	BD Biosciences	RRID:
	antibodies	647 Mouse		Human	AB_2722505
		IgG2a, κ		Pluripotent
		Isotype		Stem Cell
		Control		Transcription
				Factor Analysis
				Kit
				#51-9006410
	Secondary	Goat anti-	1:1000	ThermoFisher	RRID:
	antibodies	Rabbit IgG (H		Scientific	AB_2576217
		+ L) Highly		#A-11034
		Cross-
		Adsorbed
		Secondary
		Antibody,
		Alexa Fluor
		488
	Secondary	Goat anti-	1:1000	ThermoFisher	RRID:
	antibodies	Mouse IgM		Scientific	AB_2535807
		(Heavy Chain)		#A-21238
		Secondary
		Antibody,
		Alexa Fluor
		647

Table 2

 2 (continued )

		Antibodies used for immunocytochemistry/flow-cytometry
		Antibody	Dilution	Company Cat	RRID
				#
		Primers
		Target	Size of	Forward/Reverse primer (5′-
			band	3′)
	Differentiation	FOXA2	138 bp	AGCCCGAGGGCTACTCC/
	markers			ACGACATGTTCATGGAGCCC
	Differentiation	EOMES	90 bp	TACCCAGACCCAACCTTTCC/
	markers			GGAGGTCCATGGTAGTCCAG
	Differentiation	BRACHYURY	121 bp	CCAGTGGCAGTCTCAGGTTA/
	markers			CTACTGCAGGTGTGAGCAAG
	Differentiation	MIXL1	95 bp	AGTCCAGGATCCAGGTATGGT/
	markers			GATAATCTCCGGCCTAGCCA
	Differentiation	LHX2	199 bp	GGGTCCTCCAGGTCTGGTTC/
	markers			GTGGGGCTAGTCAAGTCTGTC
	Differentiation	MAP2	157 bp	GGCCAACTTCTCCAGTGAAC/
	markers			TGCTGGAACTCAGCAGGTAA
	House-Keeping	NUBP1	88 bp	CAGAGGGGCTTCATGTCAGG/
	Genes (qPCR)			GATTTCCTCTATAGCCGTGTCC
	Targeted	PCR/		TATCACCCGCAACTGTCTGT/
	mutation	Sequencing		CGTCGTTGTTCTTGTTGGCT
	analysis/	primer for
	sequencing	DESD214-
		E245Del and
		E245D
	Targeted	PCR/		CCCGCCGCCTCCTCCGTG/
	mutation	Sequencing		GAAAGAGGAGTCTGGGGTGC
	analysis/	primer for
	sequencing	DESS46Y
	Targeted	PCR/		GGAGGTTTTGTCTCTTCCCTTT/
	mutation	Sequencing		CTCTACCTAAACTCCCCCACCT
	analysis/	primer for
	sequencing	DESP419H
	Targeted	PCR/		ACCCTGTTACAGAAACCAGCC/
	mutation	Sequencing		AACCATGTGCTAGCCAGAGAAT
	analysis/	primer for
	sequencing	DESE439K

Acknowledgements

This research was partly supported by Sorbonne Université, France, and CNRS France, and INSERM, France. OA and PJ are supported by the AFM-Téléthon, France (Association Française contre les Myopathies) (contract number: 22142 & 28597). YH was funded by a Ph.D. fellowship from the AFM-Téléthon, France during this study. We thank Dr Brigitte Onteniente and Phenocell SAS (Grasse, France) for IPSC generation from MFM1 patients; Stéphanie Bigou and iPS cell culture facility of the Paris Brain Institute who has received funding from the program "Investissements d'avenir" ANR-10-IAIHU-06 for the clonal selection of IPSC; Dr. Nathalie Lefort (Institut Imagine, France) for her assistance concerning the analysis of SNP profiles and all personnel from the IBPS photon microscopy core facility for helpful advice and technical assistance during flow cytometry.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi. org/10.1016/j.scr.2024.103338.


Clémence Guieu-Grandsire
email: clemence.guieu-grandsire@sorbonne-nouvelle.fr

Naomi Yamaguchi
email: naomi.yamaguchi@sorbonne-nouvelle.fr

Shigeko Shinohara
email: shigeko.shinohara@sorbonne-nouvelle.fr

Acquiring two rhotics /ɾ/ and /ʁ/: complexity enlightened by French-Greek bilingualism

Introduction. Rhotic sounds are segments that are difficult for children to master, typically acquired last in several languages such as Greek, where the rhotic is an alveolar tap /ɾ/ [START_REF] Rizopoulou | An acoustic analysis of rhotic production by children and adults[END_REF] or French, which has the uvular fricative /ʁ/ as its rhotic [START_REF] Macleod | The acquisition of consonants in Québécois French: A cross-sectional study of pre-school aged children[END_REF]. This late acquisition can be explained by considering the acquisition of a phonological system as the addition of contrasts gradually enriching the internal structure of a phoneme [START_REF] Rice | Aspects of variability in child language acquisition[END_REF], which may have different specifications: an unspecified feature (absence of feature), or specified as [+] or [-]. Thus, the most specified segments, such as rhotics, will be acquired last. However, according to [START_REF] Clements | Input elaboration, head faithfulness and evidence for representation in the acquisition of left-edge clusters in West Germanic[END_REF], rhotics are universally unspecified, meaning that the phonological representations of these segments are empty. In the context of French-Greek bilingualism, this raises a major question: how will bilingual children acquire rhotics in each language if their phonological representations are not specified? Furthermore, the group of rhotics is the only one for which the uvular is acquired before the alveolar in production among French-Portuguese bilinguals (Almeida, 2012) and German-Spanish bilinguals [START_REF] Kehoe | The development of rhotics: A comparison of monolingual and bilingual children[END_REF]. How can this acquisition order be explained?

Participants. Twenty-seven children aged 2;8 to 6;2 years old were divided into four groups: eight monolingual French-speaking children, six French-Greek bilinguals living in France, five living in Greece and seven monolingual Greek-speaking children were recorded over three sessions every six months. All bilingual children were exposed to simultaneous language input from one parent speaking Greek and the other speaking French.

Method. The corpus (Bernhardt et al., 2017) used in an image naming task contained rhotics in simple onset positions (word-initial or syllable-initial). The French rhotic was produced 341 times, and the Greek rhotic 498 times. All sessions were phonetically transcribed using PHON (Hedlung & Rose, 2022), and the results were categorized into three age groups: 2;8-3;11 years (Group 1; mean=3;3 years), 4;0-4;11 years (Group 2; mean=4;4 years), and 5;0-6;2 years (Group 3; mean=5;4 years).

Results.

In Figure 1, when the target was the alveolar rhotic /ɾ/, we observed that accurate realisation by bilinguals in Greece followed that of monolingual Greek speakers, unlike the low proportion observed in bilinguals in France (up to 2.68%). When the target was the uvular rhotic, bilinguals in France in Group 1 produced it systematically as expected, while those in Greece in the same group produced it less often (83.33%). The uvular rhotic is therefore mastered much earlier by all children and is generally substituted with velar stops and affricates such as [g ͡ ʁ] or [k ͡ ʁ̥ ] for Groups 1 and 2. Bilinguals in Greece replace it with the Greek velar fricatives [ɣ] and [x] in 9.52% of cases from 3;8 years onwards. Substitutions of the alveolar rhotic are much more variable, particularly for bilinguals in France in Group 2, who also delete it in 17.39% of cases. /ɾ/ is mainly substituted with the lateral [l], but also with alveolar stops, labial, alveolar, or post-alveolar fricatives, the nasal [n], the approximants [w] or [j], and finally [ʁ], particularly among bilinguals in France with a proportion of 18.12% of cases between 4;0 and 4;11 years and 19.05% between 5;0 and 6;2 years. Discussion and conclusion. Although language dominance influences the accurate realisation of rhotics in each language, it seems that monolinguals and bilinguals behave similarly regarding the inherent structural complexity of the rhotic type. To explain these results, it is necessary to specify rhotics in French and Greek (Table 1), which differs from the proposal of Goad and Rose ( 2004) but aligns with Almeida (2012). In French, the rhotic would be specified as [DORSAL][+continuant] as well as [+sonorant], a feature particularly difficult for bilinguals in Greece who substitute it with velar fricatives. In Greek, /ɾ/ is specified as [+sonorant], explaining the substitutions with approximants and nasals. Its ambiguous manner of articulation, [+/-continuant], explains substitutions with stops. Finally, [-lateral], distinguishing it from /l/, is acquired later. The nonspecification of place [CORONAL] explains substitutions with other places. The difficulty of /ɾ/ thus lies in the combination of these feature values which can render certain phonological structures more complex than others according to [START_REF] Calabrese | A constraint-Based Theory of Phonological Markedness and Simplification Procedures[END_REF]. The specification of /ʁ/ in French is therefore simpler [START_REF] Webb | Minimalism and French /ʀ/: Phonological representations in phonetically based phonology[END_REF], which would explain the earlier mastery of /ʁ/ by all French-speaking children and its transfer into Greek by French-dominant bilinguals.

Figure 1 :

 1 Figure 1: Percentages of accurate realisations of the target /ɾ/ (left) and /ʁ/ (right).

Table 1 :

 1 Proposal of necessary specification in French and Greek, based onClements & Hume's features (1995).

		/ʁ/	/ɣ/	/d/	/l/	/ɾ/
	[sonorant]	+	-	-	+	+
	[continuant]	+	+	-	+/-	+/-
	[lateral]				+	-
	PLACE	[DORSAL] [DORSAL] [CORONAL] [CORONAL] [CORONAL]