Gaussian modeling-based multichannel audio source separation exploiting generic source spectral model Thanh T. H. Duong, Ngoc Q. K. Duong Senior Member, IEEE, Phuong Cong Nguyen, and Cuong Quoc Nguyen Abstract-As blind audio source separation has remained very challenging in real-world scenarios, some existing works, including ours, have investigated the use of a weakly-informed approach where generic source spectral models (GSSM) can be learned a priori based on nonnegative matrix factorization (NMF). Such approach was derived for single-channel audio mixtures and shown to be efficient in different settings. This paper proposes a multichannel source separation approach where the GSSM is combined with the source spatial covariance model within a unified Gaussian modeling framework. We present the generalized expectation-minimization (EM) algorithm for the parameter estimation. Especially, for guiding the estimation of the intermediate source variances in each EM iteration, we investigate the use of two criteria: [START_REF] Makino | Blind Speech Separation[END_REF] the estimated variances of each source are constrained by NMF, and (2) the total variances of all sources are constrained by NMF altogether. While the former can be seen as a source variance denoising step, the latter is viewed as an additional separation step applied to the source variance. We demonstrate the speech separation performance, together with its convergence and stability with respect to parameter setting, of the proposed approach using a benchmark dataset provided within the 2016 Signal Separation Evaluation Campaign.

KEYWORDS

Multichannel audio source separation, local Gaussian model, nonnegative matrix factorization, generic spectral model, group sparsity constraint.

I. INTRODUCTION

Real-world recordings are often mixtures of several audio sources, which usually deteriorate the target one. Thus many practical applications such as speech enhancement, sound post-production, and robotics use audio source separation technique [START_REF] Makino | Blind Speech Separation[END_REF], [START_REF] Gannot | A consolidated perspective on multimicrophone speech enhancement and source separation[END_REF] to extract individual sound sources from their mixture. However, despite numerous effort in the past decades, blind source separation performance in reverberant recording conditions is still far from perfect [START_REF] Vincent | The Signal Separation Campaign (2007-2010): Achievements and remaining challenges[END_REF], [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF]. To improve the separation performance, informed approaches have been proposed and emerged recently in the literature [START_REF] Liutkus | An overview of informed audio source separation[END_REF], [START_REF] Vincent | From Blind to Guided Audio Source Separation: How models and side information can improve the separation of sound[END_REF]. Such approaches exploit side information about either the sources themselves or the mixing condition in order to guide the separation process. Examples of the investigated Thanh T. H. Duong is with Hanoi University of Mining and Geology and International Research Institute MICA, Vietnam, e-mail: (duongthihienthanh@humg.edu.vn).

Ngoc Q. K. Duong is with Technicolor R&I, France, e-mail: (quang-khanhngoc.duong@technicolor.com).

Phuong Cong Nguyen is with MICA and Hanoi University of Science and Technology, Vietnam, e-mail: (phuong.nguyencong@hust.edu.vn).

Cuong Quoc Nguyen is with Hanoi University of Science and Technology, Vietnam, e-mail: (cuong.nguyenquoc@hust.edu.vn). side information include deformed or hummed references of one (or more) source(s) in a given mixture [START_REF] Souviraà-Labastie | Multichannel audio source separation using multiple deformed references[END_REF], [START_REF] Smaragdis | Separation by humming: User-guided sound extraction from monophonic mixtures[END_REF], text associated with spoken speeches [START_REF] Magoarou | Text-informed audio source separation. example-based approach using non-negative matrix partial co-factorization[END_REF], temporal annotation of the source activity along the mixtures [START_REF] Duong | Temporal annotationbased audio source separation using weighted nonnegative matrix factorization[END_REF], core associated with musical sources [START_REF] Hennequin | Score informed audio source separation using a parametric model of non-negative spectrogram[END_REF], [START_REF] Ewert | Score-informed source separation for musical audio recordings: An overview[END_REF], and motion associated with audio-visual objects in a video [START_REF] Parekh | Motion informed audio source separation[END_REF]. Following this trend, some recent works including ours have proposed to use a very abstract semantic information just about the types of audio sources existing in the mixture to guide the source separation. If one source in the mixture is known as "speech", then several speaker-independent speech examples can be used to create a universal speech model as presented in [START_REF] Sun | Universal speech models for speaker independent single channel source separation[END_REF]; if several types of sound sources in the mixture are known (e.g., birdsong, piano, waterfall), their audio examples found by internet search can be used to learn the corresponding universal sound class models as presented in [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF]. Such universal models were shown to be effective in guiding the source separation algorithm and resulted in promising performance. Inspired by this idea, we have further investigated the use of generic speech and noise model for single-channel speech separation in [START_REF] Duong | Single-channel speaker-dependent speech enhancement exploiting generic noise model learned by non-negative matrix factorization[END_REF] and shown its promising result in (a) the supervised case, where both speech GSSM and noise GSSM are learned during training phase, and (b) the semi-supervised case, where only the speech GSSM is pre-learned. Furthermore, we have proposed to combine the block sparsity constraint investigated in [START_REF] Sun | Universal speech models for speaker independent single channel source separation[END_REF] with the component sparsity constraint presented in [START_REF] Badawy | On-the-fly audio source separation[END_REF] in a common formulation in order to take into account the advantage of both of them [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF].

It should be noted that the works cited above [START_REF] Magoarou | Text-informed audio source separation. example-based approach using non-negative matrix partial co-factorization[END_REF], [START_REF] Ewert | Score-informed source separation for musical audio recordings: An overview[END_REF], [START_REF] Duong | Single-channel speaker-dependent speech enhancement exploiting generic noise model learned by non-negative matrix factorization[END_REF], [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF] considered only a single channel case, where the mixtures are mono, and exploited non-negative matrix factorization (NMF) [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF] to model the spectral characteristics of audio sources. Some recent works have investigated the use of the deep neural networks (DNN) to model the source spectra, where basically the types of sources in the mixture also need to be known as a side information in order to collect training data. Such DNN-based approaches were shown to offer very promising results in single-channel speech and music separation [START_REF] Huang | Joint optimization of masks and deep recurrent neural networks for monaural source separation[END_REF]- [START_REF] Hershey | Deep clustering: Discriminative embeddings for segmentation and separation[END_REF], multichannel speech separation [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF], [START_REF] Wang | Multi-channel deep clustering: Discriminative spectral and spatial embeddings for speakerindependent speech separation[END_REF]. However, they require a large amount of labeled data for training, which may not always be available and the training is usually computationally expensive.

When more recording channels are available thanks to the use of multiple microphones, a multichannel source separation algorithm should be considered as it allows to exploit important information about the spatial locations of audio sources. Such spatial information is reflected in the mixing process (usually with reverberation), and can be modeled by e.g., the interchannel time difference (ITD) and interchannel intensity difference (IID) [START_REF] Michael | A survey of convolutive blind source separation methods[END_REF]- [START_REF] Chami | A new model based underdetermined source separation[END_REF], the rank-1 time-invariant mixing vector in the frequency domain when following the narrowband assumption [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF]- [START_REF] Kitamura | Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization[END_REF], or the full-rank spatial covariance matrix in local Gaussian model (LGM) where the narrowband assumption is relaxed [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]- [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF].

In this paper, we present an extension of the previous works [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF], [START_REF] Duong | Single-channel speaker-dependent speech enhancement exploiting generic noise model learned by non-negative matrix factorization[END_REF], [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF] to the multichannel case where the NMF-based GSSM is combined with the full-rank spatial covariance model in a Gaussian modeling paradigm. Around this LGM, existing works have investigated several source spectral models such as Gaussian mixture model (GMM) [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], NMF as a linear model with nonnegativity constraints [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF], [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF], continuity model [START_REF] Duong | Multichannel harmonic and percussive component separation by joint modeling of spatial and spectral continuity[END_REF], kernel additive model [START_REF] Liutkus | Scalable audio separation with light kernel additive modelling[END_REF], heavy-tailed distributionsbased model [START_REF] Leglaive | Alphastable multichannel audio source separation[END_REF], [START_REF] Magron | Lévy NMF for robust nonnegative source separation[END_REF], and recently DNN [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF]. Focusing on NMF in this study, our work is most closely related to [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF] and [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF] as both of them use NMF within the LGM to constrain the source spectra in each EM iteration. However, our work is different from [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF] in the sense that we use the pre-trained GSSM, so that potentially the algorithm is less sensitive to the parameter initialization, and it does not suffer from the wellknown permutation problem. Our work is also different from [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF] as we exploit the mixed group sparsity constraint to guide the optimization, which allows the algorithm to automatically select the most representative spectral components in the GSSM. In addition, instead of constraining the variances of each source by NMF as done in [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF], [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF], we propose to constrain the total variances of all sources altogether by NMF and show that this novel optimization criterion offers better source separation performance. While part of the work was presented in [START_REF] Duong | Multichannel audio source separation exploiting NMF-based generic source spectral model in Gaussian modeling framework[END_REF], this paper provides more details regarding the algorithm derivation and the parameter settings. Furthermore, the source separation performance analysis and the comparison with existing approaches are extended.

The rest of the paper is organized as follows. We discuss the problem formulation and the background in Section II. We present the proposed GSSM-based multichannel source separation approach in Section III. In this section, we first present two ways of constructing the GSSM based on NMF. Then, to constrain the intermediate source variance estimates, two optimization criteria are introduced, which can be seen as either performing source variance denoising or source variance separation. The generalized EM algorithm is derived for the parameter estimation. We finally validate the effectiveness of the proposed approach in speech enhancement scenario using a benchmark dataset from the 2016 Signal Separation Evaluation Campaign (SiSEC 2016) in Section IV. For this purpose, we first analyze the convergence of the derived algorithm and investigate its sensitivity to the parameter settings in terms of source separation performance. We then show that the proposed algorithm outperforms most state-of-the-art methods in terms of the energy-based criteria.

