Simon Arberet
email: simon.arberet@irisa.fr

Senior Member, IEEE Rémi Gribonval
email: remi.gribonval@irisa.fr

Frédéric Bimbot
email: frederic.bimbot@irisa.fr

A robust method to count and locate audio sources in a multichannel underdetermined mixture

Keywords: Blind source separation, multichannel audio, delay estimation, sparse component analysis, direction of arrival

We propose a method to count and estimate the mixing directions in an underdetermined multichannel mixture. The approach is based on the hypothesis that in the neighbourhood of some time-frequency points, only one source essentially contributes to the mixture: such time-frequency points can provide robust local estimates of the corresponding source direction. At the core of our contribution is a statistical model to exploit a local confidence measure which detects the timefrequency regions where such robust information is available. A clustering algorithm called DEMIX is proposed to merge the information from all time-frequency regions according to their confidence level. So as to estimate the delays of anechoic mixtures and overcome the intrinsic ambiguities of phase unwrapping as met with DUET, we propose a technique similar to GCC-PHAT which is able to estimate delays that can largely exceed one sample. We propose an extensive experimental study which shows that the resulting method is more robust in conditions where all DUET-like comparable methods fail, that is in particular : a) when time-delays largely exceed one sample; b) when the source directions are very close.

I. INTRODUCTION

In many situations like medical imaging, musical or meeting recording, the observed data is a measurement of several signals (called sources) which have been mixed together, and it is sometimes very useful to know how many sources are in the mixture, how these sources have been mixed together, and what are the original source signals. In the context of audio sources, the measured signals are often on two channels only, that is the well-known stereophonic case, and the number of sources is often higher than the number of channels.

In this article, we consider the problem of counting and estimating the directions of the audio sources from two or more channels when there may be more sources than available channels, with an emphasis on stereophonic audio mixtures. From the mixing directions, the source signals can be estimated using the classical time-frequency masking framework [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF], [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] or using more recent approaches [START_REF] Vincent | Complex nonconvex lp norm minimisation for underdetermined source separation[END_REF]- [START_REF] Arberet | Blind spectral-GMM estimation for underdetermined instantaneous audio source separation[END_REF]. The problem we consider is related to the direction of arrival (DOA) problem, which consists in estimating the physical directions of wave propagation given a controlled sensor arrangement, where the relative positions of the sensors is known. Here, we are interested in ad hoc sensor arrangements, thus: a) the estimated "directions" cannot necessarily be related to physical directions and should rather be considered as a parameterization of the mixture; b) it can happen that some sensors are either very close or very far apart from one another. If the sensors are very close to each other, the signals can be very similar from one sensor to another (the intensity differences are very small). On the other hand, if the sensors are far from each others, the delays between sensors can be high. In these both situations, it is difficult to estimate the source directions.

Our main contribution is a new technique to count and estimate the mixing directions which is robust in the sense that it can deal with the difficult settings where the underlying sensors maybe either very close to one another or very far apart. The proposed method, called DEMIX (Direction Estimation of Mixing matrIX) [START_REF] Arberet | A robust method to count and locate audio sources in a stereophonic linear instantaneous mixture[END_REF], [START_REF]A robust method to count and locate audio sources in a stereophonic linear anechoic mixture[END_REF] is based on a clustering algorithm which gives more weight to more reliable timefrequency regions, according to a local confidence measure similar to others proposed in the litterature [START_REF] Belouchrani | Blind source separation based on timefrequency signal representations[END_REF]- [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF]. An important contribution of our work is that this confidence measure is exploited via a statistical model, which is used to weight time-frequency regions in an statistically founded way, and to detect directions via a chi-square test. We also propose a new technique, which extends GGC-PHAT [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF] to the multisource case and which, unlike DUET [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], is able to estimate time-delays that can largely exceed one sample.

We demonstrate with extensive experimental studies 1 the ability of our approach to count and estimate the source directions, by varying the mixing conditions (number of sources, distance between directions, level of reverberation), and by comparing our method with the classical DUET. The experiments show that the proposed DEMIX approach is able to a) blindly estimate the number of sources; up to 8 sources in the instantaneous case; and up to 5 sources in the anechoic case; b) estimate time-delays up to 100 samples; d) count and estimate nearby source directions, up to a distance of 10 -3 degrees, with a constant relative precision better than 10 -3 .

A. Instantaneous and Anechoic mixture model

The mixture of N audio sources on M channels can be formulated by the anechoic mixture model :

x m (t) = N n=1 a mn s n (t -δ mn )+n m (t), 1 ≤ m ≤ M (1)
In this model, each source contributes to each microphone only through the direct acoustic path, that is to say with no reflection on walls or obstacles. The parameters a mn ∈ R represent the gain (or the attenuation) and δ mn the time-delay corresponding to the path between the n-th source and the m-th microphone. The problem we address in this paper is the estimation of the number of sources N and the mixture parameters a mn and δ mn , from the only observation of the (possibly noisy) mixture signals x m (t).

Without loss of generality, we assume that M m=1 a 2 mn = 1 and that δ 1n = 0 and a 1n ≥ 0 for 1 ≤ n ≤ N , which means that we fix the gain and sign indeterminacy of the problem. If, in addition, δ mn = 0, ∀m, n, the mixture model is so-called instantaneous. Taking the Short Time Fourier Transform (STFT) X m (t, f ) of each channel x m (t) of the mixture, the mixing process can be modeled in the time-frequency domain as X(t, f ) ≈ A(f )S(t, f ) + N(t, f ) for each time frame t and normalized frequency 0 ≤ f ≤ 1/2, where bold letters such as X(t, f ) or S(t, f ) denote column vectors [X 1 (t, f ), . . . , X M (t, f )] T or [S 1 (t, f ), . . . , S N (t, f )] T , and

A(f ) = [a 1 (f ), . . . , a N (f )] is a M × N complex valued mixing matrix which columns a n (f ) = [a 1n
, a 2n e -2iπδ2nf , . . . , a M n e -2iπδ M n f ] T are related to the source locations and are called steering vectors (SV) at frequency f .

In the stereophonic case (M = 2), each column of A(f ) can be written as a two-dimensional SV:

a n (f ) = cos θ n sin θ n • e -2iπδnf ∈ C 2 . (2)
The intensity parameter (IP) θ n ∈ (-π/2, π/2] characterizes intensity difference between channels, and a possible sign difference if θ n < 0; the parameter δ n ∈ R characterizes the time delay between channels. We will generally refer to the pair (θ n , δ n ) as the (mixing) direction of the n-th source.

For the case M > 2 channels, we can generalize this pair by splitting the direction of the n-th source into its intensity profile defined by abs (a n (f )) ∈ R M , where the function abs(•) from C M to R M calculates the magnitude of each element of a vector, with abs(a n (f )) 2 = 1, and delays

∆ n = [δ 1n , δ 2n , . . . , δ M n ] T ∈ R M with δ 1n = 0.

B. Related work 1) Hypothesis:

Several existing methods attempt to estimate the mixing directions of the sources from a timefrequency representation X(t, f ) of the mixture. DUET-like methods [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF], [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] rely on the assumption (called W-disjoint orthogonality) that the mixed sources have essentially disjoint time-frequency supports, that is to say in most time-frequency points, only one source has a non-negligible contribution. This is related to the sparsity assumption on the time-frequency representation of the sources. TIFROM [START_REF] Abrard | Blind separation of dependent sources using the "time-frequency ratio of mixtures" approach[END_REF] exploits the weaker assumption that for each source, there is at least one time-frequency region where this source is dominant.

When at most one source actively contributes to a timefrequency point (t, f ), there is an index

1 ≤ n(t, f ) ≤ N such that |S n(t,f ) (t, f )| ≫ |S n (t, f )|, n = n(t, f ), so the mixing model indicates that X(t, f ) ≈ A(f )S(t, f ) ≈ S n(t,f ) (t, f ) • a n(t,f ) (f ) and the ratio R 21 (t, f ) := X 2 (t, f )/X 1 (t, f ) satis- fies : R 21 (t, f ) ≈ tan θ n(t,f ) • e -2iπδ n(t,f ) f .
So, if the sources have disjoint time-frequency supports, then all data points X(t, f ) will be aligned along the SV a n(t,f ) (f ). Also, if the sources are sparse, the data points X(t, f ) show a clear tendency to cluster along a n(t,f ) (f ) [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]. This can be seen on the scatter plot of points X(t, f ), which is a simple tool we will use later in this paper (see for example Figures 1 and2). A common approach to estimate the mixing directions is thus based on a clustering algorihm applied to the points of the scatter plots.

2) Time-delay estimation: In DUET [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], the ratios R 21 (t, f ) are computed for each time-frequency point and used to compute a local estimate of the IP θ(t, f ) := tan -1 |R 21 (t, f )| and the delay δ(t, f ) := -1 2πf ∠R 21 (t, f ) where ∠z ∈ (π, π] is the argument of the complex number z. This approach is perfectly valid if the true delay is below one sample and the gains a mn are all positive, but it may fail otherwise because of phase unwrapping ambiguities. In a sense, the problem is that a single time-frequency point does not carry alone enough information to recover the corresponding delay δ n(t,f ) . Some approaches have been proposed recently to solve these issues [START_REF] Bofill | Underdetermined blind separation of delayed sound sources in the frequency domain[END_REF]- [START_REF]A time-frequency correlation-based blind source separation method for time-delayed mixtures[END_REF] and we propose a new one in Section VI. The proposed technique is able to estimate time delays which can largely exceed one sample, as illustrated in Section VII-F.

3) Clustering: When it comes to actually clustering local estimates of the directions to get a global estimate of the directions, many authors have chosen to use a weighted smoothed histogram [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], where the amount of smoothing is determined by the shape of a "potential function" [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]. One of the difficulties with this approach consists in adjusting how much smoothing must be performed on the weighted histogram to resolve close directions without introducing spurious peaks. Moreover, the choice of the weights is also of importance. The classical approach, which consists in giving more weight to time-frequency points with more energy, might prevent the clustering step from discovering the direction of a source of weak energy. Instead of using a fixed potential function and weights based on the local energy, we introduce in Section V a new clustering algorithm which relies on the local confidence measure introduced in Section II and a statistical model described in the Section III which is used in section IV to define a proximity measure between local estimates. The proposed method is also able to count the number of sources, which is not a widely addressed task [START_REF] Vincent | The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation[END_REF], [START_REF] Luo | Active source selection using gap statistics for underdetermined blind source separation[END_REF].

II. PRINCIPLE OF THE APPROACH

The proposed approach to estimate the mixing directions is based on the same assumption as TIFROM [START_REF] Abrard | Blind separation of dependent sources using the "time-frequency ratio of mixtures" approach[END_REF] and relies on the use of time-frequency regions. The first step is a feature extraction step. The second step of our method is the clustering algorithm which is defined in section V.

A. Feature extraction

For each time-frequency region Ω t,f "in the neighborhood" of the time-frequency point (t, f ), the principle is to estimate two values:

1) an estimation û(Ω t,f ) of the SV of the most dominant source. Thus û(Ω t,f ) is called an estimated steering vector (ESV); 2) a local confidence measure, denoted T (Ω t,f ), which gets larger when the scatter plot of vectors X(τ , ω) in the region Ω t,f points more strongly in the direction of the ESV û(Ω t,f ), that is when essentially one source is active in this region. As the confidence measure T (Ω t,f ) can discriminate the regions where essentially one source is active from the other ones, it also discriminate regions where the ESV û(Ω t,f ) points in a direction of a SV from the regions where û(Ω t,f ) is unlikely to point in one of the SV direction.

B. Time-frequency regions

We consider two kinds of time-frequency regions around each time-frequency point (t, f ) : a temporal neighborhood Ω T t,f and a frequency neighborhood Ω F t,f defined by :

Ω T t,f = {(t + kL/2, f ) | |k| ≤ K} (3)
Ω F t,f = {(t, f + k/L) | |k| ≤ K}. ( 4
)
where L is the STFT window size and k ∈ Z.

C. Real-valued and complex-valued local scatter plots

Each region Ω provides a complex-valued local scatter plot X(Ω). It is a M × (2K + 1) matrix with columns X(τ , ω) ∈ C M , (τ , ω) ∈ Ω which will be used to analyse anechoic mixtures. For linear instantaneous mixtures, since the SV a n (f ) of the sources are real-valued, a real-valued local scatter plot will be used instead. It corresponds to a M × (4K + 2) matrix denoted X R (Ω) with columns ℜX(τ , ω) ∈ R M and ℑX(τ , ω) ∈ R M , (τ , ω) ∈ Ω.

