Gaël Mahé

Mériem Jaïdane

Perceptually Controlled Reshaping of Sound Histograms

Keywords: sound histogram equalization, histogram global reshaping, histogram local reshaping, noise audibility control, sparsification, histogram low-pass filtering AUD-SEN, AUD-AMCT, AUD-

Many audio processing algorithms have optimal performance for specific signal statistical distributions that may not be fulfilled for all signals. When the original signal is available, we propose to add an inaudible noise so that the distribution of the signal-plus-noise mixture is as close as possible to a given target distribution. The proposed generic algorithm (independent from the application) adds iteratively a low-power white noise to a flat-spectrum version of the signal, until the target distribution or the noise audibility is reached. The latter is assessed through a frequency masking model. Two implementations of this sound reshaping are described, according to the level of the targeted transformation and to the foreseen application: Histogram Global Reshaping (HGR) to change the global shape of the histogram and Histogram Local Reshaping (HLR) to locally "chisel" the histogram, but keeping the global shape unchanged. These two variants are illustrated by two applications where the inaudibility of the noise generated by the algorithm is required: "sparsification" for source separation, and low-pass filtering of the histogram for application of the quantization theorem, respectively. In both cases, the target histogram is reached or almost reached and the transformation is inaudible. The experiments show that the source separation performs better with HGR and that the HLR allows a better application of the quantization theorem.

I. INTRODUCTION

H ISTOGRAM equalization (HE) is well-known in image processing [START_REF]Digital Image Processing[END_REF], and is mainly used to enhance the contrast of an image. This kind of processing does not belong to the traditional toolbox of sound processing, because the "contrast" of a sound is not clearly defined and may not be enhanced through simple HE.

However, some particular audio applications are improved with a histogram equalization. A correction of non-linear distortions based on HE was proposed in [START_REF] White | Restoration of nonlinearly distorted audio by histogram equalization[END_REF]. HE was also used as a pre-processing step in automatic speech recognition, in order to increase the robustness of the latter against noise (see for example [START_REF] Hilger | Quantile based histogram equalization for noise robust large vocabulary speech recognition[END_REF], [START_REF] Al-Wakeel | Stereo-based histogram equalization for robust speech recognition[END_REF]). In all these works, HE was conceived as an enhancement processing, where the processed sound is aimed at being perceptually different from the original one.

Gaël Mahé is with the Laboratory of Informatics Paris Descartes (LIPADE) in Paris Descartes University, 45 rue des Saints Pères 75270 Paris cedex 06, France. e-mail: gael.mahe@mi.parisdescartes.fr -tel: +33 1 76 53 02 82 Mériem Jaïdane is with the Signals and Systems Lab (U2S) in the National School of Engineering of Tunis (ENIT), Université de Tunis El Manar, BP 37, 1002 Tunis-Belvédère, Tunisia. e-mail: meriem.jaidane@enit.utm.tn -tel: +216 71 874 700

In a different approach, the goal of HE is to enhance a specific later processing step that needs particular statistical properties, while it is essential to keep the perceptual properties of the signal unchanged. In [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], the audio signal is "Gaussianized" in order to better identify the non-linear system that will convey it. In [START_REF] Balazs | Time-frequency sparsity by removing perceptually irrelevant components using a simple model of simultaneous masking[END_REF]- [START_REF] Chardon | Perceptual matching pursuit with Gabor dictionaries and time-frequency masking[END_REF], the time-frequency coefficients of audio signals are "sparsified"i.e. the amount of zeros is increased [START_REF] Balazs | Time-frequency sparsity by removing perceptually irrelevant components using a simple model of simultaneous masking[END_REF], [START_REF] Pinel | Sparsification' of audio signals using the MDCT/IntMDCT and a psychoacoustic model -application to informed audio source separation[END_REF], [START_REF] Chardon | Perceptual matching pursuit with Gabor dictionaries and time-frequency masking[END_REF] or the shape parameters of their Generalized Gaussian distributions are reduced [START_REF] Mahé | Perceptually controlled doping for audio source separation[END_REF] -in order to enhance source separation [START_REF] Pinel | Sparsification' of audio signals using the MDCT/IntMDCT and a psychoacoustic model -application to informed audio source separation[END_REF], [START_REF] Mahé | Perceptually controlled doping for audio source separation[END_REF] or audio-coding [START_REF] Balazs | Time-frequency sparsity by removing perceptually irrelevant components using a simple model of simultaneous masking[END_REF], [START_REF] Chardon | Perceptual matching pursuit with Gabor dictionaries and time-frequency masking[END_REF]. In [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF], the histogram is low-pass filtered in order to fulfill the conditions of the quantization theorem [START_REF] Widrow | A Study of Rough Amplitude Quantization by Means of Nyquist Sampling Theory[END_REF] and restore the histogram of the signal from that of the sub-quantized signal. In all these applications, the transformation must be perceptually transparent.

While the algorithms operating on time-frequency distributions [START_REF] Balazs | Time-frequency sparsity by removing perceptually irrelevant components using a simple model of simultaneous masking[END_REF]- [START_REF] Chardon | Perceptual matching pursuit with Gabor dictionaries and time-frequency masking[END_REF] control the audibility of the transformation through perceptual models, none of the algorithms operating on time-domain samples [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF] provides a satisfactory perceptual control of the transformation.

Hence, our goal is to propose a generic algorithm that reshapes the histogram of the temporal samples of an audio sequence into a target histogram, under an explicit constraint of inaudibility of the additive noise generated by the transformation. The proposed algorithm is intended for any application where the original signals to be processed are available and should have a specific distribution for further optimal processing. This distribution is for example Gaussian for non-linear system identification [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], sparse for source separation [START_REF] Comon | Handbook of Blind Source Separation[END_REF], band-limited for application of the quantization theorem [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF], adapted to the error minimization in the optimal quantization [START_REF] Gersho | Principles of quantization[END_REF].

Two levels of reshaping are targeted: changing the global shape of the distribution, as in [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF]- [START_REF] Chardon | Perceptual matching pursuit with Gabor dictionaries and time-frequency masking[END_REF], or reshaping locally the histogram, as in [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF], assuming that the global shape is satisfactory.

In Section II, we will present our algorithm and discuss its features. Two versions of this algorithm will be described in Sections III-A and III-B, according to the goal of the reshaping: Histogram Global Reshaping (HGR) and Histogram Local Reshaping (HLR), respectively. Finally, some experimental results with speech and music will be presented in Section IV, related to "sparsification" for a generic algorithm of blind source separation in the time-domain, and low-pass filtering of the histogram for application of the quantization theorem.

II. PRINCIPLES OF THE ALGORITHM

A. Framework and goals

Let x be a discrete signal of finite length N , with integer values x(0) . . . x(N -1). Considering histograms whose classes are integer numbers, the histogram and the cumulative histogram of x are defined by:

∀ k ∈ Z, f x (k) = |{x(n) | x(n) = k}| F x (k) = k i=-∞ f x (i) = |{x(n) | x(n) ≤ k}|, (1)
respectively, where |S| denotes the cardinality of a set S.

Let f target be the target histogram. Our goal is to find a transformation of x into z (with integer values, too) so that:

f z ≃ f target (2) w = z -x is inaudible with reference to x (3)
where w denotes the equivalent additive transformation noise. We assess its inaudibility through a frequency masking constraint based on an frequency analysis on successive frames of duration around 20 ms. Hence, reaching (2) and (3) corresponds to solving:

min d(f z , f target ) such that : (4) γ w (m, ν) < γ mask (m, ν) ∀ frame m, frequency ν (5)
where d denotes a given distance (d TV for total variation, d KS for Kolmogorov-Smirnov...), γ w (m, ν) denotes the power spectral density of the transformation noise w in the m th frame, and γ mask (m, ν) the masking threshold of x in the m th frame, in the frequency domain. In this paper, we will consider a classical frequency-masking model [START_REF] Zwicker | The ear as a communication receiver[END_REF], [START_REF]Coding of moving pictures and associated audio for digital storage media at up to about 1,5 Mbit/s -Part 3: Audio[END_REF].

The proposed goal expressed by Eq. ( 4) and ( 5) exhibits two challenging constraints:

• the histogram optimization has to be performed on the whole signal, while the constraint is local and different for each frame; • the histogram is the histogram of the time samples, while the constraints are expressed in the frequency domain.

B. State of the art

The principle of HE in [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF] is to move the samples between neighboring classes of the histogram, from classes in excess to deficient classes as follows. Initially, z = x. For j varying from the minimum to the maximum values of x,

• If f z (j) -f target (j) = M > 0, M
samples of z of value j are randomly selected. Each of those samples gets the value j + 1, so that f z (j) = f target (j) and

f z (j + 1) = f z (j + 1) + M . • If f z (j) -f target (j) = -M < 0, M
samples of z of value j + 1 are randomly selected. Each of those samples gets the value j. If f z (j + 1) < M , the missing samples are randomly selected among those of value j + 2, and so on, until f z (j) = f target (j). At the end of this algorithm, the target histogram is exactly reached, but the algorithm does not provide any control on the audibility of the transformation. In the application that [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF] dealt with, i.e. low-pass filtering histograms, the lower the cut-off frequency was, the higher the transformation noise w was. It was only shown on some experimental examples that w was inaudible for cut-off frequencies above a threshold depending on the signal.

The first version of HE in [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], for "Gaussianization", is based on the same algorithm as in image processing. In the time-domain the principle is to find for each sample x(n) the sample z(n) so that the target cumulative distribution function F target in z(n) matches the empirical cumulative distribution function in x(n), i.e.:

F target z(n) = F x x(n) (6)
Again, this basic algorithm makes the histogram f z match exactly the target histogram f target , but does not control the audibility of the transformation of x into z.

In order to make the transformation noise w inaudible, [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] added a constraint to [START_REF] Balazs | Time-frequency sparsity by removing perceptually irrelevant components using a simple model of simultaneous masking[END_REF], leading to the following optimization problem:

∀ n , min |F target z(n) -F x x(n) | (7) |w(n)| < w max (8)
where w max was set so that the variance of w estimated on the whole sequence was equal to a target variance empirically chosen to ensure inaudibility. The drawback of this method is that it does not explicitly shape w according to a perceptual model, as expressed by Eq. ( 5) for example, which may lead to a sub-optimal trade-off between reaching the target histogram and ensuring the inaudibility of the transformation.

A perceptual control of these algorithms would require to fulfill the condition (5) on γ w (m, ν). In the time domain, this spectral shaping implies introducing a particular correlation between successive values w(n). This is not possible in Algorithm [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF] (described at the beginning of this subsection) where the samples are not processed sequentially. In Algorithm [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], described above, Equations ( 7) and (8) constitute a local optimization problem, sample by sample. Replacing (8) by ( 5) would involve neighboring samples, which makes it difficult to keep [START_REF] Pinel | Sparsification' of audio signals using the MDCT/IntMDCT and a psychoacoustic model -application to informed audio source separation[END_REF] local. The proposition that follows can be seen as a generalization of [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] that takes into account the interdependence between neighboring samples induced by the frequency-domain constraint [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF].

C. Proposed perceptually controlled histogram transformation

To circumvent the difficulty of taking into account the dependencies between each sample w(n) and its neighbors, we propose to formulate the problem in a domain of representation where the transformation noise is white.

For this purpose, we propose a process based on the flattening/recoloring scheme1 of Fig. 1, where T is a transformation (possibly non-linear) controlled by the distance between f z and f target , and assumed to be equivalent to the addition of a white noise. Hence, the choice of H -1 allows to shape in the frequency domain the transformation noise w = zx according to (5), while the specification (4) is obtained through the control of T by d(f z , f target ). We propose to implement this general scheme through the iterative algorithm illustrated by Fig. 2.

Initially,

• the original signal x is filtered by a frame-varying flattening filter H m controlled by a psychoacoustic model of x:

e x = h m * x, (9)
for each m th frame, where h m is the impulse response of H m and * means discrete convolution. While the psychoacoustic model is computed on 50%-overlapping frames, the filter coefficients are defined on nonoverlapping frames. Hence, each processed frame m corresponds to the central half of the samples of the m th analysis frame. • z = x and e z = e x . Then, for each iteration i (one cycle) and each discrete time n, we generate a random value δw e (i, n) ∼ N (0, σ 2 i ). Adding δw e (i, n) to e z (n), the n th sample of the flattened version of z, can modify z(n) . . . z(n + L) (due to the coloring filter H -1 m ), where L + 1 is the length of the impulse response of H -1 m , denoted by h -1 m , possibly infinite 2 . Denoting by e z ′ and z ′ the modified versions of e z and z,

∀ k ∈ [0, L], z ′ (n + k) = Round L j=0 h -1 m (j)e z ′ (n + k -j) (10)
If this modification of z reduces d(f z , f target ), then we add δw e (i, n) to e z (n) and we therefore modify z(n) . . . z(n + L) into z ′ (n) . . . z ′ (n + L). Otherwise we do not add δw e (i, n) to e z (n) and z remains unchanged.

Fig. 3 presents an example of H -1 m , e z and z in the frequency domain for some frame m of a piano signal, to illustrate the flattening/coloring scheme.

This procedure is repeated several times on the whole sequence. Hence, for q iterations, z = x + Round(w) [START_REF] Widrow | A Study of Rough Amplitude Quantization by Means of Nyquist Sampling Theory[END_REF] with

w = h -1 m * w e (12)
where

w e (n) = q i=1 δ(i, n)δw e (i, n) (13)
where δ(i, n) = 0 or 1 according to the decision of adding δw e (i, n) or not. All δw e (i, n) are independent. Assuming

