Manuel Moussallam
email: manuel.moussallam@telecom-paristech.fr

Laurent Daudet

Gaël Richard

Matching Pursuits with Random Sequential Subdictionaries

Keywords: Pursuits, Random Dictionaries, Sparse Approximation, Audio Signal Compression

Matching pursuits are a class of greedy algorithms commonly used in signal processing, for solving the sparse approximation problem. They rely on an atom selection step that requires the calculation of numerous projections, which can be computationally costly for large dictionaries and burdens their competitiveness in coding applications. We propose using a non adaptive random sequence of subdictionaries in the decomposition process, thus parsing a large dictionary in a probabilistic fashion with no additional projection cost nor parameter estimation. A theoretical modeling based on order statistics is provided, along with experimental evidence showing that the novel algorithm can be efficiently used on sparse approximation problems. An application to audio signal compression with multiscale time-frequency dictionaries is presented, along with a discussion of the complexity and practical implementations.

Introduction

The ability to describe a complex process as a combination of few simpler ones is often critical in engineering. Signal processing is no exception to this paradigm, also known in this context as the sparse representation problem. Yet rather simple to expose, its combinatorial nature has led researchers to develop so many bypassing strategies that it can be considered a self-sufficient topic. Throughout these years of work, many additional benefits were found arising from the dimensionality reduction. Firstly, sparsity allows faster processing, and biologically inspired models [START_REF] Smith | Efficient auditory coding[END_REF] suggest that the mammalian brain takes advantage of it. Secondly, it helps reducing both storage and broadcasting costs [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]. Finally, it enables semantic characterization, since few significant objects carry most of the signals information.

The central problem in sparse approximation theory may be written as such: Given a signal f in a Hilbert space H and a finite-size dictionary Φ = {φ γ } of M unit norm vectors (∀γ ∈ [1..M ] , φ γ 2 = 1) in H, called atoms, find the smallest expansion of f in Φ up to a reconstruction error ǫ :

min α 0 s.t. f - M γ=1 α γ φ γ 2 ≤ ǫ (1)
where α 0 is the number of non-zero elements in the sequence of weights {α γ }. Equivalently, one seeks the subset of indexes Γ m = {γ n } n=1..m of atoms of Φ and the corresponding m non zero weights {α γn } solving:

min m s.t. f -f m 2 ≤ ǫ (2)
where f m = m n=1 α γn φ γn is called a m-term approximant of f . A corollary problem is also defined, the sparse recovery problem, where one assumes that f has a sparse support Γ m in Φ (m being much smaller than the space dimension) and tries to recover it (from noisy or fewer measurement than the signals dimension). In the sparse approximation problem, no such assumption is made and one just seeks to retrieve the best m-term approximant of f in the sense of minimizing the reconstruction error (or any adapted divergence measure). Many practical problems benefit from this modeling, from denoising [START_REF] Fevotte | Denoising of musical audio using sparse linear regression and structured priors[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] to feature extraction and signal compression [START_REF] Figueras I Ventura | Low-rate and flexible image coding with redundant representations[END_REF][START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF].

Unfortunately, problem (1) is NP-hard and an alternate strategy needs to be adopted. A first class of existing methods are based on a relaxation of the sparsity constraint using the l 1 -norm instead of the nonconvex l 0 [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], thus solving a quadratic program. Alternatively, greedy algorithms can produce (potentially suboptimal) solutions to problem (1) by iteratively selecting atoms in Φ.

Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and its variants [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF][START_REF] Jaggi | High resolution pursuit for feature extraction[END_REF][START_REF] Daudet | Sparse and structured decompositions of signals with the molecular matching pursuit[END_REF][START_REF] Christensen | The cyclic matching pursuit and its application to audio modeling and coding[END_REF][START_REF] Blumensath | Gradient pursuits[END_REF] are instances of such greedy algorithms. A comparative study of many existing algorithms can be found in [START_REF] Dymarski | Greedy sparse decompositions: A comparative study[END_REF]. MP-like algorithms have simple underlying principles allowing intuitive understanding, and efficient implementations can be found [START_REF] Krstulovic | Mptk: matching pursuit made tractable[END_REF][START_REF] Mailhe | A low complexity orthogonal matching pursuit for sparse signal approximation with shift-invariant dictionaries[END_REF]. Such algorithms construct an approximant f n after n iterations by alternating two steps:

Step 1 Select an atom φ γn ∈ Φ and append it to the support Γ n = Γ n-1 γ n

Step 2 Update the approximation f n accordingly. until a stopping criterion is met or a perfect decomposition is found (i.e f n = f ). This two-step generic description defines the class of General MP algorithms [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasiincoherent dictionaries[END_REF].

When designing the dictionary, it is important that atoms locally ressemble the signal that is to be represented. This correlation between the signal and the dictionary atoms ensures that greedy algorithms select atoms that remove a lot of energy from the residual). For some classes of signals (e.g. audio), the variety of encountered waveforms implies that a large number of atoms must be considered. Eventually, one is interested in finding sparse expansions of signals in large dictionaries at reasonable complexity.

Methods addressing sparse recovery problems have benefited from random dictionary design properties (e.g with MP-like algorithms [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF][START_REF] Blumensath | Gradient pursuits[END_REF]). Randomness in this context is used for spreading information that is localized (on a few non-zero coefficients in a large dictionary) among fewer measurement vectors.

In the sparse approximation context though, the use for randomness is less obvious. Much effort has been done on designing structured dictionaries that exhibit good convergence pattern at reasonable costs [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF][START_REF] Daudet | Sparse and structured decompositions of signals with the molecular matching pursuit[END_REF][START_REF] Jost | Tree-based pursuit: Algorithm and properties[END_REF]. Typical dictionaries with limited size are based on Time-Frequency transforms such as window Fourier transforms (also referred to as Gabor dictionaries) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF][START_REF] Gribonval | Analysis of sound signals with high resolution matching pursuit[END_REF], Discrete Cosine [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF], wavelets [START_REF] Smith | Efficient auditory coding[END_REF] or unions of both [START_REF] Daudet | Sparse and structured decompositions of signals with the molecular matching pursuit[END_REF]. Such dictionaries are made of atoms that are localized in the time-frequency plance, which gives them enough flexibility to represent complex non stationary signals. However the set of considered localization reflects, in practice, a compromise.

Considering all possible localizations is, in theory, possible for discrete signals and dictionaries, however it yields very large and redundant dictionaries which raises computational issues. The standard dictionary design [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] considers a fixed subset of localizations is arguably a good option but it can never be optimal for all possible signals. Adapting the subset to the signal is then possible, but it introduces additional complexity. Finally, one may consider using multiple random subsets, perform multiple decompositions and average them in a post-processing step [START_REF] Ferrando | Probabilistic matching pursuit with gabor dictionaries[END_REF].

In this paper, an other option is considered: the use of randomly varying subdictionaries at each iteration during a single decomposition. The idea is to avoid the additional complexity of adapting the dictionary to the signal, or of having to perform multiple decompositions, while still improving on the fixed dictionary strategy.

A parallel can be traced with quantization theory. Somehow, limiting dictionary size amounts to discretizing its atoms parameter space. This quantization introduces noise in the decomposition (i.e error due to dictionary-related model mismatch). Adaptive quantization reduces the noise and multiple quantization allow to average the noise out. But the proposed technique is closer to the dithering technique [START_REF] Zamir | On universal quantization by randomized uniform/lattice quantizers[END_REF] and more generally to the stochastic resonance theory, which shows how a moderate amount of added noise can increase the behavior of many non linear systems [START_REF] Gammaitoni | Stochastic resonance[END_REF]. The proposed method only relies on a slight modification of step 1 and can therefore be applied to many pursuit algorithms. In this work, we focus on the approximation problem. Experiments on real audio data demonstrate the potential benefits in signal compression and particularly in low-bit-rate audio coding. It should be emphasized that, as such, it is not well suited for sparse recovery problems such as in the compressed sensing scheme [START_REF] Donoho | Compressed sensing[END_REF]. Indeed, searching in random subsets of Φ may burden the algorithm ability to retrieve true true sparse support of a signal in Φ.

The rest of this paper is organized as follows. Section 2 recalls the standard Matching Pursuit algorithm and some of its variants. Then, we introduce the proposed modification of Step 1 in Section 3 and expose some theoretical justifications based on a probabilistic modeling of the decomposition process. A practical application is presented in Section 4 as a proof of concept with audio signals. Extended discussion on various aspects of the new method is provided in Section 5.

Matching Pursuit for Sparse Signal Representations

Matching Pursuit Framework

