Bruno Galerne
email: bruno.galerne@parisdescartes.fr

Ares Lagae
email: ares.lagae@cs.kuleuven.be

Sylvain Lefebvre
email: sylvain.lefebvre@inria.fr

Gabor Noise by Example

Keywords: procedural noise, Gaussian texture, Gaussian random field, power spectrum estimation, non-negative basis pursuit denoising, decorrelated color space

(1) Gauss. tex.

(2) Proc. noise Figure 1: We present Gabor noise by example, a method to estimate the parameters of bandwidth-quantized Gabor noise, a procedural noise function that can generate noise with an arbitrary power spectrum, from exemplar Gaussian textures, a class of textures that is completely characterized by their power spectrum. (row 1) Gaussian texture. (row 2) Procedural noise. (insets) Estimated power spectrum.

Introduction

Since its introduction by [START_REF] Perlin | An image synthesizer[END_REF] almost three decades ago, procedural noise has become a fundamental tool in Computer Graphics. However, designing noise patterns is still hard. Implementing procedural shaders and tweaking their parameters in order to achieve a desired visual effect requires significant effort. We believe this difficulty stems, in part, from two reasons: (i) the expressiveness of current procedural noise functions is rather limited in practice (see Sec. 2.1), and (ii) selecting the appropriate noise parameters to obtain a desired noise pattern or texture is difficult, even more so when the noise function has more parameters. The evident way to address this problem is estimating the parameters of a procedural noise function from an exemplar, i.e., noise by example. However, it is not immediately clear which textures are appropriate exemplars for noise by example, since noise obviously cannot reproduce any texture. We therefore introduce a new texture class, called Gaussian textures, that is defined, similarly to noise [Lagae et al. 2010a], as what is known in statistics as a Gaussian random field [START_REF] Papoulis | Probability, Random Variables and Stochastic Processes[END_REF]. Both Gaussian textures and noise are completely characterized by their power spectrum. Ideally, a procedural noise function should thus be able to reproduce any Gaussian texture. However, current procedural noise functions have not demonstrated a diversity in noise patterns as rich as the Gaussian textures shown in Fig. 1.

In this paper, we present Gabor noise by example, a generalization of Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]] that (i) can generate a wide variety of Gaussian textures, and (ii) can estimate the noise parameters from a Gaussian exemplar. Our approach can reproduce the rich diversity of Gaussian textures shown in Fig. 1, is mathematically well-founded, and preserves the traditional advantages of procedural noise, including a compact representation and a fast on-the-fly evaluation. More specifically, the contributions of our work are:

• bandwidth-quantized Gabor noise, a generalization of Gabor noise to arbitrary power spectra that enables robust parameter estimation and efficient procedural evaluation (Sec. 4);

• a robust parameter estimation technique for bandwidthquantized Gabor noise, that automatically decomposes the noisy power spectrum estimate of an exemplar into a sparse sum of Gaussians using non-negative basis pursuit denoising [START_REF] Kim | An interior-point method for largescale l1-regularized least squares[END_REF]] (Sec. 5); and

• an efficient procedural evaluation scheme for bandwidthquantized Gabor noise, that uses multi-grid evaluation and importance sampling of the kernel parameters (Sec. 6).

Additionally, we present a maximally independent color space, a new color space for independent channel synthesis that makes the color channels of a texture maximally independent using approximate joint diagonalization [Hyvärinen et al. 2001, 11.3] (Sec. 7).

We also demonstrate an interactive editor for bandwidth-quantized Gabor noise, inspired by the Wold decomposition [START_REF] Francos | A unified texture model based on a 2-D Wold-like decomposition[END_REF]] (Sec. 8).

2 Related Work

Procedural Noise Functions

Several procedural noise functions have been proposed since the seminal work of [START_REF] Perlin | An image synthesizer[END_REF]. We refer to the survey of Lagae et al. [2010a] for a recent overview. However, the expressiveness of these noise functions is rather limited in practice. Perlin noise [START_REF] Perlin | An image synthesizer[END_REF]], sparse convolution noise as proposed by [START_REF] Lewis | Algorithms for solid noise synthesis[END_REF], and wavelet noise [START_REF] Cook | Wavelet noise[END_REF] are limited to isotropic power spectra. Anisotropic noise [START_REF] Goldberg | Anisotropic noise[END_REF]] is limited to power spectra that correspond to a coarse tiling of the frequency domain into oriented subbands. Although more complex power spectra could be obtained by using a much finer tiling, this would probably be impractical, since the method requires precomputation and storage of a noise texture for each oriented subband. Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF][START_REF] Lagae | Filtering solid Gabor noise[END_REF] has only been demonstrated on power spectra that consist of a small number of Gaussians, where the parameters of the noise corresponding to each Gaussian are set manually. We extend Gabor noise to arbitrary power spectra, where the parameters are estimated automatically from exemplars.

Example-Based Texture Synthesis

Many non-procedural example-based texture synthesis methods have been proposed. We refer to the survey of [START_REF] Wei | State of the art in example-based texture synthesis[END_REF] for a recent overview. However, these methods do not produce procedural textures. Nevertheless, for the class of Gaussian textures, our method is competitive with methods such as pyramid-based texture analysis/synthesis [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF] and parallel controllable texture synthesis [START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF] (see Sec. 9.2).

Procedural Texture Synthesis by Example

Ghazanfarpour [START_REF] Ghazanfarpour | Generation of 3d texture using multiple 2d models analysis[END_REF] presented a method for the automatic generation of 3D textures using spectral analysis of two or three 2D slices of a 3D texture, based on an earlier method [START_REF] Ghazanfarpour | Spectral analysis for automatic 3-d texture generation[END_REF] that only takes into account a single 2D texture. They generate a discrete 3D model from the 2D slices, select a low number of high-magnitude peaks in the Fourier transform of the model, and approximate it by a cosine-like summation. However, as stated by the authors, their method is not well adapted to high-bandwidth or noisy textures such as sand, and requires the user to manually set several parameters. [START_REF] Dischler | A procedural description of geometric textures by spectral and spatial analysis of profiles[END_REF] also presented a method for automatically generating a procedural description of a geometric texture from a 1D exemplar profile. They iteratively apply spectral and histogram analysis. Their method is also applicable to isotropic color textures. However, it cannot handle anisotropic textures. Lagae et al. [2010b] presented a method for procedural isotropic stochastic textures by example. They use a multi-octave wavelet noise model where the weights are computed from an exemplar, independent channel synthesis in a PCA color space, and histogram matching. [START_REF] Xue | Stone weathering in a photograph[END_REF] presented a similar method for simulating rough appearance for stone weathering in a photograph. However, both methods cannot handle anisotropic textures. Gilet et al. [2010] presented a method to generate procedural descriptions of anisotropic noisy textures by example. They use a Gabor noise model with a small number of noises and several cosines, independent channel synthesis in a PCA color space, and histogram matching. They decompose the 2D spectral domain into ellipses centered at the origin, by discretizing the power spectrum into a low number of values and fitting one or two ellipses to each binary area resulting from discretization, and associate a Gabor noise with every ellipse. Gilet and Dischler [2010] later used this method in an image-based approach for stochastic volumetric and procedural details. However, their method requires significant manual intervention and trial and error for each exemplar (see Sec. 9.2), while our approach is automatic. [START_REF] Jeschke | Estimating color and texture parameters for vector graphics[END_REF] presented a diffusion curve coloring algorithm that includes texture details based on spatially varying Gabor noise. However, their method records only a single dominant frequency, which captures only a rough impression of the texture. [START_REF] Yoon | Editing noise[END_REF]2008] presented a method to edit noise values to allow user constraints to be applied to a noise function. Their method is designed to edit individual noise values, while our approach is designed to edit the overall appearance of the noise.

Background and Motivation

Noise is defined as what is known in statistics as a Gaussian random field, which is completely characterized by its power spectrum [Lagae et al. 2010a;[START_REF] Papoulis | Probability, Random Variables and Stochastic Processes[END_REF]. The problem of procedural noise by example can thus be solved by choosing a procedural noise function and setting the parameters of the noise function such that the power spectrum of the noise is similar to that of the exemplar. This raises three important questions: (i) What are the appropriate exemplars? (ii) How is the power spectrum of the exemplar obtained? (iii) What is the appropriate procedural noise function?

Gaussian Textures

We use Gaussian textures as exemplars, since these are defined, similarly to noise, as what is known in statistics as a Gaussian random field [START_REF] Papoulis | Probability, Random Variables and Stochastic Processes[END_REF]. Both Gaussian textures and noise are completely characterized by their power spectrum. The class of Gaussian textures roughly corresponds to the class of stochastic textures. [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF] recently presented methods that can be used to generate the Gaussian version of a texture. These methods explicitly destroy the information contained in the phase spectrum, only preserving the information contained in the power spectrum. We use the methods of Galerne et al. to illustrate the class of Gaussian textures, to determine how Gaussian an exemplar is, and to provide a ground truth for our approach.

Power Spectrum Estimation

We estimate the power spectrum of the exemplar using the periodogram, since this requires only a single exemplar. The periodogram, defined as the magnitude squared of the Fourier transform, is a simple estimator for the power spectrum [Press et al. 2002, 13.4]. Unfortunately, the periodogram is inherently noisy, since it is actually a white noise process with the power spectrum as mean [Papoulis and Pillai 2002, 12.2]. Note that a less noisy estimate could be obtained by averaging multiple peridograms, but that would require several exemplars. The periodogram is based on the Fourier transform. However, the Fourier transform of a nonperiodic textures suffers from a boundary problem, which causes a cross-shaped artifact in the estimated power spectrum. We avoid this problem by using the Fourier transform of the periodic component instead [START_REF] Moisan | Periodic plus smooth image decomposition[END_REF][START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]].