d(f z , f target ) ez(n) ez(n + 1) ez(n -1) ez(N -1) ez(0) e z Q x z Initialization z(n) z(n -1) z(n + 1) z(0) z(N -1)
m th frame ON only if this will reduce d(f z , f target ) frame-based update Fig. 2. Perceptually controlled sound histogram reshaping. The initialization part is performed once on the whole signal x. Then, after copying x to z and ex to ez, the remaining part is performed with several iterations i ( ) on the whole signals ez and z. The "Q" box is a quantizer. that the δ(i, n) are independent, δ(i, .) can be modeled by a Bernoulli process of parameter p i , and then w e is a white noise of variance:

σ 2 we = (i1...iq )∈{0;1} q Pr Λ(n, i 1 . . . i q ) q k=1 i k σ 2 k ( 14
)
where Λ(n, i 1 . . . i q ) is defined in (40) (see Appendix A for proof). In the case where ∀ i, p i = p and

σ i = σ, σ 2 we = pqσ 2 (15)
Hence, the constraint γ w (m, ν) < γ mask (m, ν) (5) means:

|H -1 m (ν)| 2 σ 2 we < γ mask (m, ν) (16)
which can be implemented as:

|H -1 m (ν)| 2 = γ mask (m, ν) (17) σ we < 1 ( 18
)
where the satisfaction of the constraint ( 18) can be controlled through the choice of q and (σ i ) 1≤i≤q .

Concerning the δ independence assumption, remind that δ(i, n) depends on the effect on d(f z , f target ) of the modification of z(n) . . . z(n + L) caused by δw e (i, n). Since the δw e (i, n) are independent and the dependency of successive samples of z has only an indirect effect on the dependency of d(f z , f target ) modifications -particularly if L is highit seems reasonable to assume the Bernoullicity of δ(i, .). The latter will be checked in the experimental part (IV) and in Appendix C. This is however not a necessary condition to get w e stationary and white. Referring to Eq. ( 13), since each δw e (i, .) is stationary, w e is stationary if each δ(i, .) is stationary. Concerning the whiteness, let:

w (i) e (n) = δ(i, n)δw e (i, n) (19)
The auto-correlation of w

(i)
e is defined by:

Γ w (i) e (k) (20) = E[δ(i, n)δw e (i, n)δ(i, n + k)δw e (i, n + k)] = Pr δ(i, n)δ(i, n + k) = 1 ×E[δw e (i, n)δw e (i, n + k) | δ(i, n)δ(i, n + k) = 1]
Hence, w e is white if for each i and for k = 0:

E[δw e (i, n) | δ(i, n) = 1] = 0 (21) E[δw e (i, n)δw e (i, n + k) | δ(i, n)δ(i, n + k) = 1] = 0 (22)
The latter properties [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF][START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF] will be checked in the experimental part (IV) and in Appendix C.

The algorithm is repeated until

f z is close enough to f target or the estimated σ we (σ we = 1 N N -1 n=0 w e (n) 2 ) reaches 1.
The detailed algorithm is presented in Fig. 4. Running several iterations with small values of σ i instead of few iterations with σ i ≃ 1 allows to control a slow increase of σ we until 1, avoiding reaching 1 if not necessary. Moreover, a modification of a sample e z (n) that was not possible at an iteration may become possible at the next iteration, since the modifications of the other samples that occurred meanwhile changed the histogram. Hence, small variations distributed among all the samples can achieve the same histogram reshaping as larger variations concentrated on fewer samples, with a lower variance of the transformation noise w e .

To decide if a sample modification will make the histogram f z closer to the target histogram f target , we use a distance D (see Fig. 4) that can be different from the distance d. In particular, a distance d based on the L ∞ -norm is unlikely modified by a local modification of f z , although the latter can contribute to make f z closer to f target and thus reduce d(f z , f target ) in the long term. The distance D will be specified in Section III.

z ← x, fz ← fx, Fz ← Fx and ez ← ex Fix ε, ε ′ , M AX IT , ∆ min d , i = 0, (σ ℓ ) 1≤ℓ≤M AX IT < 1 repeat i ← i + 1 for n = 0 → N -1 do z ′ ← z, f z ′ ← fz, F z ′ ← Fz and e z ′ ← ez Generate δwe ∼ N (0, σ 2 i ) e z ′ (n) ← e z ′ (n) + δwe ∆D ← 0 for k = 0 → L do z ′ (n + k) ← Round L j=0 h -1 m (j)e z ′ (n + k -j) if z ′ (n + k) = z(n + k) then Update f z ′ Compute ∆D = D(f z ′ , ftarget) -D(fz, ftarget) end if end for if ∆D ≤ 0 then z ← z ′ , fz ← f z ′ , Fz ← F z ′ and ez ← e z ′ we(n) ← ez(n) -ex(n) end if end for until d(fz, ftarget) < ε ∨ 1 N N-1 n=0 we(n) 2 > 1 -ε ′ ∨ i > M AX IT ∨ ∆D < ∆ min d central block
Each step (i, n) of the algorithm causes either nothing, or simultaneously a decrease of D(f z , f target ) and an increase of σwe . If the convergence of D is equivalent to that of d, this makes the algorithm converge as long as it is possible to modify one sample e z (n) so that D(f z , f target ) decreases. The stop conditions of the loop however make the algorithm stop if the decrease of d during one iteration on the whole signal is too small.

The coefficients of the filter H m are updated at the same rate as the masking threshold of the signal x. Since the amount of updates of z(n) . . . z(n + L) and f z for each n depends on L, H -1 m should preferably be an FIR filter. The shorter its impulse response, the rougher the approximation of the masking threshold by |H -1 m (ν)| 2 (Eq. ( 17)). Consequently, a trade-off must be done between complexity and perceptual accuracy. The global complexity of the algorithm will we calculated in Section III.

Within this general frame, one must specify:

• the psychoacoustic model and its approximation through the frequency response |H -1 m (ν)|, which condition the perceptual control and the complexity of the algorithm;

• the number of iterations q and the standard deviations (σ i ) 1≤i≤q , which control the convergence; • the distances d and D, which define how the similarity between histograms is measured.

The latter will be specified in Section III, while the former two will be set in the experimental section (IV).

D. Perceptual issues

The coefficients of H -1 m are updated every frame, in order to match the psychoacoustic model. The phase of this filter may change, which may cause discontinuities in the output signal at some of the frame transitions. Hence, the transformation noise w may occasionally have audible discontinuities (clicks) that are not taken into account by the psychoacoustic model, since they originate from the model update itself.

To avoid this, at the end of each iteration of the algorithm, we (i) extract the transformation noise w = zx; (ii) lowpass-filter its discontinuities that occur at frame transitions and that are greater than a locally adapted threshold; (iii) add this smoothed version of w to the original signal x, thus yielding a new z. The length of the impulse response of the low-pass filter must be short enough to ensure a limited impact on P z . The details of the smoothing process are specified in Appendix B.