II. PROBLEM FORMULATION AND MODELING

In this section, we review the formulation and the Gaussian modeling framework for multichannel audio source separation. Let us formulate the problem in a general setting, where J sources are observed by an array of I microphones. The contribution of each source, indexed by j, to the microphone array is denoted by a vector c j (t) ∈ R I×1 and the I-channel mixture signal is the sum of all source images as

x(t) = J j=1 c j (t). (1)
The objective of source separation is to estimate the source images c j (t) given x(t). As the considered algorithm operates in the frequency domain, we denote by c j (n, f ) and x(n, f ) the complex-valued short-term Fourier transforms (STFT) of c j (t) and x(t), respectively, where n = 1, 2, .., N is time frame index and f = 1, 2, ..., F the frequency bin index. Equation ( 1) can be written in the frequency domain as

x(n, f ) = J j=1 c j (n, f ). (2)

A. Local Gaussian model

We consider the existing nonstationary LGM as it has been known to be robust in modeling reverberant mixing conditions and flexible in handling prior information [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. In this framework, c j (n, f ) is modeled as a zeromean complex Gaussian random vector with covariance matrix

Σ j (n, f ) = E(c j (n, f )c H j (n, f )): c j (n, f ) ∼ N c (0, Σ j (n, f )), (3)
where 0 is an I × 1 vector of zeros and H indicates the conjugate transposition. Furthermore, the covariance matrix is factorized as

Σ j (n, f ) = v j (n, f ) R j (f ), (4)
where v j (n, f ) are scalar time-dependent variances encoding the spectro-temporal power of the sources and R j (f ) are time-independent I × I spatial covariance matrices encoding their spatial characteristics when sources and microphones are assumed to be static. Under the assumption that the source images are statistically independent, the mixture vector x(n, f ) also follows a zero-mean multivariate complex Gaussian distribution with the covariance matrix computed as

Σ x (n, f ) = J j=1 v j (n, f ) R j (f ). (5)
Assuming that the mixture STFT coefficients at all timefrequency (T-F) bins are independent, the likelihood of the set of observed mixture vectors x = {x(n, f )} n,f given the set of variance and spatial covariance parameters θ = {v j (n, f ), R j (f )} j,n,f is given by

P (x|θ) = n,f 1 det (πΣ x (n, f )) e -tr(Σ -1 x (n,f ) Ψx(n,f )) , (6)
where det represents determinant of a matrix, tr() stands for matrix trace, and

Ψ x (n, f ) = E(x(n, f )x H (n, f ))
is the empirical covariance matrix. It can be numerically computed by local averaging over neighborhood of each T-F bin (n , f ) as [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF], [START_REF] Duong | Under-determined reverberant audio source separation using local observed covariance and auditory-motivated time-frequency representation[END_REF]:

Ψ x (n, f ) = n ,f w 2 nf (n , f )x(n , f )x H (n , f ), (7)
where w nf is a bi-dimensional window specifying the shape of the neighborhood such that n ,f w 2 nf (n , f ) = 1. We use Hanning window in our implementation. The quadratic T-F presentation as Ψ x (n, f ) aims to improve the robustness of the parameter estimation as it exploits the observed data in several T-F points instead of a single one. The negative log-likelihood derived from ( 6) is

L(θ) = n,f tr Σ -1 x (n, f ) Ψ x (n, f ) + log det πΣ x (n, f ) , (8)
Under this model, once the parameters θ are estimated, the STFT coefficients of the source images are obtained in the minimum mean square error (MMSE) sense by multichannel Wiener filtering as

ĉj (n, f ) = v j (n, f ) R j (f )Σ -1 x (n, f )x(n, f ). (9)
Finally, the expected time-domain source images ĉj (t) are obtained by the inverse STFT of ĉj (n, f ).

B. NMF-based source variance model

NMF has been a well-known technique for latent matrix factorization [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and shown to be powerful in modeling audio spectra [START_REF] Vincent | From Blind to Guided Audio Source Separation: How models and side information can improve the separation of sound[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF]. It has been widely applied to single channel audio source separation where the mixture spectrogram is usually factorized into two latent matrices characterizing the spectral basis and the time activation [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF]. When adapting NMF to the considered LGM summarized in Section II-A, the nonnegative source variances v j (n, f ) can be approximated as

v j (n, f ) = Kj k=1 w jf k h jkn , (10)
where w jf k is an entry of the spectral basis matrix W j ∈ R

F ×Kj +

, h jkn is an entry of the activation matrix

H j ∈ R Kj ×N +
, and K j the number of latent components in the NMF model.

To our best knowledge, this NMF formulation for the source variances within the LGM was first presented in [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF], and then further discussed in [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. However, in those works, the basis matrix W j is not a GSSM as proposed in this article (presented in Section III-A), and thus the parameters {W j , H j } were estimated differently.

C. Estimation of the model parameters

The set of parameters θ is estimated by minimizing the criterion (8) using a generalized EM algorithm (GEM) [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. This algorithm consists in alternating between E step and M step. In the E step, given the observed empirical covariance matrix Ψ x (n, f ) and the current estimate of θ, the conditional expectation of the natural statistics is computed as [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] Σ

j (n, f ) = G j (n, f ) Ψ x (n, f )G H j (n, f ) + (I -G j (n, f ))Σ j (n, f ), (11)
where

G j (n, f ) = Σ j (n, f )Σ -1 x (n, f
) is the Wiener gain, I is an I ×I identity matrix. Then in the M step, given Σ j (n, f ) the parameters θ j = {v j (n, f ), R j (f )} n,f associated to each j-th source are updated in the maximum likelihood sense by optimizing the following criterion [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]:

L(θ j ) = n,f tr Σ -1 j (n, f ) Σ j (n, f ) + log det πΣ j (n, f ) .
(12) By computing the derivatives of L(θ j ) with respect to v j (n, f ) and each entry of R j (f ) and equating them to zero, the iterative updates for these parameters are found as

R j (f ) = 1 N N n=1 1 v j (n, f ) Σ j (n, f ) (13)
v j (n, f ) = 1 I tr(R -1 j (f ) Σ j (n, f )) (14)
At each EM iteration, once v j (n, f ) is updated in the M step by ( 14), it will be further constrained by NMF as [START_REF] Duong | Temporal annotationbased audio source separation using weighted nonnegative matrix factorization[END_REF]. For this purpose, given the matrix of the current source variance estimate V j ∈ R F ×N + whose entries are v j (n, f ), the corresponding NMF parameters are estimated by minimizing the Itakura-Saito divergence, which offers scale-invariant property, as

min Hj ≥0,Wj ≥0 D(V j W j H j ), (15)
where

D(V j W j H j ) = N n=1 F f =1 d IS v j (n, f ) w jf k h jkn , and
d IS (x y) = x y -log x y -1. (16)
The parameters {W j , H j } are usually initialized with random non-negative values and are iteratively updated via the wellknown multiplicative update (MU) rules [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF].

III. PROPOSED GSSM-BASED MULTICHANNEL APPROACH

The global workflow of the proposed approach is depicted in Fig. 1. In the following, we will first review a training phase for the GSSM construction based on NMF in Section III-A. We then propose the NMF-based source variance model fitting with sparsity constraint in Section III-B. Finally, we derive the generalized EM algorithm for the parameter estimation in Section III-C. Note that we focus on NMF as the spectral model in this paper, however, the whole idea of the proposed approach can potentially be used for other spectral models than NMF such as GMM or DNN.

A. GSSM construction

In this section, we review the GSSM construction, which was introduced in [START_REF] Sun | Universal speech models for speaker independent single channel source separation[END_REF], [START_REF] Badawy | On-the-fly audio source separation[END_REF]. We assume that the types of sources in the mixture are known and some recorded examples of such sounds are available. This is actually feasible in practice. For instance, in the speech enhancement, one target source is speech and another is noise and one can easily find speech and noise recordings. We need several examples for each type of source as one recording is usually not fully representative of the others and a source like "noise" is poorly defined. Let us denote by s l j (t) a l-th single-channel learning example of j-th source and its corresponding spectrogram obtained by STFT S l j . First, S l j is used to learn the corresponding NMF spectral dictionary, denoted by W l j , by optimizing the similar criterion as [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF]:

min H l j ≥0,W l j ≥0 D(S l j W l j H l j ) (17)
where H l j is the time activation matrix. Given W l j for all examples l = 1, ..., L j of the j-th source, the GSSM for the j-t source is constructed as

U j = [W 1 j , . . . , W Lj j ], (18)
then the GSSM for all the sources is computed by

U = [U 1 , . . . , U J ]. (19)
As an example for speech and noise separation, in the practical implementation, we may need several speech examples for different male voices and female voices (e.g.

We will show in the experiment that this way of constructing the GSSM does not provide as good source separation performance as the one presented before by [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF].

B. Proposed source variance fitting with GSSM and mixed group sparsity constraint

As the GSSM is constructed to guide the NMF-based source variance constraint, we propose two fitting strategies as follows:

1) Source variance denoising: Motivated by the source variance model [START_REF] Duong | Temporal annotationbased audio source separation using weighted nonnegative matrix factorization[END_REF], when exploiting the GSSM model we propose a variant as

v j (n, f ) = Pj k=1 u jf k hjkn , (21)
where u jf k is an entry of U j , hjkn is an entry of the corresponding activation matrix

H j ∈ R Pj ×N +
. This leads to a straightforward extension of the conventional optimization criterion described by [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF] where H j is now estimated by optimizing the criterion:

min Hj ≥0 D(V j U j H j ) + λΩ( H j ), (22)
where U j is constructed by [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF] or [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF] and fixed, Ω( H j ) presents a penalty function imposing sparsity on H j , and λ is a trade-off parameter determining the contribution of the penalty. Note that as the GSSM U j constructed in (18) becomes a large matrix when the number of examples L j for each source increases, and it is actually a redundant dictionary since different examples may share similar spectral patterns. Thus to fit the source variances with the GSSM, sparsity constraint is naturally needed in order to activate only a subset of U j which represents the spectral characteristics of the sources in the mixture [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparse criteria[END_REF]- [START_REF] Hurmalainen | Group sparsity for speaker identity discrimination in factorisation-based speech recognition[END_REF].

2) Source variance separation: We propose another source variance model as

v(n, f ) = K k=1 u f k hkn , (23)
where v(n, f ) = J j=1 v j (n, f ), u f k is an entry of the GSSM model U constructed as [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and fixed, K = J j=1 P j . Under this model, let V = J j=1 V j be the matrix of the total source variance estimate, it is then decomposed by solving the following optimization problem

min H≥0 D( V U H) + λΩ( H) (24)
where Ω( H) presents a penalty function imposing sparsity on the activation matrix H = [ H 1 , ..., H J ] ∈ R K×N + the total number of rows in H. This criterion can be seen as an additional NMF-based separation step applied on the source variances, while criterion [START_REF] Uhlich | Deep neural network based instrument extraction from music[END_REF] and other existing works [START_REF] Fakhry | Audio source separation in reverberant environments using beta-divergence based nonnegative factorization[END_REF]- [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF] do not perform any additional separation of the variances, but more like denoising of the already separated variances. For the sake of simplicity, in the following, we only present the algorithm derivation for the criterion [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF], but a strong synergy can be found for the criterion [START_REF] Uhlich | Deep neural network based instrument extraction from music[END_REF].

Recent works in audio source separation have considered two penalty functions, namely block sparsity-inducing penalty [START_REF] Sun | Universal speech models for speaker independent single channel source separation[END_REF] and component sparsity-inducing penalty [START_REF] Badawy | On-the-fly audio source separation[END_REF]. The former one enforces the activation of relevant examples only while omitting irrelevant ones since their corresponding activation block in H will likely converge to zero. The latter one, on the other hand, enforces the activation of relevant components in U only. It is motivated by the fact that only a part of the spectral model learned from an example may fit well with the targeted source in the mixture, while the remaining components in the model do not. Thus instead of activating the whole block, the component sparsity-inducing penalty allows selecting only the more likely relevant spectral components from U. Inspired by the advantage of these penalty functions, in our recent work we proposed to combine them in a more general form as [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF] Ω( H) = γ

P p=1 log( + H p 1 ) + (1 -γ) K k=1 log( + h k 1 ), (25)
where the first term on the right hand side of the equation presents the block sparsity-inducing penalty, the second term presents the component sparsity-inducing penalty, and γ ∈ [0, 1] weights the contribution of each term. In [START_REF] Wang | Multi-channel deep clustering: Discriminative spectral and spatial embeddings for speakerindependent speech separation[END_REF],

h k ∈ R 1×N
+ is a row (or component) of H, H p is a subset of H representing the activation coefficients for p-th block, P is the total number of blocks, is a non-zero constant, and . 1 denotes 1 -norm operator. In the considered setting, a block represents one training example for a source and P is the total number of used examples (i.e., P = J j=1 L j ). By putting [START_REF] Wang | Multi-channel deep clustering: Discriminative spectral and spatial embeddings for speakerindependent speech separation[END_REF] into [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF], we now have a complete criterion for estimating the activation matrix H given V and the pretrained spectral model U. The derivation of MU rule for updating H is presented in the Appendix.

C. Proposed multichannel algorithm

Within the LGM, a generalized EM algorithm used to estimate the parameters {v j (n, f ), R j (f )} j,n,f by considering the set of hidden STFT coeffients of all the source images {c j (n, f )} n,f as the complete data. The overview for the GEM derivation are presented in Section II-C, and more details can be found in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF].

For the proposed approach as far as the GSSM concerned, the E-step of the algorithm remains the same as in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF]. In the M-step, we additionally perform the optimization defined either by [START_REF] Uhlich | Deep neural network based instrument extraction from music[END_REF] (for source variance denoising) or by [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF] (for source variance separation). This is done by the MU rules so that the estimated intermediate source variances v j (n, f ) are further updated with the supervision of the GSSM. The detail of overall proposed algorithm with source variance separation is summarized in Algorithm 1.

Note that this generalized EM algorithm requires the same order of computation compared to the existing method [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF] as sparsity constraint and bigger GSSM size does not significantly affect the overall computational time. As an example, for separating a 10-second long mixture presented in our experiment, both [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF] and our proposed method (when non-optimally implemented in Matlab) take about 400 seconds when running in a laptop with Intel Core i5 Processor, 2.2 GHz, and 8 GB RAM.

IV. EXPERIMENTS

A. Dataset and parameter settings

We validated the performance of the proposed approach in an important speech enhancement use case where we know already two types of sources in the mixture: speech and noise. For a better comparison with the state of the art, we used the benchmark development dataset of the "Two-channel mixtures of speech and real-world background noise" (BGN) task1 within the SiSEC 2016 [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF]. This devset contains stereo mixtures of 10-second duration and 16 kHz sampling rate. They were the mixture of male/female speeches and real-world noises recorded from different public environments: cafeteria (Ca), square (Sq), and subway (Su). Overall there were nine mixtures: three with Ca noise, four with Sq noise, and two with Su noise. The signal-to-noise ratio was drawn randomly between -17 and +12 dB by the dataset creators. Our works in single-channel case [START_REF] Duong | Single-channel speaker-dependent speech enhancement exploiting generic noise model learned by non-negative matrix factorization[END_REF], [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF] and preliminary tests on multichannel case show that only a few examples for each source could be enough to train an efficient GSSM. Thus, for training the generic speech spectral model, we took only one male voice and two female voices from the SiSEC 2015 2 . These three speech examples are also 10-second length. We performed the listening check to confirm that these examples used for the speech and noise model training are different from those in the devset, which were used for testing. For training the generic noise spectral model, we extracted five noise examples from the Diverse Environments Multichannel Acoustic Noise Database (DEMAND) 3 . Again they were 10second length and contained three types of environmental noise: cafeteria, square, metro. The STFT window length was 1024 for all train and test files. The number of NMF components in W l j for each speech example was set to 32, while that for noise example was 16. These values were found to be reasonable in [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF] and our work on single-channel case [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF]. Each W l j were obtained by optimizing (17) with 20 MU iterations.

Initialization of the spatial covariance matrices: As suggested in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF], we firstly tried to initialize the spatial covariance matrix R j (f ) by performing hierarchical clustering on the mixture STFT coefficients x(n, f ). But this strategy did not give us a good separation performance as the noise source in the considered mixtures is diffuse (i.e., it does not 2 https://sisec.inria.fr/sisec-2015/2015-underdetermined-speech-and-musicmixtures/.

3 http://parole.loria.fr/DEMAND/.

come from a single direction). Thus we initialized the noise spatial covariance matrix based on the diffuse model where noise is assumed to come uniformly from all spatial directions. With this assumption, the diagonal entries of the noise spatial covariance matrix are one and the off-diagonal entries are realvalued computed as in [START_REF] Kuttruff | Room Acoustics[END_REF]

r 1,2 (f ) = r 2,1 (f ) = sin(2πf d/v) 2πf d/v , ( 26
)
where d is the distance between two microphones and v = 334 m/s the sound velocity. The spatial covariance matrix for the speech source was initialized by the full-rank di-rect+diffuse model detailed in [START_REF] Duong | Under-determined reverberant audio source separation using a full-rank spatial covariance model[END_REF] where the speech's direction of arrival (DoA) was set to 90 degrees. This DoA initialization was chosen for balancing the fact that the speech direction can vary between 0 degree and 180 degrees in each mixture and we did not have access to the ground truth information while performing the test. The source separation performance for all approaches was evaluated by two sets of criteria. The four power-based criteria: the signal to distortion ratio (SDR), the signal to interference ratio (SIR), the signal to artifacts ratio (SAR), and the source image to spatial distortion ratio (ISR), measured in dB where the higher the better [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. The four perceptually-motivated criteria: the overall perceptual score (OPS), the target-related perceptual score (TPS), the artifact-related perceptual score (APS), and the interference-related perceptual score (IPS) [START_REF] Emiya | Subjective and objective quality assessment of audio source separation[END_REF], where a higher score is better. As power-based criteria are more widely used in source separation community, the hyper-parameters for each algorithm were chosen in order to

B. Algorithm analysis 1) Algorithm convergence: separation results as functions of EM and MU iterations:

We first investigate the convergence in term of separation performance of the derived Algorithm 1 by varying the number of EM and MU iterations and computing the separation results obtained on the benchmark BGN dataset. In this experiment, we set λ = 10 and γ = 0.2 as we will show in next section that these values offer both the stability and the good separation performance. The speech and noise separation results, measured by the SDR, SIR, SAR, and ISR, averaged over all mixtures in the dataset, illustrated as functions of the EM and MU iterations, are shown in Fig. 2.

As it can be seen, generally the SDR increases when the number of EM and MU iterations increases. With 10 or 25 MU iterations, the algorithm converges nicely and saturates after about 10 EM iterations. The best separation performance was obtained with 10 MU iterations and 15 EM iterations. It is also interesting to see that with a small number of MU iterations like 1, 2, or 3, the separation results are quite poor and the algorithm is less stable as it varies significantly even with a large number of EM iterations. This reveals the effectiveness of the proposed NMF constraint [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF].

2) Separation results with different choices of λ and γ: We further investigate the sensitivity of the proposed algorithm to two parameters λ and γ, which determine the contribution of sparsity penalty to the NMF constraint in [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF]. For this purpose, we varied the values of these parameters, λ = {1, 10, 25, 50, 100, 200, 500}, γ = {0, 0.2, 0.4, 0.6, 0.8, 1}, and applied the corresponding source separation algorithm presented in the Algorithm 1 on the benchmark BGN dataset. The number of EM and MU iterations are set to 15 and 10, respectively, as these values guarantee the algorithm's convergence shown in Fig. 2. The speech and noise separation results, measured by the SDR, SIR, SAR, and ISR, averaged over all mixtures in the dataset, represented as functions of λ and γ, are shown in Fig. 3.