Gabor Noise

We use Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF], since it has a power spectrum that can be used as a basis to approximate arbitrary power spectra. Two-dimensional anisotropic random-phase Gabor noise [START_REF] Lagae | Filtering solid Gabor noise[END_REF] is defined as

n (x; K, a, ω) = 1 √ λ i g (x -xi; K, a, ω, φi) , (1)
where K, a and ω are the amplitude, bandwidth and frequency of the noise, g is the phase-augmented Gabor kernel, the random positions {xi} are distributed according to a Poisson process with mean λ, and the random phases {φi} are distributed according to a uniform distribution on the interval [0, 2π). The phase-augmented Gabor kernel is defined as

g (x; K, a, ω, φ) = Ke -πa 2 |x| 2 cos (2πx • ω + φ) , (2)
where K, a, ω and φ are the amplitude, bandwidth, frequency and phase of the kernel. The variance of the noise is K 2 /4a 2 . The power spectrum of the noise is

Sn (ξ; K, a, ω) = K 2 8a 2 G ξ; ±ω, a 2 √ π , (3)
where G (x; µ, σ) is the 2D normalized Gaussian function with mean µ and standard deviation σ, and G (x; ±µ, σ) is a shorthand notation for G (x; µ, σ) + G (x; -µ, σ). The power spectrum of Gabor noise is thus a symmetric pair of Gaussians in the frequency domain with magnitude K, frequency ±ω and bandwidth a.

Gabor Noise by Example

We solve the problem of Gabor noise by example by decomposing the noisy power spectrum estimate of the exemplar into a sum of Gaussians, and evaluating the corresponding sum of Gabor noises. However, a straightforward implementation of this idea is problematic. A first problem is that evaluating a potentially large number of Gabor noises is inefficient (see Sec. 6.1). A second problem is that allowing arbitrary values for the bandwidth parameter of the Gabor kernels is problematic, both for the decomposition of the power spectrum (see Sec. 5.2), as well as for the evaluation (see Sec. 6.1). To overcome these problems we introduce bandwidthquantized Gabor noise (Sec. 4). We solve the first problem by grouping the Gabor kernels of the Gabor noises into a single Gabor noise. We solve the second problem by quantizing the bandwidths of the Gabor kernels. This, in turn, allows us to introduce a robust parameter estimation technique (Sec. 5) and an efficient procedural evaluation scheme (Sec. 6) for bandwidth-quantized Gabor noise.

Bandwidth-Quantized Gabor Noise

In this section, we introduce bandwidth-quantized Gabor noise, a generalization of Gabor noise to arbitrary power spectra, that quantizes the bandwidths of the Gabor kernels, in order to enable robust parameter estimation and efficient procedural evaluation.

Arbitrary Power Spectra

Our key insight for generalizing Gabor noise to arbitrary power spectra is that an arbitrary symmetric power spectrum S (ξ) can be decomposed into a sum of G symmetric pairs of Gaussians S (ξ) = g Sn g (ξ; Kg, ag, ωg) of the form of Eqn. 3. This implies that the corresponding noise n (x) can be seen as the sum of the G corresponding independent Gabor noises n (x) = g ng (x; Kg, ag, ωg) of the form of Eqn. 1. Note that we use the subscript g to index the noise functions corresponding to the Gaussians. This decomposition is similar in spirit to how an arbitrary function can be approximated to any prescribed tolerance by a sum of Gaussians [START_REF] Ferreira | A comment on the approximation of signals by gaussian functions[END_REF]]. However, in our case, all involved functions are symmetric, since they correspond to power spectra of real signals, and the weights of the Gaussians are positive, since power spectra are positive by definition. This is why we use the random-phase variant of Gabor noise [START_REF] Lagae | Filtering solid Gabor noise[END_REF] rather than the original one [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]]: The power spectrum of the original variant has an additional Gaussian centered at zero, which would make this decomposition impossible, since excess power cannot be subtracted away.

Kernel Combination

Evaluating n (x) by evaluating and summing the G corresponding Gabor noises is inefficient, especially when G is relatively large (see Sec. 6.1). We therefore combine all the Gabor kernels of the G noises into a single but more complex Gabor noise,

n (x) = 1 √ λ i 1 √ pi g (x -xi; Ki, ai, ωi, φi) , (4)
where λ = g λg and the parameters of each Gabor kernel {(Ki, ai, ωi)} are randomly chosen from {(Kg, ag, ωg)}, the parameters of the G noises, with probability pg = λg/λ.

Bandwidth Quantization

Allowing arbitrary values for the bandwidth ag is problematic, both for the parameter estimation (see Sec. 5.2) as well as for the procedural evaluation (see Sec. 6.1). We therefore quantize the bandwidth to a small set of B bandwidths B = {a b }, and regroup all the Gabor kernels which share the same bandwidth. We use the bandwidth discretization a b = 2 -b π/2 ln 2, motivated by Gabor filter bank design for texture analysis [START_REF] Bovik | Multichannel texture analysis using localized spatial filters[END_REF]], which corresponds to a set of Gaussians with a full width at half maximum (FWHM) of powers of two.

Bandwidth-Quantized Gabor Noise

We define 2D bandwidth-quantized Gabor noise as

n (x) = b∈B 1 √ λ b i 1 √ p b,i g (x -x b,i ; K b,i , a b , ω b,i , φ b,i ) ,
(5) where λ b = g λ b,g and the parameters of each Gabor kernel {(K b,i , ω b,i )} are randomly chosen from {(K b,g , ω b,g )}, the parameters of the G b noises that share the same bandwidth a b , with probability p b,g = λ b,g /λ b . Note that b G b = G, and that we use the subscript b,g to index the noise functions corresponding to the Gaussians in each bandwidth. The power spectrum of the noise is

Sn (ξ) = b∈B G b -1 g=0 K 2 b,g 8a 2 b G ξ; ±ω b,g , a b 2 √ π , (6)
i.e., a sum of Gaussians partitioned into a discrete set of bandwidths. This is illustrated in 2. The variance of the noise is

σ 2 n b g K 2 b,g 4a 2 b
, where σ 2 b,g = K 2 b,g /4a 2 b is the variance of the g-th Gaussian of bandwidth b, and σ 2 b = g σ 2 b,g is the variance of bandwidth b.

Robust Parameter Estimation

this section, we introduce a robust parameter estimation technique for bandwidth-quantized Gabor noise, which automatically decomposes the noisy power spectrum estimate of an exemplar into a sparse sum of Gaussians, ensuring a compact procedural representation.

The Parameter Estimation Problem

The goal of the parameter estimation is to determine the parameters of a bandwidth-quantized Gabor noise such that its power spectrum is close to that of the exemplar. Intuitively, this corresponds to fitting the power spectrum of the exemplar with a sum of Gaussians.

We denote the M × M discrete power spectrum estimate of the exemplar as Sex (ξ0, ξ1) (ξ0, ξ1 ∈ -M/2, . . . , M/2 -1). We discretize the power spectrum of bandwidth-quantized Gabor noise (Eqn. 6) at the same resolution of the exemplar by placing for each bandwidth b at each discrete frequency (m0, m1) (m0, m1 ∈ -M/2, . . . , M/2 -1) a Gaussian with magnitude K b,(m 0 ,m 1 ) ,

Sn (ξ0, ξ1) = b∈B (m 0 ,m 1 ) K 2 b,(m 0 ,m 1 ) 8a 2 b G (ξ0, ξ1) ; ± m0 M , m1 M , a b 2 √ π . (7)
Note that this discretization restricts the values for the frequency ω to the discrete frequencies of the exemplar, but does not restrict the values for the magnitude K. We group the parameters of Eqn. 7 into an B×M ×M parameter vector α (b, m0, m1) = K 2 b,(m 0 ,m 1 ) . We denote the discrete power spectrum of bandwidth-quantized Gabor noise with parameters α as Sn (α).