Additionally, as shown in [START_REF] Baras | Controlling the inaudibility and maximizing the robustness in an audio annotation watermarking system[END_REF], fulfilling the inaudibility constraint [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], where γ mask is provided by a frequencymasking model, ensures the inaudibility of the transformation noise, but only in average for a given set of audio signals.

Actually, the performance of the model depend on the signal, so that γ mask should be multiplied by an attenuation factor depending on the audio content.

Consequently, for each specific audio signal, all filters H m are multiplied by the same attenuation factor λ chosen in such a way that a white Gaussian noise of variance 1 filtered by the successive λH -1 m and added to the audio signal is inaudible.

III. IMPLEMENTATIONS OF THE HISTOGRAM RESHAPING

ALGORITHM

Two different goals may be assigned to this general algorithm: changing the global shape of the histogram (as in [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] and [START_REF] Mahé | Perceptually controlled doping for audio source separation[END_REF]) or, assuming the global shape is satisfying, matching accurately a target histogram (as in [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF]). Fig. 5 illustrates these two transformations: in the first case (Fig. 5a), we focus on the global shape of the histogram, while in the second case (Fig. 5b), the histogram has to be locally "chiseled", keeping the global shape unchanged. We call these two tasks Histogram Global Reshaping (HGR) and Histogram Local Reshaping (HLR), respectively. The behavior of the previously described algorithm will be oriented through the choice of the distance used to measure the dissimilarity between two histograms (respectively the Kolmogorov-Smirnov distance and the total variation distance, as will be detailed later), in two specific implementations.

A. Histogram Global Reshaping (HGR)

To change the global shape of the histogram, each step of the algorithm must reduce the difference between F z and F target , respectively the cumulative histogram of z and the target cumulative histogram. The interest of using the cumulative histogram is the underlying integration of the histogram, which smooths the local differences between f z and f target . Hence, sample modifications directed by F z -F target help focusing on the global shape of the histogram.

For this purpose, the convergence will be assessed through the Kolmogorov-Smirnov distance 3 and we define the distance D by:

D f z , f target = Z F z -F target 2 . ( 23
)
Since D is sensitive to any sample modification, it is used to decide to add (or not) δ we (i, n) to e z (n), while the Kolmogorov-Smirnov distance, sensitive to global modifications (L ∞ norm), is used only at the end of each cycle of the algorithm.

For each iteration i and each discrete time n, adding δw e (i, n) ∼ N (0, σ 2 ) to e z (n) modifies z(n) . . . z(n + L) into z ′ (n) . . . z ′ (n + L). Referring to [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF]

, each z(n + k) with 0 ≤ k ≤ L is modified by approximately h -1 (k)δw e (i, n). If F z (z(n + k) -1) < F target (z(n + k) -1) resp. F z (z(n + k) -1) > F target (z(n + k) -1) a negative (resp. positive) value h -1 (k)δw e (i, n) tends to make F z closer to F target .
Let us define the integer interval I k :