It can be seen that the proposed algorithm is less sensitive to the choice of γ, while more sensitive to the choice of λ, and the separation performance greatly decreases with λ > 10. The best choice for these parameters in term of the SDR are λ = 10, γ = 0.2. With the small value of λ (e.g., λ = 1), varying γ does not really affect the separation performance as the evaluation criteria are quite stable. We noted that with γ = 0.2, the algorithm offers 0.2 dB and 1.0 dB SDR, which are higher than when γ = 0 and γ = 1, respectively. This confirms the effectiveness of the mixed sparsity penalty [START_REF] Wang | Multi-channel deep clustering: Discriminative spectral and spatial embeddings for speakerindependent speech separation[END_REF] in the multichannel setting.

C. Comparison with the state of the art

We compare the speech separation performance obtained on the BGN dataset of the proposed approach with its close prior art (i.e. Arberet's algorithm [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF]) and other state-of-theart methods presented at the SiSEC campaign over different years since 2013. The results of these methods were submitted by the authors and evaluated by the SiSEC organizers [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF], [START_REF] Ono | The 2013 Signal Separation Evaluation Campaign[END_REF], [START_REF] Ono | The 2015 Signal Separation Evaluation Campaign[END_REF]. All comparing methods are summarized as follows:

• Martinez-Munoz's method (SiSEC 2013) [START_REF] Ono | The 2013 Signal Separation Evaluation Campaign[END_REF]: this algorithm exploits source-filter model for the speech source and the noise source is modeled as a combination of pseudo-stationary broadband noise, impulsive noise, and

Algorithm 1 Proposed GSSM + SV separation algorithm Require: Mixture signal x(t)
List of examples of each source in the mixture {s l j (t)} j=1:J,l=1:Lj Hyper-parameters λ, γ, MU-iteration Ensure: Source images ĉj (t) separated from x(t) -Compute the mixture STFT coefficients x(n, f ) ∈ C F ×N and then Ψ x (n, f ) ∈ C I×I by ( 7) -Construct the GSSM model U j by [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF], then U ∈ R F ×K + by ( 19) -Initialize the spatial covariance matrices R j (f ), ∀j, f (see Section IV-A) -Initialize the non-negative time activation matrix for each source H j randomly, then

H = [ H 1 , . . . , H J ] v ∈ R K×N + -Initialize the source variance v j (n, f ) = [U j H j ] n,f
// Generalized EM algorithm for the parameter estimation: repeat // E step (perform calculation for all j, n, f ): 11) // M step: updating spatial covariance matrix and unconstrained source spectra

Σ j (n, f ) = v j (n, f )R j (f ) // eq. (4) Σ x (n, f ) = J j=1 v j (n, f ) R j (f ) // eq. (5) G j (n, f ) = Σ j (n, f )Σ -1 x (n, f ) // Wiener gain Σ j (n, f ) = G j (n, f ) Ψ x (n, f )G H j (n, f ) + (I - G j (n, f ))Σ j (n, f ) // eq. (
R j (f ) = 1 N N n=1 1 vj (n,f ) Σ j (n, f ) // eq. (13) v j (n, f ) = 1 I tr(R -1 j (f ) Σ j (n, f
)) // eq. ( 14)