Using the above notation, the goal of the parameter estimation is to find a parameter vector α satisfying three constraints: (i) Sn (α) is close to Sex, i.e., the power spectrum of the noise is close to that of the exemplar, (ii) α is non-negative, i.e., all Gaussians are positive, and (iii) α is sparse, i.e., the representation is compact. Once such a parameter vector α is found, the parameters of the bandwidthquantized Gabor noise are fully determined. Every non-zero entry of α (b, m0, m1) determines a Gaussian with bandwidth b , magnitude K = α (b, m0, m1) and frequency ω = (m0/M, m1/M ). The number of non-zero entries of α also determines the number of Gaussians G, the number of bandwidths B, and the number of in each bandwidth G b . We define the sparseness ratio of α as G/M 2 , the number of Gaussians relative to the number of pixels of the exemplar.

Solving the Parameter Estimation Problem using Non-Negative Basis Pursuit Denoising

We solve the parameter estimation problem as a non-negative basis pursuit denoising (NNBPDN) problem (e.g., [START_REF] Kim | An interior-point method for largescale l1-regularized least squares[END_REF]

]), minimize Sn(α) -Sex 2 2 + ν α 1 subject to α ≥ 0 , (8)
where • 1 is the ℓ 1 -norm and ν is a regularization parameter. We use non-negative basis pursuit denoising because it satisfies the three constraints above, in particular, because it is designed to produce sparse solutions. This guarantees a compact procedural representation, and can deal with the noisy power spectrum estimate of the exemplar. Note that with an ℓ 2 rather than an ℓ 1 -norm regularization, nearly none of the components of α be zero.

The regularization parameter ν determines the sparseness of the solution. However, the relation between ν and the sparseness depends on the exemplar. We therefore introduce a relative regularization parameter κ, and define ν as κνmax, where νmax is the smallest value of ν for which α = 0 is a solution of Eqn. 8, i.e., νmax = 2 S T n (Sex) ∞ [START_REF] Kim | An interior-point method for largescale l1-regularized least squares[END_REF]]. We expose the relative regularization parameter κ to the user in order to allow the user to trade sparseness for accuracy (and vice versa). Additionally, we solve Eqn. 8 following the homotopy strategy, i.e., for a sequence of kappa values 1/2 i , initializing each step with the result of the previous one. This is more efficient, and provides the user with intermediate solutions. This is illustrated in Fig. 3.

One of the reasons why we have quantized the bandwidth in Sec. 4.3 is to guarantee a robust optimization: A fixed set of bandwidths allows us to formulate our optimization as a convex optimization, while a sufficiently spaced set of bandwidths ensures a well-conditioned optimization.

Implementation

We numerically solve Eqn. 8 using the fast iterative shrinkagethresholding (FISTA) algorithm of [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF], combined with the duality gap stopping criterion associated to the corresponding convex problem [START_REF] Kim | An interior-point method for largescale l1-regularized least squares[END_REF][START_REF] Mairal | Convex and network flow optimization for structured sparsity[END_REF]]. We take symmetry into account by restricting (m0, m1) to the positive half-plane. For large Gaussians, we restrict the possible values for the frequency ω to the discrete frequencies of a sub-sampled version of the exemplar, thus ensuring that Eqn. 8 does not become ill-conditioned. We implemented our solver in Matlab.

Efficient Procedural Evaluation

In this section, we introduce an efficient procedural evaluation for bandwidth-quantized Gabor noise, that uses multi-grid evaluation and importance sampling of the kernel parameters.

Multi-Grid Evaluation

The procedural evaluation of noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]] uses a grid with a cell size equal to the kernel radius r to efficiently genthe random positions of the kernels that overlap the point of evaluation. The kernel radius r, and thus the cell size of the grid, is determined by the bandwidth a of the Gabor kernels, and is defined as r = 1/a. However, this does not work for bandwidthquantized Gabor noise, which consists of kernels with different bandwidths. One of the reasons why we have quantized the bandwidth in Sec. 4.3 is to be able to evaluate bandwidth-quantized Gabor noise: It allows us to interpret bandwidth-quantized Gabor noise as a sum of B independent Gabor noises that each consist only of Gabor kernels with bandwidth b , which, in turn, allows us to efficiently evaluate bandwidth-quantized Gabor noise using a multi-grid evaluation, where a grid with a cell size of r b = 1/a b is associated to each of the B noises. Note that evaluating bandwidth-Gabor noise by independently evaluating and summing the G noises instead would be inefficient, because each grid has an associated setup cost and G is typically much larger than B. Also note that for the multi-grid evaluation the actual value of the bandwidths {a b } does not matter.

Importance Sampling of Kernel Parameters

The procedural evaluation of Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]] allows the user to trade quality for speed (and vice versa) by exposing a parameter called the impulse budget. The impulse budget N is defined as the expected number of kernels that overlaps the point of evaluation, and determines the impulse density λ, which is in turn defined as λ = N/πr 2 , where r is the kernel radius. However, for bandwidth-quantized Gabor noise, which can be seen as a sum of B sums of G b Gabor noises n b,g , the impulse density λ b,g of each of the noises must be determined.

We distribute the impulse budget N over the G noises proportion-ally to their variance σ 2 b,g , since the contribution of each of these noises can differ significantly. This is done by setting N b,g to σ 2 b,g /σ 2 N . This implies that N b = σ 2 b /σ 2 N , the impulse budget for each bandwidth, is proportional to the variance of that bandwidth, and that p b,g = σ 2 b,g /σ 2 b , the kernel parameters, are selected with a probability proportional to the variance of the corresponding Gaussian.

In contrast to Gabor noise, bandwidth-quantized Gabor noise thus requires sampling discrete probability distributions to determine the parameters of each kernel. The efficiency of the sampling procedure is essential for performance. Additionally, the efficiency of the preprocess that is required for the sampling procedure is essential as well, since these distributions change whenever the noise parameters change (e.g., during editing).

We therefore use the alias method [START_REF] Walker | An efficient method for generating discrete random variables with general distributions[END_REF]], which features a constant-time sampling procedure and a linear-time preprocess. To sample a discrete probability distribution, the alias method builds two tables, the alias table A and the probability table F , after which sampling is done by generating a random index i, and returning the i-th value with probability F [i], or the A[i]-th value with probability 1 -F [i]. Several linear-time algorithms for building the tables required by the alias method are available [START_REF] Walker | An efficient method for generating discrete random variables with general distributions[END_REF][START_REF] Vose | A linear algorithm for generating random numbers with a given distribution[END_REF]]. However, these algorithms require linear auxiliary storage (typically two linked lists), which can be problematic (e.g., when using CUDA or GLSL) and inefficient (due to dynamic memory allocations). Following a hint in the work of Vose [1991, Sec. VI], we have implemented an in-place linear-time algorithm that does not require auxiliary storage.

Implementation

In contrast to Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]], bandwidth-quantized Gabor noise operates at much higher impulse densities. We therefore generate random numbers distributed according to a Poisson distribution with mean λ using the Gaussian approximation P = ⌊λ + √ λX + 1/2⌋, where X is distributed according to the standard normal distribution. We generate random numbers distributed according to the standard normal distribution using the basic form of Box-Muller transform X = -2 ln (U1) cos (2πU2), where U1 and U2 are distributed according to the standard uniform distribution. We seed the random number generator of each cell using an approach based on random number tables, seed(b, x, y) =

P b [b%L] ⊕ Px[x%L] ⊕ Py[y%L]
, where P b , Px and Py are random number tables of size L (typically 256), and % and ⊕ denote modulo and XOR. Please see the supplemental material for a more detailed discussion. We implemented our procedural evaluation in CUDA.

Maximally Independent Color Space

In this section, we introduce the maximally independent color space, a new color space for independent channel synthesis, which we use to synthesize color Gaussian textures with our method. Note that this is only one specific way to obtain colored noise-based textures. Our noise patterns can be used in any way conventional noise patterns are used.

Previous work, inspired by the work of [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF], typically performs independent channels synthesis in a PCA color space, obtained by diagonalizing the covariance matrix using principal component analysis (PCA). However, the use of a PCA color space in this context is somewhat contrived, both in theory, since color channels with zero covariance are not necessarily independent, and in practice, since some previous work reports it works well (e.g., [START_REF] Qin | Aura 3d textures[END_REF]), while some reports the contrary (e.g., [START_REF] Kopf | Solid texture synthesis from 2D exemplars[END_REF]). We therefore introduce the maximally independent color space, a new color space for independent channel synthesis where the color channels are maximally independent.

We explain our new maximally independent color space, as well as its relation to the PCA color space used in previous work, using techniques from independent component analysis (ICA) [Hyvärinen et al. 2001].

We define the lagged correlation matrices [Hyvärinen et al. 2001, 18.1]