D. Principal Component Analysis and confidence measure

To compute the ESV û(Ω t,f ) and their corresponding local confidence measure T (Ω t,f ), one can simply rely on Principal Component Analysis (PCA) applied to vectors X(τ , ω) in the region Ω t,f . Performing a PCA on the local scatter plot X(Ω) (resp. X R (Ω)), we obtain a principal component (PC) as a unit vector û(Ω) ∈ C M (resp. û(Ω) ∈ R M ) defined up to a multiplicative factor e iψ (resp. up to a sign ±1), as well as the real-valued positive eigenvalues in decreasing order λ1 (Ω) ≥ . . . ≥ λM (Ω) ≥ 0 of the M ×M complex Hermitian positive semi-definite matrix X(Ω)X H (Ω) (resp. the real symmetric positive semi-definite matrix X R (Ω)(X R (Ω)) T ). We define the (empirical) confidence measure as :

T (Ω) := λ 1 (Ω) 1 M -1 M m=2 λ m (Ω) .
(5)

We will discuss in Section III why this measure can also be viewed as a signal-to-noise ratio between the dominant source and the contribution of the other sources (plus noise). Thus it is useful to express it in the deciBel (dB) scale : 10 log 10 ( T (Ω)).

Figure 1 shows the local scatter plot of X R (Ω) in two timefrequency regions of an audio mixture : as expected from the theoretical results of Section III, the confidence measure is high when essentially one source is active, and low when many sources are simultaneously active. Figure 2(a) displays the real-valued global scatter plot for all time-frequency points weighted by their energy, which is used in standard approaches to determine the mixing directions. One can see that there are many "outliers", i.e. points with high energy that are not close to the mixing directions. Conversely, Figure 2(b) displays the collection of vectors ±10 log 10 T (Ω t,f ) • û(Ω t,f ) obtained by PCA for all timefrequency regions of the signal. One can observe that points of Figure 2(b) are better concentrated along the mixing directions and thus should be better candidates to estimate the mixing direction using a clustering algorithm. This will be confirmed experimentally.

III. LOCAL STATISTICAL MODEL

In this section, we analyze the relation between the (empirical) local confidence measure and the reliability of the ESV, based on a simple statistical model of the mixing process in a time-frequency region.

In the instantaneous mixing model (δ n = 0, ∀n), the mixing matrix A(f ) is a real-valued matrix A which does not depend on the frequency. By taking the real or imaginary part of the complex-valued mixture model

X(t, f ) = AS(t, f ) + N(t, f ), an equivalent real-valued model is ex- pressed as : [ℜX(t, f ), ℑX(t, f )] = A[ℜS(t, f ), ℑS(t, f )] + [ℜN(t, f ), ℑN(t, f )].
In the proposed model, we assume that a source s n , with a SV a n , is the most active source in the region Ω. The values of the real and imaginary parts of the STFT of this source in the region Ω are modeled as independent centered random variables with (large) variance σ 2 s . The contribution of the other sources, including possibly noise, are modeled by an independent, isotropic, M -dimensional centered distribution with covariance matrix σ 2 n I M . Therefore, the entries

ℜX(τ, ω), ℑX(τ, ω), (τ, ω) ∈ Ω of the scatter plot X R (Ω) = a n • S R n (Ω) + N R (Ω) have covariance : Σ X = σ 2 s a n a T n + σ 2 n I M . (6)
The largest eigenvalue of Σ X is λ 1 = σ 2 s + σ 2 n associated to the PC u 1 and the remaining eigenvalues are λ 2 = . . . = λ M = σ 2 n . It follows that the direction of the SV a n coincides with the PC u 1 . The PC u 1 is defined up to a multiplicative factor. Nevertheless, as a n is normalised, we can impose that u 1 be normalised too and then we have u 1 = ±a n . Thus the ESV is defined up to a sign. The "true" confidence measure defined as :

T := λ 1 1 M -1 M m=2 λ m = σ 2 s /σ 2 n + 1 (7)
can be viewed as a signal-to-noise ratio between the dominant source and the contribution of the other sources (plus noise).

A. Asymptotic distributions

If the observation of the scatter plot X R (Ω) were sufficient to get a perfect estimate of the covariance matrix Σ X , the analysis would be over. However, in practice, the PC û(Ω) and the local confidence measure T (Ω) are computed on samples of only |Ω| points. Hence, û(Ω) and T (Ω) only provide an estimation of a n and T . It is nonetheless possible to rely on random matrix theory so as to quantify the precision of these estimates as a function of the sample size |Ω|.

Let

Σ X = UΛU T , with Λ = diag([λ 1 , • • • , λ M ]
), be the spectral decomposition of Σ X , and ΣX = ÛΛ ÛT , with [START_REF] Anderson | An Introduction to Multivariate Statistical Analysis[END_REF]Theorem 13.5.1] that the following properties are satisfied :

Λ = diag([ λ1 , • • • , λM ]), the one of the empirical covariance matrix ΣX := |Ω| -1 X R (Ω)(X R (Ω)) T . When X R (Ω) is Gaussian, it is known from [19, Corollary 7.2.3] that the empirical covariance matrix ΣX follows a Wishart distribution |Ω| -1 W M (Σ X , |Ω| -1) of dimension M . If M = 2, since the eigenvalues in Λ are pairwise distinct, it follows from
|Ω| -1 • ( λ1 -λ 1 ) L -→ N (0, 2λ 2 1 ) (H 1 ) |Ω| -1 • M m=2 ( λm -λ m ) L -→ N 0, 2 M m=2 λ 2 m (H 2 ) |Ω| -1 • (û 1 -u 1 ) L -→ N (0, V) (H 3 )
where L -→ denotes convergence in law when the sample size |Ω| gets large. The M × M covariance matrix V is given by :

V = λ 1 m≥2 λ m (λ m -λ 1 ) 2 u m u T m = σ 2 s σ 2 n + 1 • σ 2 n σ 2 s 2 (I M -u 1 u T 1 ) = T (T -1) 2 • (I M -u 1 u T 1 ).
The same properties can be proved [START_REF]Asymptotic theory for principal component analysis[END_REF] when M > 2 and

λ 1 > λ 2 ≥ • • • ≥ λ M .
From now, we will assume that H 1 , H 2 , H 3 hold, which is likely to hold for a much wider class of distributions of X R (Ω) than the only Gaussian distribution. Discussing the full extent of the validity of H 1 , H 2 , H 3 is beyond the scope of this paper.

Although the convergences in H 1 , H 2 , H 3 are asymptotic (for a large sample size |Ω|), we will use a pretty small sample size (|Ω| = 10). Experimental results in section VII show that the method is nonetheless working well.

B. Robust confidence measure

So as to evaluate the quality of the ESV using the asymptotic relation of hypothesis H 3 , we define a measure called robust empirical confidence, which is a bound that has a high probability to lie below the true (unknown) confidence measure.

Definition (Robust empirical confidence) Using the empirical confidence measure T (Ω), we define the robust empirical confidence T (Ω) of level 1α by :

T (Ω) := T (Ω)e -q(α) q 2M (|Ω|-1)(M -1) (8)
where q(α) is the quantile of order 1α of the Gaussian law with zero mean and variance one.

The quantile q(α) = F -1 (1 -α) is defined by the inverse function F -1 of the cumulative distribution function (CDF) F (q) = q -∞ 1 √
2π e -y 2 /2 dy of law N (0, 1). Thanks to Theorem 3.1 below, for a fixed and sufficiently large sample size |Ω|, the value of the "true" confidence measure T (Ω) is guaranteed to be higher than the one of the robust empirical confidence T (Ω) with a probability 1α. In our experiments, we chose q(α) = 6.3.

Theorem 3.1: (Confidence interval for the confidence measure) Assume that H 1 , H 2 hold, and let T be defined by Equation ( 7) and T :=

λ1 1 M -1 P M m=2 λm
. Then :

P T ≥ T e -q(α) q 2M (|Ω|-1)(M -1) -----→ |Ω|→∞ 1 -α (9)
where q(α) is the quantile of level 1α of law N (0, 1).

Proof: From H 1 and H 2 , denoting μ := ( λ1 , 1 M -1 M m=2 λm ) and µ := (σ 2 s + σ 2 n , σ 2 n ), we have : |Ω| -1 • (μ -µ) L -→ N 0, 2 • diag µ 2 1 , µ 2 2 M -1 . Writing 1 2 ln T = f (μ) with f (x 1 , x 2 ) = 1 2 ln x 1 -1 2 ln x 2 , from [19, Theorem 4.2.3] we have with d = ( ∂f ∂xi µ ) i=1,2 : |Ω| -1 1 2 ln T - 1 2 ln T L -→ N 0, 2 • d T diag µ 2 1 , µ 2 2 / (M -1) d . One can easily check that d T diag µ 2 1 , µ 2 2 / (M -1) d = M 4(M -1)
, and we obtain :

2(M -1) M • |Ω| -1 1 2 ln T - 1 2 ln T L -→ N (0, 1).
To conclude, for sufficiently large sample size |Ω|, we have :

P T ≤ T e -q q 2M (|Ω|-1)(M -1) = P 1 2 ln T - 1 2 ln T ≥ q M 2(|Ω| -1)(M -1) = F (q).
where F (q) is the CDF of law N (0, 1). Then, there is a quan-

tile q(α) = F -1 (1 -α) such that T ≥ T e -q(α) q 2M (|Ω|-1)(M -1)
with probability exceeding 1α.

C. Precision of the direction estimation

Let us now come back to the relation between the above defined T (Ω) and the covariance matrix V of the asymptotic distribution of û(Ω) around the SV a n . Proposition 3.2: If the mixture covariance matrix Σ X is defined as in Equation ( 6) and if hypothesis (H 3 ) holds, then the distribution of the PC û1 converges in law, when the sample size |Ω| is large, to :

û1 L -→ N (a n , σ 2 (T (Ω)) • R) (10)
with :

σ 2 (T ) := T (|Ω| -1) • (T -1) 2 (11)
R := I M -a n a T n . ( 12
)
and also :

E û1 -a n 2 = (M -1) • σ 2 (T ). ( 13
)
Proof: According to H 3 , the asymptotic distribution of û1 is :

û1 L -→ N (a n , (|Ω| -1) -1 V)
Also, from Equation ( 10), the square of the asymptotic distance between û1 and a n is :

û1 -a n 2 = σ 2 (T (Ω)) • Ξ
where the random variable Ξ ∼ χ 2 (M -1) is distributed according to a χ 2 distribution with M -1 degrees of freedom.

D. χ 2 test to decide if two regions correspond to the same direction

We assume two independent random variables û(Ω 1 ), û(Ω 2 ) distributed according to the following laws :

û(Ω 1 ) ∼ N (a n , σ 2 (T (Ω 1 )) • R) (14) û(Ω 2 ) ∼ N (a n ′ , σ 2 (T (Ω 2 )) • R) (15)
We want to test the hypothesis

H 0 = {a n = a n ′ } against the hypothesis H 1 = {a n = a n ′ }. If a n = a n ′
, the square distance between û(Ω 1 ) and û(Ω 2 ) is given by :

û(Ω 1 ) -û(Ω 2 ) 2 = σ 2 (T (Ω 1 )) + σ 2 (T (Ω 2 )) Ξ (16)
where Ξ ∼ χ 2 (M -1). Thus :

P û(Ω 1 ) -û(Ω 2 ) 2 σ 2 (T (Ω 1 )) + σ 2 (T (Ω 2 )) ≥ q χ 2 M -1 (α) = α where q χ 2 M -1 (α)
is the quantile of level α of the χ 2 law with M -1 degrees of freedom. Thus, the χ 2 test of level α consists in rejecting hypothesis H 0 if :

û(Ω 1 ) -û(Ω 2 ) 2 ≥ q χ 2 M -1 (α) • σ 2 (T (Ω 1 )) + σ 2 (T (Ω 2 )) (17)
IV. PROXIMITY MEASURE Before using a clustering algorithm for estimating the mixture directions, we need to define a proximity measure so as to quantify the degree of similarity between pairs û(Ω), T (Ω) . As these pairs are indexed by a time-frequency region Ω, we use the term region to indicate these pairs. We first define the distance measure between two SV u 1 and u 2 , and then we use the asymptotic results from section III to define the distance measure between two regions û(Ω 1 ), T (Ω 1 ) and û(Ω 2 ), T (Ω 2 ) .

A. Distance measure between two steering vector (SV)