I k = [z(n + k), z ′ (n + k)[ if z(n + k) ≤ z ′ (n + k) [z ′ (n + k), z(n + k)[ if z ′ (n + k) < z(n + k) (24)
For each k from 0 to L, D(f z , f target ) varies by: where

∆ k D = I k F (k) z ′ -F target 2 -F (k-1) z ′ -F target 2 , ( 25
)
f z ′ (z ′ (n + k)) ← f z ′ (z ′ (n + k)) + 1 f z ′ (z(n + k)) ← f z ′ (z(n + k)) -1 if z ′ (n + k) < z(n + k) then I k = [z ′ (n + k), z(n + k)[ ∆ k D -= I k |F z ′ -Ftarget| 2 In I k , F z ′ ← F z ′ + 1 else I k = [z(n + k), z ′ (n + k)[ ∆ k D -= I k |F z ′ -Ftarget| 2 In I k , F z ′ ← F z ′ -1 end if ∆D ← ∆D + I k |F z ′ -Ftarget| 2 -∆ k D -
F (k) z ′ is the new cumulative histogram after modifying z(n + k) (with F (-1) z ′ = F z ). If the total variation of D(f z , f target ) resulting from the L + 1 modifications of z is negative, i.e. ∆D = L k=0 ∆ k D ≤ 0, (26)
we carry out these modifications of e z (n) and z(n) . . . z(n + k), otherwise we cancel them (see Fig. 4 with d = d KS and central block detailed in Fig. 6).

Note that a convergence according to

D(f z , f target ) (L 2 norm) means a convergence according to d KS (f z , f target ) (L ∞ norm))
. Hence, reducing D(f z , f target ) contributes, in the long term, to a decrease of d KS (f z , f target ).

B. Histogram Local Reshaping (HLR)

For a fine adjustment of a histogram having the same global shape as f target , samples of z have to move locally from values in excess in the histogram to deficient values. We assess the convergence through the total variation distance 4 (d TV ), which is sensitive to local differences between histograms.

For each i, n, adding δw e (i,

n) ∼ N (0, σ 2 ) to e z (n) mod- ifies z(n) . . . z(n + L) into z ′ (n) . . . z ′ (n + L). Each of these L + 1 modifications causes a variation of d TV (f z , f target ) in {0; ± 1
N }, according to the relative values of f target (z ′ (n+k)) and f z (z ′ (n + k)) on the one hand, of f target (z(n + k)) and f z (z(n + k)) on the other (see Table I). If the global variation of d TV resulting from the L + 1 modifications of z is negative then we carry out these modifications, otherwise we cancel them.

The detailed algorithm is presented in Fig. 4 with D = d = d TV and central block detailed in Fig. 7. This algorithm ensures the decrease of d TV (f z , f target ).

C. Computational complexity

From Fig. 4, the mean algorithmic complexity C after the initialization phase, measured by the number of 4 Let f and g be two histograms as defined in Section II. The total variation distance is defined by: multiplications-accumulations (MAC), is approximated as follows:

d TV (f, g) = 1 2N k∈Z |f (k) -g(k)|. TABLE I POSSIBLE VARIATIONS OF d TV (fz, ftarget). n k = n + k ftarget(z ′ (n k )) ftarget(z ′ (n k )) > fz(z ′ (n k )) ≤ fz(z ′ (n k )) ftarget(z(n k )) ≥ fz(z(n k )) 0 +1/N ftarget(z(n k )) < fz(z(n k )) -1/N 0 if f z ′ (z(n + k)) > ftarget(z(n + k)) then ∆D ← ∆D -1 2N else ∆D ← ∆D + 1 2N end if if f z ′ (z ′ (n + k)) < ftarget(z ′ (n + k)) then ∆D ← ∆D -1 2N else ∆D ← ∆D + 1 2N end if f z ′ (z ′ (n + k)) ← f z ′ (z ′ (n + k)) + 1 f z ′ (z(n + k)) ← f z ′ (z(n + k)) -1
C ≃ N it N (L + 1) L + 1 + Pr z ′ (n) = z(n) C ′ , (27)
where N it and C ′ denote the number of iterations and the mean complexity of the central block, respectively. Referring to Fig. 6 and7,

C ′ = 3E[|z ′ (n) -z(n)||z ′ (n) = z(n)] + 4
for HGR (28) 2 ≤ C ′ ≤ 4 for HLR [START_REF] Paillard | Perceval: perceptual evaluation of the quality of audio signals[END_REF] We empirically found the expectation in Eq. 28 equal to 1.5 to 3, and Pr z ′ (n) = z(n) ≃ 0.5, so that Eq. 27 becomes:

C ≃ N it N (L + 1)(L + κ), (30)
with 2 ≤ κ ≤ 7.5. Although the complexity is linear in N , the multiplicative factor is quadratic in L, so that flattening filters with long impulse responses can slow down the execution.

IV. EXPERIMENTAL RESULTS

The algorithm (with its two variants HGR and HLR) was tested on music (instruments and singing voice) and on speech. In the following experiments, the filters H -1 m approximating the masking threshold [START_REF] Roberts | Digital signal processing[END_REF] are FIR filters of orders L. The coefficients of H -1 m are those of the auto-regressive model of the inverse of the masking threshold. Indeed, considering 1/γ mask (m, ν) as a power spectral density, one can compute the corresponding auto-correlation coefficients by inverse DFT. Then we derive from these coefficients an auto-regressive model of order L corresponding to a transfer function of form σ m /(1 [START_REF] Roberts | Digital signal processing[END_REF], Chap. 11). The corresponding frequency response is a smooth version of 1/γ mask (m, ν). Hence, the frequency response of

+ L i=1 a (m) i z -i ) (see
H -1 m (z) = (1 + L i=1 a (m) i z -i )/σ m approximates the masking threshold.
This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.

The final version of record is available at http://dx.doi.org/10.1109/TASLP.2018.2836143

Copyright (c) 2018 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

The choice of L depends on the sampling frequency and the shape of the masking threshold. In Subsection IV-A, where music signals sampled at 32kHz are processed, we chose L = 50. In Subsection IV-B, where narrow-band speech is processed, the masking threshold is smoother than for music at higher sampling frequency, so that L = 10 is sufficient (reminding that increasing L increases the complexity).

To set the attenuation factor λ applied to all H -1 m , mentioned in Subsection II-D, the masking was assessed through:

• the Objective Difference Grade (ODG) for music, predicted by PQevalAudio 5 according to PEAQ [START_REF]Recommendation BS.1387: Method for objective measurement of perceived audio quality[END_REF], which provides values between -4 and 0. A value greater than -1 indicates an inaudible distortion. We fixed λ for the ODG to be -0.5 ± 0.1. • the Mean Opinion Score (MOS) for speech, estimated by PESQ [START_REF]Perceptual evaluation of speech quality (PESQ), an objective method for end-to-end speech quality assessment of narrow-band telephone networks and speech codecs[END_REF], which provides scores between -0.5 and 4.5.

An estimated MOS greater than 4 indicates an inaudible distortion. We fixed λ for the MOS to be 4.1 ± 0.1. We present one application for each algorithm. Video and audio demonstrations, as well as the Scilab codes of the algorithms, are available at http://up5.fr/SoundHistograms

A. Histogram Global Reshaping for "Sparsification"

We ran the HGR algorithm on various musical tracks of the QUASI database [START_REF] Vincent | The signal separation evaluation campaign (2007-2010): Achievements and remaining challenges[END_REF], [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF] sampled at 32 kHz. The masking threshold was updated every 8ms according to the MPEG-1 first model [START_REF]Coding of moving pictures and associated audio for digital storage media at up to about 1,5 Mbit/s -Part 3: Audio[END_REF], computed on 50%-overlapping frames of length 512 (16ms), apodized by a Hann window. We set for σ i a constant value of 0.25.

Here, we present the results for an 11s sequence of bass and a 10 s sequence of piano. These sequences have generalized Gaussian (GG) distributions 6 of shape parameters 1.6 and 2.1, respectively (estimation by the moment method [START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF]). We propose to divide their shape parameters by 2 by running the Histogram Global Reshaping algorithm (Fig. 6) with target histograms f target matching GG distributions of shape parameters 0.8 and 1.05, respectively. Fig. 8 presents the convergence of the algorithm in terms of Kolmogorov-Smirnov distance, in parallel with the estimation of σ we , for both the bass sequence (8a) and the piano sequence (8b). For the bass, the algorithm reaches the bound of the inaudibility constraint (σ we < 1) just at the convergence point, while in the case of the piano, the bound of the inaudibility constraint is reached several iterations before reaching the target distribution. In each iteration, each second of signal needs 90 to 120 s to be processed, with our Scilab implementation on a standard computer.

The assumption that each δ(i, .) can be modeled by a 5 http://www-mmsp.ece.mcgill.ca/Documents/Software/ Packages/AFsp/PQevalAudio.html Note that PQevalAudio is calibrated to signals sampled at 48 kHz, while our signals are sampled at 32 kHz. Consequently, we up-sampled them to 48 kHz before measurement 6 The probability density function of a generalized Gaussian distribution is of the form: f (x) = β/ 2αΓ(1/β) exp -(|x -µ|/α) β , where µ, α, and β denote the mean, scale and shape parameters, respectively. β = 2 for a Gaussian distribution. Γ denotes the gamma function, defined for each z ∈ C \ Z by: Γ Bernoulli process was validated by statistical tests detailed in Appendix C. Fig. 9 compares the cumulative histograms F x , F z and F target at the stop-point of the algorithm (end of the 10 th iteration) for bass and piano (9a and 9b, respectively). For the bass (Fig. 9a), the global shape of f z matches that of f target , which leads to very close cumulative histograms. The histogram f z has a final estimated shape parameter of 1.2 7 . As indicated by Fig. 9b, the final histogram for piano is intermediate between the original one and the target one. The transformed signal z has an estimated shape parameter of 1.7.

(z) = ∞ 0 t z-1 e -t dt
The inaudibility of the transformation was assessed through the Objective Difference Grade (ODG), predicted by PQe-valAudio. We obtained an ODG of -0.9 for the bass, and -1.0 for the piano (inaudible distortions).

To compare the performance of this algorithm to that of [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], we plotted d KS (f z , f target ) as a function of the ODG for both algorithms, for the previous bass and piano signals. As indicated by Fig 10, the HGR algorithm reaches a much better trade-off between audio quality and histogram reshaping.

A possible application of this sparsification is source separation in time domain, since blind source separation algorithms (BSS) are known to reach better performance with sparse ODG bass [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] bass (HGR) piano [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] piano (HGR) Fig. 10. Kolmogorov-Smirnov distance to the target histogram, d KS (fz, ftarget), as a function of the ODG, for the bass and piano signals, both processed by the algorithm [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF] and by the HGR algorithm. For [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF], each point corresponds to a given value of wmax. For HGR, each point corresponds to the result of an iteration. signalsi.e. with lower shape parameters in the case of generalized Gaussian distribution -and to behave badly if the distribution is close to a Gaussian (see for example [START_REF] Comon | Handbook of Blind Source Separation[END_REF], [START_REF] Vincent | From blind to guided audio source separation: How models and side information can improve the separation of sound[END_REF]). Our proposal is applicable to contexts where the original sources are available before mixing, like informed source separation (ISS8 ) [START_REF] Liutkus | Informed source separation: a comparative study[END_REF].

To illustrate this proposal, we compared the performance of a classic BSS algorithm on stereo mixtures of pairs of instruments (among singing voice, piano, acoustic guitar, electric guitar, and keyboards), sparsified according to our algorithm, to the performance obtained without sparsification. The experiment was run as illustrated by Fig. 11. We sparsified each source signal according to the HGR algorithm, with a target shape parameter equal to half of that of the original distribution of the signal. We chose as mixture matrix:

A = 1/3 2/3 2/3 1/3 (31)
We used the Scilab implementation 9 of the FastICA algorithm [START_REF] Hyvärinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] to separate the sources. Since the Independent Component Analysis (ICA) of two mixtures of two sources is the easiest task of BSS, the sources are often well separated without any pre-processing. Consequently we focused the experiment on mixtures badly separated by FastICA. We ran this experiment on four duos from the QUASI database. For each of them, we evaluated:

• the sparsity of the original and sparsified signals, through the shape parameter β of the distribution; • the quality of the sparsified mix, compared to the original one, through the ODG (predicted by PQevalAudio); • the performance of the source separation, through the Source to Distortion Ratio (SDR), the Source to Interference Ratio (SIR) and the Source to Artifact Ratio (SAR) as defined by [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. The results are summarized in Table II and Fig 12 10 and can be heard on the aforementioned web page. For mixtures A, B and D, the sparsification clearly enhances the performance of source separation. Note that the SIR and SDR increase even if the shape parameter β is not much reduced (see Mixtures C and D). For mixtures that are already fairly separated without sparsification, like Mixture C, the sparsification leads to a lower enhancement of separation. In Experiment B, since the shape parameter of the guitar signal is around 3, sparsifying this signal could make its distribution close to a Gaussian one. Consequently, we let it unsparsified.

The ODG of the mixes (Table II) and the SxRs (Fig. 12) measured separately the perceptual impairments of the mixes due to sparsification and the separation enhancement due to sparsification, respectively. To measure the global gain due to sparsification, encompassing both effects, we computed the perceptual scores provided by the PEASS toolkit [START_REF] Emiya | Subjective and objective quality assessment of audio source separation[END_REF], using the original signals (x 1 and x 2 in Fig. 11) as reference signals for both couples of separated signals: (x 1 , x2 ) and (ẑ 1 , ẑ2 ). PEASS provides four perceptual scores (PS): overall (OPS), target-related (TPS), interference-related (IPS), and artifactsrelated (APS). These scores are summarized in Fig. 13 for four considered mixes. The same observations as for Fig. 12 can be made. Hence, when taking as reference the unsparsified signals, sparsification-mixing-unmixing process leads to better separation quality than mixing-unmixing.

B. Low-pass filtering of the histogram

We will illustrate the Histogram Local Reshaping algorithm in the case of low-pass filtering of the histogram for the application of the quantization theorem.

Let x be a discrete-valued signal with (discrete) probability density function (PDF) f x . The characteristic function of x is defined by DFT -1 (f x ). According to the sub-quantization theorem [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF], adapted from the quantization theorem [START_REF] Widrow | A Study of Rough Amplitude Quantization by Means of Nyquist Sampling Theory[END_REF] to discrete-valued signals, if the characteristic function of x is equal to zero for frequencies |ν| > 1 2K in [-1 2 ; 1 2 ], then the probability density function of x can be derived from that of the signal x Q resulting from the sub-quantization of x with a factor K.

The original PDF is recovered through filtering the PDF of x Q by a filter of frequency response G(ν):

G(ν) = 1 R(ν) if |ν| < 1 2K 0 otherwise (32)
where:

R(ν) =      K if ν ∈ Z sin(πKν) sin(πν) exp(jπν) ν / ∈ Z and K even sin(πKν) sin(πν) ν / ∈ Z and K odd (33)
The same results hold for histogram instead of PDF.

Hence, to fulfill the condition of the sub-quantization theorem, the histogram must be low-passed, with a cutoff frequency 1/2K. Since this may not change the global shape of the histogram, the HLR algorithm is appropriate for this purpose. The HLR algorithm was run on speech signals from the TIMIT database [START_REF] Garofolo | Timit acoustic-phonetic continuous speech corpus[END_REF], sampled at 8kHz. The masking threshold was approximated by the power spectral density of the signal, minus a given offset [START_REF] Paillard | Perceval: perceptual evaluation of the quality of audio signals[END_REF], computed on 50%-overlapping frames of length 256 (32ms), apodized by a Hamming window. We set K = 16, so that the target histogram is the result of the low-pass filtering of f x with a normalized cut-off frequency 1/32.

We present the results for a 3s sentence pronounced by a female and a male speakers. Fig. 14 illustrates the convergence of the algorithm in terms of total variation distance, in parallel with the estimation of σ we , for constant σ i = 0.25, M AX IT = ∞ and ∆ min d = 0 (see Fig. 7). In each iteration, each second of signal needs 3s to be processed, with our Scilab implementation on a standard computer. After a strong decrease during the first iteration, the distance to the target histogram decreases more and more slowly, while the estimated σ we reaches 1. The irregular evolution of d TV (f z , f target ) for a large number of iterations is caused by the discontinuities smoothing procedure at the end of iteration, which effect on d TV may not be negligible anymore compared to the low decrease of d TV during the iteration. If we stop the algorithm after 100 iterations, the histogram is really lowpassed, as illustrated by Fig. 15 for the female speaker (a similar result was obtained for the male speaker), whereas z remains perceptually identical to x. The Mean Opinion Score (MOS) of z compared to x is estimated by PESQ [START_REF]Perceptual evaluation of speech quality (PESQ), an objective method for end-to-end speech quality assessment of narrow-band telephone networks and speech codecs[END_REF]. The estimated MOS is 4.4 for the female speaker (4.3 for male speaker), which indicates an inaudible distortion. Comparatively, low-pass the histogram as in [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF] leads to a signal y of which histogram f y is equal to f target , but with an estimated MOS of 3.9 for the female speaker and 3.7 for the male speaker (slightly audible distortion).

Note that δ(i, .) here is far from a Bernoulli process, but is stationary, and the conditions [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF][START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF] for whiteness of w e are fulfilled for most iterations (see Appendix C).

We sub-quantized with a factor K = 16 the original signal x and its transformed versions z and y into x Q , z Q , and y Q , respectively. Then we retrieved each original PDF from that of the quantized signal, using the filter G (32). The results are summarized in Table III. The proposed algorithm leads to a recovery error which is much better than without low-pass filtering the histogram, higher than that of [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF], but with a guaranteed inaudible transformation noise. Hence, HLR is an efficient approach for inaudibly lowpass-filtering sound histograms, which makes possible to at last take advantage from the powerful Widrow quantization theorem for audio signals.

|DF T -1 (f z )| |DF T -1 (f x )| |DF T -1 (f target )|

V. CONCLUSION

We have proposed a sound histogram reshaping method that makes any histogram closer to a target histogram, while controlling the inaudibility of the transformation. The principle is to add iteratively a low-power white noise to a flat-spectrum version of the signal, until the target distribution or the noise audibility is reached. This scheme was applied through two variants, Histogram Global Reshaping (to change the global shape of the histogram) and Histogram Local Reshaping (to locally "chisel" the histogram, keeping its global shape unchanged).

This algorithm for perceptually controlled reshaping of sound histograms opens new perspectives in various audio processing applications where the original signal is available and its distribution preferably fulfills specific requirements: bandlimited characteristic function for histogram restoration [START_REF] Halalchi | Revisiting quantization theorem through audiowatermarking[END_REF]; Gaussianity for non-linear audio system identification [START_REF] Mezghani-Marrakchi | Nonlinear audio systems identification through audio input Gaussianization[END_REF]; sparsity for informed audio source separation [START_REF] Comon | Handbook of Blind Source Separation[END_REF]; specific distribution for optimal audio quantization [START_REF] Gersho | Principles of quantization[END_REF].

Future works will focus on reducing the complexity factor due to the length of the coloring filter impulse response. The solution could be based on a perception-based filter-bank (auditory filters mimicking the auditory system [START_REF] Strahl | Analysis and design of gammatone signal models[END_REF], [START_REF] Necciari | Audlet filter banks: A versatile analysis/synthesis framework using auditory frequency scales[END_REF]) and white noise addition in each frequency channel, which would avoid the memory effect responsible for complexity.

Audio signals have a malleability bounded by the properties of human audition, which was extensively exploited by audiocoding [START_REF] Spanias | Audio Signal Processing and Coding[END_REF], [START_REF] Moreau | Outils pour la compression des signaux -Applications aux signaux audio[END_REF] and watermarking [START_REF] Hua | Twenty years of digital audio watermarking -a comprehensive review[END_REF]. This perceptually controlled reshaping of sounds histogram is another way of exploring this malleability. checking whether the T i (j) are independent and match a geometric law of parameter p i .

We assessed the independence as follows [START_REF] Boudec | Performance Evaluation of Computer and Communication Systems[END_REF]. First, a Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test checks the stationarity 11 . Then, if the process can be considered as stationary, the independence can be assessed through the autocorrelation coefficients ρ(k): if the latter are inside the 95 % confidence interval around 0, the hypothesis of independence can be validated.

Assuming that T i follows a geometric law, the estimation of the parameter p i according to the maximum likelihood is given by pi = 1/(1 + T i ), where T i is the empirical mean of T i . We assessed the goodness of fit of the geometric law of parameter pi graphically and through the chi-squared test for fit. For each i, we computed the test statistics for 1000 draws of 1000 individuals of T i and observed the percentage of statistics χ 2 above the 5 % critical value.

The results are summarized in Table IV and Fig. 17, for the experiments of subsection IV-A. The assumption of independence for T i is valid for both bass and piano. As illustrated by Fig. 17, the histogram of T i matches well a geometric law, except in the tail of the distribution for the piano (for probabilities lower than 10 -3 , emphasized here by the log-scale), which leads to poor results in the chi-squared tests.

For the experiments of subsection IV-B, the stationarity of

| δ(i, n) = 1] | δ(i, n)δ(i, n + k) = 1]
out of the 95 % confidence interval around 0 female 9 31 0 male 1 45 0 article.

Fig. 1 .

 1 Fig. 1. Basic principle of the perceptually controlled sound histogram transformation.

Fig. 3 .

 3 Fig. 3. For a frame m of a piano signal sampled at 32 kHz, magnitude frequency response |H -1 m (ν)| of the coloring filter and magnitude spectra (in decibels) |Ez(ν)| and |Z(ν)| of signals ez and z. The spectra are computed by Discrete Fourier Transform on 512 samples multiplied by a Hann window. The filter H -1 m corresponds to the 256 central samples of this frame.

Fig. 4 .

 4 Fig. 4. Perceptually controlled histogram transformation. D is a distance that can be different from d. The central block depends on the application case: HGR (see Fig. 6) or HLR (see Fig 7).

Fig. 5 .

 5 Fig. 5. Examples of Histogram Global Reshaping (HGR, a) and Histogram Local Reshaping (HLR, b). The total variation distance (d TV ) is more sensitive to local differences than the Kolmogorov-Smirnov distance (d KS ).

Fig. 6 .

 6 Fig. 6. Central block of perceptually controlled histogram transformation (see Fig 4) for Histogram Global Reshaping.

Fig. 7 .

 7 Fig. 7. Central block of perceptually controlled histogram transformation (see Fig 4) for Histogram Local Reshaping.

  The signals have integer values in [-2 15 ; 2 15 -1].

Fig. 8 .

 8 Fig.8. Histogram Global Reshaping for an 11s bass sequence (a) and a 10s piano sequence (b). Initial shape parameters 1.6 and 2.1 respectively. Target shape parameters 0.8 and 1.05 respectively. Evolution of the estimated standard deviation of we (σw e ) and of the Kolmogorov-Smirnov distance to the target distribution (d KS (fz, ftarget)). The dotted lines correspond to the prolongation of the algorithm beyond its normal stop-point (when σwe reaches 1).

Fig. 9 .

 9 Fig. 9. Target cumulative histogram Ftarget, difference between the cumulative histogram Fx and Ftarget, difference between cumulative histogram "sparsified" z and Ftarget, for an 11s bass sequence (a) and a 10s piano sequence (b).

Fig. 12 .

 12 Fig. 12. Performance of source separation with and without sparsification for 4 mixtures A, B, C, and D. Vertical bars represent confidence intervals (when applicable). In Mixture B, the guitar signal was not sparsified. SIR = Source to Interference Ratio; SDR = Source to Distortion SAR = Source to Artifact Ratio.

Fig. 13 .

 13 Fig. 13. Global perceptual impact of sparsification on source separation for 4 mixtures A, B, C, and D. OPS = Overall Perceptual Score, TPS = Targetrelated Perceptual Score, IPS = Interference-related Perceptual Score, APS = Artifacts-related Perceptual Score. Vertical bars represent confidence intervals (when applicable).

Fig. 14 .

 14 Fig. 14. Histogram Local Reshaping for a 3s speech sequence from a female (a) and a male (b) speakers, with target histogram defined by low-pass filtering of the original histogram, with cut-off frequency 1/32. Evolution of the Total Variation distance d TV (fz, ftarget) and of the estimated variance of we (σw e ). According to the constraint (18), the algorithm stops when σwe reaches 1.

Fig. 15 .

 15 Fig. 15. After 100 iterations of the HLR algorithm, characteristic functions of the original signal (|DFT -1 (fx)|) and of the signal with low-pass-filtered histogram (|DFT -1 (fz)|), compared to the target characteristic function |DFT -1 (ftarget)|.

Fig. 17 .

 17 Fig. 17. How T i fits a geometric law of parameter pi , for the first and last iterations, in the bass (a) and piano (b) experiments.

TABLE II RESULTS

 II OF SPARSIFICATION OF THE MIXTURES.

	Mixture	instruments	shape parameter	ODG mix
	A	voice piano		1.8 → 1.6 2.1 → 1.7		-0.8
	B	guitar keyboards	3 1.8 → 1.3		-0.4
	C	voice 1 voice 2		1.1 → 1.0 1.0 → 0.9		-0.6
	D	guitar solo guitar acoustic	1.5 → 1.5 1.4 → 1.1		-0.8
	x 1					x1
	x 2			A	BSS	x2
		sp a rs ifi c a ti o n	z 1 z 2	A	BSS	ẑ1 ẑ2
	Fig. 11. Blind source separation (BSS) of a stereo mixture (mixture matrix
	A), without (top) and with (bottom) sparsification of the source signals.
	SxR (dB)
				Mixture

TABLE III PDF

 III RECOVERY ERROR, AS THE TOTAL VARIATION DISTANCE BETWEEN THE ORIGINAL AND THE RECOVERED PDF.

		Recovery error
	Signal	Female speaker	Male speaker
	Original x	1.3 × 10 -1	1.3 × 10 -1
	HLR-modified z	1.9 × 10 -2	2.7 × 10 -2
	Transformed as in [10] y	2.4 × 10 -3	2.2 × 10 -3

TABLE V NUMBER

 V OF ITERATIONS (AMONG 100) HAVING EMPIRICAL E[δwe(i, n) | δ(i, n) = 1] AND MORE THAN 5 OR 10 (AMONG 100) EMPIRICAL E[δwe(i, n)δwe(i, n + k) | δ(i, n)δ(i, n + k) = 1] OUT OF THE 95 % CONFIDENCE INTERVAL AROUND ZERO.

	Number of iterations having
	more than 5	more than 10

E[δwe(i, n) E[δwe(i, n)δwe(i, n + k)

Note that the flattening filter H actually flattens the masking threshold of the signal, not the signal spectrum itself, as will be seen later (see Eq. (17))

Let f and g be two histograms and F and G the respective corresponding cumulative histograms, as defined in Section II. The Kolmogorov-Smirnov distance is defined by:d KS (f, g) = 1 N sup k∈Z |F (k) -G(k)|.

Since the final Kolmogorov-Smirnov distance,

× 10 -3 , is very weak, the shape parameter estimation method may not be accurate

The typical application case of ISS is the following: musical sources are recorded separately in a studio on different tracks; the tracks are watermarked before mixing; then, any final user can separate the sources from the mixture by using the information conveyed by the watermark, assuming the latter has been specifically designed for this purpose.

http://research.ics.aalto.fi/ica/fastica/

In some cases, the separation performed by FastICA is variable. Hence, we ran the separation algorithm 20 times for each mixture and displayed the means and the confidence intervals (if applicable) of the SxR values.

We used the "grocer" toolbox of Scilab: https://atoms.scilab.org/toolboxes/grocer

δ(i, n)δw e (i, n), [START_REF] Hua | Twenty years of digital audio watermarking -a comprehensive review[END_REF] where ∀ i, δw e (i, n) ∼ N (0, σ 2 i ) and all δw e (i, n) are independent, δ(i, .) is a Bernoulli process of parameter p i .

A. Case ∀ i, p i = p and σ i = σ

The probability density function of w e (n) is given by:

where

Hence, w e (n) has zero mean and variance:

B. General case

For (i 1 . . . i q ) ∈ {0; 1} q , let us define the event:

probability density function of w e (n) is given by: where

Hence, w e (n) has zero mean and variance:

APPENDIX B SMOOTHING THE DISCONTINUITIES OF THE TRANSFORMATION NOISE w

The smoothing process is illustrated by Fig. 16. The general principle is to decompose w into two components, a continuous part w cont and a piece-wise-constant signal w disc (with jumps), and then to low-pass filter w disc .

The transformation noise w is obtained by subtracting z to x and then differentiated. The difference signal dw is thresholded only at block transitions and then integrated, yielding w cont . Subtracting the resulting signal to w yields a piece-wiseconstant signal w disc , containing the discontinuous component of w. Low-pass filtering w disc and adding the result to w cont provides a signal w, corresponding to w with smoothed blocktransitions. Adding it to x provides the new z.

The thresholding is performed as follows. Let dw(n) = w(n)w(n -1). We define dw med (n) as the median value of dw on the interval n ± M 1 , where M 1 is a small fixed integer. Let

The low-pass filter has a triangular impulse response and is non-causal in such a way as to introduce no delay.

In the experiments described in Section IV, M 1 = 2, M 2 corresponds to 4ms and the low-pass filter has an impulse response of length 1ms.

APPENDIX C CHECKING WHETHER δ(i, .) IS A BERNOULLI PROCESS

For each iteration i, we want to check whether the decision δ(i, n) of adding δ we (n) to e z (n) actually follows a Bernoulli process of parameter p i . Denoting by T i the random variable "number of zeros between two ones in δ(i, .)", this means T i is rejected for i = 1 (p 1 decreases during the iteration), and accepted at 1 to 10% risk for i ≥ 2. For both speakers and all iterations, the majority of the auto-correlation coefficients are out of the 95% confidence interval around 0. Hence, the T i (j) cannot be considered as independent. Moreover, as illustrated by Fig. 18, the histogram of T i does not match a geometric law: there are too much zeros and high values, i.e. δ(i, .) contains too long sequences of zeros or ones. Consequently, δ(i, .) here cannot be considered as a Bernoulli process.

For each of the 100 iterations, we computed the empirical values of E[δw e (i, n) | δ(i, n) = 1] and E[δw e (i, n)δw e (i, n + k) | δ(i, n)δ(i, n + k) = 1] for k = 1 to 100 (see conditions [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF][START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF]), and we counted how many of them are out of the 95 % confidence interval around zero. The results reported in Table V indicate that the conditions [START_REF] Liutkus | Gaussian processes for underdetermined source separation[END_REF][START_REF] Varanasi | Parametric generalized Gaussian density estimation[END_REF] are fulfilled or almost fulfilled for most iterations.

ACKNOWLEDGMENT

The authors thank the reviewers for their comments and suggestions, which were helpful to clarify and enrich this