C I (ξ) of the M ×N RGB image I (x) as C I (ξ) = 1 M 2 x (I (x) -µ I ) (I (x + ξ) -µ I ) T =   C I R ,I R (ξ) C I R ,I G (ξ) C I R ,I B (ξ) C I G ,I R (ξ) C I G ,I G (ξ) C I G ,I B (ξ) C I B ,I R (ξ) C I B ,I G (ξ) C I B ,I B (ξ)   , (9)
where ξ is the lag and µ I is the mean of I, and C I i ,I j (ξ) is the cross-correlation of i-th channel with the j-th channel shifted over ξ. The lagged correlation matrices of the image I after applying a color transform T are TC I (ξ) T T .

Using this notation, the PCA color space used in previous work is obtained by finding a matrix TP CA that diagonalizes the correlation matrix C I (0) using PCA. Although this implies that the color channels have zero covariance, the matrix TP CA does not diagonalize the lagged correlation matrices C I (ξ) for ξ = 0, which implies that the channels are not independent.

We instead obtain our new maximally independent color space by finding the matrix TAJD that approximately diagonalizes all lagged correlation matrices C I (ξ), and not just C I (0). This is motivated by the fact that diagonalizing the lagged correlation matrices transforms the mixed channels into independent ones in the ICA framework [Hyvärinen et al. 2001, 18.1]. We find the matrix TAJD approximate joint diagonalization (AJD) [Hyvärinen et al. 2001, 11.3]. More specifically, we minimize the objective function J ,

J (T) = ξ off TC I (ξ) T T , (10)
where off (A) = i =j a 2 ij is a measure for the diagonality of a matrix, using the algorithm of [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF].

The value of J (T) is a measure for the independence of the color channels in the color space determined by T. This allows us to compare the effectiveness of the RGB, PCA and maximally independent color space by comparing the values of J (TRGB = I), J (TP CA) and J (TAJD). Please refer to the supplemental materials for more details.

Interactive Noise Editing

In this section, we present an interactive editor for bandwidthquantized Gabor noise, inspired by the Wold decomposition [START_REF] Francos | A unified texture model based on a 2-D Wold-like decomposition[END_REF]]. Although our approach is exemplar-based, it still allows editing for maximum flexibility.

Our editor is inspired by the Gabor noise widgets of Lagae et [2009, 3.3], which allow the user to manipulate the parameters of a Gabor noise by directly manipulating the corresponding Gaussians in the spectral domain. However, in contrast to Lagae et al., we have to deal with tens or even hundreds of Gaussians. We therefore group the Gaussians, and allow the user to manipulate the groups rather than the individual Gaussians. We group the Gaussians automatically inspired by the Wold decomposition [START_REF] Francos | A unified texture model based on a 2-D Wold-like decomposition[END_REF], which states that a Gaussian random field can be decomposed into a stochastic component, a set of strongly anisotropic components, and a set of periodic components. We use guidelines from the Wold decomposition [START_REF] Liu | Modeling spatial and temporal textures[END_REF]] to perform a basic decomposition of the sparse set of Gaussians. Please refer to the supplemental materials for more details. Note that, in contrast to previous work [START_REF] Francos | A unified texture model based on a 2-D Wold-like decomposition[END_REF][START_REF] Liu | Modeling spatial and temporal textures[END_REF]] our goal is not to fully decompose a discrete power spectrum, but rather to provide meaningful groups of Gaussians to the user, which is less challenging. We allow the user to manipulate groups by rotating, scaling and translating groups, which affects the frequency ω of the Gaussians in the group, and by changing the magnitude and bandwidth of groups, which affects the magnitude K and bandwidth a of the in the group. We also allow the user to manually group Gaussians, and to cut, copy and paste groups.

Results, Comparisons and Discussion

Results

In Fig. 4, we show detailed results of our approach. Col. 1 shows the exemplar. Col. 2 shows the Gaussian version of the exemplar. Differences between col. 1 and col. 2, if any, are due to the the fact that the exemplar is not necessarily Gaussian. Col. 2 also serves as ground truth for our approach. Col. 3 shows the Gaussian texture corresponding to the sparse power spectrum estimated using our parameter estimation. Differences between col. 2 and col. 3, if any, are due to the approximate nature of the parameter estimation. Col. 4 shows the noise generated by our procedural evaluation. Differences between col. 3 and col. 4, if any, are due to the various approximations in the procedural evaluation.

In Fig. 1, we show more results of our approach for the grayscale case. This figure shows the Gaussian version of the exemplar (row 1) and the procedural noise (row 2), and illustrates that our approach works very well for Gaussian textures. In Fig. 5 (and in Fig. 7(col. 1-2)), we show more results of our approach for the color case. This figure shows the exemplar (row 1) and the procedural noise (row 2), and illustrates that our approach works well for many nearly-Gaussian textures (a-e), but less well for non-Gaussian texures (f-h), including structured textures and textures with non-Gaussian color distributions. Please refer to the supplemental material for detailed results for a set of 44 textures, including all examples in the paper. The two parameters of our approach are (i) the relative regularization parameter κ, which trades sparseness for accuracy in the parameter estimation, and (ii) the impulse budget N , which trades quality for speed in the procedural evaluation. We used the default parameters κ = 1/256 and N = 1024 for all examples, which in practice makes our method automatic. Note that for many exemplars, better results could be obtained by tweaking κ and N .

The performance of our parameter estimation is roughly two minutes per texture (unparallelized Matlab, 2.67GHz Intel Xeon X5650). However, we believe this could be significantly improved if needed. The performance of our procedural evaluation is roughly 30 FPS (CUDA, NVIDIA Quadro 5000). Both numbers are for the grayscale case, an image resolution of 128 × 128, and the default parameters. The performance is roughly linear in the image resolution and N . The performance is only marginally lower in the color case, since typically one channel dominates in the maximally independent color space. Please refer to the supplementary materials for exact values.

The sparseness of our representation of the estimated power spectrum is on the average roughly 3.5%, for the default value of κ.

Please refer to the supplementary materials for exact values. These results suggest that the power spectrum of Gaussian textures can be sparsely represented.

Our maximally independent color space improves results over the PCA color space, but in practice this difference is barely noticeable. For nearly-Gaussian exemplars, J (TRGB)>>J (TP CA)> J (TAJD), but J (TAJD) is relatively close to J (TP CA), e.g., for Fig. 4(a), the values are respectively 0.9962, 0.0067 and 0.0028. Please refer to the supplementary materials for all values. Nevertheless, the fact that J (TP CA) is relatively close to J (TAJD), the optimal value, explains why the PCA color space has been successful for stochastic textures in previous work. We have not found any Gaussian texture where independent channel synthesis in the maximally independent color space fails, which seems to suggest that the color channels of Gaussian textures are largely independent. For non-Gaussian exemplars, such as [Kopf et al. 2007, Fig. 2], the maximally independent channels might still be too correlated to allow successful independent channel synthesis. Interestingly, for this example, both J (TRGB) and J (TP CA) are relatively large, which seems to suggest that neither color space will allow successful independent channel synthesis.

In Fig. 6 and in the video sequences in the supplemental materials we show results of our interactive noise editor. These demonstrate several editing operations. Note that none of the other examples has been edited, and that the editor adapts the impulse budget to guarantee interactive performance.

Our approach inherits all advantages of Gabor noise, including surface mapping without a texture parameterization (surface Gabor noise) and analytic texture filtering [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF]].

Comparisons

In Fig. 7, we compare our approach to the methods of [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF], [START_REF] Lefebvre | Parallel controllable texture synthesis[END_REF], Lagae et al. [2010b], and Gilet et al. [2010]. Please refer to the supplemental material for a more detailed version of this comparison for a set of 19 textures. Col. 1 shows the exemplar. Col. 2 shows our results. Col. 3 shows the results of Gilet et al. These results were generated by the authors of [Gilet et al. 2010] themselves, who reported that significant manual intervention and trial and error were needed to achieve these results, while our approach is automatic. Additionally, their results typically exhibit a lack of high-frequency detail, increased regularity, and/or subtle color shifts. Col. 4 shows the results of Lagae et al. Their method cannot handle anisotropy. Col. 5 shows the results of Heeger and Bergen. Their method cannot handle anisotropy well. Note that similar results would probably be obtained with anisotropic noise [START_REF] Goldberg | Anisotropic noise[END_REF], which uses the same steerable filters. Col. 6 shows the results of Lefebvre and Hoppe. Their method is not procedural, and cannot preserve the anisotropic detail well in all cases, although this could be alleviated by manually constraining the jitter.

The method of [START_REF] Ghazanfarpour | Spectral analysis for automatic 3-d texture generation[END_REF] cannot handle noisy textures such as sand, while our method can, and requires the to manually set parameters (i.e., Erra and Err φ ), while our is automatic. The method of [START_REF] Dischler | A procedural description of geometric textures by spectral and spatial analysis of profiles[END_REF] cannot handle anisotropic textures, while our method can.

Discussion