Because of the sign indeterminacy (or more generally, in the complex case, the phase indeterminacy) of the ESV, we define the square of the direction distance measure between two unit SV u 1 and u 2 by :

d 2 (u 1 , u 2 ) := min |z|=1,z∈C u 1 -zu 2 2 = 2 (1 -| u 1 , u 2 |) . ( 18
)
As an exemple, let u 1 and u 2 be two real 2-dimensional SV, with IP θ 1 and θ 2 and such that θ 1θ 2 ≈ 0. Then

d 2 (u 1 , u 2 ) = 2 (1 -| cos (θ 1 -θ 2 ) |) ≈ (θ 1 -θ 2 ) 2 .

B. Distance measure between two regions

According to the model developped in section III, two regions Ω 1 and Ω 2 belong to the same cluster if the asymptotic Gaussian distribution of their respective ESV û(Ω 1 ) and û(Ω 2 ) have the same mean, that is u(Ω 1 ) = u(Ω 2 ). The χ 2 test of Equation ( 17) can be used to test this hypothesis.

Define D T as :

D T ((u 1 , T 1 ) , (u 2 , T 2 )) := d(u 1 , u 2 ) σ 2 (T 1 ) + σ 2 (T 2 ) , (19)
and let ζ = q χ 2 M -1 (α). The hypothesis u(Ω 1 ) = u(Ω 2 ) is accepted (with confidence level α) if :

D T û(Ω 1 ), T (Ω 1 ) , û(Ω 2 ), T (Ω 2 ) ≤ ζ (20)
In other words, if [START_REF]Asymptotic theory for principal component analysis[END_REF] holds then û(Ω 1 ) and û(Ω 2 ) are considered as sufficiently close to one another to be merged. This distance measure D T is similar to the known Fisher criterion function [START_REF] Duda | Pattern classification ans Scene Analysis[END_REF]. In practice the occurrence of the unknown "true" confidence values T (Ω) are not accessible. Thus, we replace them by their empirical estimate T (Ω), or a more pessimistic estimate T (Ω) < T (Ω) defined in Equation ( 8) and depending on quantile q(α).

We will see in the next section how the distance measure D T , obtained by our statistical model, is used to cluster the directions, and also (on Figure 3) how the shape of the clusters, determined by measure D T , is matching with the data.

V. DIRECTION ESTIMATION WITH A CLUSTERING ALGORITHM

In this section we describe the two proposed DEMIX (Direction Estimation of Mixing matrIX) algorithms used to estimate the mixture directions by combining and clustering the local estimates û(Ω), T (Ω) .

According to the classification of Theodoridis and Koutroumbas [START_REF] Theodoridis | Pattern recognition[END_REF], there are three families of clustering algorithms : 1) partitional clustering based on a cost function, 2) hierarchical clustering, 3) sequential clustering. We begin by a brief review of the principles and pros and cons of these families, before describing our algorithms which are related to the BSAS algorithm [START_REF] Theodoridis | Pattern recognition[END_REF], of the third family.

A. Classification of clustering algorithms 1) Partitional clustering based on a cost function:

We could do the clustering with a standard algorithm like the Kmeans (also known as LBG [START_REF] Linde | An algorithm for vector quantizer design[END_REF]), which minimises iteratively the within-class variance. But this approach has two main drawbacks :

• There is no guarantee to converge to a global optimum, and the obtained centroids can depend on the initial partitions. Thus the general strategy for the problem, as it is actually done by DUET, is to run the algorithm many times with random initial partitions. Another solution adopted by the ELBG [START_REF] Patanè | The enhanced LBG algorithm[END_REF] is to call for a "roulette" mechanism typical of genetic algorithms so as to avoid convergence to a local optimum. We will evaluate the performance of this last approach in section VII.

• The number N of clusters must be provided as an input parameter. A large number of attempts have been made to estimate the appropriate N , but their performance are often data dependent [START_REF] Xu | Survey of clustering algorithms[END_REF] and/or computationally intensive [START_REF] Luo | Active source selection using gap statistics for underdetermined blind source separation[END_REF].

2) Hierarchical clustering: Hierarchical clustering (HC) algorithms organise data into a hierarchical structure according to the proximity matrix [START_REF] Xu | Survey of clustering algorithms[END_REF]. There are some heuristics to estimate the number N of clusters, but the main problem of HC algorithms is that the computational complexity is generally at least O(Q 2 ), which is a limit for our application where the number Q of data points is generally larger than one million.

3) Sequential Clustering: These algorithms produce a single clustering [START_REF] Theodoridis | Pattern recognition[END_REF] and are thus faster. Algorithm like the Basic Sequential Algorithmic Scheme (BSAS) [START_REF] Theodoridis | Pattern recognition[END_REF], does not know a priori the number of clusters and thus can be used to count the number N of sources.

B. DEMIX-Instantaneous

First we present DEMIX-Instantaneous, which is designed for instantaneous mixtures, then we present in section V-C DEMIX-Anechoic for anechoic mixtures.

The first step of the algorithm consists in iteratively creating clusters by selecting regions Ω k with highest empirical confidence T (Ω k ) and aggregating to them other regions which ESV are sufficiently close to û(Ω k ), in the sense of Eq. ( 20)). The number K of clusters thus created depends on the structure of the scatter plot of all the regions û(Ω), T (Ω) (see Figure 2). The second step of the algorithm is to estimate the centroid ûc k of each cluster by first selecting a subset of confident regions and then weighting these regions according to their confidence value. As some clusters may be created by some outliers, we finally use a statistical test to eliminate unreliable clusters and keep N ≤ K clusters which centroids provide the estimated directions of the mixing matrix. Below, we detail each step of the algorithm.

1) Cluster creation : The first step of the algorithm iteratively creates K clusters C k ⊂ P where P is the set of all regions Ω of the scatter plot. As each cluster corresponds to an estimated mixture direction, the aim of this step is also to get a first estimation of the number of sources.

1.1) initialize : K = 0, P K = P 0 = P ; 1.2) find the region Ω K ∈ P K with highest confidence:

Ω K := arg max Ω∈P K T (Ω); 1.
3) create a cluster C K with all regions Ω ∈ P such that û(Ω) is sufficiently close to û(Ω k ); 1.4) update P K+1 = P K \ C K by removing regions of cluster C K which are still in P K ; 1.5) stop if P K = ∅, otherwise increment K ← K + 1 and go back to 1.2.

Note that, in step 1.3, the newly created cluster may contain regions already contained in previous clusters.

2) Direction estimation : after creating K clusters {C k } K k=1 , we estimate their centroids u(C k ). As we can see on Figure 3, the distribution of the points/regions around a mixture direction k is "symmetrical" only when the confidence measure is large enough. Thus, so as to have non biased estimation for direction k, the estimation is based on a subset C ′ k ⊂ C k of confident regions which have a confidence measure larger than an adaptative threshold. To define this threshold, we first set the condition that a region cannot be used to estimate two different directions; secondly, so as to have a maximum number of regions to estimate a direction, we force the threshold to be as small as possible.

The estimation is finally done in the following steps : 2.1) determine the confidence threshold :

η k := max Ω∈C k ∩[∪ j =k Cj ] T (Ω) (21)
2.2) keep only regions with sufficiently high empirical confidence values (the other ones are no longer used in the rest of the algorithm) :

C ′ k := Ω ∈ C k T (Ω) ≥ η k .
T " (u, T ) , " û(Ω k ), b T (Ω k ) "" = ζ
In light of the statistical model developed in Section III, Eq. ( 10)-( 12), each ESV û(Ω) of the thresholded cluster C ′ k is distributed as N (u k , σ 2 (T ) • R) . The minimum variance unbiased estimator of the "true direction" u k is given by :

v k := Ω∈C ′ k σ -2 (T (Ω)) • û(Ω) Ω∈C ′ k σ -2 (T (Ω)) (22)
In practice, since the ESV û(Ω) are only defined up to a sign, we multiply each direction with the sign ε(Ω) such that the correlation ε(Ω) • û(Ω), û(Ω K ) with the direction û(Ω K ) is positive. Moreover, the true confidence levels must be replaced with empirical estimates :

ṽk := Ω∈C ′ k σ -2 ( T (Ω)) • sign ( û(Ω), û(Ω K ) ) • û(Ω) Ω∈C ′ k σ -2 ( T (Ω)
) .

(23) Finally, the direction estimation formula for cluster C k after normalisation is :

u(C k ) := ṽk / ṽk (24)
3) Cluster Elimination : the last step of the algorithm consists in eliminating unreliable clusters. When two clusters are very close, it may be the consequence of two nearby directions, or one of these clusters may have been created because of some outliers. However, if we can define a confidence measure for a cluster, we can assume that two clusters which are close to each other are associated with two distinct directions, if and only if these two clusters have a large confidence measure.

Under the model developed in section III, the minimum variance unbiased estimator defined by ( 22) is distributed as :

v k ∼ N u k , Ω∈C ′ k σ -2 (T ) -1 R . ( 25
)
Thanks to Proposition 3.2, the error of the estimation v k of direction u k is characterized by :

σ 2 (C k ) := E v k -u k 2 = (M -1)• Ω∈C ′ k σ -2 (T (Ω)) -1 (26
) In practice, since the "true" confidence measure T (Ω) is not known, we replace it by the value of the robust empirical confidence measure T (Ω). The variance measure σ 2 (C k ) can be converted into a confidence measure using the inverse function of Equation [START_REF] Févotte | Two contributions to blind source separation using time-frequency distributions[END_REF].

The cluster elimination step consists in reiterating the Cluster creation step of section V-B1, with cluster centroids u(C k ) and their associated confidence measure T (C k ) as input instead of regions û(Ω), T (Ω) .

To summarize, cluster C j will merge with cluster

C o = C j if : D T û(C j ), T (C j ) , û(C o ), T (C o ) ≤ ζ c ( 27
)
The value of ζ c has to be larger than the threshlold ζ of Eq. ( 20) to efficiently eliminate unreliable clusters. Notice that if we knew (a priori) the number N of directions, we could implement this cluster elimination step by keeping the N directions which have the largest confidence value T (C j ).

C. DEMIX-Anechoic

We now detail the DEMIX-Anechoic algorithm. The main difference with the DEMIX-Instantaneous algorithm lies in the cluster creation step, as each mixture direction, in addition to be characterized by an intensity profile abs (a n (f )) ∈ R M , is also characterised by interchannel delays ∆ n . Thus the centroid of a cluster is now defined by a frequency dependent steering vector u C k (f ) parameterized by both :

• a frequency independent intensity profile abs(u C k (f )).

• frequency dependent phases on each channel determined by the delays ∆ k . The main changes in the algorithm are: a) the incorporation of the time-delay estimation step, b) the test to determine when the complex-valued ESV û(Ω) of a region is sufficiently close to a cluster.

1) Cluster creation and delay estimation : This step follows the same iterative procedure as for DEMIX-Instantaneous as described in section V-B1, except for step 1.3 divided now in 2 steps : 1.3.a) create a temporary cluster C K with all regions Ω ∈ P K for which abs(û(Ω)) is sufficiently close to abs(û(Ω K )), that is to say regions Ω such that :

d(abs(û(Ω)), abs(û(Ω K ))) ≤ ζ 2 • σ( T (Ω K )),
where T (Ω K ) is defined in Equation ( 8) and ζ 2 is a threshold. 1.3.b) estimate the interchannel delays ∆ K for C K ;

if ∆ K is considered as well identified (cf Section VI) : define the centroid u C K (f ) using the intensity profile abs(û(Ω K )) and the delays ∆ K ; create the cluster C K with all regions Ω ∈ P sufficiently close to u C K (f ); otherwise : reject the cluster C K := CK ; In Step 1.3.b, we need to compute the distance between an ESV û(Ω) and the centroid SV u C K (f ), which is frequency dependent. Therefore, we consider as sufficiently close all regions Ω ∈ P such that :

D T û(Ω), T (Ω) , u C K (f (Ω)), T (Ω K ) ≤ ζ, ( 28
)
where f (Ω) is the central frequency of the time-frequency region Ω.

2) Direction estimation and cluster elimination:

The direction estimation step and the cluster elimination step are similar to the last two steps of DEMIX-Instantaneous. The main difference is that, instead of using the distance d(., .) in the definition of D T (see Equation ( 19)), we use the distance :

d c (u Ci (•), u Cj (•)) = d(u Ci (f ), u Cj (f ))df. ( 29
)
so as to take into account the fact that the steering vectors depend on the frequency f .

VI. TIME-DELAY ESTIMATION