V j = {v j (n, f )} n,f V = J j=1 V j
// MU rules for NMF inside M step to further constrain source spectra by the GSSM for iter = 1, ..., MU-iteration do for p = 1, ..., P do pitched interferences. The parameter estimation is based on the MU rules employed in non-negative matrix factorization.

Y p ← 1 + Hp 1 end for Y = [Y 1 , . . . , Y P ] for k = 1, ..., K do z k ← 1 + h k 1 end for Z = [z 1 , . . . , z K ] // Updating activation matrix V = U H H ← H U ( V V .-2 ) U ( V .-1 )+λ(γY+(1-γ)Z) . 1 2 // eq. (31) end for v j (n, f ) = [U j H j ] n,f //
• Wang's method [START_REF] Wang | A region-growing permutation alignment approach in frequency-domain blind source separation of speech mixtures[END_REF] (SiSEC 2013): this algorithm performs well-known frequency domain independent component analysis (ICA). The associated permutation problem is solved by a novel region-growing permutation alignment technique. • Le Magoarou's method [START_REF] Magoarou | Text-informed audio source separation. example-based approach using non-negative matrix partial co-factorization[END_REF] (SiSEC 2013): this approach uses text transcript of the speech source in the mixture as prior information to guide the source separation process.

The algorithm is based on the nonnegative matrix partial co-factorization. • Bryan's method [START_REF] Bryan | An efficient posterior regularized latent variable model for interactive sound source separation[END_REF] (SiSEC 2013): this interactive approach exploits human annotation on the mixture spectrogram to guide and refine the source separation process.

The modeling is based on the probabilistic latent component analysis (PLCA), which is equivalent to NMF. • Rafii's method [START_REF] Rafii | Online REPET-SIM for real-time speech enhancement[END_REF] (SiSEC 2013): this technique uses a similarity matrix to separate the repeating background from the non-repeating foreground in a mixture. The underlying assumption is that the background is dense and low-ranked, while the foreground is sparse and varied. • Ito's method [START_REF] Ito | Permutation-free convolutive blind source separation via full-band clustering based on frequencyindependent source presence priors[END_REF] (SiSEC 2015): this is a permutationfree frequency-domain blind source separation algorithm via full-band clustering of the time-frequency (T-F) components. The clustering is performed via MAP estimation of the parameters with EM algorithm. • Liu's method [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF] (SiSEC 2016): the algorithm performs Time Difference of Arrival (TDOA) clustering based on GCC-PHAT. • Wood's method [START_REF] Wood | Blind Speech Separation and Enhancement With GCC-NMF[END_REF] (SiSEC 2016): this recently proposed algorithm first applies NMF to the magnitude spectrograms of the mixtures with channels concatenated in time. Each dictionary atom is clustered to either the speech or the noise according to its spatial origin. • Arberet's method [START_REF] Arberet | Nonnegative matrix factorization and spatial covariance model for under-determined reverberant audio source separation[END_REF]: using the similar local Gaussian model, the algorithm further constrains the intermediate source variances by unsupervised NMF with criterion [START_REF] Badawy | On-the-fly audio source separation -a novel user-friendly framework[END_REF]. Such algorithm is implemented by Ozerov et. al. in [37]. This method is actually the most relevant prior art to compare with as it falls in the same LGM framework.

The proposed approach with different variants are summarized as:

• GSSM + SV denoising: The proposed GSSM + full-rank spatial covariance approach where the estimated variances of each sources V j are further constrained by criterion [START_REF] Uhlich | Deep neural network based instrument extraction from music[END_REF]. We submitted results obtained by this method to the SiSEC 2016 BGN task and obtained the best performance over the actual test set in term of SDR [START_REF] Liutkus | The 2016 signal separation evaluation campaign[END_REF]. • GSSM + SV separation: The proposed approach with source variance separation by optimizing criterion [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF]. In order to investigate the benefit of the sparsity constraint, we further report the results obtained by this method when λ = 0. Finally, to confirm the effectiveness of the GSSM construction by [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF], we report the results obtained when the GSSM of the same size is learned jointly by concatenating all example's spectrograms S l j as [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF]. In this case, only the component sparsity is applied (i.e., γ = 0) as block does not exist. This setting is named "GSSM'+component sparsity" in Table 1. The separation results obtained by different methods for each noisy environment (Ca, Sq, Su), and the average overall mixtures are summarized in Table 1. The boxplot to illustrate the variance of the results obtained by the two proposed approaches is shown in Fig. 4. It is interesting to see that the results obtained by the proposed approach without sparsity constraint were lower than that of Arberet's method for all noisy environments, even the former used the pre-trained GSSM while the latter is completely unsupervised. It reveals that the GSSM itself is redundant and contains some irrelevant spectral patterns with the actual sources in the mixture. Thus constraining the source variances by the GSSM without a relevant spectral pattern selection guided by the sparsity penalty is even worse than the unsupervised NMF case where the spectral patterns are randomly initialized and then updated by MU rules. The importance of such sparsity penalty is explicitly confirmed by the fact that the results obtained by the proposed approach with sparsity constraint are far better than the setting without the sparsity constraint. Also, it is not surprising to see that the "GSSM + SV denoising" clearly outperforms Arberet's method (except for the ISR and the TPS) in all noisy environments as the former exploits additional information about the types of sources in the mixtures in order to learn the GSSM in advance. The "GSSM + SV separation" offers better separation performance in terms of SDR, SIR, OPS, IPS, on square and subway environments as well as on average compared to the "GSSM + SV denoising" and the "GSSM' + component sparsity". This confirms the effectiveness of the proposed source variance separation criterion [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF] and the GSSM construction [START_REF] Duong | Speech enhancement based on nonnegative matrix factorization with mixed group sparsity constraint[END_REF].