A simpler alternative to our parameter estimation would be to obtain the K and ω parameter of each Gabor kernel by importance sampling the power spectrum of the exemplar during the procedural evaluation. However, it is not clear how to obtain the a parameter in this case, and fixing it to a large or small bandwidth would respectively be inaccurate or inefficient. Additionally, this alternative is not compact.

An alternative to non-negative basis pursuit denoising are methods based on Gaussian mixture models (GMM) and expectation maximization (EM), such as the one by [START_REF] Papas | Goal-based caustics[END_REF]. However, non-negative basis pursuit denoising has several advantages over these methods for our problem: (i) it allows the user to specify the desired accuracy of the approximation rather than the desired number of Gaussians; (ii) it does not require an initialization; and (iii) it provides a solution in terms of quantized rather than continuous bandwidths.

Compared to our approach, Gilet et al. [2010] additionally use cosines and histogram matching, which can improve the results in some cases. Although our approach is compatible with these additions, we have chosen not use them, since the usage of cosines and histogram matching is not compatible with surface Gabor noise and analytic filtering of Gabor noise [START_REF] Lagae | Procedural noise using sparse Gabor convolution[END_REF].

Conclusion

We have presented Gabor noise by example, a generalization of Gabor noise that can generate a wide variety of Gaussian textures, and can estimate the noise parameters from an exemplar Gaussian texture. Gabor noise by example preserves the traditional advantages of procedural noise, including a compact representation and fast on-the-fly evaluation, and is mathematically well-founded.

For all of our examples, the Gaussian version of the exemplar and the procedural noise are virtually indistinguishable. This means that all remaining information that makes the exemplar different from its Gaussian version is in the phase spectrum of the exemplar. We therefore believe that understanding the phase spectrum of textures is an important direction for future research.

Figure 2 :

 2 Figure 2: Bandwidth-quantized Gabor noise. (a) A power spectrum composed of Gaussians with discrete bandwidths. (b-g) The power spectrum split out into its discrete bandwidths. The bandwidth quantization enables robust parameter estimation and efficient procedural evaluation.

Figure 3 :

 3 Figure 3: Robust parameter estimation for bandwidth-quantized Gabor noise. (a) A Gaussian exemplar and its power spectrum estimate. (b-h) Sparse decompositions of the power spectrum estimate into Gaussians with discrete bandwidths for decreasing values of the relative regularization parameter κ. The sparse decomposition ensures a compact procedural representation.

Figure 4 :

 4 Figure 4: Detailed results. (col. 1) Exemplar. (col. 2) Gaussian version of the exemplar. (col. 3) Gaussian texture corresponding to the sparse power spectrum estimated using our robust parameter estimation. (col. 4) Procedural noise obtained using our efficient procedural evaluation. (insets) The power spectrum estimate or power spectrum. (row 1) Grayscale results. (row 2) Color results obtained using our maximally independent color space.

Figure 5 :

 5 Figure 5: More results. (row 1) Exemplar. (row 2) Procedural noise. (insets) Sparse power spectrum representation obtained using our robust parameter estimation. (a-e) Nearly-Gaussian textures. (f-h) Non-Gaussian textures.

Figure 6 :

 6 Figure 6: Interactive noise editing. Two screen captures of an interactive editing session. Please refer to the supplemental materials for the corresponding video sequence.

Acknowledgements

We would like to thank the anonymous reviewers and Jean-Michel Dischler, Fredo Durand, Guillaume Gilet and Gabriel Peyré.

Ares Lagae is a Postdoctoral Fellow of the Research Foundation -Flanders (FWO). This work was partly supported by ANR project MATAIM (ANR-09-BLAN-0029-01), KU Leuven project CREA/08/017, ANR project SIMILAR-CITIES (2008-COORD-021-01), EU project VERVE (http://www.verveconsortium.eu/), Autodesk (research and software donations), NVIDIA (Academic Partnership Program) and Adobe (research donation).


Laurent Girin
email: laurent.girin@gipsa-lab.grenoble-inp.fr

Jonathan Pinel

Informed Audio Source Separation from Compressed Linear Stereo Mixtures

In this paper, new developments concerning a system for informed source separation (ISS) of music signals are presented. Such system enables to separate I > 2 musical instruments and singing voices from linear instantaneous stationary stereo (2-channel) mixtures, based on audio signal natural sparsity, pre-mix source signal analysis, and side-information embedding (within the mix signal). The foundations of the system have been presented in previous papers, within the framework of uncompressed (16-bit PCM) mix signals.

In the present paper, we extend the study to compressed mix signals. For instance, we use a MPEG-AAC codec and we show that the ISS process is quite robust to compression, opening the way for "real-world" karaoke/soloing/remixing applications for downloadable music.

INTRODUCTION

Nowadays, the public is more and more eager to interact with musical pieces during restitution. For example, more and more commercial hits are distributed with several versions, some being instrumental versions (for karaoke), other being remixes. The karaoke phenomenon gets generalized from voice to instruments in musical video games such as Rock Band and Guitar Hero. In this case, the interaction with the musical piece is based on the inclusion of the multitrack recording in the video game. Indeed, multitrack formats (i.e. concatenation of the separate instruments/voices recordings in a single large file) offer the possibility to manipulate separately the different sound sources while remixing the musical piece in real-time, e.g. controling their volume and spatialization, or adding individual audio effects, leading to a new "active" mode of music listening.

Several multitrack formats are present on the music market, but their use remains marginal for the general public since those formats are confronted to the predominance of the stereo (2-channel) format. Moreover, the music industry is still reluctant to release the multitrack version of musical hits. Because of this historical inertia of usage and technology, it seems quite unrealistic to think that most of the music listeners would be ready to switch rapidly from stereo to multitrack usage if they were offered new interactive possibilities during music listening. In contrast, it is likely that a significant amount of music listeners would entertain the active listening experience if they could access to it from usual stereo format. This is the goal of the DReaM project. 1 This project aims at developping source separation methods that enable the efficient separation of many different sound sources from only 2-channel mixture signals. Hence we enter the socalled underdetermined difficult (because non-invertible) configuration of the source separation problem. Moreover, the mixture signals should be realistic music pieces, approaching professional music production, and the separation should be processed in real-time with reasonable computation cost, so that real-time source manipulation and remixing can follow. Such ambitious goal is out of reach for current Blind (or semi-blind) Source Separation (BSS) techniques and related Independent Components Analysis (ICA), i.e. separation techniques that only proceed from (determined and overdetermined) mixture signals with very few information about the sources and mixture process characteristics [START_REF] Cardoso | Blind signal separation: Statistical principles[END_REF][2] [START_REF] Comon | Handbook of Blind Source Separation -Independent Component Analysis and Applications[END_REF]. It is even still a quite difficult task for informed source separation (ISS) methods, where additionnal prior "external" information about the sources is provided, for example with midi/score files [START_REF] Dubnov | Optimal filtering of an instrument sound in a mixed recording using harmonic model AES 42 ND INTERNATIONAL CONFERENCE[END_REF] [START_REF] Woodruff | Remixing stereo music with score-informed source separation[END_REF] or melody humming [START_REF] Smaragdis | Separation by "humming": User-guided sound extraction from monophonic mixtures[END_REF].

Therefore, in the DReaM project, an extreme ISS configuration is considered, in the sense that the source sig-nals and the mixture process are assumed to be perfectly known. However, the separation does not exploit those information directly: we rather propose a two-step specific coder-decoder configuration. The coder corresponds to the music signal production level (e.g. music recording/mixing in studio) where the source signals are assumed to be available (because of separate recordings) and the mixing process is controlled. The decoder corresponds to the music signal restitution level (e.g. audio-CD at home or mobile player), where only the stereo mix signal is available. Parameters that characterise the source signals and the mixing process are embedded into the mixture signal at the coder level, so that they can be retrieved at the decoder and exploited for source signals separation from the mix signal.

In [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8], we considered linear instantaneous stationary stereo (LISS) mixtures, a.k.a. constant-gain stereo panning. Basically, the embedded information consists of i) the mixture matrix, and ii) the indexes of the two predominant sources in each time-frequency (TF) region as provided by an analysis of the source signals at the coder. Thus, source separation at the decoder uses those information to perform local 2 × 2 inversion of the mixture in each TF region. The separation performances obtained by this system are quite impressive: Gains of about 20dB-SDR [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF] are obtained for all sources of 5source mixtures, enabling realistic applications such as generalized karaoke/soloing. Another implementation of the ISS principle within the DReaM project can be found in [START_REF] Liutkus | Informed source separation using latent components[END_REF] within the general framework of Non-negative Matrix Factorization (NMF).