In this section, we present a method that estimates the time-delay of directions. We begin with a presentation of the approach for stereophonic mixtures, where only one delay needs to be estimated, before extending it to more channels.

A. Principle of the method 1) Case where only one source is active:

To explain the basic idea of the method, let us assume for a moment that only one source n is active in time frame t. Then, for each frequency, the DUET ratio satisfies R 21 (t, f ) = tan( θ(t, f ))e i b φ(t,f ) ≈ a2n a1n e -2iπf δn , and the Inverse Fourier Transform (IFT) of R 21 (t, f )/|R 21 (t, f )| ≈ e -2iπf δn yields a Dirac at time δ n :

r 21 (τ ) := R 21 (t, f ) |R 21 (t, f )| e i2πf τ df ≈ δ(τ -δ n ) (30)
The GCC-PHAT method [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF] consists in detecting the peak in function r 21 (τ ) with the following estimator :

δ n := arg max τ r 21 (τ ) (31)
2) Case where more than one source are active: In practice, one rarely observes an entire time frame t where only one source is active, but as indicated in section V-C1 one can determine a set Cn of time-frequency regions, which have similar intensity profiles, and where it is likely that only one source is active. In each of the regions Ω ∈ Cn , the phase of the ESV û(Ω) is e iφ(Ω) ≈ e -2iπδnf , where f = f (Ω) is the "central frequency" of the time-frequency region. The accuracy of this approximation is related to the value of σ 2 ( T (Ω)) as given in Equation [START_REF] Févotte | Two contributions to blind source separation using time-frequency distributions[END_REF]. One can expect to obtain a more accurate estimate of the phase for a given frequency f , by weighting all estimates corresponding to timefrequency regions Ω with central frequency f according to their precision. For that purpose, we propose the following estimator :

R Cn (f ) := Ω w f (Ω)e i b φ(Ω) Ω w f (Ω) ≈ e -2iπδnf (32)
with

w f (Ω) := 1/σ 2 ( T (Ω)) if Ω ∈ C n and f = f (Ω) 0 otherwise .
(33) The IFT of R Cn (f ) is expected to be approximately a Dirac delta at time δ n .

r Cn (τ ) := R Cn (f )e i2πf τ df ≈ δ(τ -δ n ) (34)
The highest peak of this function provides the final time-delay estimate :

δ n := arg max τ r Cn (τ ) (35)
In practice, we consider that a well identified peak is found if the amplitude r Cn ( δ n ) of the main peak of function r Cn (τ ) exceeds that of all other possible peaks by at least 3dB. Notice that our time-delay estimator extends the GCC-PHAT estimator [START_REF] Knapp | The generalized correlation method for estimation of time delay[END_REF] to the case of multiple sources.

B. Delay estimation for more than two channels

For more than two channels, there are M -1 > 1 interchannel delays to estimate, and their definition depends on the channel we choose as a reference. Notice that, even for a time-frequency point where only one source is active, if the intensity on a channel m is close to zero (that is ℜeX m (t, f ) ≈ 0 and ℑmX m (t, f ) ≈ 0), then the phase estimation ∠X m (t, f ) = tan -1 ℑmXm(t,f ) ℜeXm(t,f ) on that channel is unstable, as well as the phase difference between channel m and a channel k = m. To avoid these intrinsic phase unstabilities, we propose to choose as a reference, for cluster C K , the channel with the largest intensity in the intensity profile abs(û(Ω K )) = (u m ) M m=1 of region Ω K : we let m K := arg max m |u m | and estimate the interchannel delays δ m,mK between each channel m = m K and the reference channel m K , using the stereophonic time-delay estimation method described in section VI-A. We consider that ∆ K is well identified if all delays of ∆ K are well identified.

C. Discrete time implementation

In practice, time-frequency representations are only computed with a discrete grid of frequencies. As a consequence, the estimators defined in Equations ( 34) and ( 35) only provide time-delays on a discrete time grid. If the IFT of Equation ( 34) is computed with the same frequencies as those used by the STFT X(t, f ), then the temporal resolution of the delay estimator of Equation ( 35) is one sample. It is nevertheless possible to increase this resolution by zero padding or "spectral zooming" [START_REF] Hoyer | The zoom fft using complex modulation[END_REF] the function of Equation (32).

VII. EXPERIMENTAL STUDY

In order to evaluate the DEMIX-Instantaneous and Anechoic methods, we propose in this section experiments aiming at :

• comparing the proposed DEMIX-Instantaneous algorithm with classical clustering approaches using the (K-Means like) ELBG algorithm and variants. • testing the limits of the DEMIX-Instantaneous algorithm on anechoic mixtures. • comparing the ability of DEMIX-Anechoic and DUET 2in estimating the of anechoic and convolutive mixtures. All experiments were performed in the stereophonic case. Thus the SV are defined according to Equation (2), and in the instantaneous case the ESV are only parametrized by an IP θ(t, f ).

A. Clustering algorithm variants

The purpose of these four variants of ELBG is to diagnose the success and failure of the DEMIX algorithm. In other words, we want to understand the impact on the results of : a) the "local smoothing" effect of PCA, which replaces a pointwise estimate of an IP at a given time-frequency point with a smoothed estimate averaged on a time-frequency region; b) the use of a confidence measure rather than the energy value to give more weight to the IP of specific timefrequency region.

We compare the DEMIX-Instantaneous algorithm with ELBG [START_REF] Patanè | The enhanced LBG algorithm[END_REF], which is an improvement of the classical LBG (i.e. K-means) algorithm [START_REF] Linde | An algorithm for vector quantizer design[END_REF], on instantaneous mixtures. We considered four variants of the ELBG algorithm :

• ELBG on the IP θ(t, f ) obtained from the time-frequency bins X(t, f ). That is to say the classical ELBG;

• WELBG (a weighted variant of ELBG) on the IP θ(t, f ) obtained from the time-frequency bins X(t, f ) using the amplitude ρ(t, f ) = X(t, f ) as a weight; • ELBG on the IP θ(t, f ) obtained from the ESV û(t, f ) after the PCA; • WELBG on the IP θ(t, f ), using the optimal weight 1 σ 2 T (t, f ) where σ 2 (T ) is defined in Equation [START_REF] Févotte | Two contributions to blind source separation using time-frequency distributions[END_REF]. The different algorithms tested in this study are represented on Figure 4. We also propose one variant of DEMIX-Instantaneous where, instead of estimating the mixing direction with a weighted mean using Equations ( 23) and ( 24), we take the "best region", that is the centroid û(Ω K ) as a direction estimation (we call this variant DEMIX-Instantaneous-BR).

x1(t) x2(t) X1(t, f ) X2(t, f ) θ(t, f ) T (t, f ) scale Ω T v STFT local PCA Polar θ(t, f ) ρ(t, f ) ELBG WELBG ELBG WELBG DEMIX

B. Performance measures

As the proposed DEMIX methods are able to both estimate the number of sources and the mixing directions, we propose two measures to evaluate the performance of each of them.

1) Counting accuracy: A first measure of performance is the rate of success in the estimation of the number of sources. This measure is applied only on DEMIX, because DUET, ELBG and its variants do not estimate the number of sources.

2) Accurate directions estimation: We propose a performance measure called the mean direction error (MDE) which is the mean distance between true directions and estimated ones, computed with an optimized permutation to best match directions.

For a linear instantaneous mixture, given the true directions A = [a 1 . . . a N ] and estimated ones A = [ a 1 a 2 . . . a N ] the mean direction error (MDE) is defined as :

MDE (A, Â) := min π∈S N 1 N N n=1 d(a n , âπ(n) ) (36)
where S N is the permutation group of size N .

To also measure the error in terms of relative precision, we also define the relative mean direction error (RMDE) as the MDE divided by the min-distance between true directions :

RMDE (A, Â) := MDE (A, Â) min n =n ′ d(a n , a n ′ ) . ( 37
)
The RMDE is zero if and only if the estimate is perfect, while if the RMDE is close to one, the estimation error is of the same order of magnitude as the distance between true directions, indicating a very poor estimation quality. The generalisation of the RMDE for anechoic mixtures is done by replacing in Equations ( 36) and (37), the distance d(., .) with the distance d c (., .) defined in Equation (29).

C. Evaluation signals

We evaluated the clustering algorithms over the speech data in [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF], that is to say speech sources of duration 11.9 s, sampled at 8 kHz from 30 English speakers (males and females) from as many different audio books. For each configuration of the mixing parameters, we generated T = 10 mixtures from different source signals.

D. Method parameters 1) STFT parameters:

For each method, a STFT X(t, f ) is computed as a first step. Except for DUET, we combine different scales corresponding to frame sizes of 2 n samples, ranging from 128 samples (16 ms) to 65536 samples (8.2 s). For DUET, a frame size of 512 samples (64 ms) is used. We used a Hanning window with a half-frame overlap.

2) DEMIX specific parameters: For DEMIX-Instantaneous and DEMIX-Anechoic, we used the following parameters : the size of the neighborhood is set to |Ω| = 10, which corresponds to the optimum value in the tests of the TIFROM-CF method [START_REF] Puigt | Time-frequency ratio-based blind separation methods for attenuated and time-delayed sources[END_REF]. The value of the thresholds are set to : ζ = 3.3, ζ c = 9.5, and (for DEMIX-Anechoic only) ζ 2 = 2.33.

3) ELBG specific parameters: Since ELBG and its variants are randomly initialized, we run them I = 10 times for each test mixture and focus on the smallest error (RMDE) over these 10 runs, which thus gives an optimistic estimate of their performance.

E. Evaluations on Instantaneous mixtures

In this section, mixtures are obtained using instantaneous mixtures without noise, that is using Equation (1) with δ mn = 0, ∀m, n and n m (t) = 0, ∀ t.

1) Experimental protocol: First, we study the performance of the different algorithms depending on the number of sources, and second we fix the number of sources to three, and we vary the distance between these three sources.

a) 1 st experiment : N equally spaced sources: in the first experiment, noiseless linear instantaneous mixtures are generated with mixing matrices in the most favorable shape, that is where all directions are equally spaced (as in [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]), with a number of directions going from N = 2 to N = 10.

b) 2 nd experiment : 3 sources getting closer and closer: in the second experiment, 3 sources are placed with the following IP : θ l+2 = π 4 + l ∆θ π/180, with l ∈ {-1, 0, 1}. In this experiment, we only vary the angular distance ∆θ between IP in order to test the robustness of the algorithm when directions get close to each other (∆θ small).

2) Results: We observe (Tab I) that up to N = 8 sources, DEMIX estimates correctly the number of directions in more than seven cases out of ten, but when N > 9 it fails to count the number of sources.

As can be seen on Figures 5 and6, DEMIX-Instantaneous (DEMIX-Inst) yields a better performance than the best among I = 10 instances of the ELBG algorithms. The comparison of DEMIX-Instantaneous with DEMIX-Instantaneous-BR (DEMIX-Inst-BR) shows that the weighted mean using the confidence measure to estimate the directions as opposed to taking the "best" region, significantly improves the performance.

A remarkable fact is the behavior of DEMIX with three directions getting very close, compared to all other algorithms. The RMDE of the four ELBG variants approaches 1 when the distance between true directions gets close to zero (see Figure 6). In other words, the ELBG variants essentially confuse all directions. On the opposite, DEMIX remains very robust when the directions are very close to each other: as reported in Figure 6, the RMDE of DEMIX-Instantaneous (DEMIX-Inst-M) remains below 10 -3 until the directions get closer than 10 -3 degrees (∆θ < 10 -3 ). Also we notice, by observing the results for the four variants of the ELBG, that the ELBG results are not significantly improved if the local estimation θ(t, f ) of a direction are replaced by those obtained from local PCA θ(t, f ). The use of the confidence measure to "boost" the most reliable IP has a much more significant impact on the performance.

F. Evaluations on synthetic anechoic mixtures

We propose an experiment to test the limits of DEMIX-Instantaneous and the behavior of DEMIX-Anechoic as well as DUET on anechoic mixtures, by varying smoothly the degree of "anechoism" from a near instantaneous mixture to a "strong" anechoic mixture. The larger the angular distance ∆θ between IP, the larger the relative error of DUET. This is probably due to the fact that for the clustering step, DUET uses the ratio R 21 (τ , ω) = X 2 (τ , ω)/X 1 (τ , ω) which becomes unstable when the IP θ n get closer to π/2. Thus, for any delay, the performance of DUET is worse than that of DEMIX-Anechoic. Also, for any angular distance between IP (∆θ = 5 or 25 degrees), DEMIX-Anechoic performs better than DUET when the delay is lower than 0.1 sample, but especially when the delay δ is higher than one sample 3 .