When compared to the top-performing state-of-the-art methods in the SiSEC campaigns, the proposed approach performs generally better in terms of the energy-based criteria but worse for the perceptually-motivated ones. Especially in Ca environment the OPS obtained by the proposed approach is far below those offered by other methods. This may be due to the fact that the hyper-parameters were optimized for the SDR, but not the OPS. The "GSSM + SV separation" with sparsity constraint outperforms all other methods, but Wang's approach, in terms of the SDR, the most important energybased criterion, at all noisy environment. This confirms the effectiveness of the proposed approach where the GSSM is successfully exploited in the LGM framework. It should be noted that Wang's method [START_REF] Wang | A region-growing permutation alignment approach in frequency-domain blind source separation of speech mixtures[END_REF] is based on the frequencydomain ICA so it is not applicable for under-determined mixtures where the number of sources is larger than the number of channels. Also, in this method, an additional postfiltering technique was applied to the separated speech source so as to maximize the denoising capability.

V. CONCLUSION

In this paper, we have presented a novel multichannel audio source separation algorithm weakly guided by some source examples. The considered approach exploits the use of generic source spectral model learned by NMF within the well-established local Gaussian model. In particular, we have proposed a new source variance separation criterion in order to better constrain the intermediate source variances estimated in each EM iteration. Experiments with the benchmark dataset from the SiSEC campaigns have confirmed the effectiveness of the proposed approach compared to the state of the art. Motivated by the effectiveness of the GSSM, future work can be devoted to extending the current approach in order to exploit in addition the use of a generic spatial covariance model, which remains to be defined. In addition, the theoretical grounding of the source variance separation criterion needs to be further investigated. Another promising investigation could be extending the idea of source variance separation to DNNbased models inspired by the work of Nugraha et al. [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF].

VI. ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers of this manuscript and [START_REF] Duong | Multichannel audio source separation exploiting NMF-based generic source spectral model in Gaussian modeling framework[END_REF] for their helpful and constructive comments that greatly contributed to improving the quality of the paper. We also would like to thank Professor Nobutaka Ono for providing us results of the SiSEC 2013 and the SiSEC 2015.

APPENDIX DERIVATION OF MU RULE FOR UPDATING THE ACTIVATION MATRIX IN ALGORITHM1

Let L( H) denote the minimization criterion [START_REF] Nugraha | Multichannel audio source separation with deep neural networks[END_REF] with the mixed sparsity constrained Ω( H) defined as in [START_REF] Wang | Multi-channel deep clustering: Discriminative spectral and spatial embeddings for speakerindependent speech separation[END_REF] and D(• •) being IS divergence. The partial derivative of L( H) with respect to an entry h kn is This ∇ h kn L( H) can be written as a sum of two nonnegative parts, denoted by ∇ + h kn L( H) ≥ 0 and ∇ - h kn L( H) ≥ 0, respectively, as

∇ h kn L( H) = F f =1 u f k 1 [U H] n,f - v(n, f ) [U H]
∇ h kn L( H) = ∇ + h kn L( H) -∇ - h kn L( H) (28)
with

∇ + h kn L( H) F f =1 u f k 1 [U H] n,f + λ.γ + H p 1 + λ.(1 -γ) + h k 1 , ∇ - h kn L( H) F f =1 u f k v(n, f ) [U H] 2 n,f . (29)
Following a standard approach for MU rule derivation [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Févotte | Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis[END_REF]), h kn is updated as

h kn ← h kn ∇ - h kn L( H) ∇ + h kn L( H) .η , (30)
where η = 0.5 following the derivation in [START_REF] Lefèvre | Itakura-Saito non-negative matrix factorization with group sparsity[END_REF], [START_REF] Févotte | Algorithms for nonnegative matrix factorization with the beta-divergence[END_REF], which was shown to produce an accelerated descent algorithm. Putting (29) into [START_REF] Winter | MAP-based underdetermined blind source separation of convolutive mixtures by hierarchical clustering and 1 -norm minimization[END_REF] and rewriting it in a matrix form, we obtain the updates of H as

H ← H U ( V V .-2 )
U ( V .-1 ) + λ(γY + (1 -γ)Z)

Fig. 1 .

 1 Fig. 1. General workflow of the proposed source separation approach. The top green dashed box describes the training phase for the GSSM construction. Bottom blue boxes indicate processing steps for source separation. Green dashed boxes indicate the novelty compared to the existing works [36]-[38].

  , 5 examples in total), and examples of different types of noise such as those from outdoor environment, cafeteria, waterfall, street, etc.,. (e.g., 6 examples in total). Note that as another variant investigated in this work, the GSSM U j can be constructed differently by first concatenating all examples for each source (S j = [S 1 j , . . . , S Lj j ]), and then performing NMF on the concatenated spectrogram only once by optimizing the criterion min Hj ≥0,Uj ≥0 D(S j U j H j ).

Fig. 2 .

 2 Fig. 2. Average separation performance obtained by the proposed method over stereo mixtures of speech and noise as functions of EM and MU iterations. (a): speech SDR, (b): speech SIR, (c): speech SAR, (d): speech ISR, (e): noise SDR, (f): noise SIR, (g): noise SAR, (h): noise ISR

Fig. 3 .

 3 Fig. 3. Average separation performance obtained by the proposed method over stereo mixtures of speech and noise as functions of λ and γ. (a): speech SDR, (b): speech SIR, (c): speech SAR, (d): speech ISR, (e): noise SDR, (f): noise SIR, (g): noise SAR, (h): noise ISR

  updating constrained spectra until convergence -Source separation by multichannel Wiener filtering (9) -Time domain source images ĉj (t) are obtained by the inverse STFT of ĉj (n, f ).

Fig. 4 .

 4 Fig.4. Boxplot for the speech separation performance obtained by the proposed "GSSM + SV denoising" (P1) and "GSSM + SV separation" (P2) methods.

. 1 2 , 1 + 1 +

 211 [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] where V = U H, Y = [Y 1 , . . . , Y P ] with Y p , p = 1, . . . P an uniform matrix of the same size as H p whose entries are Hp 1 , and Z = [z 1 , . . . , z K ] with z k , k = 1, . . . K a uniform vector of the same size as h k whose entries are h k 1 .

https://sisec.inria.fr/sisec-2016/bgn-2016/


Cornel Theodor Hangan
email: theodor.hangan@uha.fr

Marius Murea

Tewfik Sari
email: tewfik.sari@uha.fr

Poleni curves on surfaces of constant curvature *

Keywords:

. The above equation is yet meaningful on a surface if one interprets k as the geodesic curvature of the curve. In this paper we solve the above equation on a surface of constant curvature.

Elastic Poleni curves

In [START_REF] Langer | The total squared curvature of closed curves[END_REF] the authors show that on surfaces of Gaussian curvature G, elastic curves γ are solutions of the intrinsic equation

k ′′ g + 1 2 k 3 g + k g (G -λ) = 0,
where k g denotes the geodesic curvature of γ, λ is a constant and the primes indicate derivation with respect to the arc-length parameter s on γ. As the function f (s) = 2 cosh(s) satisfies the equation

f ′′ + 1 2 f 3 -f = 0,
there exist on a surface S of constant curvature elastic curves with intrinsic equation

(1) k g (s) = 2 cosh (s) .

If S is the euclidean plane, the solution of ( 1) is Giovanni Poleni's curve introduced in 1729 in relation with the tractrix [START_REF] Poleni | Epistolarum Mathematicarum fasciculus. Patavii 1729. See also Revue du palais de la Découverte, courbes mathématiques[END_REF]. Its parametric equations are [START_REF] Brunnett | Elastic curves on the sphere[END_REF] x (s) = s -2 tanh (s) , y (s) = 2 cosh (s) . satisfies the differential equation

f ′′ + 1 2 f 3 -α 2 f = 0.
Therefore there exist on surface S elastic curves with intrinsic equation

(3) k g (s) = 2α cosh (αs) .

If S is the euclidean plane, equation (3) writes

x ′′ = -2αy ′ cosh(αs) , y ′′ = 2αx ′ cosh (αs) .

Integrating these equations, with initial condition

x(0) = 0, x ′ (0) = 1, y(0) = 2 α , y ′ (0) = 0, one finds x(s) = s -2 α tanh(αs), y(s) = 2 α cosh (αs) .

Notice that the curve (αx (s/α) , αy (s/α)) is simply the Poleni curve [START_REF] Brunnett | Elastic curves on the sphere[END_REF], so that the study of equation [START_REF] Euler | Additamentum De curvis elasticis, in the Methods inveniendi lineas curvas maximi minimivi proprietate gaudentes[END_REF] in the euclidean plane reduces to the study of equation [START_REF] Bryant | Reduction for constrained variational problems and 1 2 k 2 ds[END_REF]. For the sphere and the hyperbolic space, there are no dilatations and equation ( 3) deserves a special study. In the following, we call "Poleni curves" the solutions of the intrinsic equation [START_REF] Euler | Additamentum De curvis elasticis, in the Methods inveniendi lineas curvas maximi minimivi proprietate gaudentes[END_REF]. Geodesics on cylinders having as directrix the curve (2) are extremals of Sadowsky's functional