The mixture signals considered in [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8] were 16-bit PCM uncompressed signals, and the side-information was embedded using a high-capacity watermarking technique [START_REF] Pinel | A highcapacity watermarking technique for audio signals based on MDCT-domain quantization[END_REF] well suited for such signals: straightly stated, the watermarking is shown to be inaudible and have negligible effects on separation performances. Therefore, the resulting ISS system can be applied to audio-CD and .wav music, and a real-time software demonstrator of the decoder has already been achieved and presented in [START_REF] Marchand | Interactive music with active audio CDs[END_REF]. Since audio compression [START_REF] Brandenburg | Overview of MPEG audio: Current and future standards for low bit-rate audio coding[END_REF] has become a very popular tool for ubiquitous listening of "dematerialized" music, audio applications must consider compressed audio formats. The goal of the present study is to demonstrate the feasability of ISS within the audio compression framework. Therefore, in the present paper we still consider LISS mixtures, but the mix signal is compressed, for instance using MPEG-AAC [START_REF]11 MPEG. Information technology-Generic coding of moving pictures and associated audio information[END_REF]. Because audio compres-sion and watermarking are very difficult to conciliate, the side-information is here transmitted using metadata segments of the compressed binary stream, or adjunct channels, such as in MPEG-SAC/SAOC [START_REF] Herre | MPEG surround -the ISO/MPEG standard for efficient and compatible multi-channel audio coding[END_REF] [START_REF] Herre | New concepts in parametric coding of spatial audio: from SAC to SAOC[END_REF]. We show that, although the separation performances are degraded compared to the use of PCM mix signals, the sparsenessbased index-embedding ISS technique is remarkably robust to audio compression of the mix signal at reasonable bitrates. We also show that using a refined coding of the side-information enables a quite reasonable additional transmission cost compared to the bitrate of the compressed mix. This paper is organized as follows. In Section 2, we present an overview of the ISS system based on audio signals sparsity in the TF domain and local inversion of the mixture. The general principle of the system is common for both uncompressed and compressed cases, but some implementation details in Section 2 rather concern the uncompressed case. The adaptation of the system to the new case of compressed mix signals is provided in Section 3. Results obtained for music mixtures are presented in Section 4, before concluding in Section 5.

ISS BASED ON SIGNAL T-F SPARSITY AND LOCAL MIX INVERSION

General Principle

The functional schema of the proposed ISS method is represented on Fig. 1 for the case of uncompressed signals (as previously presented in [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF] [START_REF] Parvaix | Informed source separation of linear instantaneous under-determined audio mixtures by source index embedding[END_REF]). The decoder is the "active player", that can process separation only on mix signals that have been generated by the corresponding coder. At the coder, the mix signal is generated as a linear instantaneous stationary stereo (LISS) mixture, i.e. summation of individual source signals with constantgain panning coefficients. Then, the system looks for the two sources that better "explain" the mixture (i.e. the two source signals that are predominant in the mix signal) at different time intervals and frequency channels, and the corresponding source indexes are encoded and embedded into the mixture signal as side-information. At the decoder, the only available signal is the (embedded) mix signal. The side-information is extracted from the transmitted mix signal, decoded and used to separate the source signals by a local time-frequency mixture inversion process. In [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF] [START_REF] Parvaix | Informed source separation of linear instantaneous under-determined audio mixtures by source index embedding[END_REF], the embedding of the sideinformation is done by using a high-capacity watermarking method similar to the one presented in [START_REF] Pinel | A highcapacity watermarking technique for audio signals based on MDCT-domain quantization[END_REF], and the AES 42 ND INTERNATIONAL CONFERENCE, Ilmenau, Germany, 2011 July 22-24 (watermarked) mix signal is quantized to 16-bits PCM (for audio-CD / wav file compatibility).