G. Comparison between DEMIX-Anechoic and DUET on room simulated anechoic mixtures

As a third experiment, we compared the performance of the proposed DEMIX-Anechoic algorithm with the classical 3 Note that other methods like TIFROM [START_REF] Puigt | Time-frequency ratio-based blind separation methods for attenuated and time-delayed sources[END_REF] and TIFCORR [START_REF]A time-frequency correlation-based blind source separation method for time-delayed mixtures[END_REF] are able to estimate delays longer than only one sample. Thus, it will be interesting to evaluate all these methods within the framework of an evaluation campaign. DUET algorithm on anechoic mixtures obtained by an anechoic room simulation with the RoomSim MATLAB toolbox [START_REF] Campbell | A matlab simulation of "shoebox" room acoustics for use in research and teaching[END_REF].

1) Experimental protocol: the experimental protocol is the same as in [START_REF]A robust method to count and locate audio sources in a stereophonic linear anechoic mixture[END_REF], but applied to the signals described in section VII-C. So we invite the reader to refer to [START_REF]A robust method to count and locate audio sources in a stereophonic linear anechoic mixture[END_REF] for a detailed description of the protocol.

The experience consists in estimating the performance of algorithms by changing the number of sources from N = 2, to N = 7.

2) Results: a) counting the sources: DEMIX-Anechoic estimates the number of sources with more than 80% of success until N = 5 (see Tab II).

b) direction estimation: Tab III shows that the average RMDE of DEMIX-Anechoic is consistently better than that of DUET by a factor of at least 10.

Since the RMDE for DEMIX can only be measured when a correct number of sources is estimated, it was not computed when N > 6 with DEMIX-Anechoic.

H. Evaluation in reverberant situation

In this section, we evaluate the robustness of DEMIX-Anechoic and DUET when there is some reverberation, which corresponds to a more realistic recording situation than the instantaneous and anechoic situation.

1) Experimental protocol: similarly to the experimental protocol of section VII-F, we generated a sterephonic mixture of 3 sources with IP : θ l+2 = π 4 + l ∆θ π/180, with ∆θ ∈ {5, 25} degrees and with delays δ n ∈ {-1, 0, +1}, to enable DUET to correctly estimate the delay as follows.

The mixture was obtained as the sum of an anechoic part and a reverberant part :

x m (t) = 2) Results: Figure 8 shows that the RMDE between the estimated anechoic mixture directions and the anechoic part of the true mixture directions decreases when the SN R in gets lower than 0 dB, and the lower the SN R in , the higher the RMDE for the two methods. This is due to the fact that the convolutive model of Equation (38) differs more and more from the anechoic one of Equation (1) but probably also because adding reverberation decreases the sparsity of the source images in the TF domain. Also, whatever is the SN R in DEMIX-Anechoic obtains better results than DUET, whereas the absolute delay δ of the two side sources is set to 1, which is the most favorable configuration for the DUET method according to Figures 7(a

VIII. SUMMARY AND CONCLUSIONS

We have proposed a new approach to estimate the spatial directions of an unknown number of sources from a multichannel mixture in a possibly underdetermined and anechoic setting. The experiments on stereophonic recordings have illustrated the ability of the proposed method : a) to count the number of sources until 8 sources in the instantaneous case and 5 sources in the anechoic case; b) to robustly estimate very close directions that classical clustering algorithms like K-Means or ELBG failed to estimate; c) to estimate delays as large as 100 samples in simulated anechoic mixtures.

The proposed method relies on a statistical model of the mixture and exploits a certain level of sparsity of the timefrequency representations to extract local estimates of the directions. Our main contribution is the use of a confidence value to robustly estimate the mixing directions as well as the number of sources, together with a method similar to GCC-PHAT to estimate the time delays of anechoic mixtures. The method seems essentially limited by the fact that it relies on the assumption that each target source significantly "emerges" from the others in sufficiently many time-frequency regions. This condition is likely to fail when the mixture is made of "too many" sources or when the sources representations are not sparse enough. One way to deal with these cases for mixtures with M > 2 channels would be to replace the confidence measure by a measure which indicates the likelihood that at most M -1 sources are active. This would require adequate modifications of the clustering algorithm which may become significantly more complex. Another interesting perspective is to extend the present method to the convolutive case by considering sparse filters, the anechoic case being the special case of a 1-sparse filter. Yet, in the general case, the intensity difference of a direction at different frequencies would no longer be constant, and other techniques must be found to cluster directions estimated at different frequencies.

Fig. 1 .

 1 Fig. 1. Local scatter plots in two time-frequency regions. Lines indicate true source directions. STFT window size is L = 4096 and the size of the neighborhood is |Ω| = 99.

  (a) Real-valued global scatter plot of points X R (t, f )/ X R (t, f ) weighted by their energy ρ 2 (t, f ) = X R (t, f ) 2 . (b) Scatter plot of points û(Ω t,f ) weighted by their confidence measure 10 log 10 b T (Ω t,f ). Neighborhood size is |Ω t,f | = 10.

Fig. 2 .

 2 Fig.2. Comparison of the scatter plots of points used in the standard approach and the ones used by the PCA approach. STFT window size is L = 4096.

Figure 3 Fig. 3 .

 33 Figure 3 illustrates this process. 2.3) estimate the centroid u(C k ) using the regions from cluster C ′ k .

Fig. 4 .

 4 Fig. 4. Block diagram of the different tested algorithms and data flow

Fig. 5 .

 5 Fig. 5. Relative mean direction error (RMDE) as a function of the number of sources for DEMIX-Instantaneous (DEMIX-Inst), its variant DEMIX-Instantaneous-BR (DEMIX-Inst-BR) and the best instance (over 10) of the four variants of the ELBG

Fig. 6 .

 6 Fig. 6. Relative mean direction error (RMDE) as a function of the angular distance in degrees (∆θ) between the 3 directions for DEMIX-Instantaneous (DEMIX-Inst), its variant DEMIX-Instantaneous-BR (DEMIX-Inst-BR) and the best instance (over 10) of the four variants of the ELBG

  angular distance between IP is ∆θ = 25 degrees.

  angular distance between IP is ∆θ = 5 degrees.

Fig. 7 .

 7 Fig. 7. Relative mean direction error (RMDE) as a function of the absolute delay δ of the two side sources. nb of sources 2 3 4 5 6 7 DEMIX Anec 80 100 100 90 10 0

2 N 10 PFig. 8 .

 2108 Fig. 8. Relative mean direction error (RMDE) between the estimated anechoic mixture directions and the anechoic part of the true mixture directions, as a function of the SN R in , for values ∆θ = 25 degrees and ∆θ = 5 degrees of the angular distance between IP, and for an absolute delay δ = 1 of the two side sources.

  ) and 7(b).

TABLE II FREQUENCY

 II OF CORRECT COUNTS OF THE NUMBER OF SOURCES (IN %)

The proposed method has been submitted to the SISEC

evaluation campaign[START_REF] Vincent | The 2008 signal separation evaluation campaign: A community-based approach to large-scale evaluation[END_REF] for underdetermined instantaneous speech and music mixtures and obtained the best results for the mixing system estimation task.

We thank S. Rickard and C. Fearon for graciously providing their implementation of DUET[START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF].


M Bergounioux
email: maitine.bergounioux@univ-orleans.fr

L Piffet
email: loic.piffet@univ-orleans.fr

A second-order model for image denoising

Keywords: Second order total variation, image reconstruction, denoising, texture, variational method Mathematics Subject Classification (2000) 65D18, 68U10, 65K10

We present a variational model for image denoising and/or texture identification. Noise and textures may be modelled as oscillating components of images. The model involves a L 2 -data fitting term and a Tychonov-like regularization term. We choose the BV 2 norm instead of the classical BV norm. Here BV 2 is the bounded hessian function space that we define and describe. The main improvement is that we do not observe staircasing effects any longer, during denoising process. Moreover, texture extraction can be performed with the same method. We give existence results and present a discretized problem. An algorithm close to the one set by Chambolle [9] is used: we prove convergence and present numerical tests.

Introduction

Variational models in image processing have been extensively studied during the past decade. There are used for segmentation processes (geodesic or geometric contours) and restoration purpose as well. We are mainly interested in the last item which involves denoising or deblurring methods and textures extraction as well. Roughly speaking image restoration problems are severely ill posed and a Tychonov-like regularization process is needed. The general form of such models consists in the minimization of an "energy" functional :