S = γ κ 2 1 + τ 2 κ 2 2 ds
where κ, τ represent the curvature and the torsion of the curve γ, in the threedimensional euclidean space E 3 , see [START_REF] Th | Elastic helices[END_REF]. An analogous result is valid for geodesics on conuses having as directrix elastic Poleni curves on the sphere S 2 , see [START_REF] Th | Poleni curves, elastic curves and elastic strips[END_REF]. This example motivates our research. Notice that the Poleni curve (2) in the euclidean plane has an axe of symmetry and it has an asymptote as s → ±∞. In fact, on any surface, a solution of equation ( 3) is asymptotically geodesic and has an axis of symmetry. If S is the sphere, the Poleni curve corresponding to α = 1, is (see Proposition 1):

x = 2 cos(s) tanh(s/2), y = 2 sin(s) tanh(s/2),
where (x, y) are the stereographical coordinates. It is a symmetric double spiral which tends to the circle of radius 2 centered at the origin. Hence (Theorem 2), on the sphere, the situation is similar to the euclidean case: there exist two geodesics (here circles) which are characterized by the fact that the first one is an axis of symmetry of the curve and the second one is asymptotic to the curve. If S is the hyperbolic half plane, the Poleni curve corresponding to α = 1 is (see Proposition 3):

x = s, y = cosh(s).

Notice that this is the catenary (or chainette) curve of the euclidean space. Hence (Theorem 4), on the hyperbolic space, the Poleni curve is asymptotic to all geodesics that intersect it transversally at a unique point and among these geodesics there is one which is an axis of symmetry of the curve. The study of elastic curves goes back to Euler in 1744 [START_REF] Euler | Additamentum De curvis elasticis, in the Methods inveniendi lineas curvas maximi minimivi proprietate gaudentes[END_REF]. His classification of plane elastic curves can be found in Love's book (see [START_REF] Love | A treatise on the mathematical theory of elesticity[END_REF], Section 263) or in the historical survey by Truesdell [START_REF] Truesdell | The influence of elasticity on analysis : the classical heritage[END_REF]. This study was the birth of the Euler Lagrange calculus of variation. In Euler's classification of plane elastic curves, Poleni's curve, called also "la courbe des forçats", occupies a special place: discarding the circle and the straight line, it is the only elastic plane curve that does not need elliptic functions for its parametrization.

The study of euclidean elastic curves (together with the extension to the non-euclidean case) was revisited by Langer and Singer [START_REF] Langer | The total squared curvature of closed curves[END_REF] and also by Bryant and Griffiths [START_REF] Bryant | Reduction for constrained variational problems and 1 2 k 2 ds[END_REF], using Griffiths's formalism based on the theory of exterior differential systems, and by Jurdjevic [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], using Pontryagin's Maximum Principle of optimal control theory. The papers [START_REF] Bryant | Reduction for constrained variational problems and 1 2 k 2 ds[END_REF][START_REF] Jurdjevic | Non-Euclidean elastica[END_REF][START_REF] Langer | The total squared curvature of closed curves[END_REF] contain important results on the classification of elastic curves in the non euclidean case. Even though they mentioned the special case of asymptotically geodesic elastic curves (see especially Figure 8 of [START_REF] Bryant | Reduction for constrained variational problems and 1 2 k 2 ds[END_REF] and Figure 2.c of [START_REF] Langer | The total squared curvature of closed curves[END_REF]), they all don't give explicit formula for the Poleni curve. On the sphere S 2 , the Poleni curve is explicitly mentioned by Brunnett and Crouch, see Figure 7 in [START_REF] Brunnett | Elastic curves on the sphere[END_REF], but its double spiral behavior is not apparent in their figure .  This paper is organized as follows. In Section 2 we recall some facts on the riemannian geometry of surfaces. In Section 3, equation (3) of Poleni curves is expressed in the geographical coordinates of the sphere S 2 and numerical solutions of the obtained equation are given. In Section 4, equation (3) is expressed in the stereographical projection of the sphere and the obtained equation is explicitly solved in the case α = 1. In Section 5, equation ( 3) is expressed in the hyperbolic half plane and the obtained equation is explicitly solved in the case α = 1. In Section 6 we compare our results with the parameterizations given by Jurdjevic in [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF]. In Section 7, the Poleni curves of geodesic curvature f (s) = 2α cosh(αs) are considered. Discarding the constant curvature curves, these curves are the only elastic curves that can be parametrized by elementary functions.

Geodesic curvatures and intrinsic equations

On the surface S, the geodesic curvature k g of the curve γ represents the modulus of the normal component of the acceleration vector A := ∇ γ ′ γ ′ . Clearly, s being an arc-length parameters on γ, let γ ′ = dγ/ds be the unit tangent vector field along γ. Denoting by g (•, •) the scalar product on the tangent bundle of the surface S, it follows that

γ ′ 2 := g (γ ′ , γ ′ ) = 1.
Let ∇ denotes the operator of covariant (absolute) differentiation on the surface. Then

g (∇ γ ′ γ ′ , γ ′ ) = 1 2 ∇ γ ′ g (γ ′ , γ ′ ) = 0, so the acceleration vector A = ∇ γ ′ γ ′ is orthogonal to γ ′ along the curve. Equa- tion (3) of Poleni curves becomes therefore (4) ∇ γ ′ γ ′ = 2α cosh (αs)
.

If along γ, an orthonormal frame field {T = γ ′ , N } is chosen, in order to solve (4) one has to integrate the system (5)

∇ γ ′ γ ′ = 2α cosh (αs) N.
Following [START_REF] Vranceanu | Leçons de Géométrie Différentielle[END_REF], we recall some general formulas from the theory of surfaces when the local coordinates x 1 , x 2 are used. Let (6)

ds 2 = 2 i,j=1 g ij dx i dx j = g 11 dx 2 1 + 2g 12 dx 1 dx 2 + g 22 dx 2 2
be the riemannian metric on the surface S and let

(7) L = 1 2 2 i,j=1
g ij dx i dx j be the Lagrange function associated to [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF]. The Euler-Lagrange system corre-

sponding to L is (8) d ds ∂L ∂x ′ i - ∂L ∂x i = 0, i = 1, 2
and we use it to deduce the components of the acceleration vector A = ∇ γ ′ γ ′ . This is done by isolating second order derivatives of the coordinates x i in (8) i.e. transforming it in (9)

x ′′ i + 2 j,k=1 Γ i jk x ′ j x ′ k = 0, i = 1, 2.
Now the left side part of (9) represents the component A i of the acceleration vector A with respect to x i , for i = 1, 2. In the sequel, we apply this procedure for several curvilinear systems of coordinates on surfaces of constant curvature.

3 Geographical coordinates on the sphere S 2 Denote φ and θ the longitude and the colatitude of a point on the unit sphere S 2 . The cartesian coordinates of the point are

x = sin θ cos φ, y = sin θ sin φ, z = cos θ, 0 < θ < π
and the riemannian metric of S 2 becomes ( 10)

ds 2 = dθ 2 + sin 2 θdφ 2 .
The Lagrange function ( 7) associated to ds 2 is

(11) L = 1 2 (θ ′ ) 2 + sin 2 θ (φ ′ ) 2
and the Euler-Lagrange equations write

d ds θ ′ -sin θ cos θ (φ ′ ) 2 = 0, d ds sin 2 θφ ′ = 0
so that, transforming in the form [START_REF] Poleni | Epistolarum Mathematicarum fasciculus. Patavii 1729. See also Revue du palais de la Découverte, courbes mathématiques[END_REF] gives

θ ′′ -sin θ cos θ (φ ′ ) 2 = 0, φ ′′ + 2 cos θ sin θ θ ′ φ ′ = 0.
The left side parts of these equations represent the components A θ , A φ of the acceleration vector A of a curve γ along which θ and φ are functions of arc-length parameter s. The vector fields along the curve, T (θ ′ , φ ′ ), N (-φ ′ sin θ, θ ′ / sin θ) are unitary and orthonormal as

(12) g (T, T ) = γ ′ 2 = (θ ′ ) 2 + sin 2 θ (φ ′ ) 2 = 1 and g (N, N ) = (-φ ′ sin θ) 2 + sin 2 θ (θ ′ / sin θ) 2 = sin 2 θ (φ ′ ) 2 + (θ ′ ) 2 = 1 g (T, N ) = θ ′ (-φ ′ sin θ) + sin 2 θφ ′ (θ ′ / sin θ) = 0.
It follows that the system of differential equations

θ ′′ -sin (θ) cos (θ) (φ ′ ) 2 = - 2α cosh (αs) φ ′ sin (θ) (13) φ ′′ + 2 cos (θ) sin (θ) θ ′ φ ′ = 2α cosh(αs) θ ′ sin (θ) (14)
translates the vector equation [START_REF] Th | Elastic helices[END_REF].

V V O Figure 2:
The curve of geographic coordinates (φ i , θ i ) (left) and its projection on the symmetry plane (right). The asymptotic plane is orthogonal to OV We have solved numerically (13)-( 14), in the case α = 1, together with the conditions θ (0) = π/4, θ ′ (0) = 0, φ (0) = π/4 and φ ′ (0) = -√ 2 by a Runge-Kutta algorithm of fourth order using the symbolic computing system MAPLE. We have obtained θ i , θ ′ i , φ i , φ ′ i approximation of θ(i∆s), θ ′ (i∆s), φ(i∆s), φ ′ (i∆s) respectively.