MDCT x L [n] x R [n] X R [f,t] s 1 [n] s 2 [n] s I [n] Oracle estimator S I [f,t] S 1 [f,t] I ft Water- marking Mixing (LISS) Mixing matrix A Water- mark extraction MDCT MDCT Decoding 2 × 2 inver- sion A IMDCT IMDCT IMDCT ] , [ ˆ1 t f S ] , [ ˆt f S I ] [ ˆ1 n s ] [ ˆ2 n s ] [ ˆn s I MDCT MDCT MDCT MDCT Coding IMDCT IMDCT ] , [ t f X W L ] , [ t f X W R ] [ ~n x W L 16-bits PCM 16-bits PCM ] [ ~n x W R ] , [ ~t f X W L ] , [ ~t f X W R ] [ ~n x W L ] [ ~n x W R X L [f,t] I ft

Time-Frequency Decomposition

The source signals of interest are voice/instrument signals playing a same piece of music (but recorded separately). They are non-stationary, with possibly large temporal and spectral variability, and they generally strongly overlap in the time domain. Decomposing the signals in the time-frequency (TF) domain leads to a sparse representation, i.e. few TF coefficients have a high energy and the overlapping of signals is much lower in the TF domain than in the time domain [START_REF] Plumbley | Sparse representations in audio and music: from coding to source separation[END_REF]. Therefore, the separation of source signals can be carried out more efficiently in the TF domain [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF][19] [START_REF] Araki | Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors[END_REF]. The Modified Discrete Cosine Transform (MDCT) [START_REF] Princen | Analysis/synthesis filter bank design based on time domain aliasing cancellation[END_REF] is used here as the TF decomposition since it presents several properties very suitable for the present problem: good energy concentration (hence emphasizing audio signals sparsity), very good robustness to quantization (hence robustness to quantization-based watermarking), orthogonality and perfect reconstruction. We will see in the following that the MDCT is also particularly related to the problem of mix compression.

On Fig. 1, the MDCT is applied on the source signals and on the mixture signal at the input of the coder to enable the selection of predominant sources in the TF domain. In the case of uncompressed signals, watermarking of the resulting side-information is applied on the MDCT coefficients of the mix signal (see [START_REF] Pinel | A highcapacity watermarking technique for audio signals based on MDCT-domain quantization[END_REF]) and the time samples of the watermarked mix signal are provided by inverse MDCT (IMDCT). At the decoder, the (PCM-quantized) mix signal is MDCT-transformed and the side-information is extracted from the resulting coefficients. We will see that in the case of compressed signals, no IMDCT reconstruction is necessary at the coder since the compressed binary stream is transmitted (together with the side-information; see Section 3.2). In both uncompressed and compressed cases, source separation is carried out at the decoder also in the MDCT domain, and the resulting separated MDCT coefficients are used to reconstruct the corresponding time-domain separated source signals by IMDCT.

We do not provide here the detailed MDCT equations, since they can be found in many dedicated papers, e.g. [START_REF] Princen | Analysis/synthesis filter bank design based on time domain aliasing cancellation[END_REF]. Let us just mention here that the MDCT / IMDCT is applied on signal time frames of W = 2048 samples (46.5ms for 44.1kHz-sampling), with a 50%-overlap between consecutive frames (of 1024 frequency bins). The frame length W is chosen to follow the dynamics of music signals while providing a frequency resolution suitable for the separation. Appropriate windowing is applied at both analysis and synthesis to ensure the "perfect reconstruction" property [START_REF] Princen | Analysis/synthesis filter bank design based on time domain aliasing cancellation[END_REF].

Local mix inversion and sources selection

Since the MDCT is a linear transform, the LISS source separation problem remains LISS in the transformed domain. For each frequency bin f and time bin t, we thus have:

X[ f ,t] = A • S[ f ,t] (1)
where

X[ f ,t] = [X 1 [ f ,t], X 2 [ f ,t]] T denotes the stereo mixture coefficients vector, S[ f ,t] = [S 1 [ f ,t], • • • , S I [ f ,t]]
T denotes the I-source coefficients vector, and A denotes the mixing matrix that gathers the stereo panning coefficients. Because of audio signal sparsity in the TF domain, only at most 2 sources are assumed to be relevant, i.e. of significant energy, at each TF bin ( f ,t). Therefore, the mixture is locally given by:

X[ f ,t] ≈ A I f t S I f t [ f ,t] (2)
where I f t denotes the set of 2 relevant sources at TF bin ( f ,t), i.e. the two source signals that locally "better explain" the mixture. A I f t represents the corresponding 2×2 mixing sub-matrix. If I f t denotes the complementary set of non-active (or at least poorly active) sources at TF bin ( f ,t), the source signals are estimated by [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]:

ŜI f t [ f ,t] = A -1 I f t X[ f ,t] ŜI f t [ f ,t] = 0 ( 3
)
where A -1 I f t denotes the inverse of A I f t (which is assumed to be invertible and well-conditioned). Note that such separation technique exploits the 2-channel spatial information of the mixture signal and relaxes the restrictive assumption of a single active source at each TF bin, as made in [START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] [START_REF] Araki | Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors[END_REF].

The side-information that is transmitted between coder and decoder (in addition to the mix signal) consists of the coefficients of the mixing matrix A and the combination of indexes I f t that identifies the predominant sources in each TF bin. This contrasts with classic blind and semiblind separation methods where those information have to be estimated from the mix signal only, generally in two steps which can both be a very challenging task and source of significant errors.

As for the mixing matrix, the number of coefficients to be transmitted is quite low in the present LISS configuration (only I fixed coefficients for each piece of music, if A is made of normalized column vectors depending on source azimuths). Therefore, the transmission cost of A is negligible compared to the transmission cost of I f t , and in the following we do not detail the encoding and transmission of A.

As for the source indexes, the optimal I f t is estimated at the ISS coder using the source signals, the matrix A, and the mixture signals. This is done using an Oracle estimator, as introduced in [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF] for providing upper bounds for the performances of source separation algorithms. Exploiting the reconstruction properties of the MDCT, the overall best separation in the time domain in the mean squared error (MSE) sense is obtained by finding the optimal combination of source signals in the MSE sense at each TF bin separately [START_REF] Vincent | Oracle estimators for the benchmarking of source separation algorithms[END_REF]:

Ĩ f t = arg min I f t ∈P I ∑ i=1 Ŝi [ f ,t] -S i [ f ,t] 2 ( 4
)
where P represents the set of all possible combinations I f t and the I estimated source signals Ŝi ( f ,t) are provided by (3). In the case I is limited to a small number of sources (typically about 5 for standard western popular music), Ĩ f t can be found by exhaustive search, and coded with a very limited number of bits before being embedded into the mixture signal (see Section 3.3).

Source separation at the decoder

At the coder, Ĩ f t is determined from the "original" mix signal (Equation ( 4)). In contrast, at the decoder, only the transmitted mix signal is available, and the source separation is obtained by applying Equation (3) using Ĩ f t and using the MDCT coefficients of the transmitted mix signal instead of the MDCT coefficients of the original mix signal. In the uncompressed case, these MDCT coefficients are the watermarked coefficients recalculated from the watermaked (and PCM quantized) time-domain signal. It has been shown in [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8] that the influence of the watermarking (and PCM quantization) on separation performance is negligible (for 5-source mixtures). This is because the side-information bitrate that is used to encode the optimal combination for each TF bin (64kbps per channel; see [START_REF] Parvaix | Informed source separation of linear instantaneous under-determined audio mixtures by source index embedding[END_REF] and below) is significantly lower than the maximal bitrate allowed by the high-capacity watermarking technique (approx. 250kbps per channel [START_REF] Pinel | A highcapacity watermarking technique for audio signals based on MDCT-domain quantization[END_REF]). For such bitrate the watermarked MDCT coefficients remain very close to the original MDCT coefficients, hence the separation performances are almost identical for unwatermarked and watermarked mix signals. In fact, the degradation is mainly due to the "failure" of the sparsity assumption, i.e. the fact that nonpredominant but non-null sources interfere in the local inversion process.

Those performances are described in details in [START_REF] Parvaix | Informed source separation of linear instantaneous under-determined audio mixtures by source index embedding[END_REF] for "real-world" 5-source LISS music mixtures of different musical styles. Basically, source enhancement from input (mix) to output (separated) ranges from 17dB to 25dB depending on sources and mixture, which is remarkable given the difficulty of such underdetermined mixtures. The rejection of competing sources is very efficient and the source signals are clearly isolated, as confirmed by listening tests. Artefacts (musical noise) are present but are quite limited. The quality of the isolated source signals makes them usable for generalized karaoke/soloing/remix applications.

ISS FOR COMPRESSED MIX SIGNAL

In the present study, MPEG-AAC compression is applied to the mix signal at the output of the coder, and corresponding decoding is applied at the input of the decoder, as illustrated on Fig. 2. This combination of ISS with compression has two main consequences: First, it induces a perturbation of the mixture signal MDCT coefficients that is likely to influence the separation performances; and second it modifies the side-information embedding process. We develop those two points in the next subsections, and we propose a new coding scheme for the side-information that ensures both lower side-info bitrate and more efficient source separation as compared to our previous works [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8].

Separation from compressed signal

As for the uncompressed case, the separation at the decoder consists of applying (3) using the MDCT coefficients of the transmitted mix signal instead of the coefficients of the original mix signal. In the new compressed case, the MDCT coefficients X C [ f ,t] are recalculated from the compressed (decoded) mix signal xC . Therefore, we have:

Ŝ Ĩ f t [ f ,t] = A -1 Ĩ f t • X C [ f ,t] Ŝ Ĩ f t [ f ,t] = 0 (5)
As opposed to the watermarking process of the uncompressed case, the compression of the mix signal is here expected to modify the MDCT coefficients X C [ f ,t] in a significant way, i.e. with potential perturbation of the separation process (compared to the use of the original MDCT coefficients X[ f ,t]). The degree of degradation is expected to depend on the compression bitrate, and this dependancy may not be straightforward: on the one hand, a poor-quality low-bitrate mix signal is expected to lead to degraded source separation, but on the other hand the good quality of a "reasonable-rate" mix ensured by the "global" perceptual masking of quantization error does not guarantee to offer good quality of separated sources. The purpose of this paper is precisely to report the results of an experimental study of the robustness of separation performances to compression, for realworld music signals mixed with the LISS configuration and compressed at different rates (see Section 4).

An important remark can be made here. In the present study we use an MPEG-AAC codec [START_REF]11 MPEG. Information technology-Generic coding of moving pictures and associated audio information[END_REF]. It happens that MPEG-AAC compression is based on the quantization of MDCT coefficients. Therefore, in a more "compact" structure, the quantization, source selection and separation processes could be applied on the same MDCT coefficients, without recalculation of these coefficients from the decoded signal at the decoder. In the present study we do not achieve such integrated system. Rather, we use an "external" AAC codec, because this approach is the "easy way" to carry first experiments on the merging of ISS and compression without having to investigate any codec source code. This suboptimal approach has the advantage to be independent of the codec: in the diagram of Fig. 2, the AAC codec can be replaced by any other codec without changing the global ISS structure. The development of a (computationaly more efficient) fully integrated AAC-ISS system where the separation is directly carried out from the quantized coefficients is part of our current works. The present study can be seen as a preliminary feasability study for such advanced AAC-ISS merged technology, assuming that the degradation of the MDCT coefficients recalculated from the decoded mix is similar to the degradation of the MDCT coefficients quantized in the compression process.

Side-information embedding

As mentioned before, in [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8] a high-rate watermarking technique [START_REF] Pinel | A highcapacity watermarking technique for audio signals based on MDCT-domain quantization[END_REF] was used to embed the side-information (A and Ĩ f t ) within the mix signal waveform. In the present case, audio compression and watermarking are quite antinomic: it is difficult to simultaneously encode a AES 42 ND INTERNATIONAL CONFERENCE, Ilmenau, Germany, 2011 July 22-24 2: Functional diagram of the ISS system for compressed stereo mix. signal or a set of coefficients with as few binary resource as possible, and embed additional information within the same signal or set of coefficients. Therefore, we consider that the side-information is here transmitted using metadata segments of the compressed binary stream. It can also be transmitted using a dedicated channel. 2 In the present paper, we do not detail the metadata segment format and the corresponding embedding of the side-information. However, we provide in Section 4 the side-info bitrates obtained for our test signals, as derived from the variable-rate encoding process described below.

MDCT x L [n] x R [n] X L [f,t] X R [f,t] s 1 [n] s 2 [n] s I [n] Oracle estimator AAC Compr- ession S I [f,t] S 1 [f,t] I[f,t]
I[f,t] 2 × 2 inver -sion Mixing matrix A IMDCT IMDCT IMDCT ] , [ ˆ1 t f S ] , [ ˆt f S I ] [ ˆ1 n s ] [ ˆ2 n s ] [ ˆn s I MDCT MDCT MDCT MDCT ] [n x C L ] [n x C R Fig.
Overall, the spirit of the present ISS-AAC system is connected to the spirit of MPEG-SAC/SAOC [START_REF] Herre | MPEG surround -the ISO/MPEG standard for efficient and compatible multi-channel audio coding[END_REF][16] and related pioneering works on parametric multi-channel joint audio coding (e.g. [START_REF] Faller | Parametric joint-coding of audio sources[END_REF]) where spatialization parameters that enable to respatialize a multi-channel audio scene from a mono/stereo downmix signal are transmitted with the latter. However, it must be noted that in contrast these systems, the goal of the proposed ISS system is to completely separate the source signals and not only to resynthesize/respatialize the audio scene from down-2 Note that the combination of compression and watermarking has been tackled in several prospective papers, e.g., [23][24], but this remains a difficult task, especially for "high" data rates. Moreover, in this specific ISS framework, joint watermarking/compression of MDCT coefficients is likely to provide a higher degradation of the separation than compression alone. This point is part of our future investigations. mix signals. As a result, the nature of transmitted sideinformation and the way it is exploited are quite different from spatial audio coding systems. In particular, the spatialization parameters in SAC/SAOC are used to "redistribute" the content of spectral bands of the downmix signal across the different output channels, but they cannot separate the contribution of two different sources that are present within the same band (hence the sources are "respatialized together" and not clearly separated; see, e.g., [START_REF] Faller | Parametric joint-coding of audio sources[END_REF]). In contrast, the separation of two overlapping sources is precisely the goal of our local inversion process, that is carried out for each TF bin.

Side-information coding

In [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8], the total number of sources I was fixed for all MDCT frames and frequencies of a given piece of music. Therefore, the number of source combinations was also fixed and Ĩ f t was coded with a fixed number of bits (e.g., 10 combinations of 2 active sources out of 5, coded with 4 bits). In the present study, we refine the coding of Ĩ f t by exploiting the fact that musical sources generally have some temporal sparsity, i.e. they sometimes are silent. Therefore, before making the mix, each source signal is labeled into non-silent / silent sections, and for each MDCT frame (i.e. each time bin t), a Ibit code c 1 is transmitted to provide the combination of AES 42 ND INTERNATIONAL CONFERENCE, Ilmenau, Germany, 2011 July 22-24 non-silent sources. The estimation of Ĩ f t by ( 3) is then carried out for each frequency bin f only among the nonsilent sources, leading to a shorter Ĩ f t -code c 2 compared to the fixed-size code of our previous work. Let us give the following example: if only 4 sources out of 5 are nonsilent on a given frame, the combination of active sources at each frequency bin only concerns 2 sources out of 4, hence 6 combinations instead of 10, and a 3-bit code can be used instead of a 4-bit code for each frequency bin of the considered frame. Note that i) when only two sources are non-silent on a given frame, the analysis of ( 3) is useless, the 2 active sources at each frequency bin being the 2 non-silent sources, ii) when one source is non silent, the inversion process for this source is reduced to dividing either the left or right mixture coefficient by the corresponding mixing coefficient, and iii) obviously when all sources are silent, all estimated source coefficients are set to zero. In all those three cases, c 2 is not transmitted, since the decoding of active sources can be made from c 1 alone, leading to a particularly important bit saving.