F(u) = u -u d X + R(u) , u ∈ Y ⊂ X ,
where X, Y are (real) Banach spaces, R is a regularization operator, u d is the observed (or measured) image and u is the image to recover or to denoise. The first term is the fitting data term and the second one permits to get a problem which is no longer ill posed via a regularization process. The most famous model is the Rudin-Osher-Fatemi denoising model (see [START_REF] Acar | Analysis of bounded variation penalty methods for ill-posed problems[END_REF], [START_REF] Osher | Nonlinear total variation based noise removal algorithms[END_REF]). This model involves a regularization term that preserves discontinuities, what a classical H 1 -Tychonov regularization method does not. The observed image to recover is split in two parts u d = u + v where u represents the oscillating component (noise or texture) and v is the smooth part (often called the cartoon component). So we look for the solution as u + v with v ∈ BV (Ω) and u ∈ L 2 (Ω), where BV (Ω) is the functions of bounded variation space defined on an open set Ω ( [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization[END_REF]). The regularization term involves only the cartoon component v, while the remainder term u d -v represents the noise to be minimized. We get min v∈BV (Ω)

µ|v| BV (Ω) + 1 2 u d -v 2 L 2 (Ω) , (P ROF )
where µ > 0. This problem has a unique solution in BV (Ω). This functional space is the good one to deal with discontinuous functions that imply that the derivative (in the distribution sense) may be a measure (for example a Dirac measure). This model is used for denoising purpose. However, the use of the BV norm implies numerical perturbations. The computed solution turns to be piecewise constant which is called the "staircasing effect". Therefore, though noise can be successfully removed, the solution is not satisfactory. This variational model has been improved using different functional spaces, for the data fitting term or the regularizing term.

Recently people considered that an image can be decomposed into many components, each component describing a particular property of the image ( [START_REF] Aubert | Modeling very oscillating signals, application to image processing[END_REF][START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Osher | Image decomposition and restoration using total variation minimization and the H 1 norm[END_REF][START_REF] Osher | Modeling textures with total variation minimization and oscillating patterns in image processing[END_REF][START_REF] Osher | Image denoising and decomposition with total variation minimization and oscillatory functions. Special issue on mathematics and image analysis[END_REF]] and references therein for example). We shall assume as well that the image we want to recover from the data u d can be decomposed as f = u + v where u and v are functions that characterize different parts of f (see [START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Osher | Image denoising and decomposition with total variation minimization and oscillatory functions. Special issue on mathematics and image analysis[END_REF][START_REF] Yin | A comparison of three total variation based texture extraction models[END_REF] for example).

Components u and v belong to different functional spaces: u is the non regular part and belongs to L 2 (Ω) while v is a more regular part and belongs to BV 2 (Ω) (that we define in the sequel). The remainder term u = f -v involves the oscillating component (as noise and/or texture) and possibly contours. Such decompositions have been already performed [START_REF] Aubert | Modeling very oscillating signals, application to image processing[END_REF][START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF] using the so called Meyer-space of oscillating functions G [START_REF] Meyer | Oscillating patterns in image processing and nonlinear evolution equations[END_REF] instead of BV 2 (Ω). So far, the modelling we propose is not the same: the oscillating component will be a priori included in the non regular remainder term part u := f -v while v involves the cartoon. Our philosophy is different than the one used in Meyer approach. Here we perform a second order analysis to get sharper result so that we look for the smooth (cartoon) part. The oscillating part, texture and/or noise and possible contours will be part of the remainder term and are not modelled a priori.

The paper is organized as follows : in the next section, we present the functional framework and the space BV 2 (Ω) with useful properties. Section 3 is devoted to the variational model. In section 4, we focus on the discretization process. We present a Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] like algorithm in section 5 and numerical tests are reported in the last section.

2 The space BV 2 (Ω)

Let Ω be an open bounded subset of R n , n ≥ 2 (practically n = 2) smooth enough (with the cone property and C 1 for example). Following Demengel [START_REF] Demengel | Fonctions à hessien borné[END_REF], we define the space of bounded hessian functions that we call BV 2 (Ω). We first recall the definition and the main properties of the space BV (Ω) of functions of bounded variation (see [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization[END_REF][START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF] for example), defined by

BV (Ω) = {u ∈ L 1 (Ω) | Φ 1 (u) < +∞},
where

Φ 1 (u) := sup Ω u(x) div ξ(x) dx | ξ ∈ C 1 c (Ω), ξ ∞ 1 . (1)
The space BV (Ω), endowed with the norm u BV (Ω) = u L 1 + Φ 1 (u), is a Banach space. The derivative in the sense of distributions of every u ∈ BV (Ω) is a bounded Radon measure, denoted Du, and Φ 1 (u) = Ω |Du| is the total variation of u. We next recall standard properties of functions of bounded variation [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF][START_REF] Attouch | Variational analysis in Sobolev and BV spaces : applications to PDEs and optimization[END_REF].

u → Φ 1 (u) is lower semi-continuous (denoted in short lsc) from BV (Ω) to R + for the L 1 (Ω) topology. 3. BV (Ω) ⊂ L 2 (Ω) with continuous embedding, if n = 2. 4. BV (Ω) ⊂ L p (Ω) with compact embedding, for every p ∈ [1, 2), if n = 2.
Now we extend this definition to the second derivative (in the distributional sense). Recall that the Sobolev space is defined as

W 1,1 (Ω) = { u ∈ L 1 (Ω) | ∇u ∈ L 1 (Ω) }
where ∇u stands for the first order derivative of u (in the sense of distributions).

Definition 1 A function u ∈ W 1,1 (Ω) is Hessian bounded if Φ 2 (u) := sup Ω ∇u, div(ξ) R n | ξ ∈ C 2 c (Ω, R n×n ), ξ ∞ 1 < ∞,
where div(ξ) = (div(ξ 1 ), div(ξ 2 ), . . . , div(ξn)),

with ∀i,

ξ i = (ξ 1 i , ξ 2 i , . . . , ξ n i ) ∈ R n and div(ξ i ) = n k=1 ∂ξ k i ∂x k . BV 2 (Ω) is defined as the space of W 1,1 (Ω) functions such that Φ 2 (u) < +∞. Remark 1 If V = R n×n , ξ ∞ = sup x∈Ω n i,j=1 ξ j i (x) 2 .
We give thereafter many useful properties of BV 2 (Ω) (proofs can be found in [START_REF] Demengel | Fonctions à hessien borné[END_REF][START_REF] Piffet | Modèles variationnels pour l'extraction de textures 2D[END_REF]).

Theorem 1 The space BV 2 (Ω) endowed with the following norm

f BV 2 (Ω) := f W 1,1 (Ω) + Φ 2 (f ) = f L 1 + ∇f L 1 + Φ 2 (f ), (3)
where Φ 2 is given by [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF], is a Banach space.

Proposition 2 A function u belongs to BV 2 (Ω) if and only if u ∈ W 1,1 (Ω) and ∂u ∂x i ∈ BV (Ω) for i ∈ {1, . . . , n}. In particular

Φ 2 (u) n i=1 Φ 1 ∂u ∂x i n Φ 2 (u).
Remark 2 The previous result shows that

BV 2 (Ω) = u ∈ W 1,1 (Ω) | ∀i ∈ {1, . . . , n}, ∂u ∂x i ∈ BV (Ω) .
We get a lower semi-continuity result for the semi-norm Φ 2 as well.

Theorem 2 The operator Φ 2 is lower semi-continuous from BV 2 (Ω) endowed with the strong topology of

W 1,1 (Ω) to R. More precisely, if {u k } k∈N is a sequence of BV 2 (Ω) that strongly converges to u in W 1,1 (Ω) then Φ 2 (u) lim inf k→∞ Φ 2 (u k ). Remark 3 In particular, if lim inf k→∞ Φ 2 (u k ) < ∞, then u ∈ BV 2 (Ω).
We have embedding results as well:

Theorem 3 (Demengel [10]) Assume n ≥ 2. Then BV 2 (Ω) ֒→ W 1,q (Ω) with q n n -1 , (4)
with continuous embedding. Moreover the embedding is compact if q < n n-1 . In particular

BV 2 (Ω) ֒→ L q (Ω) for q n n -2 if n > 2 (5) BV 2 (Ω) ֒→ L q (Ω), ∀q ∈ [1, ∞[, if n = 2. ( 6
)
In the sequel, we set n = 2 and Ω is a subset of R 2 , so that BV 2 (Ω) ⊂ H 1 (Ω) with continuous embedding and BV 2 (Ω) ⊂ W 1,1 (Ω) with compact embedding.

The variational model

We now assume that the image we want to recover from the data u d can be decomposed as f = u + v where u and v are functions that characterize different parts of f (see [START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Osher | Image denoising and decomposition with total variation minimization and oscillatory functions. Special issue on mathematics and image analysis[END_REF][START_REF] Yin | A comparison of three total variation based texture extraction models[END_REF] for example). Components u and v belong to different functional spaces: v is the (smooth) second order part and belongs to BV 2 (Ω) while u is the remainder term that should involve noise and/or texture. We consider the following cost functional defined on BV 2 (Ω) :

F (v) = 1 2 u d -v 2 L 2 (Ω) + λΦ 2 (v) + δ v W 1,1 (Ω) , (7)
where λ, δ ≥ 0. We are looking for a solution to the optimization problem inf{

F (v) | v ∈ BV 2 (Ω) } (P)
The first term

u d -v 2 L 2 (Ω)
of F is the fitting data term. Here we have chosen the L 2norm for simplicity but any L p norm can be used (p ∈ [2, +∞)). We shall investigate in a future work the very case where p = 1; indeed we need to develop an approximation process to deal with this additional non differentiable term. Other terms are Tychonovlike regularization terms. The term λΦ 2 (v) + δ v W 1,1 (Ω) is nothing else that the BV 2 (Ω) norm of v. However, we have split it because the δ-part is not useful from the modelling point of view. It is only necessary to prove existence of solutions. We shall choose δ = 0 for numerical tests.

If the image is noisy, the noise is considered as a texture and will be involved in the remainder term u := u d -v: more precisely v will be the part of the image without the oscillating component, that is the denoised part and u is expected to involve noise, contours (and part of texture). In the sequel we shall focus on the denoising process. Such an approach has already been used by Hinterberger and Scherzer [START_REF] Hinterberger | Variational methods on the space of functions of bounded Hessian for convexification and denoising[END_REF] with the BV 2 (Ω) space. Their algorithm is different from the one we use here. Note that if we decide to split the function we look for in more than two components (f = u + v + w, where u is the texture part, v the cartoon and w the noise for example, with appropriate models ) then the minimization problem is a structured optimization problem to which standard decomposition methods can be applied (for example alternating minimization (Gauss Seidel) or parallel methods). Therefore, the study of (P) is quite significant from this point of view. For parallel methods one may consult [START_REF] Attouch | A Parallel Splitting Method for Coupled Monotone Inclusions[END_REF]. First we give a general existence and uniqueness result for problem (P).

Theorem 4 Assume that λ > 0 and δ > 0. Problem (P) has a unique solution v.

Proof.-We first prove existence by using the direct method in calculus of variations. Let vn ∈ BV 2 (Ω) be a minimizing sequence, i.e.

lim n→+∞ F (vn) = inf{ F (v) | v ∈ BV 2 (Ω) } < +∞.
The sequence (vn) n∈N is bounded in BV 2 (Ω). With the compactness result of Theorem 3, this yields that (vn) n∈N strongly converges (up to a subsequence) in W 1,1 (Ω) to v * ∈ BV 2 (Ω). Theorem 2 gives the following:

Φ 2 (v * ) lim inf n→+∞ Φ 2 (vn). So F (v * ) lim inf n→+∞ F (vn) = min v∈BV 2 (Ω) F (v),
and v * is a solution to (P). Uniqueness is straightforward with the strict convexity of F due to the term

u d -v 2 L 2 (Ω) .

Remark 4

It is an open question to know if we really need δ > 0. We could expect a Poincaré-Wirtinger inequality in the BV 2 (Ω)-space which is not very difficult to prove using appropriate "density" results. However, the boundary conditions we have to add are not clear (we have to deal with "second order" Neumann boundary conditions) .

It is not straightforward to get the suitable boundary conditions that allow to assert that the semi-norm Φ 2 is equivalent to the BV 2 (Ω)-norm. If we assume δ = 0 we are not sure to prove existence of a solution without additional assumptions and we lose uniqueness.

The discretized problem

Problem (P) can be written equivalently as inf

v∈BV 2 (Ω) u d -v 2 L 2 (Ω) 2λ + Φ 2 (v) + δ v W 1,1 (Ω) , ( P)
where δ has been replaced by δ λ . We are going to compute the solution numerically. We first present the discretization process.

Discretization of problem P

We assume for simplicity that the image is squared with size N × N . We note X := R N ×N ≃ R N 2 endowed with the usual inner product and the associated Euclidean norm

u, v X := 1 i,j N u i,j v i,j , u X := 1 i,j N u 2 i,j . (8)
We set Y = X × X. It is classical to define the discrete total variation as following (see for example [START_REF] Aubert | Mathematical Problems in Image Processing, Partial Differential Equations and the Calculus of Variations[END_REF]): the discrete gradient of the numerical image u ∈ X is ∇u ∈ Y :

(∇u) i,j = (∇u) 1 i,j , (∇u) 2 i,j , (9)
where (∇u

) 1 i,j = u i+1,j -u i,j if i < N 0 if i = N, and (∇u) 2 i,j = u i,j+1 -u i,j if j < N 0 if j = N. The (discrete) total variation corresponding to Φ 1 (u) is given by J 1 (u) = 1 i,j N (∇u) i,j R 2 , (10)
where (∇u

) i,j R 2 = (∇u) 1 i,j , (∇u) 2 i,j R 2 = (∇u) 1 i,j 2 + (∇u) 2 i,j 2 .
The discrete divergence operator div is the opposite of the adjoint operator of the gradient operator ∇:

∀(p, u) ∈ Y × X, -div p, u X = p, ∇u Y , so that (div p) i,j =      p 1 i,j -p 1 i-1,j if 1 < i < N p 1 i,j if i = 1 -p 1 i-1,j if i = N +      p 1 i,j -p 2 i,j-1 if 1 < j < N p 2 i,j if j = 1 -p 1 i,j-1 if i = N. (11)
To define a discrete version of the second order total variation Φ 2 we have to introduce the discrete Hessian operator. For any v ∈ X, the Hessian matrix of v, denoted Hv is identified to a X 4 vector: (Hv) i,j = (Hv) 11 i,j , (Hv) 12 i,j , (Hv) 21 i,j , (Hv) 22 i,j , with, for every i, j = 1, . . . , N,

(Hv) 11 i,j =    v i+1,j -2v i,j + v i-1,j if 1 < i < N, v i+1,j -v i,j if i = 1, v i-1,j -v i,j if i = N, (Hv) 12 i,j =    v i,j+1 -v i,j -v i-1,j+1 + v i-1,j if 1 < i N, 1 j < N, 0 if i = 1, 0 if i = N, (Hv) 21 i,j =    v i+1,j -v i,j -v i+1,j-1 + v i,j-1 if 1 i < N, 1 < j N, 0 if i = 1, 0 if i = N, (Hv) 22 i,j =    v i,j+1 -2v i,j + v i,j-1 if 1 < j < N, v i,j+1 -v i,j if j = 1, v i,j-1 -v i,j if j = N.
The discrete second order total variation corresponding to Φ 2 (v) is defined as

J 2 (v) = 1 i,j N (Hv) i,j R 4 . (12)
The discretized problem stands

inf v∈X u d -v 2 X 2λ + J 2 (v) + δ(|v| + J 1 (v)), (P d ) where |v| := 1 i,j N |v i,j | .
In the finite dimensional case we have an existence result with δ = 0.

Theorem 5 Problem P d has a unique solution for every δ ≥ 0 and λ > 0.

Proof.-The cost functional

F δ := u d -v 2 X 2λ + J 2 (v) + δ(|v| + J 1 (v)) ,
is continuous and coercive because of the term u d -v 2 X . In addition it is strictly convex so that we get the result.

For numerical purpose we shall set δ = 0. In fact, we have performed tests with δ = 0 and very small δ = 0 (as required by the theory to get a solution to problem P): results were identical.

In the sequel, we adapt the method by Chambolle in the BV (Ω)-case, to the second order framework. We briefly recall the original result of [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Consider the finite dimensional optimization problem derived from the discretization of the ROF model :

inf u∈X J 1 (u) + 1 2λ u d -u 2 X . (P 1 d )
The following result holds :

Proposition 3 The solution to (P 1 d ) is given by u = u d -P λK1 (u d ), ( 13
)
where

K 1 = {div g | g ∈ Y, g i,j R 2 1 ∀i, j = 1, • • • , N },
and P λK1 is the orthogonal projector operator on λK 1 .

We have a very similar result that we present in next sections. We first describe the conjugate function of J 2 4.2 The J 2 Legendre-Fenchel conjugate function

As the function J 2 is positively homogeneous, the Legendre-Fenchel conjugate

J * 2 (v) = sup u u, v X -J 2 (u),
is the characteristic function of a closed, convex set K

J * 2 (v) = 1 K (v) = 0 if v ∈ K +∞ else. As J * * 2 = J 2 , we get J 2 (v) = sup u∈K v, u X .
We use the inner scalar product of X 4 :

p, q X 4 = 1 i,j N p 1 i,j q 1 i,j + p 2 i,j q 2 i,j + p 3 i,j q 3 i,j + p 4 i,j q 4 i,j , for every p = p 1 , p 2 , p 3 , p 4 , q = q 1 , q 2 , q 3 , q 4 ∈ X 4 . So, for every v ∈ X,

J 2 (v) = sup p∈C p, Hv X 4 , (14)
where the feasible set is

C := { p ∈ X 4 | p i,j R 4 1, ∀ 1 ≤ i, j ≤ N } .
Let us compute the adjoint operator of H (which is the discretized "second divergence" operator) :

∀p ∈ X 4 , ∀v ∈ X H * p, v X = p, Hv X 4 .
We verify that H * : X 4 → X satisfies for every p = (p 11 , p 12 , p 21 , p 22 ) ∈ X 4

(H * p) i,j =        p 11 i-1,j -2p 11 i,j + p 11 i+1,j if 1 < i < N p 11 i+1,j -p 11 i,j if i = 1,
p 11 i-1,j -p 11 i,j if i = N, +        p 22 i,j-1 -2p 22 i,j + p 22 i,j+1 if 1 < j < N, p 22 i,j+1 -p 22 i,j if j = 1,
p 22 i,j-1 -p 22 i,j if j = N, +                                            p 12 i,j-1 -p 12 i,j -p 12 i+1,j-1 + p 12 i+1,j if 1 < i, j < N, p 12
i+1,j -p 12 i+1,j-1

if i = 1, 1 < j < N, p 12 i,j-1 -p 12 i,j if i = N, 1 < j < N, p 12 i+1,j -p 12 i,j if 1 < i < N, j = 1,
p 12 i,j-1 -p 12 i+1,j-1 if 1 < i < N, j = N, p 12
i+1,j

if i = 1, j = 1, -p 12 i+1,j-1 if i = 1, j = N, -p 12 i,j if i = N, j = 1,
p 12 i,j-1 if i = N, j = N, (15)
+                                            p 21 i-1,j -p 21 i,j -p 21 i-1,j+1 + p 21 i,j+1 if 1 < i, j < N, p 21 i,j+1 -p 21 i,j if i = 1, 1 < j < N, p 21 i-1,j -p 21 i-1,j+1 if i = N, 1 < j < N, p 21 i,j+1 -p 21 i-1,j+1 if 1 < i < N, j = 1,
p 21 i-1,j -p 21 i,j if 1 < i < N, j = N, p 21
i,j+1

if i = 1, j = 1, -p 21 i,j if i = 1, j = N, -p 21 i-1,j+1 if i = N, j = 1,
p 21 i-1,j if i = N, j = N, Theorem 6 The Legendre-Fenchel conjugate of J 2 is J * 2 = 1 K2
where

K 2 := {H * p | p ∈ X 4 , p i,j R 4 1, ∀i, j = 1, . . . , N } ⊂ X. (16)
Proof.-We have already mentioned that J * 2 = 1 K where K is of a non empty, closed, convex set subset of X. We first show that K 2 ⊂ K. Let be v ∈ K 2 . The discretized version of the definition of the second order total variation gives

J 2 (v) = sup ξ∈K2 ξ, v X ,
and ξ, v -J 2 (v) 0 for every ξ ∈ K 2 and v ∈ X. This gives for every

v * ∈ K 2 J * 2 (v * ) = sup v∈K2 v * , v -J 2 (v) 0. ( 17
)
As J * 2 takes only one finite value, we get J * 2 (v * ) = 0, which yields that v * ∈ K. Therefore, K 2 ⊂ K and as K is closed we finally obtain K2 ⊂ K.

In particular

J 2 (v) = sup ξ∈K2 v, ξ sup ξ∈ K2 v, ξ sup ξ∈K v, ξ = sup ξ∈K v, ξ -J * 2 (ξ) = J * * 2 (v).
As J * * 2 = J 2 , we have

sup ξ∈K2 v, ξ = sup ξ∈ K2 v, ξ = sup ξ∈K v, ξ . (18)
Now, let us assume there exists v * ∈ K such that v * / ∈ K2 . On may strictly separate v * and the closed convex set K2 . There exists α ∈ R and v 0 such that

v 0 , v * > α sup u∈ K2 v 0 , u .
We obtain sup

ξ∈K v 0 , ξ v 0 , v * > α sup u∈ K2 v 0 , u = sup u∈K v 0 , u ,
that gives a contradiction. We have proved that

K = K2 . As K 2 is closed we get K = K 2 .
Remark 5 We proved the previous result for the convenience of the reader. Indeed, one may note that J 2 is the support function of K 2 which is the conjugate function of the indicator function 1 K2 of K 2 (see [START_REF] Ekeland | Convex Analysis and Variational problems[END_REF] p. [START_REF] Meyer | Oscillating patterns in image processing and nonlinear evolution equations[END_REF]). Therefore, as K 2 is closed and convex and J 2 is continuous we get J * 2 = 1 * * K2 = 1 K2 .

Case where δ = 0

Now we focus on the case where δ = 0. Indeed, the δ-term in definition 7 was needed as a tool in the infinite dimensional framework to prove existence of solutions to (P). However, we have seen that the finite dimensional problem has a solution even if δ = 0, which is the most interesting case. The problem we have to solve turns to be

min v∈X u d -v 2 X 2λ + J 2 (v). (P 2 d )
As in the BV-case (proposition 3) we have a characterization of the solution.

Theorem 7 The solution v of (P 2 d ) verifies

v = u d -P λK2 (u d ),
where P λK2 is the orthogonal projector operator on λK 2 .

Proof.-The proof is similar to the one by Chambolle [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] but we give it for convenience. The solution v to (P 2 d ) is characterized by

0 ∈ ∂ J 2 (v) + 1 2λ v -u d 2 2 = v -u d λ + ∂J 2 (v), that is u d -v λ ∈ ∂J 2 (v).
As J 2 is proper, convex and continuous, then

v * ∈ ∂J 2 (v) ⇐⇒ v ∈ ∂J * 2 (v * ). So v ∈ ∂J * 2 ( u d -v λ ) ⇐⇒ 0 ∈ -v + ∂J * 2 ( u d -v λ ) ⇐⇒ 0 ∈ u d -v λ - u d λ + 1 λ ∂J * 2 ( u d -v λ ).
This means that w

= u d -v λ is the solution to min w 1 2 w - u d λ 2 X + 1 λ ∂J * 2 (w). As J * 2 = 1 K2 this implies that u d -v λ is the orthogonal projection of u d λ on K 2 : u d -v λ = P K2 ( u d λ ) . As P K2 ( u d λ ) = 1 λ P λK2 (u d )
, this ends the proof.

A fixed-point algorithm to compute ∂J 2

In [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] Chambolle proposed a fixed-point algorithm to compute P λK1 (f ) (and ∂J 1 (f )).

Let us briefly recall the main idea that we use again in the BV 2 context. To compute P λK1 (f ) we have to solve

min{ λdiv p -u d 2 X | p i,j R 2 1 ∀i, j = 1, . . . , N }. (19)
that can be solved with a fixed-point method :

p 0 = 0 (20a) p n+1 i,j = p n i,j + τ (∇(div p n -u d /λ)) i,j 1 + τ (∇(div p n -u d /λ)) i,j R 2 (20b)
In addition, a convergence result is provided :

Theorem 8 ([9], Theorem 3.1) Assume that τ satisfies τ 1/8. Then λdiv p n converges to P λK1 (u d ) as n → +∞.

Remark 6 In [8] J.F. Aujol proves that the modified algorithm :

p 0 = 0 (21a) p n+1 i,j = p n i,j + τ (∇(div p n -u d /λ)) i,j max (1, p n i,j + τ (∇(div p n -u d /λ)) i,j R 2 ) (21b) converges if τ 1/4.
We extend this result to the second-order case. To compute P λK2 (u d ) we have to solve

min λH * p -u d 2 X | p ∈ X 4 , p i,j 2 R 4 -1 0, i, j = 1, . . . , N . (P ′ ) Let us denote F (p) = λH * p -u d 2 X and g i,j (p) = p i,j 2
R 4 -1 = (p 11 i,j ) 2 + (p 12 i,j ) 2 + (p 21 i,j ) 2 + (p 22 i,j ) 2 -1.
First order optimality conditions give the existence of Lagrange multipliers α i,j , (i, j) ∈ {1, . . . , N } 2 , such that

∇F (p) + N i,j=1 α i,j ∇g i,j (p) = 0, ( 22a
)
α i,j 0 and α i,j g i,j (p) = 0, (i,

j) ∈ {1, . . . , N } 2 . ( 22b
)
It is easy to see that

∇F (p) = 2λH [λH * p -u d ] and that N i,j=1
α i,j ∇g i,j (p) = 2α i,j (p 11 i,j , p 22 i,j , p 12 i,j , p 21 i,j )

1 i,j N .
Therefore relations ( 22) are equivalent to

∀1 ≤ i, j ≤ N H λH * p -u d i,j + α i,j p i,j = 0, ( 23a
)
∀1 ≤ i, j ≤ N α i,j 0 and α i,j g i,j (p) = 0. (23b)

Let us compute the multipliers α i,j more precisely :

-If α i,j > 0 then p i,j R 4 = 1. -If α i,j = 0 then (H [λH * p -u d ]) i,j = 0.
In both cases we get

∀1 ≤ i, j ≤ N α i,j = H λH * p -u d i,j R 4
and we finally obtain the following equality :

∀(i, j) ∈ {1, . . . , N } 2 , H λH * p -u d i,j + H λH * p -u d i,j R 4 p i,j = 0. ( 24
)
We use a semi-implicit gradient method to solve these equations, namely : Choose τ > 0, and 1. Set p 0 = 0, n = 0 2. p n is known. For (i, j) ∈ {1, . . . , N } 2 , compute p n+1 i,j with

p n i,j = p n+1 i,j + τ H H * p n -u d /λ i,j + H H * p n -u d /λ i,j R 4 p n+1 i,j
. This is equivalent to

p n+1 i,j = p n i,j -τ (H [H * p n -u d /λ]) i,j 1 + τ (H [H * p n -u d /λ]) i,j R 4 . ( 25
)
The algorithm step τ is related to the adjoint operator H * norm that we call κ in the sequel and we first give a κ estimate:

Lemma 1 The adjoint operator H * norm, κ satisfies κ ≤ 8 .
Proof.-The definition of κ reads κ = sup

p X 4 1 H * p . As H * p X = sup q∈ BX (0,1) H * p, q X and ∀q ∈ X 4 H * p, q X = p, Hq X 4 Hq X 4 p X 4 ,
we get

H * p X |||H||| p X 4 . (26)
For any q ∈ X 4

Hq 2 X 4 = 1 i,j N q i+1,j -2q i,j + q i-1,j 2 + q i,j+1 -q i,j -q i-1,j+1 + q i-1,j 2 + q i+1,j -q i,j -q i+1,j-1 + q i,j-1 2 + q i,j+1 -2q i,j + q i,j-1

2 4 1 i,j N q 2 i+1,j + q 2 i,j + q 2 i,j + q 2 i-1,j + q 2 i,j+1 + q 2 i,j + q 2 i-1,j+1 + q 2 i-1,j
+q 2 i+1,j + q 2 i,j + q 2 i+1,j-1 + q 2 i,j-1 + q 2 i,j+1 + q 2 i,j + q 2 i,j + q 2 i,j-1 4 × 16 q 2 X .

Finally |||H||| 8, and with relation ( 26), H * p X 8 p X 4 .

We deduce that κ 8 (and κ 2 64).

Theorem 9 Let be τ 1/64. Then λ (H * p n ) n converges to P λK2 (u d ) as n → +∞.

Proof.-It is easy to check (by induction) that for every n 0 and i, j, p n i,j R 4

1. Set n > 0 and let be η n = (p n -p n+1 )/τ . Denote by (Un)n the sequence defined by :

Un = H * p n -u d /λ 2 X , ∀n 0.
The outline of the proof is the following : we first prove that (Un)n is a (nonnegative) non increasing sequence and therefore it is convergent. The proof provides an equality (29) that allows to show that the sequences (p n )n and (p n+1 )n have the same cluster points which are the unique solution to a projection problem. Therefore the whole sequence (p n )n converges.

• We first prove that (Un)n is a non increasing sequence if τ 1/κ 2 where κ is the norm of H * . We get

U n+1 = H * (p n -τ η n ) -u d /λ 2 X = -τ H * η n + H * p n -u d /λ , -τ H * η n + H * p n -u d /λ X = Un -2τ H * η n , H * p n -u d /λ X + τ 2 H * η n 2 X = Un -2τ η n , H H * p n -u d /λ X 4 + τ 2 H * η n 2 X Un -τ 2 η n , H H * p n -u d /λ X 4 -κ 2 τ η n 2 X 4 Un -τ   N i,j=1 2η n i,j , H H * p n -u d /λ i,j R 4 -κ 2 τ η n i,j 2 R 4   As η n i,j = H H * p n -u d /λ i,j + ρ n i,j , with ρ n i,j = H λH * p n -u d /λ i,j R 4 p n+1 i,j ;
we obtain

2η n i,j , H H * p n -u d /λ i,j R 4 -κ 2 τ η n i,j 2 R 4 = η n i,j , H H * p n -u d /λ i,j R 4 + η n i,j , H H * p n -u d /λ i,j R 4 -κ 2 τ η n i,j 2 R 4 = η n i,j , η n i,j -ρ n i,j R 4 + H H * p n -u d /λ i,j + ρ n i,j , Hl H * p n -u d /λ i,j R 4 -κ 2 τ η n i,j 2 R 4 = η n i,j 2 R 4 -η n i,j , ρ n i,j R 4 + H H * p n -u d /λ i,j 2 R 4 + ρ n i,j , H H * p n -u d /λ i,j R 4 -κ 2 τ η n i,j 2 R 4 = η n i,j 2 R 4 -η n i,j , ρ n i,j R 4 + H H * p n -u d /λ i,j 2 R 4 + ρ n i,j , η n i,j -ρ n i,j R 4 -κ 2 τ η n i,j 2 R 4 = (1 -κ 2 τ ) η n i,j 2 R 4 + H H * p n -u d /λ i,j 2 R 4 -ρ n i,j 2 R 4 .
Finally

U n+1 Un-τ   N i,j=1 (1-κ 2 τ ) η n i,j 2 R 4 + H H * p n -u d /λ i,j 2 R 4 -ρ n i,j 2 R 4   (27)
As

p n+1 i,j 1, ρ n i,j R 4 (H [λH * p n -u d /λ]) i,j R 4 this yields that if τ 1/κ 2 ,
then the sequence (Un) is non increasing. Thus, the sequence (Un)n is convergent to some m. Moreover, relation (27

) gives 0 ≤ N i,j=1 (1-κ 2 τ ) η n i,j 2 R 4 + H H * p n -u d /λ i,j 2 R 4 -ρ n i,j 2 R 4 ≤ Un -U n+1 τ . (28)
Passing to the limit in relation (28) gives

lim n→+∞ N i,j=1 (1-κ 2 τ ) η n i,j 2 R 4 + H H * p n -u d /λ i,j 2 R 4 -ρ n i,j 2 R 4 = 0. (29)
• Let us prove that we may extract a subsequence (p n k ) k of (p n )n such that (p n k ) k and (p n k +1 ) k have the same limit.

As (p n )n is bounded, we call p one cluster point and (p n k ) k the corresponding subsequence of (p n )n such that p n k → p, as k → +∞. Let us call p′ the limit (up to a subsequence) of p n k +1 . We obtain with relation ( 25)

p′ i,j = pi,j -τ (H [H * p -u d /λ]) i,j 1 + τ (H [H * p -u d /λ]) i,j R 4 . ( 30
)
We note ρ et η the respective limits of ρ n k et η n k when k tends to +∞. With (29), we obtain

N i,j=1 (1-κ 2 τ ) ηi,j 2 R 4 + H H * p-u d /λ i,j 2 R 4 -ρi,j 2 R 4 = 0.
As the terms of the sum are nonnegative, we get

∀i, j (1 -κ 2 τ ) ηi,j 2 R 4 = 0 and H H * p -u d /λ i,j 2 R 4 -ρi,j 2 R 4 = 0.
-If κ 2 τ < 1, then ηi,j = 0 for all i, j, and so p′ = p.

-If κ 2 τ = 1, then, for all i, j, ρi,j

R 4 = (H [H * p -u d /λ]) i,j R 4 , that is to say H H * p -u d /λ i,j R 4 p′ i,j R 4 = H H * p -u d /λ i,j R 4 .
This implies p′ i,j R 4 = 1 or H H * p -u d /λ i,j R 4 = 0.

-If (H [H * p -u d /λ]) i,j R 4 = 0, then relation (30) gives p′ i,j = pi,j .

-If p′ i,j R 4 = 1, then 1 = pi,j -τ (H [H * p -u d /λ]) i,j R 4 1 + τ (H [H * p -u d /λ]) i,j R 4 pi,j R 4 + τ (H [H * p -u d /λ]) i,j R 4 1 + τ (H [H * p -u d /λ]) i,j R 4
, therefore pi,j R 4 1 and (together with pi,j R 4 1) pi,j R 4 = 1. We deduce that pi

,j -τ H H * p -u d /λ i,j R 4 = pi,j R 4 + τ H H * p -u d /λ i,j R 4 .
Since triangular inequality turns to be an equality, there exists β ∈ R * so that τ H H * p -u d /λ i,j = β pi,j .

As p′

i,j = 1, relation (30) implies p′ i,j = 1 -β 1 + |β| pi,j , so that 1 -β 1 + |β| = 1; this yields β 0 and p′ i,j = pi,j .

Finally, p = p′ , and ∀i, j, H λH * p -u d i,j + H λH * p -u d i,j R 4 pi,j = 0. This is the Euler equation for (P ′ ). Therefore p is a solution to (P ′ ). With uniqueness of the projection, we deduce that the sequence (λH * p n )n converges to P λK2 (u d ).

We conclude with Lemma 1 since τ 1/64 then τ ≤ 1/κ 2 .

Numerical results

In this section we briefly present numerical tests for the denoising process. A full comparison with existing methods will be performed in a forthcoming paper.

Examples

Throughout this section, we consider the following images that are degraded with a white Gaussian noise with standard deviation σ: We perform numerical tests with different values of σ. In any case, we observe that the model is quite efficient for image restoration. Moreover, we note that we lose details information when parameter λ increases, what was expected. However, especially when the "observed" image is very noisy, we have a blurriness (subjective) feeling, that we do not have when restoration is performed with the standard ROF model. Checking what happens precisely on slices (lines) of the image (Figure 8 for example), we remark that the BV 2 -model keeps contour information pretty well, anyway better than expected watching the image.

Numerical tests have been performed on different machines so that we cannot report precisely on the CPU time. However, the result is quite satisfactory after few iterations so that the process is not too slow. In what follows, the stopping criterion has been set to a maximal number of iterations itmax that can be chosen arbitrary large. The algorithms have been implemented with MATLAB c software. We give al so the Signal to Noise Ratio (SNR) 1

: SN R(v) = 20 log 10 u L 2 u -v L 2
, where u is the expected image and v is the restored one. SN R(u d ) gives the observed SNR (with the noisy input image).

We have performed tests for two σ values. In the first example σ = 0.15 and we stopped after itmax=5000 iterations (Figure 3) and in the second case σ = 0.25 (the noise is more important). The noisy images are below: 1 There are different ways to compute this quantity with MATLAB c . We used the following syntax : SNR = 20*log10(norm(u(:))/norm(u(:)-v(:))) As expected, we see on Figure 3 that the smoothing process is more efficient when λ is large. For both images, the result is satisfactory for λ ≃ 10.

Sensibility with respect to λ parameter

Sensitivity with respect to iterations number itmax

We fix λ = 15 and σ = 0.15. Figures 6 and7 give the behavior of a slice (line) during iterations (we can see more easily how noise is removed). The algorithm converges well: the quality of restoration is improved as the number of iterations grows. Noise is removed and contours are preserved. It is well known that the ROF model makes staircasing effect appear, since the resulting image is piecewise constant on large areas. We first compare the two models on test images that are not very noisy. In Figure 9 we see that piecewise constant areas appear with ROF, which is not the case with the BV 2 model. To focus on this fact, we have selected a line of the image that meets contours. precisely what happens: the image restored with ROF is clearly piecewise constant, and the BV 2 restored one seems to be blurred (Figure 12). However, this is an optical effect: considering a slice shows that the BV 2 model removes noise significantly and contours are better preserved: the amplitude of high peaks that correspond to contours is not changed, which is not the case in ROF-model (Figure 11).

Texture extraction

We do not report much on texture extraction process. The parameter tuning is slightly different but the algorithm behaves similarly (see [START_REF] Echegut | A variational model for image texture identification[END_REF]). Many variational methods have been developed for texture extraction (see [START_REF] Aubert | Modeling very oscillating signals, application to image processing[END_REF][START_REF] Aubert | Image decomposition into a bounded variation component and an oscillating component[END_REF] and the references therein). We shall precisely compare the BV 2 method to the existing ones in a forthcoming paper. We present an example in Figures 13 and14. The second order approach via the BV 2 space seems promising. Many improvements have to be performed. The algorithm is still slow (though we get acceptable results for quite few iterations ≃ 30). We currently investigate the Nesterov-techniques to speed up the method (see Weiss et al. [START_REF] Weiss | Efficient schemes for total variation minimization under constraints in image processing[END_REF] or Fadili-Peyré [START_REF] Fadili | Total Variation Projection with First Order Schemes[END_REF]). Moreover, we have to look for modifications of the variational model using different norms (for example the L 1 norm) for the fitting data term. Furthermore, coupling existing techniques for texture extraction with the second order approach should give quite performing results. This will be done in future works.

Fig. 1

 1 Fig. 1 Test images ("Shapes" (size : 740 × 740) and "Lena" (size : 512 × 512) )

Fig. 2

 2 Fig. 2 Noisy images -white Gaussian noise with standard deviation σ.

Fig. 3 Fig. 4

 34 Fig. 3 Solution v -Standard deviation σ = 0.15 -5 000 iterations

Fig. 5 15 (Fig. 6 Fig. 7

 51567 Fig. 5 Noisy images with lineσ = 0.15

  (a) Noisy image (SNR 16.9) (b) Noisy slice (c) ROF model solution (SNR 23.82) (d) ROF model slice (e) BV 2 model solution -(SNR 25.5) (f) BV 2 model slice

Fig. 8 Fig. 9

 89 Fig. 8 Comparison between ROF and BV 2 modelsσ = 0.25, λ = 25, 5 000 iterations.

Figures 8 and 11

 11 Figures 8 and 11 are obtained for λ = 25 and λ = 50 respectively and 5 000 iterations. Figures 10 and 12 are obtained for 10 000 iterations and λ = 50 and showprecisely what happens: the image restored with ROF is clearly piecewise constant, and the BV 2 restored one seems to be blurred (Figure12). However, this is an optical effect: considering a slice shows that the BV 2 model removes noise significantly and contours are better preserved: the amplitude of high peaks that correspond to contours is not changed, which is not the case in ROF-model (Figure11).

  (a) Noisy slice (b) Original slice (c) BV 2 model -50 iterations. (d) BV 2 model -5 000 iterations.(e) ROF model -5 000 iterations.

Fig. 11

 11 Fig. 11 Zoom on "Lena" slicesσ = 0.25, λ = 50

( a ) 5 Fig. 14

 a514 Fig. 13 Texture extraction -Rescaled texture image u d -v at 10 and 15 iterations

  Proposition 1 Let Ω be an open subset of R n with Lipschitz boundary. 1. For every u ∈ BV (Ω), the Radon measure Du can be decomposed into Du = ∇u dx + D s u, where ∇u dx is the absolutely continuous part of Du with respect of the Lebesgue measure and D s u is the singular part. 2. The mapping