Let V be the peak of the curve i.e. the point with maximal curvature 2 (for s = 0) and denote by t V the tangent at the solution curve in V . We observe on Figure 2 that the plane orthogonal to t V in V is a plane of symmetry for the curve. It contains a first self-intersecting point. Despite the curve seems to be closed, in fact it is composed by two symmetrical spirals which wind around a diameter of S 2 and these spirals tend asymptotically to the equatorial circle corresponding to V , see Theorem 2.

At this stage, the computed φ i , θ i , φ ′ i and θ ′ i verify the equations ( 13) and ( 14). The equation ( 12) has not been used yet. The computed quantities satisfy

-2 • 10 -6 ≤ (θ ′ i ) 2 + sin 2 θ i (φ ′ i ) 2 -1 ≤ 6 • 10 -6 ,
in other words, the computed φ i , θ i , φ ′ i and θ ′ i verify the three equations ( 12), ( 13) and ( 14).

Stereographical projection

If one projects the unit sphere S 2 from its north pole on the plane tangent to S 2 at its south pole and one transports its metric on that plane, one gets on the plane the Riemann metric

ds 2 = dx 2 + dy 2 1 + 1 4 (x 2 + y 2 ) 2 .
The geodesics of this metric are the images of the great circles of S 2 , these are the circles or straight lines of the plane which contain two diametrically opposed points of the circle of radius 2 centered on the origin. Starting as before with the Lagrange function

L = 1 2 x ′2 + y ′2 1 + 1 4 (x 2 + y 2 ) 2 ,
the Euler-Lagrange equations are

d ds x ′ 1 + 1 4 (x 2 + y 2 ) 2 + 1 2 x x ′2 + y ′2 1 + 1 4 (x 2 + y 2 ) 3 = 0, d ds y ′ 1 + 1 4 (x 2 + y 2 ) 2 + 1 2 y x ′2 + y ′2 1 + 1 4 (x 2 + y 2 ) 3 = 0.
Isolating the second order derivatives x ′′ , y ′′ one gets the components of the acceleration

A x = x ′′ - x x ′2 -y ′2 + 2yx ′ y ′ 2 + 1 2 (x 2 + y 2 ) , A y = y ′′ - y y ′2 -x ′2 + 2xx ′ y ′ 2 + 1 2 (x 2 + y 2 )
.

The vector fields T (x ′ , y ′ ), N (-y ′ , x ′ ) along the curve γ are unitary and orthogonal as

(15) γ ′ (s) 2 = (x ′ (s)) 2 + (y ′ (s)) 2 1 + 1 4 (x(s)) 2 + (y(s)) 2 2 = 1
and equation ( 5) is expressed by the system (16)

x ′′ - x x ′2 -y ′2 + 2yx ′ y ′ 2 + 1 2 (x 2 + y 2 ) = -2αy ′ cosh(αs) , ( 17
) y ′′ - y y ′2 -x ′2 + 2xx ′ y ′ 2 + 1 2 (x 2 + y 2 ) = 2αx ′ cosh(αs)
.

Proposition 1 Up to isometry, the parametric equations of the Poleni curve with α = 1 on the sphere are (18)

x(s) = 2 cos(s) tanh(s/2), y(s) = 2 sin(s) tanh(s/2).

Proof-Up to isometry, one can restrict to the initial condition

x (0) = 0, x ′ (0) = 1, y (0) = 0, y ′ (0) = 0.
We use the symbolic computing system MAPLE to write Taylor series of the solution of (16,17), with α = 1, satisfying the previous initial conditions. At order 7 one obtains

x(s) = s - 7 12 s 3 + 11 120 s 5 + 0(s 7 ), y(s) = -s 2 + 1 4 s 4 - 11 360 s 6 + 0(s 7 ).
Then we compute

x(s)

x 2 (s) + y 2 (s) = 1 - s 2 2 + s 4 24 + 0(s 6 ) = cos(s), y (s)
x 2 (s) + y 2 (s) = s -s 3 6 + s 5 120 + 0(s 6 ) = sin(s),

x 2 (s) + y 2 (s) 2 = s 2 - s 3 24 + s 5 240 + 0(s 6 ) = tanh(s/2).
Hence formulas (18) are conjectured. By straightforward calculations we verify that ( 18) is actually the solution of (16,17), with α = 1, satisfying the specified initial conditions.

Theorem 2 Let V be the peak of the Poleni curve, with α = 1, on the sphere. Then the plane orthogonal to the Poleni curve at V is a symmetry plane for the curve. The curve has an infinite number of self-intersecting points of arc-length ±s n given by s n = π 2 + nπ, n ∈ N contained in the symmetry plane. The Poleni curve is composed of two symmetrical spirals which wind around the diameter determined by V and tend asymptotically to a circle of radius 1.

Proof-Let us project the unit sphere S 2 from point -V on the plane tangent to S 2 at V . We chose cartesian coordinates x and y on this plane such that t V projects on the vector (1, 0). Then the Poleni curve is given by equations (18). Thus x 2 (s) + y 2 (s) tends to 4. The circle of radius 2 correspond to the circle of radius 1 on the sphere, orthogonal to V . The arc-lengths ±s n of the self-intersecting points are given by equation cos(s n ) = 0. Thus s n = π 2 + nπ.

Figure 3 shows the Poleni curve (18) and its image via the transformation

(19) (x, y) → xf x 2 + y 2 , yf x 2 + y 2 ,
where f (r) = -ln(1-r/2) r = 1 2 + O(r). This transformation is a magnifying glass of the neighborhood of the circle of radius 2. It shows many self-intersection Figure 3: The Poleni curve (18) and its image via the transformation (19), showing the self-intersection points on the y axis. points of the Poleni curve. In the geographic coordinates (φ, θ) the curve (18) is given by φ(s) = s, θ(s) = arccos 1 cosh(s) .

In the cartesian coordinates (x, y, z) such that x 2 + y 2 + z 2 = 1, the Poleni curve has the following parametric equations (20)

x(s) = cos(s) tanh(s), y(s) = sin(s) tanh(s), z(s) = 1 cosh(s) .

The Poincaré model of hyperbolic geometry

In the half-plane P = {(x, y) ∈ R 2 : y > 0} the riemannian metric and the Lagrange function are

ds 2 = dx 2 + dy 2 y 2 , L = 1 2 x ′2 + y ′2 y 2 .
The Euler-Lagrange equations are

d ds x ′ y 2 = 0, d ds y ′ y 2 + x ′2 + y ′2 y 3 = 0.
The x, y components of the acceleration are

A x = x ′′ -2 x ′ y ′ y , A y = y ′′ + x ′2 -y ′2 y
and equation ( 5) is expressed by the system

(21) x ′′ -2 x ′ y ′ y = -2αy ′ cosh(αs) , y ′′ + x ′2 -y ′2 y = 2αx ′ cosh(αs)
.

The Poleni curves (x(s), y(s)) are the solutions of ( 21) satisfying

(22) (x ′ (s)) 2 + (y ′ (s)) 2 (y(s)) 2 = 1.
Proposition 3 Up to isometry, the parametric equations of the Poleni curve, with α = 1, on the hyperbolic half-plane are

(23) x(s) = s, y(s) = cosh(s).
Proof-Up to isometry, one can restrict to the initial condition

x (0) = 0, x ′ (0) = 1, y (0) = 1, y ′ (0) = 0.
As in the proof of Proposition 1, we use the symbolic computing system MAPLE to write Taylor series of the solution and we conjecture formulas (23). By straightforward calculations we verify that ( 23) is actually the solution of ( 21), with α = 1, satisfying the previous initial conditions. Theorem 4 Let V be the peak of the Poleni curve. The geodesic which is orthogonal to the Poleni curve at point V is an axis of geodesical symmetry of the curve and it is asymptotic to the curve as s → ±∞. A geodesic is asymptotic to the Poleni curve as s → +∞ or s → -∞ if and only if it intersects transversally the Poleni curve at an unique point.

Proof-Up to isometry, the Poleni curve P is the graph y = cosh(x), with V = (0, 1). The geodesic which is orthogonal to P at V is the y axis. It is an axis of symmetry of P. A vertical line is the only geodesic that can be asymptotic to P. Let us prove that any vertical line x = a is asymptotic to P.

Let P = (s, cosh(s)) be a point of P. The geodesic passing by P and orthogonal to the vertical line x = a intersects this line at point Q (see Figure 4). Its equation is

(x -a) 2 + y 2 = R 2 , where R 2 = (s -a) 2 + cosh 2 (s).
Hence the distance from point P to the geodesic x = a is

d P (P, Q) = a s 1 + (f ′ (x)) 2 f (x) dx , f (x) = R 2 -(x -a) 2 .
Thus

d P (P, Q) = argtanh |s -a| (s -a) 2 + cosh 2 (s)
tends to 0 as s → ±∞. x 2 + y 2 < 1} via the transformation

z = x + iy ∈ P → z -i z + i = x 2 + y 2 -1 -2ix x 2 + (y + 1) 2 ∈ D.
In particular, the isometric image of (23) under the isometry z → -1 z of P is the curve.

x(s) = -s s 2 + cosh 2 (s) , y(s) = cosh(s) s 2 + cosh 2 (s)
.

The geodesics which are asymptotic to this curve are the y axis and all the circles with center on the x axis and tangent to the y-axis. All these geodesics intersect the Poleni curve transversally at a unique point.