Finally, this new improved coding scheme has a noticeable consequence on the quality of the separated sources: in our previous works, silent source signals happened to be estimated as non-silent in TF bin where all sources had very low energy and the "two predominant sources" assumption was not verified (a situation that happens more frequently in the higher region of the spectrum) . This resulted in musical noise for the reconstructed sources in sections where the original sources were silent. With the new coding scheme, silent sections are ensured to be reconstructed as silent sections. This provides more natural reconstructed sources and this is a convenient property for karaoke/soloing applications.

EXPERIMENTS

Data and performance measures

Tests have been processed with 44.1kHz-sampled music signals. We report the results for 3 complete songs made of LISS mixtures of 5 sources each (recorded separately in studio conditions), from different artists and different musical styles (see Fig. 3). The mixture matrices A were set arbitrarily for each piece of music, so that the resulting mixture seems (to us) well musically balanced and the 2 × 2 submatrices are well-conditioned, e.g.: A = 0.93 0.80 0.71 0.60 0.50 0.37 0.60 0.71 0.80 0.87 [START_REF] Smaragdis | Separation by "humming": User-guided sound extraction from monophonic mixtures[END_REF] The quality of separated sources has been assessed by both informal listening tests with high-quality headphones, and log power ratios, as defined in [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. For clarity, we only report here the source-to-distortion ratio (SDR) that provides an overall separation performance criterion (that gathers the influence of interfering sources and artefacts [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]). We also provide in the legends of Fig. 3 typical input signal-to-interference ratios (SIR IN ) to underline that musical sources do not contribute with the same power in a well musically balanced mix (note that those input SIR are averaged on sections where all 5 sources are active to avoid under-representation biais). Hence, gains of the form SDR-SIR IN are good indicators of separation efficiency.

Separation results

Fig. 3 complete suppression of any given source). This has been assessed in our previous studies [7][8]. This result also stands for compressed mix at reasonable bitrates (say, over 128kbps): although the quality of the separated signals from the compressed mix is lower than from the uncompressed mix, this quality is sufficient for individual source enhancement or subtraction. For lower bitrates, a more or less annoying level of musical noise corrupts the estimated sources. In some cases, this noise is correlated with one of the sources, for example it can be correlated with the beat of the drums. The degree of pollution of this noise depends on the considered source signal and on the application. For example, in a remix application (where all sources are remixed with different spatialization and possible addition of audio effects), the noise on a given separated source can be masked by the other remixed sources. In the case of generalized karaoke however, at low compression rates, the noise can prevent the complete suppression of a given source. Demo samples can be downloaded at http://www.gipsa-lab.fr/∼laurent.girin/demo /LISS-AAC-demo.rar.

Side-information bitrates

The distribution of the number of non-silent sources and the corresponding side-information bitrates (average coding cost of c 1 and c 2 ) are provided in Table 1 for the 3 test signals. The additional bitrate is about 100kbps for the Rock and Electro-Jazz examples and it is 77kbps for the Funk piece (revealing its sparser musical arrangement). Even if it not much lower than the mix bitrate, this is quite moderate since it enables to recover 5 source signals from the mix. Note that the side-information bitrate was 128kbps in the basic configuration of our previous works [7][8]. Therefore, the two-step c 1 -c 2 coding leads to about 1/3 bit saving. A deepened investigation of the distribution of code c 2 in the TF plan reveals that this distribution is highly structured: c 2 is often identical for large subblocks of the TF plan, revealing the "clustering" organization of the spectral content of the sources in the TF plane. Therefore, the side-information can be coded even more efficiently. This is part of our future works.

CONCLUSION

The sparse ISS method described in [START_REF] Parvaix | Informed source separation of underdetermined instantaneous stereo mixtures using source index embedding[END_REF][8] has been [START_REF] Marchand | Interactive music with active audio CDs[END_REF] are thus likely to be extended to dematerialized compressed music. The simplicity of the side-information is compensated by an efficient exploitation of signals sparsity and spatial information, as long as the MDCT coefficients are not too much degraded by the compression. This study proves the feasability of a smarter and more computationaly efficient "integrated" ISS-AAC system, where the quantization and separation processes would be carried out on the same MDCT coefficients, as discussed in Section 3.1. In the present study, we did not provide a detailed characterization of the quantization effects on the separation performance (e.g. in different subbands, at different rates, etc.) We plan to do that within the framework of this integrated system, and derive strategies for the control of the bit allocation process by target separation performance. Future extensions of this work will also deal with more complex types of mixture such as binaural/convolutive mixtures, a larger number of sources, and a more efficient coding of the side-information.

Fig. 1 :

 1 Fig. 1: Functional diagram of the ISS system for uncompressed stereo mix.

  provides the average separation SDR obtained with the 3 test signals, for each source signal (averaged over the total duration of the piece of music), and for AAC bitrates ranging from 192 to 64kbps. The results for PCM mix signal are also plotted for comparison (with no watermarking). It can be seen that, as could be expected, the compression of the mix signal leads to a degrada-

	obtained with PCM mix signals enables major appli-
	cations such as remixing or generalized karaoke (i.e.
	tion of the output SDR, and this degradation increases as
	the compression bitrate decreases (since compressed and
	PCM mix signals differ more and more). However, the
	degradation is quite regular and moderate untill 128kbps:
	for example, only 0.5-2dB SDR are lost from PCM to
	192kbps for test signals 1 and 3, depending on the source
	signals, and another 0.5-2dB interval is lost from 192 to
	160kbps for all 3 test signals (for test signal 2, the tran-
	sition from PCM to 192kbps is more abrupt, indicating a
	possible dependancy of the compression effects on mu-
	sical style and source composition; this point will be in-
	vestigated in details in future studies). Since the PCM-
	reference separation performances range within 15-24dB
	gains (SDR-SIR IN ), source separation can be extended
	to compressed source signal at reasonable bitrates, i.e.
	from 192 to 128kbps, with performances that remain
	remarkable, given the low quantity of transmitted side-
	information (see below for additional side-information
	rates). Below 128kbps, the degradation gets worse, but
	this is not surprising, since it is in line with the degrada-
	tion of compressed signal quality as binary resource gets
	sparse (going from 192 to 160kbps is easier than going
	from 96 to 64kbps).
	Those measures are confirmed by informal listening
	tests. Although not perfect, the separation quality

Table 1 :

 1 % of frames with i non-silent sources, i = 1 to 5, and corresponding side-information bitrates. audio signals, with quite robust separation performances at reasonable compression bitrates, and a moderate side-information bitrate (e.g. about 100kbps to separate 5 sources). Specific applications such as active-listening from audio-CD

	Signal	1	2	3	4	5	bitrate (kbps)
	Rock	1.9	0.3 22.3 28.7 46.8 102
	Funk	5.1 11.6 27.5 38.2 17.6 77
	Elec.-Jazz 3.8	3.9 13.2 31.0 48.0 100
	extended to compressed

Disque Repensé pour l'ecoute active de la Musique -This project is supported by the French National Research Agency (ANR) -Grant CONTINT

CORD 006