The Matching Pursuit (MP) algorithm [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and its variants in the General MP family [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasiincoherent dictionaries[END_REF] all share the same underlying iterative two-step structure, as described by Figure 1. Specific algorithms differ in the way they perform the atom selection criteria C (Step 1), and the approximant construction A (Step 2). Standard MP as originally defined in [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] is given by:

C(Φ, R n f ) = arg max φγ ∈Φ | R n f, φ γ | (3) A(f, Φ Γ n ) = n-1 i=0 R i f, φ γi φ γi (4)
where R n f denotes the residual at iteration n (i.e R n f = f -f n ). Orthogonal Matching Pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] is based on the same selection criteria, but the approximant update step ensures orthogonality between the residual and the subspace spanned by the already selected atoms

V Γ n = span{φ γn } γn∈Γ n A OMP (f, Φ Γ n ) = P V Γ n f (5)
where P V is the orthonormal projector onto V.

The algorithm stops when a predefined criterion is met, either an approximation threshold in the form of a Signal-to-Residual Ratio (SRR):

SRR(n) = 10 log f n 2 f -f n 2 (6)
or a bit budget (equivalently a number of iterations) in coding applications.

Weak MP

Often, the size of Φ prevents the use of the criterion C as defined in 3, especially for infinite-size dictionaries. Instead, one can satisfy a weak selection rule C weak (Φ, R n f ) = φ γ weak such that:

| R n f, φ γ weak | ≥ t n sup φγ ∈Φ | R n f, φ γ | (7)
with 0 < t n ≤ 1. In practice, a weak selection step is implemented by limiting the number of atomic projections in which the maximum is searched. This is equivalently seen as a subsampling of the large dictionary Φ. Temlyakov [START_REF] Temlyakov | A criterion for convergence of weak greedy algorithms[END_REF] proved the convergence when tn n < ∞ and stability properties have been studied by Gribonval et al in [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasiincoherent dictionaries[END_REF].

Let us denote Φ = {φ γ } a large multiscale time-frequency dictionary with atom parameter γ living in S × U × Ξ a finite subset of R + × R 2 . S denotes the set of analysis scales, U is the set of time indexes and Ξ the set of frequency indexes (see for example [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]). Computing C(Φ, R n-1 f ) requires, in the general case, the computation of all atom projections R n-1 f, φ γ which can be prohibitive. A weak strategy consists in computing only a subset of these projections, thus using a coarse subdictionary Φ I ⊂ Φ, which can be seen as equivalent to a subsampling of the parameter space S × U × Ξ. For instance, time indexes can be limited to fractions of the analysis window size.

Reducing the size of the dictionary is interesting for coding purposes. Actually, the smaller the atom parameter space, the cheaper the encoding of atom indexes. Figure 2 shows an example of such dictionary subsampling in a time-frequency plane. Each point represents the centroid of a time-frequency atom φ γ , black circles corresponds to a coarse subsampling of the parameter space and constitute the subdictionary Φ I 1 ⊂ Φ. The full dictionary Φ figured by the small red crosses, is 16 times bigger than Φ I 1 .

However, the choice of subsampling has consequences on the decomposition. Figure 2 pictures an example of a different coarse subdictionary defined by another index subset Φ I 2 ⊂ Φ. For most signals, the decomposition will be different when using Φ I 1 or Φ I 2 . Moreover, there is no easy way to guess which subdictionary would provide the fastest and/or most significant decomposition. Using one of these subdictionaries instead of Φ implies that available atoms may be less suited to locally fit the signals. MP would thus create energy in the time frequency plane where there is none in the signal, an artifact known as dark energy [START_REF] Sturm | Dark energy in sparse atomic decompositions[END_REF].

Higher Resolution and Locally-Adaptive Matching Pursuits

To tackle this issue, high resolution methods have been studied. First, a modification of the atom selection criterion has been proposed by Jaggi et al [START_REF] Jaggi | High resolution pursuit for feature extraction[END_REF] that not only takes energy into account but also local fit to the signal. Other artifact preventing methods have also been introduced (e.g for pre-echo [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF] or dark energy [START_REF] Sturm | Sparse approximation and the pursuit of meaningful signal models with interference adaptation[END_REF]). A second class of bypassing strategies are locally-adaptive pursuits. An atom is first selected in the coarse subdictionary Φ I , then a local optimization is performed to find a neighboring maximum in Φ. Mallat et al proposed a Newton method [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], Goodwin et al focused on phase tuning [START_REF] Goodwin | Matching pursuit and atomic signal models based on recursive filter banks[END_REF], Gribonval tuned a chirp parameter [START_REF]Fast matching pursuit with a multiscale dictionary of gaussian chirps[END_REF] and Christensen et al addressed both phase and frequency [START_REF] Christensen | The cyclic matching pursuit and its application to audio modeling and coding[END_REF][START_REF] Sturm | Cyclic matching pursuits with multiscale time-frequency dictionaries[END_REF]. An equivalent strategy is defined for gammachirps in [START_REF] Lewicki | Coding time-varying signals using sparse, shift-invariant representations[END_REF]. Locally-adaptive methods in time-frequency dictionaries implement a form of time and / or frequency reassignment of an atom selected in a coarse subdictionary.

These strategies yield representations in the large dictionary at a reduced cost. In many applications, the faster residuals energy decay compensates the slight computational overhead of the local optimization. In a previous work [START_REF] Moussallam | Audio signal representation for factorisation in the sparse domain[END_REF], we emphasized the usefulness of locally-adaptive methods for shift-invariant representation and similarity detection in the transform domain. Though, the resulting atoms are from the large dictionary, i.e their indexes are more costly to encode. Moreover, these method require an additional parameter estimation step that may increase the overall complexity.

Statistics and MP

Different approaches have been proposed to enhance MP algorithms using statistics. Ferrando et al [START_REF] Ferrando | Probabilistic matching pursuit with gabor dictionaries[END_REF] were (to our knowledge) the first to propose to run multiple sub-optimal pursuits and retrieve meaningful atoms in an a posteriori averaging step. A somehow similar approach is introduced in [START_REF] Durka | Stochastic time-frequency dictionaries for matching pursuit[END_REF] where the authors emphasize a lack of precision of their representation due to the structure of the applied dictionary. To avoid these decomposition artifacts, they follow a Monte-Carlo like method, where the set of atom parameter is randomly chosen before each decomposition followed by an averaging step.

Elad et al [START_REF] Elad | A plurality of sparse representations is better than the sparsest one alone[END_REF] used the same paradigm for support recovery, taking advantage of many suboptimal representations instead of a single one of higher precision. Similar approaches are adopted within a Bayesian framework in [START_REF] Schniter | Fast bayesian matching pursuit: model uncertainty and parameter estimation for sparse linear models[END_REF][START_REF] Zayyani | Bayesian pursuit algorithm for sparse representation[END_REF], although the subset selection is a byproduct of the choice of the prior. Their work is related to the branch of compressed sensing that uses pursuits and random measurement matrices to retrieve sparse supports of various classes of high dimensional signals [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit[END_REF]. A recent work by Divekar et al [START_REF] Divekar | Probabilistic matching pursuit for compressive sensing[END_REF] proposed a form of probabilistic pursuit. Such work differs from ours in the sense that it is signal-adaptive and tuned for the support recovery problem.

Matching Pursuits with Sequences of Subdictionaries

Proposed new algorithm

This paper presents a modification of MP which consists in changing the subdictionary at each iteration. Let us define a sequence I = {I n } n∈[0..m-1] of length m. Each I n is a set of indexes of atoms from Φ. At iteration n, the algorithm can only select an atom in the subdictionary Φ I n ⊂ Φ. We call such algorithms Sequential Subdictionaries Matching Pursuits (SSMP).

This method is illustrated by Figure 3. When the sequence I is random, we call such algorithms Random Sequential Subdictionaries (RSS) Pursuits. The proposed modification can be applied to any algorithm from the General MP family (RSS MP, RSS OMP, etc..). After m iterations, the algorithm outputs an approximant of the form ( 8)

f m = m-1 n=0 α n φ I n γn (8)
where φ I n γn ∈ Φ I n . SSMP can be seen as a Weak General MP algorithm as described in [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasiincoherent dictionaries[END_REF]. Indeed, it amounts to reducing the atom selection choice to a subset of Φ, which defines the weak selection rule [START_REF] Divekar | Probabilistic matching pursuit for compressive sensing[END_REF]. The originality of this work is the fact that we change the subdictionary at each iteration while it is usually fixed during the whole decomposition. The sequence of subdictionaries is known in advance and is thus signal independent.

The sequence I might be built such that lim m→∞ m-1 n=0 Φ I n = Φ, although in this case, the modified pursuit is not equivalent to a pursuit using the large dictionary Φ. However, it allows the selection of atoms from Φ that were not in the initial subdictionary Φ I 0 . Step 1 is modified so as to have the search take place in a different subdictionary at each iteration. The dictionary subsampling is controlled by a fixed pre-defined sequence I.

A first example

Let Φ be a Gabor dictionary. Each atom in Φ is a Gabor atom that can be parameterized by a triplet (s, u, ξ) ∈ S × U × Ξ defining its scale, time and frequency localization respectively. The continuous-time version of such atom is given by:

φ s,u,ξ (t) = 1 √ s g t -u s e iξ(t-u) (9)
where g can be a Gaussian or Hann window. A discrete Gabor dictionary is obtained by sampling the parameter space S × U × Ξ (and the window g). Let N be the signal dimension and K be the number of Gabor scales (s 1 , s 2 , .., s K ). For each scale s k , u takes integer values between 0 and N and ξ from 0 to s k /2. The total size of Φ in this setup is

M = K k=1 s k N 2 atoms. Let f ∈ R N be
a m-sparse signal in Φ, meaning that f is the sum of m (randomly chosen) components from the full dictionary Φ. We are interested in minimizing the reconstruction error ǫ = f -fn 2 f 2 given the size constraint n = m/2.

Let Φ I 0 be a subdictionary of Φ defined by selecting in each scale s k only atoms at certain time localizations. Subsampling the time axis by half the window length is standard in audio processing (i.e. ∀u ∃p ∈ [0, ⌊2N/s k ⌋ -1], s.t. u = p.s k /2). Φ I 0 can be seen as a subsampling of Φ. The size of Φ I 0 in this setup is down to K × N atoms. MP using Φ I 0 during the whole decomposition is a special case of SSMP where all the subdictionaries are the same. We label this algorithm Coarse MP. Now let us define a pseudo-random sequence of subdictionaries Φ I n that are translated versions of Φ I 0 . Each Φ I n has the same size than Φ I 0 but its atoms are parameterized as such:

∀u ∃p ∈ [0, ⌊2N/s k ⌋ - 1], s.t. u = p.s k /2 + τ n k where τ n k ∈ [-s k /4, s k /4
] is a translation parameter different for each scale s k of each subdictionary Φ I n . Each Φ I n can also be seen as a different subsampling of Φ. The size of each subdictionary Φ I n is thus also K × N . MP using the sequence of subdictionaries Φ I n is labeled Random SSMP (RSS MP).

Finally, a MP using the full dictionary Φ during the whole process is labeled Full MP. We compare the reconstruction error achieved with these three algorithms (i.e Coarse MP, RSS MP and Full MP).

Figure 4 gives corresponding results with the following setting: N = 10000, m = 300, S = [32, 128, 512] and thus M = 840000, averaged over 1000 runs. Although RSS MP is constrained at each iteration to choose only within a limited subset of atoms, it exhibit an error decay close to the one of Full MP.

If f was exactly sparse in Φ I 0 , then keeping the fixed subdictionary (Coarse MP) would provide a faster decay than using the random sequence (RSS MP). However in our setup, such case is unlikely to happen. As stressed by many authors cited above, the choice of a fixed subdictionary is often suboptimal with respect to a whole class of signals, as it may not have the basic invariants one may expect from a good representation (for instance, the shift invariance). A matter of concern is now to determine when RSS MP performs better than Coarse MP in terms of approximation error decay.

In the next subsection we make use of order statistics to model the behavior of the two strategies (i.e. MP with fixed coarse subdictionary and MP with varying subdictionaries) depending on the initial distribution of projections.

Order Statistics

We need to introduce a few tools from the order statistics theory. The interested reader can refer for example to [START_REF] Nagaraja | Order statistics from discrete distributions[END_REF][START_REF] Gupta | On the order statistics from equally correlated normal random variables[END_REF] for more details. Let z 1 , z 2 , ..z n be n i.i.d samples drawn from a continuous random variable Z with probability density f Z and a distribution function F Z . Let us denote Z 1:n , Z 2:n , .., Z n:n the order statistics. The random variable Z i:n represents the i th smallest element of the n samples. The probability density of Z i:n is denoted f Z i:n and is given by:

f Z i:n (z) = n! (n -i)!(i -1)! F Z (z) i-1 f Z (z)(1 -F Z (z)) n-i (10)
The density of extremum values is easily derived from this equation, in particular the maximum value of the sequence has density:

f Z n:n (z) = nF Z (z) n-1 f Z (z) (11)
The moments of the order statistics are given by the formula:

µ (m) i:n = E(Z m i:n ) = ∞ -∞ z m f Z i:n (z)dz (12)
And the variance is denoted σ 2 i:n . For convenience we will write the expectation

µ i:n = µ (1)
i:n .

Order statistics in an MP framework

Given the greedy nature of MP, a statistical modeling of the series of projection maxima gives meaningful insights about the overall convergence properties. Let f be a signal in R N . At the first iteration, the projection of the residual R 0 f = f over a complete dictionary Φ of M unit normed atoms

{φ i } i∈[1,M] (M > N ) is given by: ∀i ∈ [1..M ], α i = R 0 f, φ i (13)
Let us denote z i = |α i | and consider them as M i.i.d samples drawn from a random variable Z living in [0, f 2 ]. Note that such an assumption holds only if one considers that atoms in Φ are almost pairwise orthogonals, i.e that Φ is quasi-incoherent. MP selects the atom φ γ0 = arg max z i whose weights absolute value is given by the M th order statistics of Z: Z M:M .

|α γ0 | = max φγ ∈Φ | R 0 f, φ γ | = Z M:M (14)
Standard MP constructs a residual by removing this atoms contribution and iterating. Hence, weights of selected atoms remains independent and at iteration n the M -n-th order statistics of Z describes the selected weight:

|α γn | = max φγ ∈Φ | R n f, φ γ | = Z M-n:M (15)
the energy conservation in MP allows to derive [START_REF] Gammaitoni | Stochastic resonance[END_REF].

f 2 = R n f 2 + n-1 i=0 |α γi | 2 (16)
Combining ( 16) and ( 15) and taking the expectation, one gets ( 17)

E( R n f 2 ) = f 2 - n-1 i=0 µ (2) M-i:M (17)
We recognize the second order moment of µ

(2)

M-i:M = σ 2 M-i:M + µ 2 M-i:M .
Similarly, from ( 16) we can derive the variance of the estimator:

V ar( R n f 2 ) = n-1 i=0 n-1 j=0 cov(Z 2 M-i:M , Z 2 M-j:M ) (18)
Given a distribution model for Z, Equations ( 17) and ( 18) provide mean and variance estimates of the residuals energy decay with a standard MP on a quasi-incoherent dictionary Φ.

Redrawing projection coefficients: changing the dictionaries

The idea at the core of the new algorithm described in section 3.1 is to draw a new set of projections at each iteration by changing the dictionary in a controlled manner. Let us now assume that we know a sequence of complete quasi-incoherent dictionaries {Φ i } i∈[0,n] . Let us make the additional assumptions that the projection of f has the same distribution in all these dictionaries.

At the first iteration, the process is similar, so the atom φ γ0 in Φ 0 is selected as arg max φγ ∈Φ0 | R 0 f, φ γ | with the same weight as [START_REF] Fevotte | Denoising of musical audio using sparse linear regression and structured priors[END_REF]. After subtracting the atom we have a new residual R 1 f . The trick is now to search for the most correlated atom in Φ 1 . The absolute values of the corresponding inner products define a new set of M i.i.d samples

z 1 i : ∀i ∈ [1..M ], z 1 i = |α i | = | R 1 f, φ i | (19)
Let us assume that the z 1 i samples have the same distribution law than the z i , except that the sample space is now smaller because R 1 f has less energy than f . This means the new random variable Z1 is distributed like Z × R 1 f f . Let us denote Z1 M:M the M -th order statistic for Z1, its second moment is:

E(Z1 2 M:M ) = E(Z 2 M:M ).E( R 1 f 2 f 2 ) + cov(Z 2 M:M , R 1 f 2 f 2 ) ( 20
)
and we know that

R 1 f 2 = f 2 -Z 2 M:M hence: cov(Z 2 M:M , R 1 f 2 ) = -cov(Z 2 M:M , Z 2 M:M ) = (µ (2)
M:M ) 2 -µ (4) M:M (21)
which leads to:

E(Z1 2 M:M ) = µ
(2)

M:M .(1 - µ (2) M:M f 2 ) + 1 f 2 (µ (2)
M:M ) 2 -µ (4) M:M = µ
(2)

M:M - µ (4) M:M f 2 Since the new residual R 2 f has energy R 2 f 2 = R 1 f 2 -Z1 2 M:M , taking the expectation: E( R 2 f 2 ) = f 2 -2µ
(2)

M:M + µ (4) M:M f 2 (22)
and we see higher moments of the M -th order statistic appear in the residuals energy decay formula. At the n-th iteration of this modified pursuit, we end up with a rather simple formula [START_REF] Jost | Tree-based pursuit: Algorithm and properties[END_REF] where higher moments of the M -th order statistic of Z appear.

E( R n f 2 ) = f 2 + n i=1 (-1) j n i µ (2i) M:M f 2(i-1) (23)
Similarly, a variance estimator can be derived that only depends on the moments of the highest order statistic:

V ar( R n f 2 ) = n i=1 n j=1 (-1) j+i n i n j cov(Z (2i) M:M , Z (2j)
M:M ) f 2(i+j-1) (24)

Comparison with simple models

Let us now compare the behavior of the new samples strategy to the one where they are fixed. The relative error ǫ(n

) = 10 log R n f 2 f 2
after n iterations is computed for both strategies. Knowing the probability density function (pdf) f Z , Equations ( 17) to [START_REF] Krstulovic | Mptk: matching pursuit made tractable[END_REF] provide closed-form mean and variance estimate of ǫ(n) for both strategies. For example if Z ∼ U(0, 1), its order statistics follow a beta distribution:

Z k:M ∼ Beta(k, M + 1 -k) and all moments µ (m)
M-k:M can be easily found. A graphical illustration of potential of the new strategy is given by Figure 5. For three different distributions of Z, namely uniform, normal (actually half normal since no negative values are considered) and exponential. The density function of two order statistics of interest (i.e the maximum f Z M:M , and the M/2-th element f Z M/2:M ) is plotted. These two elements combined provide an insight on how fast the expected value of selected weight will decline using the fixed strategy. Some comments can be made:

• In the uniform case: one might still expect to select relatively large coefficients (i.e good atoms) after M/2 iterations with the standard strategy. In this case, redrawing the samples (i.e changing the subdictionary) appears to be detrimental to the error decay rate.

• For normally and exponentially distributed samples, the expected value of a coefficient selected in a fixed sequence drops more quickly towards small values (i.e there are relatively fewer large values).

Relative error profiles indicate that selecting the maximum of redrawn samples can prove beneficial in terms of minimizing the reconstruction error. It is promising for sparse approximation problems such as lossy compression.

In practice, the uniform model is not well suited for sparse approximation problems. It basically indicates that the chosen dictionary is poorly correlated with the signal (e.g. white noise in Dirac dictionaries). Exponential and normal distributions, on the opposite, are closer to practical situations where the signal is assumed sparse in the dictionary, with additional measurement noise such as the one presented in section 3.2. In these cases, changing the subdictionaries seems to be beneficial.

Let us again stress that we have only presented a model for error decays with two strategies and a set of strong assumptions. It is not a formal proof of a guaranteed faster convergence. However, it gives insight on when the proposed new algorithm may be useful and when it may not, along with estimators of the decay rates when a statistical modeling of signal projections on subdictionaries is available. To our knowledge, such modeling was not proposed before for MP decompositions.

Matching Pursuit with Random Time Frequency Subdictionaries for audio compression

Audio compression with union of MDCT dictionaries

Sparse decompositions using Matching Pursuit have been proven by Ravelli et al [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF] to be competitive for low bit-rate audio compression. In spite of Gabor dictionaries, they use the Modified Discrete Cosine Transform (MDCT) [START_REF] Princen | Analysis/synthesis filter bank design based on time domain aliasing cancellation[END_REF] that is at the core of most transform coders (e.g. MPEG1-Layer III [START_REF]Coding of moving pictures and associated audio for digital, storage at speed up to about 1[END_REF], MPEG 2/4 AAC). Let us recall the main advantage of RSSMP compared to Coarse MP: reconstruction error can be made smaller with the same number of iterations. If the random sequence of subdictionaries is already known by both coder and decoder then the cost of encoding each atom is the same in both cases.

For this proof of concept, a simplistic coding scheme is used as presented in Figure 6. For each atom in the approximant, the index is transmitted alongside its quantized amplitude (using a uniform mid-tread quantizer). The cost of encoding an atom index is log 2 (M ) bits where M is the size of the dictionary in which the atom was chosen. A quantized approximation fn of f is obtained after n iterations with a characteristic SN R (same as SRR but computed on fn instead of f n ) and associated bit-rate.

The signals are taken from the MPEG Sound Quality Assesment Material (SQAM) dataset, their sampling frequency is reduced to 32000 Hz. Three cases are compared in this study:

Coarse MP Matching Pursuit with a union of MDCT basis with no additional parameter. The MDCT basis have sizes from 4 to 512 ms (8 dyadic scales from 2 7 to 2 14 samples) hop size is 50% in each scale.

LoMP Locally Optimized Matching Pursuit (LoMP) with a union of MDCT Basis. This is a locallyadaptive pursuit where an atom is first selected in the coarse dictionary. Then its time localization is optimized in a neighbourhood to find the best local atom in the full dictionary. This pursuit is a slightly suboptimal equivalent of a Matching Pursuit using the full dictionary (Full MP). It is nevertheless less computationally intensive. However, it requires the computation (and the transmission) of an additional parameter (the local time shift) per atom (see [START_REF] Moussallam | Audio signal representation for factorisation in the sparse domain[END_REF] for more details).

RSS MP MP with a random sequence of subdictionaries. Each subdictionary is itself a union of MDCT basis, each basis of scale s k is randomly translated in time by a parameter

τ i k ∈ [-s k /4 : s k /4]. It is worth noticing that each scale is independently shifted. The sequence T = {τ i k } k=1..K,i=1.
.n is known in advance at both encoder and decoder side, i.e. there is no need to transmit it.

Step SSMP and Jump SSMP are two deterministic SSMP algorithms that are further described in Section 5.

Sparsity results

Figure 7 shows the number of iterations needed with the three described algorithm (and two variants discussed in 5.1) to decompose short audio signals from the MPEG SQAM dataset [2] (a glockenspiel, an orchestra, a male voice , a solo trumpet and a singing voice), at 10 dB of SRR. One can verify that RSS MP yields sparser representations (fewer atoms atoms are needed to reach a given SRR) than Coarse MP. This statement remains valid at any given SRR level in this setup.

Experiments where run 100 times with the audio signal being randomly translated in time at each run. Figure 7 shows that empirical variance remains low, which gives us confidence in the fact that RSS MP will give sparser representations in most cases. So far we have not found a natural audio example contradicting this.

Compression results

Figure 7 shows that the locally-adaptive algorithm (LoMP) is better in terms of sparsity of the achieved representation than RSS MP. However, each atom in these representations is more expensive to encode. Actually, one must transmit an additional local time shift parameter. With RSS MP no such parameter needs to be transmitted. This gives RSS MP a decisive advantage over the two other algorithms for audio compression.

As an example, let us examine the case of the last audio example labeled vega. This is the first seconds from Susan Vega's song Tom's Dinner. The three algorithms are run with a union of 3 MDCT scales (atoms have lengths of 4, 32 or 256 ms). Again these scores have been averaged over 100 runs with the signal being randomly offset at each run and showed very little empirical variance. In order to reach e.g. 20 dB of SRR (i.e before weight quantization):

• Coarse MP selected 6886 atoms. If each has a fixed index coding cost of 18 bits and weight coding cost of 16 bits, a total cost of 231 kBits would be needed.

• LoMP selected 3691 atoms. With same index and weight cost plus an additional time shift parameter whose cost depends on the length of the selected atom, the total cost would be 149 kBits.

• RSS MP selected 3759 atoms (with empirical standard deviation of 16 atoms). Each atom has the same fixed cost as in the Coarse MP case. The average total size is thus 126 kBits.

Using an entropy coder as in Figure 6 does not fundamentally change these results. Figure 8 summarizes compression results with Coarse MP, LoMP and RSS MP for various audio signals from this dataset. Although these comparisons do not use the most efficient quantization and coding tools, it is clear from this picture that the proposed randomization paradigm can be efficiently used in audio compression with greedy algorithms. In this proof of concept, the quality measure adopted was the SN R. However, such algorithms and dictionaries are known to introduce disturbing ringing artifacts and pre-echo. We have not experienced that the proposed algorithm increased nor reduced these artifacts compared to the other pursuits. Furthermore, pre-echo control and dark energy management techniques can be applied with this algorithm just as with any other. Audio examples are available online 1 .

Discussion

Having presented the novel algorithm and an application to audio coding, further questions are investigated in this section. First an experimental study is conducted to illustrate why the sequence of subdictionaries should be pseudo-randomly generated. Then, the RSS scheme is applied to orthogonal pursuit and finally a complexity study is provided.

Random vs Deterministic

In the above experiments we have constructed a sequence of randomly varying subdictionaries by designing a sequence of time shifts T = {τ i k } k=1..K,i=1..n by means of a pseudo-random generator with a uniform distribution:

∀(k, i)τ i k ∼ U(-s k /4, s k /4
). There are other ways to construct such a sequence, some of which are deterministic. However, we have found that this pseudo-random setup is interestingly the one that gives the best performances.

To illustrate this, let us recall the setup of Section 4 and compare the sparsity of representations achieved with sequences of subdictionaries built in the following manner:

RANDOM: T is randomly chosen: ∀(k, i)τ i k ∼ U(-s k /4, s k /4)
STEP: T is constructed with stepwise increasing sequences :

∀(k, i)τ i k = mod(τ i-1 k + 1, s k /2
) and ∀k, τ 0 k = 0. This sequence yields dictionaries that are quite close from one iteration to the next.

JUMP: T is constructed with jumps : ∀(k, i)τ i k = mod(τ i-1 k + s k /4 -1, s k /2
)and ∀k, τ 0 k = 0. This sequence yields dictionaries that are very dissimilar from one iteration to the next, but are quite close to dictionaries 2 iterations later.

Figure 7 shows that among all cases, the random strategy is the best one in terms of sparsity. The worst strategy is STEP. This can be explained by the fact that, from one iteration to the next, the subdictionaries are well correlated. Indeed all the analysis windows are shifted simultaneously by the same (little) offset. The JUMP strategy is already better. Here the different analysis scales are shifted in different ways. Low correlation between successive subdictionaries appears to be an important factor. Experiments were run 100 times and the hierarchical pattern remain unchanged for higher SRR levels.

Many more deterministic strategies have been tried, the random strategy remains the most interesting one when no further assumption about the signal is made. We explain this phenomenon by noticing that pseudo random sequences are more likely to yield subdictionaries that are very uncorrelated one to another not only on a iteration-by-iteration basis but also for a large number of iterations.

Orthogonal pursuits

Other SSMP family members can be created with a simple modification of their atom selection rule. In particular, one may want to apply this technique to Orthogonal Matching Pursuit. The resulting algorithm would then be called RSS OMP.

To evaluate RSS OMP, we have recreated the toy experiment of [START_REF] Blumensath | Gradient pursuits[END_REF] and use 1000 random dictionaries of size 128×256 with atoms drawn uniformly from the unit sphere. Then, 1000 random signals are created using 64 elements from the dictionaries with unit variance zero mean Gaussian amplitudes. For each signal and corresponding dictionary, the decomposition is performed using MP and OMP with a fixed subdictionary of size 128 × 64 (Coarse MP and Coarse OMP), a random sequence of subdictionaries of size 128 × 64 (RSS MP and RSS OMP), and the complete dictionary (Full MP and Full OMP). Matlab/Octave code to reproduce this experiment is available online.

Averaged results are shown Figure 9. Fixed subdictionary cases (Coarse MP/OMP) show a saturated pattern caused by incompleteness of the subdictionary. Although working on subdictionaries four times smaller, pursuits on random sequential subdictionaries exhibit a behavior close to pursuits on the complete

(M N ) O(LN ) O(M N ) O(LN ) O(LN ) O(LN ) Selection O(M ) O(L) O(M ) O(L) O(L) O(L) Step 2 Gram Matrix 0 0 O(iN ) O(iN ) 0 O(iN ) Weights 0 0 O(i 2 ) O(i 2 ) 0 O(i 2 ) Residual O(N ) O(N ) O(iN ) O(iN ) O(N ) O(iN ) Total O(M N ) O(LN ) O(i 2 + M N ) O(i 2 + LN ) O(LN ) O(i 2 + LN )
Table 1: Complexities of the different steps for various algorithm in the general case (no update trick for projection). N is the signal dimension. M is the number of atoms in Φ and L the number of atoms in subdictionaries. The iteration number i is always smaller than N .

dictionary. Potential gains are even more visible with the OMP algorithm, where the probabilistic parsing of the large dictionary allows for a much better error decay rate than the fixed subdictionary case (Coarse OMP).

Computational Complexities 5.3.1. Iteration complexities in the general case

Let f be in R N and Φ be a redundant dictionary of M atoms also in R N . Let Φ I be a subdictionary of Φ of L atoms. Complexities in the general case are given at iteration i by Table 1. When no update trick is used, Step 1 requires the projection of the N dimensional residual onto the dictionary followed by the selection of the maximum index. Step 2 is rather simple for standard MP, it only involves an update of the residual through a subtraction of the selected atom. For OMP though, a Gram Matrix computation followed by an optimal weights calculation is needed to ensure orthogonality between the residual and the subspace spanned by all previously selected atoms. Complexities of this two additional steps increases with the iteration number and has led to the development of bypassing techniques that limit complexity using iterative QR factorization trick [START_REF] Blumensath | Gradient pursuits[END_REF]. Complexities in Table 1 are given according to these tricks.

Update tricks and short atoms

Although RSS MP has the same complexity in the general case as Coarse MP, there are practical limitations to its use. First, Mallat et al [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] proposed a simple update trick that effectively accelerates a decomposition. Projection step is performed using previous projection values and the pre-computed inner products between atoms of the dictionary. Changing the subdictionary on an iteration-by-iteration basis prevents from using this trick. An even greater optimization is available when using fast transforms and more specifically when only local updates are possible. One of the main accelerating factor proposed for MP in the MPTK framework [START_REF] Krstulovic | Mptk: matching pursuit made tractable[END_REF] and for OMP in [START_REF] Mailhe | A low complexity orthogonal matching pursuit for sparse signal approximation with shift-invariant dictionaries[END_REF] is to limit the number of projections that require an update to match the support of the selected atom. Again, when changing dictionaries at each iteration, this trick can no longer be used to accelerate Step 1.

Using fast transforms, the complexity of the projection step at first iteration reduces to O(M log N ) for Full MP [START_REF] Mailhe | A low complexity orthogonal matching pursuit for sparse signal approximation with shift-invariant dictionaries[END_REF]. Additionally if one uses short atoms of size P << N , it further reduces to O(M log P ). When local update can be used instead of full correlation, the projection step after iteration 1 has an even more reduced complexity of O(αP log P ) in Φ I where α = L N is the redundancy factor and of βP log P in Φ where β = M N = α M P . Regardless of the dictionary, the residual only needs a local update around the chosen atom. Table 2 shows that in this context, the proposed modification becomes less competitive with respect to the standard algorithm on fixed subdictionary and can even be slower than the full dictionary case if L > P M L . However, it can still be considered as a viable alternative since memory requirements for full dictionary projections can get prohibitive.

In order to limit the computational burden, it is possible to change the subdictionary only once every several iterations. If the subdictionary Φ I i remains unchanged for J consecutive iterations, then the usual tricks for speeding up the projection can be applied. When Φ I i changes at iteration i + J, all projections need be recomputed. Although it can accelerate the algorithm, one may expect a resulting loss of quality.

Sparsity/Complexity tradeoff

Another way to tackle the complexity issue is by choosing an appropriate size for the subdictionaries. Since RSS MP yields much sparser solutions with dictionaries of the same size than Coarse MP, we can further reduce the size of the subdictionaries in order to accelerate the computation while still hoping to get compact representations. Indeed, in the discussion above, we have only compared complexities per iteration, but the total number of iterations is also important.

To illustrate this idea, we have used subsampled subdictionaries with a subsampling factor varying from 1 (no subsampling) to 32. The subsampling is performed by limiting the frame indexes in which we calculate the MDCT projections. It is important to notice that the subdictionaries are no longer redundant, and not even complete when the subsampling factor is greater than 2. Figure 10 shows for two signals how the subsampling impacts the sparsity of the representation (here expressed as a bit-rate calculated as in section 4) and the computational complexity relative to the reference Coarse MP algorithm for which we have implemented the local update trick.

Dictionaries used in this setup are multiscale MDCT with 8 different scales (from 4 to 512 ms), the experiment is run 100 times on a dual core computer up to a SRR of 10 dB. Using random subdictionaries 10 times smaller than the fixed one of Coarse MP, we can still have good compression ( cost 30 % less than the reference) with the computation being slightly faster.

More interestingly, this subsampling factor, associated with the RSS MP algorithm, controls a sparsity/complexity tradeoff that allows multiple use cases depending on needs and resources.

Conclusion

In this work, we proposed a modification of greedy algorithms of the MP family. Using a pseudo random sequence of subdictionaries instead of a fixed subdictionary can yield sparser approximations of signals. This sparsity basically comes at no additional coding cost if the random sequence is known in advance, thus giving our algorithm a clear advantage for compression purposes as shown here for audio signals. On the downside, the modification may increase complexity, but we have proposed a tradeoff strategy that helps reducing computation times.

The idea presented in this work can be linked to existing techniques from other domains, where randomness is used to enhance signal processing tasks. Dithering for quantization is one such technique. Another example is spread spectrum in communication [START_REF] Dixon | Spread Spectrum Systems[END_REF], where a signal is multiplied by a random sequence before being transmitted on a bandwidth that is much larger than the one of the original content. Here, encryption is more important a goal than compression. However, performing RSS MP of a signal somehow requires the definition of a key (the pseudo-random sequence) that must be known at the reception side in order to decode the representation. Moreover, the representation may also live in a much larger space than the original content. It is worth noticing for instance the work of Puy et al [START_REF] Puy | Universal and efficient compressed sensing by spread spectrum and application to realistic Fourier imaging techniques[END_REF] where the authors apply spread spectrum techniques to the compressed sensing problem.

Other experiments should be conducted for feature extraction tasks, since the good convergence properties indicates that the selected components carry meaningful information. For instance, music indexing [START_REF] Ravelli | Audio signal representations for indexing in the transform domain[END_REF] and face pattern recognition [START_REF] Phillips | Matching pursuit filters applied to face identification[END_REF] use Matching Pursuits to derive the features that serve for their processing. Further work will have to investigate whether RSS MP can improve on the performance of other algorithms on such tasks.

Figure 1 :

 1 Figure 1: Block diagram for General MP algorithms. An input signal f is decomposed onto a dictionary Φ by iteratively selecting atoms φγ n (Step 1) and updating an approximant fn (Step 2).

Figure 2 :

 2 Figure 2: Time Frequency grids defined by Φ a redundant monoscale time-frequency dictionary, Φ I 1 and Φ I 2 , two subsets of Φ. Atoms from Φ I 2 are almost the same as Φ I 1 but with an additional time and/or frequency offset.

Figure 3 :

 3 Figure 3: Block diagram of Sequential Subdictionaries Pursuits (SSMP).Step 1 is modified so as to have the search take place in a different subdictionary at each iteration. The dictionary subsampling is controlled by a fixed pre-defined sequence I.

Figure 4 :

 4 Figure4: Evolution of normalized approximation error ǫ with Matching Pursuit using a fixed subdictionary (Coarse MP), a pseudo random sequence of subdictionaries (RSS MP) and the full original dictionary (Full MP). Results averaged over 1000 runs.

Figure 5 :Figure 6 :

 56 Figure 5: Left : Pdf of Z and pdfs for the maximum and the M/2-th order statistic of Z for uniform, normal (σ = 1) and exponential (µ = 1) distributions. Right: relative error decay (mean and variances) as a function of iteration number. M = 100

Figure 7 :

 7 Figure 7: Number of iterations needed to achieve 10 dB of SRR. Means and standard deviation for various 4 seconds length audio signals from the MPEG SQAM dataset with MP on various fixed and varying sequences of subdictionaries and random initial time offset. Step SSMP and Jump SSMP are two deterministic SSMP.

Figure 8 :

 8 Figure 8: SNR/Bitrate curves for 6s-length signals from the MPEG SQAM dataset with multi-scale MDCT dictionaries (2 7 ,2 10 and 2 13 samples per window).

Figure 9 :

 9 Figure 9: Comparison between algorithm MP (thin) and OMP (bold) with fixed coarse subdictionary Φ I 0 (dashed dotted blue), random sequential subdictionaries Φ I n and full dictionary Φ (dashed red). Input signal is a collection of N/2 vectors from Φ with normal weights. Results averaged over 1000 runs. Here N = 128, M = 256Step Full MP Coarse MP Full OMP Coarse OMP RSS MP RSS OMPStep 1 Projection O(M N ) O(LN ) O(M N ) O(LN ) O(LN ) O(LN ) Selection O(M ) O(L) O(M ) O(L) O(L) O(L) Step 2 Gram Matrix 0 0 O(iN ) O(iN ) 0 O(iN ) Weights 0 0 O(i 2 ) O(i 2 ) 0 O(i 2 ) Residual O(N ) O(N ) O(iN ) O(iN ) O(N ) O(iN ) Total O(M N ) O(LN ) O(i 2 + M N ) O(i 2 + LN ) O(LN ) O(i 2 + LN )

Figure 10 :

 10 Figure 10: Illustration of the sparsity/complexity tradeoff that can be achieved by controlling the size of the subdictionaries in RSS MP.

Table 2 :

 2 Complexities after iteration 1 of the different steps for various algorithm in structured time-frequency dictionaries with fast transform and limited atomic support P << N . In this case, the local update speed up trick can not be applied to RSS MP.

http://www.tsi.telecom-paristech.fr/aao/?p=531

Acknowledgments

We are very grateful to the anonymous reviewers for their detailed comments that helped improve the quality of this paper. We also would like to thank Dr. Bob Sturm for the interest he showed in our work and his relevant remarks and advice. Discussing with him was both an opportunity and a pleasure. This work was partly supported by the QUAERO programme, funded by OSEO, French State Agency for innovation.


Manuel Moussallam

Laurent Daudet

Gaël Richard

RANDOM TIME-FREQUENCY SUBDICTIONARY DESIGN FOR SPARSE REPRESENTATIONS WITH GREEDY ALGORITHMS

Keywords: Matching Pursuits, Random Subdictionaries, Sparse Audio Coding

Sparse signal approximation can be used to design efficient low bitrate coding schemes. It heavily relies on the ability to design appropriate dictionaries and corresponding decomposition algorithms. The size of the dictionary, and therefore its resolution, is a key parameter that handles the tradeoff between sparsity and tractability. This work proposes the use of a non adaptive random sequence of subdictionaries in a greedy decomposition process, thus browsing a larger dictionary space in a probabilistic fashion with no additional projection cost nor parameter estimation. This technique leads to very sparse decompositions, at a controlled computational complexity. Experimental evaluation is provided as proof of concept for low bit rate compression of audio signals.

INTRODUCTION

Randomness has proven surprisingly useful in a wide variety of computational and statistical fields. In communications, spread spectrum techniques, where a signal is modulated by a random binomial sequence before transmission allows for better bandwidth management. Quantization has long taken advantage of the dithering technique that uses randomness to avoid perceptually disturbing artefacts linked to the quantization noise. More generally, stochastic resonance theory has shown how a moderate amount of added noise can increase the average behavior of many non-linear systems. More recently, tremendous work has been achieved on the compressive sampling scheme, making use of random measurement matrices. Behind all these examples lies a common intuition: controlling the random part of a system is better than having to deal with colored measurement noise or transmission errors, or signal-dependent deviations. The key point is to spread the information where it can be efficiently found. Having to guess where discriminant low-level features are hidden in huge-dimensional spaces is too costly or simply not feasible. In such cases, randomness can be used as a powerful sieve by information miners.

At the opposite of uniform random distributions is the concept of sparsity. Sparse coding of digital signals has been the subject of many works in the past few years for audio [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF], images [START_REF] Figueras I Ventura | Low-rate and flexible image coding with redundant representations[END_REF] or video streams. Low bitrate coders have been designed and proved to be competitive with state of the art industrial solutions. The core idea is to decompose an original signal f as a combination of (a few) objects from a dictionary Φ of indexed elementary waveforms. In a coding framework, the coder has to transmit the set of indexes of the LD is on a joint position between Univ. Paris Diderot and Institut Universitaire de France non-zero coefficients, together with their quantized values. Here, a crucial yet often underestimated issue is the choice of the size of Φ, that always resorts to a tradeoff. If the dictionary is small (slight or no overcompleteness), computations are fast, the index coding cost per coefficient is low, but many coefficients may be needed. If it is large (i.e., dense in the parameter space), we have a fast decay of the approximation error as a function of the approximation order, but the cost of encoding the indexes is higher, and computations get cumbersome. Indeed, the computational complexity associated with these sparse techniques, as opposed to suboptimal but much simpler transform-based coders, is probably the main limitation to their widespread use in practical applications. Strategies have been proposed to lower the computational cost of using large dictionaries, based on local adaptation of the selected atoms [START_REF]Fast matching pursuit with a multiscale dictionary of gaussian chirps[END_REF] or probabilistic approaches [START_REF] Durka | Stochastic timefrequency dictionaries for matching pursuit[END_REF] where successive runs with random sub-optimal atom selection are performed, then averaging yields a robust sparse approximation. Yet, these approaches are still associated with high index coding costs, if the atoms parameter space is large.

In this work, we propose a different paradigm that mitigates the drawbacks of using a large dictionary while keeping most of the benefits. Based on the algorithm described in [START_REF] Moussallam | Matching pursuits with random sequential subdictionaries[END_REF], a single run is performed using varying subdictionaries. These subdictionaries have limited size, but are designed so as to evenly span a much larger dictionary space. In this work, we use a simplistic audio coding example as a proof of concept to demonstrate the usefulness of randomization for sparse representation problems. The key point in our technique is that the choice of subdictionary is not adaptive, but is parameterized by a fixed pseudo-random sequence, also known by the decoder. In other words, we have the (theoretical) complexity of working with a small dictionary, and the small coding costs, but the whole large dictionary is spanned. It should be emphasized that, unlike a compressive sensing framework, we our goal is to design a 'standard' coding scheme with maximal efficiency at the cost of computational complexity at the encoder, and minimal decoding complexity at the decoder.

The rest of this paper is organized as follows. Section 2 recalls the Matching Pursuit framework for compression. In Section 3, the novel approach is presented along with considerations on the random sequence design, and in Section 4 trivial audio compression serves as a proof of concept for the suitability of the proposed method.

SIGNAL COMPRESSION WITH GREEDY ALGORITHMS

Let f a discrete signal living in R N . Greedy algorithms iteratively decompose f using m elements from a dictionary Φ = {φγ} of elementary objects called atoms by alternating two steps 1: Select an atom in the dictionary and 2: Update approximant and residual.

where δ(f, fm) is a distorsion measure, usually in the form of a quadratic reconstruction error. Ideally, the dictionary is chosen to be related to the signal's characteristics. Speech can be efficiently compressed using carefully designed codebooks, 2D wavelets bases are useful for image compression. Music representation is better achieved using dictionaries of windowed cosines (MDCT) or gabor functions. The algorithm selects at iteration i the atom that maximizes a correlation function, usually an energy criterion:

φ ⋆ γ i = arg max φγ ∈Φ | R i f, φγ | (2)
where R i f = f -fi is the residual signal. The approximant update depends on the nature of the algorithm but can generally be seen as a descent in a direction defined by the newly selected atom (plain Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency IEEE Trans[END_REF]), the subspace spanned by all previously selected atoms (Orthogonal MP) or a gradient estimate ( [START_REF] Blumensath | Gradient pursuits[END_REF]). This criterion can also be modified to take perceptual models into account or dealing with pre-echo artefacts. However, there is always, as mentioned in the introduction, the central problem of choosing the size of Φ. From a continuous parameter space, choosing a dictionary Φ for practical use amounts to discretizing the parameters.

In the compressive sensing framework, this leads to the so-called basis mismatch problem [START_REF] Chi | Sensitivity to basis mismatch in compressed sensing[END_REF]: the chosen representation space is not exactly the one where f has the sparsest representation.

Pursuits in time-frequency dictionaries

Time frequency dictionaries such as Gabor Dictionaries and MDCTbased dictionaries are well suited for audio signals. Time and frequency resolutions are constrained by the scale of the chosen transform that defines the analysis window size and the overlap intervals between consecutive analysis windows. The finer this analysis grid, the better the chances of choosing well-localized atoms that remove a lot of perceptually relevant energy. However, the complexity gets higher. Both for tractability and compression purposes, the timefrequency grid that serves for inner product computations need to be quantized. Let Φ = {φ s,u,ξ } be a dictionary of localized waveforms of scale s ∈ S, time position u ∈ U and frequency index ξ ∈ Ξ.

The set of available atom indexes is denoted Γ = S × U × Ξ.

Quantization of the time frequency plane

A convenient way to model Φ that coincides with practical implementations is to see it as a union of monoscale dictionaries: Φ = s∈S Φ ♯s each of which defines a quantization of the time frequency plane with resolution defined by the set of indexes ♯s = Us × Ξs. If the quantization is uniform, the size of Φ ♯s is

Ts = N ∆su × Fs 2∆sξ ( 3
)
where ∆su is the time interval between frames for the scale s, Fs is the sampling frequency and ∆sξ is the frequency resolution. The subdictionary Φ ♯s can be an orthonormal base (e.g an MDCT with 50 % overlap [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]), in this case ∆su = s/2 and ∆sξ = Fs/s and then Ts = N It can also be overcomplete (Gabor dictionaries with more than 50% overlap: ∆su < s/2, Ts > N ) or span only a limited subspace (∆su > s/2, Ts < N ). ♯s defines a two dimensional lattice in the time-frequency plane that can be seen as a quantization of the underlying continuous time-frequency parameters. Indeed, the choice of a suboptimal atom can be understood as a quantization error artefact. Figure 1 shows how different lattices can fit different components of a signals, here transients and harmonics of a glockenspiel signal. The quantization error is greatly lowered by the concatenation of all the bases in the dictionary [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]. Nonetheless, traditional Matching Pursuit-based strategies are using the same dictionary during the whole decomposition process. By doing so, the a priori choice of lattices introduces a bias in the decomposition as explained in Durka's work [START_REF] Durka | Stochastic timefrequency dictionaries for matching pursuit[END_REF]. While their solution would be to run multiple decompositions with different lattices and averaging the results in a Monte Carlo fashion, we propose an novel approach inspired by the dithered quantization technique.

PURSUITS WITH A RANDOM SEQUENCE OF SUBDICTIONARIES

Instead of choosing an analysis grid once and for all, the lattices are chosen so as to span the largest possible space in an ergodic fashion during the decomposition. For each scale, a sequence of lattices ∆s♯ = ♯ 0 s , ♯ 1 s , ..♯ k s is used and at iteration i an atom is selected in Φi = s∈S Φ ♯ i s . doing so, we expect the equivalent of the quantization error to be evenly spread among the selected atoms, thus removing the bias. This technique is conceptually equivalent to adding a uniform noise in the time-frequency domain before quantizing it. In the overall, we hope to promote the selection of more salient features than with a fixed lattice. This technique is completely non-adaptive. The sequences ∆s♯ are known in advance and independant from the signal. They are also known by the decoder and therefore there is no need to encode them. The virtual cost of the index of atom i is then down to O(log 2 ( s∈S T i s )).

Learning the sequence

Since we want the subset sequence to be independant of the signal, we can try to estimate its desirable properties. When a signal model is available (e.g sinusoidal+transient modelling of audio, edges and textures modelling of images), one can try to manually design a sequence that will minimize the quantization error (i.e the suboptimality factor) under a dictionary size constraint. In this work however, we have no signal model and we are interested in designing a universal sequence for audio signals. Figure 2 shows a decomposition of a short glockenspiel signal with MP over a full temporal-resolution multiscale discrete Gabor dictionary. Frequency resolution is constrained in each scale, the full temporal resolution is achieved by performing Short Time Fourier Transforms with high overlapping between consecutive analysis frames (i.e ∆su = 1 sample). Atoms are clearly not uniformly distributed in the time-frequency plane, which would be the case for white noise. Most real life signals are

Fig. 2. Time-frequency joint and marginal distribution of atoms in the decomposition of a glockenspiel signal with MP on a multiscale

Gabor dictionary with full time resolution.

structured, a complete randomness in the choice of an analysis lattice would not give optimal performance.

Time resolution subsampling

In particular, for audio signals, the frequency resolution problem seems to be efficiently adressed by the use of multiple scales. The time resolution, however, presents a more interesting challenge. Let us consider a reference coarse lattice ♯s with overlap of ∆su = s/2 between consecutive frames and define the local time shift τ of an atom relative to this coarse grid. Then the distribution of τ in the interval [-s 4 , s 4 ] ressembles a uniform distribution. To verify this statement, we decomposed short audio signals from the MPEG SQAM test database up to the first 1000 atoms in a full resolution multiscale Gabor dictionary and calculated the Gini index of τ . This index quantifies how far from the uniform distribution a candidate distribution stands and has been demonstrated [START_REF] Hurley | Comparing measures of sparsity[END_REF] to be a suitable sparsity measure.

Figure 3 shows that distribution for the decomposition of the orchestra signal ressembles the one for the decomposition of white noise: atoms become more and more uniformly spread. The glockenspiel signal's atoms are less uniformly distributed, but their distribution can not be considered sparse after a few hundreds iterations. From this observation, we state that an efficient and simple way to simulate a pursuit in a large dictionary is to use orthonormal basis randomly shifted at each iteration by a random variable τ .

SCALABLE SPARSE CODING OF AUDIO SIGNALS

Performances with a simple encoding model

In order to demonstrate the potential benefits of this technique, we compared traditional approaches to the new one in a simplified audio coding task. We considered 3 cases:

MP with a union of 8 MDCT: ∀s, ∆su = s/2, Ts = N and T = SN .

MP with a union of 8 MDCT with full temporal resolution (for computational reasons, an approximate based on local optimization after a coarse grid search is used as in [START_REF] Moussallam | Matching pursuits with random sequential subdictionaries[END_REF]) Ts = N s/2, T = s∈S N s/2. MP with Random Sequence of Subdictionaries that are 8 MDCT randomly shifted in time. The subdictionary size is constant T i s =

N .

The distorsion measure that we used was the Signal to Noise Ratio (SNR), and bitrates are estimated using the upper bound C(fm) = m (log 2 (T ) + Q), where we assume a simple uniform mid-tread quantizer with Q bits per coefficient. Note that using an entropy coder instead has also been tested, with similar results. Figure 4 summarizes the coding scheme. Figure 5 shows that the novel algorithm gives better performances in all cases but at very low bit-rates. This can be explained by the fact that, as exposed in 3.2, in the first iterations, the full resolution dictionary succesfully locates the most prominent features and it actually compensates the additional costs. At bigger rates, however, atoms are more evenly spread and the lower index cost favors the randomized method. To summerize our results, this randomized greedy pursuit allows to have the costs of the small dictionary with a decay rate close to the one on the large dictionary.

Controlling computational complexity

Although the theoretical complexity of the novel algorithm is equivalent to the one of the fixed small dictionary case, in practice it can get much higher. Inner products computation, actually, are usually performed using Fast transforms, and in the fixed dictionary case, most of these products remains unchanged from one iteration to the next, thus greatly reducing the cost of all but the first iteration. In our case, however, this trick cannot be applied anymore since the central point is to change all the projections at every iteration. There is nonetheless one parameter that allows us to control the complexity of the algorithm and it is the size of the subdictionaries. So far, have considered Φi to be overcomplete, by decreasing the time resolution of each lattices ♯ i s , less inner products need to be comat each iteration and the algorithm is much faster. By doing so, however, we can severely damage the convergence rate of the algorithm.

To evaluate the performances, we used different time resolution ∀s, ∆su ∈ [s/2, s, 2 × s, .., 12 × s]. We used the same simple encoding scheme and compared execution times and bitrates between the randomized method and a reference MP with fixed dictionary of 8 MDCT with the same optimizations than in [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF]. All algorithms were given the same approximation quality target of 5 of SNR.

Figure 6 shows that the time resolution parameter can serve as a control parameter for limiting complexity while still improving the rate. A time resolution of ∆su = 5 × s which is equivalent to subdictionaries of size T i s = N/10 yields a bitstream 20% smaller than MP with the fixed grid, with a slightly faster computation time. With a speech signal, all cases took slightly longer than the reference but the bitrate was always smaller. Finally, the glockenspiel signal gave the less competitive results: its components are so well localized that the subsampling quickly penalizes the achieved bitrate. However in this case, the randomized MP with overcomplete subdictionaries took less than 3 times longer than the reference MP. This can be explained by the fact that our algorithm selects better atoms than the reference and thus, much less iterations are needed to achieve the same approximation level.

CONCLUSION

The proposed algorithm appears to be suitable for sparse approximation of complex signals. The potential benefits are in low bitrate compression, and we exhibit several sound examples were these advantages show off. The unsupervised nature of the algorithm and the randomness introduced in the atom selection makes it very easy to design worst-case scenarii for which the algorithm would converge slower than a pursuit over a fixed dictionary. However, on average, and with a small empirical variance, the proposed scheme appears to have the coding costs of the small dictionary with a decay rate close to the one on the large dictionary.

In conclusion, adding randomness to the parameter space within a greedy sparse decomposition process can be highly beneficial. Although on different paradigms, this is reminiscent of both dithered quantization techniques and compressive sampling strategies. Whether such random sub-dictionary techniques can be combined with other sparse decomposition schemes (e.g. iterated thresholding) is still an open issue.

Fig. 1 .

 1 Fig. 1. Projection of a glockenspiel signal quantized with 3 different time frequency lattices ♯s for s ∈ [128, 1024, 8192].

Fig. 3 .

 3 Fig. 3. Gini index for the time localization of the first 1000 atoms in a full resolution multiscale MDCT dictionary for glockenspiel, orchestra and white noise signals.

Fig. 5 .

 5 Fig. 5. SNR achieved at different rates for 5 signals (6 seconds) from MPEG -SQAM database. Results averaged over 20 runs.

Fig. 4 .

 4 Fig. 4. Simplistic encoding scheme with greedy decompositions: bitstream includes atom indexes from the sequence of subdictionaries Φ i and the quantized weights.


Manuel Moussallam

Laurent Daudet

Gaël Richard

AUDIO SIGNAL REPRESENTATIONS FOR FACTORIZATION IN THE SPARSE DOMAIN

Keywords: Sparse Representation, Matching Pursuit, Audio Signal Decomposition, Audio Similarity, Factorization

In this paper, a new class of audio representations is introduced, together with a corresponding fast decomposition algorithm. The main feature of these representations is that they are both sparse and approximately shift-invariant, which allows similarity search in a sparse domain. The common sparse support of detected similar patterns is then used to factorize their representations. The potential of this method for simultaneous structural analysis and compressing tasks is illustrated by preliminary experiments on simple musical data.

INTRODUCTION

The need for intelligent storage of ever-growing digital multimedia data has posed an interesting challenge on the Computer Science community. The rise of search engines, file-sharing plateforms and social networks on the internet, has proven that the ability to quickly retrieve consistent information from huge databases is now as important as managing its storage. One can see it as a joint optimization problem of space consumption and data access time but also as a concurrency between machine-oriented and human-readable systems. The compromise is set by the choice of a representation (i.e a transformation) of the data, that will hopefully be much sparser than the original data while keeping as much information as possible (ideally all of it), that is to say capturing only its very essence, and making it easily parsable by machines and/or humans.

A typical such task is archival of audio data, where indexing and compressing at the same time is crucial. Intuition (along with considerable research in the past years) tells us that a good way to proceed is to decompose the time-series that constitute the raw data in audio, as a combination of elementary sound objects or atoms. This idea is also strongly linked to the belief that human brain somehow processes sounds in a similar fashion, at least at a perception level [START_REF] Lewicki | Efficient coding of natural sounds[END_REF].

So far, audio coding schemes have essentially focused on reducing short-term redundancies (usually frame-by-frame). The possibility of using long term redundancies, as can be done on text using sequential data compression methods [START_REF] Ziv | A universal algorithm for sequential data compression[END_REF] where repeating patterns are detected and factorized, remains a challenging issue for audio signals. However, pattern repetition is an extremely common feature in music, with groups of notes called "motifs" in classical music, or "riffs" in popular music. At a larger time-scale, the verse-chorus structure is another example of (near) repetition. In a radio station, the same popular song may be played tens of times during a single day. For archival purposes, one may be interested in audio representations where these repetitions, or near-repetitions, appear clearly.

Looking for large-scale similarities directly in the PCM domain is doomed to failure, first because the data rates would make this an extremely high-dimensional search problem, but more fundamentally because most of the time the repeats are not perfect at the bit level : note levels, durations, etc. may be slightly different. A better idea is to look for similarities in a sparse domain, for instance one of the transform-domain used in audio coding, e.g. the Modified Discrete Cosine Transform (MDCT). Recently, Ravelli [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF] has enlightened the potential of a multiscale MDCT-based greedy decompositions for very low bit-rate coding : this larger dictionary allows for even sparser decompositions. Working in a sparse domain has substantial advantages: first, the amount of data is reduced, with hope that pattern classification becomes tractable again ; then, the information is sorted, with the largest coefficients accounting for the most salient sound structures : similarity search can then be progressive, first at a crude level, and then with finer and finer details.

However, these MDCT-based decompositions suffer from a major drawback : because they rely on some a priori slicing of the signal, they are not shift-invariant. In other words, the representation of a time-shifted signal x(•τ ) can be quite different from the representation of x with the time parameters shifted by τ . All the timing (i.e., phase) information is here hidden in the amplitude of the MDCT parameters.

The main contributions of this paper can be summarized as follows :

• we propose a new shift-invariant decomposition (in the sense of [START_REF] Blumensath | Sparse and shift-invariant representations of music[END_REF]), based on adding time-shifts to the multiscale MDCT of [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF],

• we provide an efficient decomposition scheme, called LOMP (for Locally Optimized Matching Pursuit), that effectively finds near shift-invariant decompositions at the cost of a small increase in computational complexity per iteration, compared to standard shift-dependent Matching Pursuit (note that this drawback is offset by a significantly faster energy decay, hence the addition of shift-invariance results in a faster algorithm per dB of SNR). • we assess the potential of these decompositions for factorizing repeated musical patterns, on preliminary tests.

The rest of the paper goes as follows: Section 2 briefly introduces the multiscale MDCT Matching Pursuit framework and its observed limitations for our purposes. The novel shift-invariant decomposition, together with an effective decomposition algorithm, that reduces these limitations is presented in Section 3. Then, in Section 4, this new representation is applied to audio factorization and preliminary compression experiments illustrate its potential.

SIGNAL MODEL

Sparse Representation problem

Let x ∈ R N be a digital signal of length N . Let Φ = {φ k } k=1..K be a dictionary of K waveforms. We are looking for the smallest combination of I atoms from Φ that minimizes a quadratic error ǫ = x-P I i=1 αiφi 2 . For very large signals such as for audio, the most widely-used approach is to use a greedy algorithm, that iteratively builds a representation by selecting at each step the best atom in the dictionary and subtracting it from the signal. Convergence is assured by the fact that the overall energy of the residual is strictly decreasing, but there is no guarantee of finding the global minimum. The basic version is called Matching Pursuit [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF], and it has been extensively studied, modified and improved [START_REF] Gribonval | Analysis of sound signals with high resolution matching pursuit[END_REF]. In practice, one has to stop the algorithm after a finite number of steps, usually when the approximation x = P ı∈I αiφi fulfills a fidelity depth criteria, expressed as a Signal-To-Residual Ratio (SRR) :

SRR = 10 log " x 2 x -x 2 « ( 1
)
The approximation is fully characterized by the vector α which has I non zero coefficients. If the dictionary is redundant enough (i.e K > N ), one might expect sparsity (i.e I << N ) to arise at an acceptable fidelity (SRR) level.

Multiscale MDCT dictionary

The Multiscale MDCT dictionary as described in [START_REF] Ravelli | Union of MDCT bases for audio coding[END_REF] consists in a union of 8 MDCT basis of different scales Φ = S Φm, distributed in a dyadic way. The multiscale nature of the dictionary provides a structural decomposition of the sounds, with harmonics (using large scales), transient parts (using short, well localized atoms) and speech or noise (using middle-sized atoms). Efficient implementations of the Matching Pursuit algorithm take advantage of the hierarchical structure of the dictionary, and the use of Fast Fourier Transform to achieve nearly real-time performances.

Similarity detection and time invariance

The issue of similarity detection in the sparse domain has been studied for example in [START_REF] Sturm | On similarity search in audio signals using adaptive sparse approximations[END_REF]. However, our main goal is quite different, since we seek to use the similarity information as a tool for factorizing similar audio segments. Thus, an additional constraint on the representation is that it should have some robustness to time shifts (in Blumensath's sense [START_REF] Blumensath | Sparse and shift-invariant representations of music[END_REF] ) that may result from non-aligned frame slicings. Representations obtained with the multi-scale dictionary above are not time-invariant, as with any MDCT-based decomposition. Actually, the phase information is embedded in the α coefficients, and thus they greatly vary with temporal offsets. This is a serious obstacle to similarity detection, since two similar audio segments may have different sparse decompositions if they are not sample-aligned.

LOCAL OPTIMIZATION OF THE MATCHING PURSUIT

Shift-invariant multiscale MDCT dictionary

An effective yet computationally intensive solution to the above issue is to use a fully shift-invariant dictionary as in [START_REF] Blumensath | Sparse and shift-invariant representations of music[END_REF]. Let Φ G be a multiscale MDCT dictionary, and let Φ G m be a single scale subdictionary of size Lm. The MDCT is a 50% overlapped transform, therefore, for a N -length signal, Φ G m is an orthonormal basis of R 2N Φ G m can be upgraded in the following manner: for each atom φm,i of Φ G m , add all possible shifted versions φm,i,τ = φm,i * δτ with τ living in ˆ-Lm 4 , Lm 4 ˜to the dictionary, thus yielding Φ F m a highly redundant dictionary of size N.Lm. Repeating this operation in all the basis, a global dictionary Φ F = S Φ F m can be built.

Hybrid approach: Locally Optimized Matching Pursuit (LOMP)

Using Φ F leads to sparser, clearer representations, but the complexity of a Matching Pursuit over Φ F (called here FullMP) would be prohibitive for real data analysis. We propose here a hybrid approach that uses atoms from Φ F at the cost of a search over Φ G . Alternatively, this can be seen as an attempt to extract the phase information out of α into a separate vector τ , thus yielding phaseinvariance properties. At iteration n let r n be the residual, the best atom

φ G kmax = arg max φ k ∈Φ G | φ k , r n | of length Lm is selected in Φ G . Then a local optimization is performed yielding φ F kmax,τmax = arg max τ ∈[-Lm 4 , Lm 4 ] | D r n , (φ G kmax * δτ ) E | (2)
where * denotes the convolution product and δτ is a discrete dirac impulse located in τ (this convolution is simply a time shift by τ samples). The modified algorithm is given in Algorithm 1, the major difference with the standard MP algorithm is summarized in step 4 in bold. It is worth noticing that contrary to the method proposed in [START_REF] Gribonval | Analysis of sound signals with high resolution matching pursuit[END_REF], where smaller sub-atoms are constructed at each iteration to fit locally the energy distribution, the LOMP algorithm performs a temporal realignement of the selected atom, which can be implemented efficiently using FFTs. Update residual : r n+1 ← r nαmax.φ kmax,τmax 7: until a stopping condition is met

Convergence discussion

As a weak (or partially) greedy algorithm, one can prove (see for instance [START_REF] Temlyakov | A criterion for convergence of weak greedy algorithms[END_REF]) the convergence of Algorithm 1. Moreover, one might easily be convinced that, at least for the first iterations, the convergence rate of the FullMP algorithm is faster than standard MP on the original (and much smaller) MDCT dictionary (called here "MP"). For LOMP algorithm, an intermediate behavior can be expected.

Figure 1 offers an insight of this for a short segment of orchestra music taken from the MPEG SQAM database. Implementation of all three algorithms using Python programming language showed an equivalent complexity for MP and LOMP, on a standard PC. Indeed, at each iteration : normalized by the SRR gain, LOMP has a smaller computation time than the standard (non shift-invariant) MP.

Time invariance discussion

LOMP atoms are more expensive to encode since an additional timeshift coefficient is to be encoded. However, since the convergence rate is often faster, one might hope to obtain sparser representation.

FACTORIZATION IN THE SPARSE DOMAIN

Factorization paradigm

Let x and ỹ be the two LOMP representation in Figure 2 (c) and (d). The set of atoms {φi} and coefficient vector α for x and ỹ are almost equals. In other words, both representation share the same sparse support and their differences are entirely characterized by the sets of time-shifts τ x and τ y . This has two important consequences: first the common sparse support denotes a similarity between x and y. Second, it means that one only needs to encode α and the set of atoms {φi} once, as a common factor to describe both x and ỹ. Alternatively, ỹ can be described using the sparse support of x and τ y . Now this paradigm can be extended to any x and y. Knowing a representation x as a reference, let ŷ be an approximation of y, which support is the same as x and has a set of time-shift τ y . A method to build such ŷ is described by Algorithm 2

Algorithm 2 Factorization

Input: y , Φ , x = P i∈I αiφi * δ(τ x i ) Output: τ y such that y = P i∈I αiφi * δ(τ y i ) + r 1: r 0 := y 2: for all i ∈ I do Update residual : r n+1 ← r nαi.(φi * τmax) 5: end for This paradigm is close to the information theory concept of distributed source coding where repetitions are considered as correlated sources that are not co-located [START_REF] Pradhan | Distributed source coding using syndromes (discus): design and construction[END_REF]. However, in audio and musical streams, a prior detection of these repeating segments is required.

Similarity Detection

If x and y are either identical or simply similar (which might be the case for repeated choruses in a pop song, chords sequences in classical music, or even whole songs and news reports on radio broadcast recordings), then chances are that atoms from x are useful to de- To illustrate this method, Figure 3 presents the SRR of ŷ when using 4 different {xi}i=1..4 as reference signal for the factorization, and at different decomposition levels. Here y is the first 8 piano notes (2 s, 1st half of the 1st bar) taken from the left channel of a piano recording of Bach's first Prelude in C Major, from the Well-Tempered Clavier. The first case is simply x1 = y, obviously the factorization is complete and the reached SRR matches the reference's one. More interestingly, x2 is the exact same musical pattern played right afterwards in the piece (2nd half of the 1st bar), thus the correlation between x2 and y is far lower. Third example, x3 is the right channel corresponding to y. Finally, for the irrelevant case, x4 is a few seconds exctracted from an orchestra recording, therefore having no similarity with y .

One can observe that a factorization at a shallow depth (3 to 5 dB) allows us to detect the three types of similarity and reject the irrelevant case. Moreover, the detection process yields a preliminary decomposition of y in the same time that is reveals the signal structure.

Compression Experiments

Knowing its reference x, a factorized approximation ŷ of y is described completely by the set of atom time shifts τ y . One can picture the factorization process as a greedy decomposition in which the sequence of atom choice in the dictionary and their amplitude is arbitrarily fixed (only local time-shifts are allowed). This strong constraint can result in a lower approximation precision compared to MP decomposition, but its cost is also fairly reduced, which means that high compression levels can be achieved.

Working on the same Bach's signal, different decomposition techniques are compared: MP and LOMP on a 8xMDCT dictionary, and factorization using LOMP decompositions of x1, x2 and x3, in a simple encoding task. For MP and LOMP decompositions, atom coefficients are encoded using a simple 7-bit uniform scalar quantizer and plain entropy coding of the quantized coefficients, atom indexes have a fixed cost of log 2 (2.N ) and for LOMP atoms, an additional cost for the time shifts is also evaluated by plain entropy coding. For factorized representations, only the set of time shifts is encoded.

Results are shown Figure 4. The quality of the compression is here measured by the Signal-To-Noise Ratio (SNR) of the decoded signal. At very low bitrates, the factorization approach is interesting, even when the reference is a different (albeit closely related) signal. When near exact redundancy is detected and used for factorization, important coding gains may be achieved with this technique.

CONCLUSION AND FUTURE WORK

We have presented a multiscale greedy decomposition of audio signals that introduces robustness to time shifts, while increasing the sparsity of the representation, at a negligible computational penalty. Using this framework, similarity detection can be performed along with an enhanced compressive scheme through factorization in the representation domain, even for non-perfect redundancies. Many extensions of this work are yet to be undertaken, among which a study of the semantic similarities that may or may not be factorized. The scalability of this approach will also be investigated so as to design a complete audio coding scheme based on this model.

Algorithm 1 repeat 3 : 4 : 5 :

 1345 Locally Optimized Matching Pursuit (LOMP) Input: x , Φ Output: α, τ and r such that x = P i∈I αi.(φi * δ(τi)) + r 1: r0 := x 2: Find index kmax = arg max φ k ∈Φ | φ k , r n | Compute local optimal time-shift τmax = arg maxτ | r n , (φ kmax * δτ ) | Update projection score : αmax ← r n , φ kmax,τmax 6:

Fig. 1 .

 1 Fig. 1. Comparison of (a) : convergence rates and (b) : Complexities in computational time , for MP , FullMP and LOMP decomposition of 4096 samples of orchestra signal over a dictionary of 3 MDCT basis (32,128,512).

Fig. 2 .

 2 Fig. 2. Energy distributions of a bell sound decomposed by MP and LOMP algorithms. Left: original signal decomposition at 20 dB of SRR with 8 MDCT basis with MP (a) and LOMP (c). Right: decomposition of the same signals shifted in time by 100 samples with MP (b) and LOMP (d).

3 :

 3 Compute local optimal time-shiftτmax = arg maxτ | r n , (αi.φi * δτ ) | 4:

Fig. 3 .

 3 Fig. 3. SRR reached by factorizing y (few piano notes) using 4 different signals xi as references, at 5 different levels of decomposition

Fig. 4 .

 4 Fig. 4. Bitrates and SNR obtained with MP, LOMP or factorization using different references

This work was partly supported by the QUAERO Programme, funded by OSEO, French State agency for innovation and by the DReaM project (ANR-09-CORD-006) of the French National Research Agency CONTINT program.