Jurdjevic's parametrization of non euclidean elastic curves

Jurdjevic gave the following parametrization (see [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], Formula (10) page 117) of an elastic curve on the sphere of geodesical curvature k(s) :

(24) x = ± cos(φ) 1 - k 2 M , y = sin(φ) 1 - k 2 M , z = k √ M ,
where φ(s) is defined by

φ ′ = √ M H 1 M -k 2 , H 1 (s) = H - 1 2 k 2 (s).
The constants H and M in these formulas are given by (see [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], Proposition 4)

(u ′ ) 2 + u 3 -4u 2 (H -1) + 4u(H 2 -M ) = 0, u(s) = k 2 (s).
In the case k(s) = 2 cosh(s) , the function

u(s) = k 2 (s) is a solution of (25) (u ′ ) 2 + u 3 -4u 2 = 0. Hence H = 2, M = H 2 = 4 and φ ′ (s) = 2 2 - 2 cosh 2 (s) 4 - 4 cosh 2 (s) = 1. Thus φ(s) = s. On the other hand 1 - k 2 M = 1 - 1 cosh 2 (s) = | tanh(s)|.
Consequently formula (24) writes :

(26) x(s) = ± cos(s)| tanh(s)|, y(s) = sin(s)| tanh(s)|, z(s) = 1 cosh(s)
.

Actually, these formulas define two symmetrical Poleni curves on the sphere, compare with (20). Notice that there is no choice of the signe ± in (26) which gives a complete Poleni curve, see Figure 6. Jurdjevic gave also the parametrization of elastic curve on the hyperbolic space

H = {(x, y, z) ∈ R 3 : x 2 + y 2 -z 2 = -1, z > 0}.
He distinguishes three cases. The first case corresponds to M < 0, see [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], Formula (10), page 117. Notice that this formula contains a misprint : M should be replaced by |M |. The parametrization of the elastic curve is where φ(s) is defined by

(27) x = ± cos(φ) k 2 -M -1, y = sin(φ) k 2 -M -1, z = k √ -M ,
φ ′ = √ -M H 1 M -k 2 .
The second case corresponds to M > 0, see [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], Formula (12) page 120. The elastic curve is parametrized by

(28) x = - k √ M , y = cosh(φ) 1 + k 2 M , y = sinh(φ) 1 + k 2 M ,
where φ(s) is defined by

φ ′ = - √ M H 1 M + k 2 .
In the third case, corresponding to M = 0, the parametrization (which appears without number in page 122 of [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF]) is

(29) x = k 2 φ -k 2 + 1 2k , y = kφ, z = k 2 φ + k 2 + 1 2k ,
where φ(s) is defined by

φ ′ = H 1 k 2 .
In the hyperbolic case H 1 (s) = H -k 2 (s)/2, as in the spherical case, but constants H and M are given now by (see [START_REF] Jurdjevic | Non-Euclidean elastica[END_REF], Proposition 4)

(u ′ ) 2 + u 3 -4u 2 (H + 1) + 4u(H 2 -M ) = 0,
The Poleni curve satisfies (25). Thus, it corresponds to H = 0, M = H 2 = 0. In this case φ ′ = -1/2. Thus φ(s) = -s/2. Consequently formula (29) writes :

(30) x(s) = s 2 + cosh 2 (s) -4 4 cosh(s) , y(s) = -s cosh(s) , z(s) = s 2 + cosh 2 (s) + 4 4 cosh(s)

The coordinates (x 1 , y 1 ) in the Poincaré hyperbolic half space P are linked to the coordinates (x, y, z) by the relations

x 1 = y x -z , y 1 = 1 z -x .
Thus, the curve in P corresponding to the curve (30) is

x 1 = s 2 , y 1 (s) = cosh(s) 2 .
This curve is isometric to the Poleni curve (23) since the transformation x 1 = x/2, y 1 = y/2 is an isometry of P.

7 Poleni curves of geodesic curvature f (s) = 2α cosh(αs)

We analyze the spherical case first. We use cartesian coordinates. Here we have

H = 1 + α 2 , M = H 2 = (1 + α 2 ) 2 .
Using Jurdjevic's parametrization (24), we obtain, after some lengthy computations and simplifications, the following formulas (the details will appear in a forthcoming paper)

Proposition 5 Up to isometry, the parametric equations of the Poleni curve on the sphere S 2 are

x(s) = (1 -α 2 ) sin(s) + 2α cos(s) tanh(αs) 1 + α 2 , y(s) = (α 2 -1) cos(s) + 2α sin(s) tanh(αs) 1 + α 2 , z(s) = 2α (1 + α 2 ) cosh(αs) .

Figures 7 and8 show the curves corresponding to α = 2 and α = 1/2. From this parametrization we deduce the following result Theorem 6 The plane which is orthogonal to the Poleni curve at its peak V is a plane of symmetry for the curve. The diameter contained in the plane which is orthogonal to the symmetry plane and makes the angle

β(α) = 2 arctan(α)
with the vector OV , is asymptotic to the curve. The Poleni curve has an infinite number of self-intersection points of arc-lengths ±s n which are given by the solutions of equation We analyze now the hyperbolic case. We use cartesian coordinates on H 2 . Here we have

s n + arctan 2α tanh(αs n ) 1 -α 2 = nπ, n ∈ N. V V
H = α 2 -1, M = H 2 = (α 2 -1) 2 .
Using Jurdjevic's parametrization (28) in the case α 2 = 1, we obtain the following formulas (the details will appear in a forthcoming paper)

Proposition 7 Up to isometry, the parametric equations of the Poleni curve, with α 2 = 1, on the hyperbolic space are Notice that this family of elastic curves does not contain the Poleni curve (23) which corresponds to α = 1. From the parametrization we deduce the following result.

x(s) = 2α |α 2 -1| cosh(αs) , y(s) = (1 + α 2 ) sinh(s) -2α cosh(s) tanh(αs) |α 2 -1| , z(s) = (1 + α 2 ) cosh(s) -2α sinh(s) tanh(αs) |α 2 -1| . 0 2 1/2 -5 -1/5 2 1/2 -5 -1/5
Theorem 8 For all α = ±1, the geodesic which is orthogonal to the Poleni curve at its peak, is an axis of symmetry of the curve and there exist a unique geodesic, which is asymptotic to the curve as s → ±∞. When |α| < 1 the Poleni curve has no self-intersection. When |α| > 1 the Poleni curve has a unique selfintersection point. The arc-lengths ±s 0 of the self-intersecting point are given by the solution of equation tanh(s 0 ) = 2α 1 + α 2 tanh(αs 0 ). The Poleni curves corresponding to α and -α respectively, are isometric. Recall that in Theorem 4 we proved that in the case α = 1, all the geodesics which intersect the Poleni curve transversally at a unique point are asymptotic to the curve as s → ±∞. Figure 9 shows the curves in P and D corresponding to the Poleni curve of H 2 defined by the parametrization given in Proposition 7, with various values of α. In P, the peak V of the curve and its tangent vector t V are

V = 0, α -1 α + 1 , t V = α -1 α + 1 , 0 ,
thus V tends to the "boundary" as α → ±1. In P the asymptotic geodesic is the circle of radius 1 centered at the origin. In D the asymptotic geodesic is the vertical axis.

Figure 10 shows various Poleni curves in the hyperbolic half-space P, and the corresponding curves in the disk D, where the initial condition is fixed at x 1 (0) = 0, x ′ 1 (0) = 1, y 1 (0) = 0, y ′ 1 (0) = 0, where (x 1 , y 1 ) are the coordinates in P. The parametrization of these curves is

Figure 1 :

 1 Figure 1: Poleni's curve It is plotted in Figure 1. Actually, for every α ∈ R, the function f (s) = 2α cosh(αs)

  a) 2 +y 2 = R 2

Figure 4 :

 4 Figure 4: Vertical lines are asymptotic to the Poleni curve y = cosh(x) in P.

Figure 5 :

 5 Figure 5: Isometric Poleni curves in P and D.

Figure 5

 5 Figure 5 shows various images of the Poleni curve (23) under the hyperbolic rotations z = x + iy → az + b -bz + a , a 2 + b 2 = 1, together with their corresponding curves in the unit disk D = {(x, y) ∈ R 2 : x 2 + y 2 < 1} via the transformation

Figure 6 :

 6 Figure 6: The curve (26) with sign -(left) and sign + (right), seen from the north pole.

Figure 7 :Figure 8 :

 78 Figure 7: Poleni curves in the sphere with α = 2 (left) and α = 0.5 (right).

Figure 9 :

 9 Figure 9: Poleni curves in P and D, corresponding to the curve given by the parametrization of Proposition 7, with α = 2, α = 1/2, α = -5 and α = -1/5.

x 1 ( 1 -

 11 s) = (α + 1)(α 2 + 1) sinh[(α-1)s] α-(α 2 -1) cosh(s) sinh(αs) 2α (cosh[(α -1)s] + 1) + (α -1) 2 cosh(s) cosh(αs) , y 1 (s) = (α + 1) 2 cosh(αs) 2α (cosh[(α -1)s] + 1) + (α -1) 2 cosh(s) cosh(αs) ,This parametrization reduces to the Poleni curve (23) when α = 1. It shows how the Poleni curve (23) is the limiting case between Poleni curves with selfintersection and those without self-intersection.

Figure 10 :

 10 Figure 10: Poleni curves in P and D, corresponding to α = 0 (a geodesic), α = 0.65, α = 1 (the Poleni curve), α = 1.35 and α = 2.