# lowrank_timefrequency_synthesis__a4d4788b.pdf

Low-Rank Time-Frequency Synthesis

C edric F evotte Laboratoire Lagrange (CNRS, OCA & Universit e de Nice) Nice, France cfevotte@unice.fr

Matthieu Kowalski Laboratoire des Signaux et Syst emes (CNRS, Sup elec & Universit e Paris-Sud) Gif-sur-Yvette, France kowalski@lss.supelec.fr

Many single-channel signal decomposition techniques rely on a low-rank factorization of a time-frequency transform. In particular, nonnegative matrix factorization (NMF) of the spectrogram the (power) magnitude of the short-time Fourier transform (STFT) has been considered in many audio applications. In this setting, NMF with the Itakura-Saito divergence was shown to underly a generative Gaussian composite model (GCM) of the STFT, a step forward from more empirical approaches based on ad-hoc transform and divergence specifications. Still, the GCM is not yet a generative model of the raw signal itself, but only of its STFT. The work presented in this paper fills in this ultimate gap by proposing a novel signal synthesis model with low-rank time-frequency structure. In particular, our new approach opens doors to multi-resolution representations, that were not possible in the traditional NMF setting. We describe two expectation-maximization algorithms for estimation in the new model and report audio signal processing results with music decomposition and speech enhancement.

1 Introduction

Matrix factorization methods currently enjoy a large popularity in machine learning and signal processing. In the latter field, the input data is usually a time-frequency transform of some original time series x(t). For example, in the audio setting, nonnegative matrix factorization (NMF) is commonly used to decompose magnitude or power spectrograms into elementary components [1]; the spectrogram, say S, is approximately factorized into WH, where W is the dictionary matrix collecting spectral patterns in its columns and H is the activation matrix. The approximate WH is generally of lower rank than S, unless additional constraints are imposed on the factors.

NMF was originally designed in a deterministic setting [2]: a measure of fit between S and WH is minimized with respect to (w.r.t) W and H. Choosing the right measure for a specific type of data and task is not straightforward. Furthermore, NMF-based spectral decompositions often arbitrarily discard phase information: only the magnitude of the complex-valued short-time Fourier transform (STFT) is considered. To remedy these limitations, a generative probabilistic latent factor model of the STFT was proposed in [3]. Denoting by {yfn} the complex-valued coefficients of the STFT of x(t), where f and n index frequencies and time frames, respectively, the so-called Gaussian Composite Model (GCM) introduced in [3] writes simply

yfn Nc(0, [WH]fn), (1)

where Nc refers to the circular complex-valued normal distribution.1 As shown by Eq. (1), in the GCM the STFT is assumed centered (reflecting an equivalent assumption in the time domain which

Authorship based on alphabetical order to reflect an equal contribution. 1A random variable x has distribution Nc(x|µ, λ) = (πλ) 1 exp (|x µ|2/λ) if and only if its real and imaginary parts are independent and with distribution N(Re(µ), λ/2) and N(Im(µ), λ/2), respectively.

is valid for many signals such as audio signals) and its variance has a low-rank structure. Under these assumptions, the negative log-likelihood log p(Y|W, H) of the STFT matrix Y and parameters W and H is equal, up to a constant, to the Itakura-Saito (IS) divergence DIS(S|WH) between the power spectrogram S = |Y|2 and WH [3].

The GCM is a step forward from traditional NMF approaches that fail to provide a valid generative model of the STFT itself other approaches have only considered probabilistic models of the magnitude spectrogram under Poisson or multinomial assumptions, see [1] for a review. Still, the GCM is not yet a generative model of the raw signal x(t) itself, but of its STFT. The work reported in this paper fills in this ultimate gap. It describes a novel signal synthesis model with low-rank time-frequency structure. Besides improved accuracy of representation thanks to modeling at lowest level, our new approach opens doors to multi-resolution representations, that were not possible in the traditional NMF setting. Because of the synthesis approach, we may represent the signal as a sum of layers with their own time resolution, and their own latent low-rank structure.

The paper is organized as follows. Section 2 introduces the new low-rank time-frequency synthesis (LRTFS) model. Section 3 addresses estimation in LRTFS. We present two maximum likelihood estimation approaches with companion EM algorithms. Section 4 describes how LRTFS can be adapted to multiple-resolution representations. Section 5 reports experiments with audio applications, namely music decomposition and speech enhancement. Section 6 concludes.

2 The LRTFS model

2.1 Generative model

The LRTFS model is defined by the following set of equations. For t = 1, . . . , T, f = 1, . . . , F, n = 1, . . . , N:

fn αfnφfn(t) + e(t) (2)

αfn Nc(0, [WH]fn) (3) e(t) Nc(0, λ) (4)

For generality and simplicity of presentation, all the variables in Eq. (2) are assumed complexvalued. In the real case, the hermitian symmetry of the time-frequency (t-f) frame can be exploited: one only needs to consider the atoms relative to positive frequencies, generate the corresponding complex signal and then generate the real signal satisfying the hermitian symmetry on the coefficients. W and H are nonnegative matrices of dimensions F K and K N, respectively.2 For a fixed t-f point (f, n), the signal φfn = {φfn(t)}t, referred to as atom, is the element of an arbitrary t-f basis, for example a Gabor frame (a collection of tapered oscillating functions with short temporal support). e(t) is an identically and independently distributed (i.i.d) Gaussian residual term. The variables {αfn} are synthesis coefficients, assumed conditionally independent. Loosely speaking, they are dual of the analysis coefficients, defined by yfn = P

t x(t)φ fn(t). The coefficients of the STFT can be interpreted as analysis coefficients obtained with a Gabor frame. The synthesis coefficients are assumed centered, ensuring that x(t) has zero expectation as well. A low-rank latent structure is imposed on their variance. This is in contrast with the GCM introduced at Eq. (1), that instead imposes a low-rank structure on the variance of the analysis coefficients.

2.2 Relation to sparse Bayesian learning

Eq. (2) may be written in matrix form as

x = Φα + e , (5)

where x and e are column vectors of dimension T with coefficients x(t) and e(t), respectively. Given an arbitrary mapping from (f, n) {1, . . . , F} {1, . . . , N} to m {1, . . . , M}, where M = FN, α is a column vector of dimension M with coefficients {αfn}fn and Φ is a matrix of size T M with columns {φfn}fn. In the following we will sometimes slightly abuse notations by

2In the general unsupervised setting where both W and H are estimated, WH must be low-rank such that K < F and K < N. However, in supervised settings where W is known, we may have K > F.

indexing the coefficients of α (and other variables) by either m or (f, n). It should be understood that m and (f, n) are in one-to-one correspondence and the notation should be clear from the context. Let us denote by v the column vector of dimension M with coefficients vfn = [WH]fn. Then, from Eq. (3), we may write that the prior distribution for α is

p(α|v) = Nc(α|0, diag(v)) . (6)

Ignoring the low-rank constraint, Eqs. (5)-(6) resemble sparse Bayesian learning (SBL), as introduced in [4, 5], where it is shown that marginal likelihood estimation of the variance induces sparse solutions of v and thus α. The essential difference between our model and SBL is that the coefficients are no longer unstructured in LRTFS. Indeed, in SBL, each coefficient αm has a free variance parameter vm. This property is fundamental to the sparsity-inducing effect of SBL [4]. In contrast, in LRTFS, the variances are now tied together and such that vm = vfn = [WH]fn .

2.3 Latent components reconstruction

As its name suggests, the GCM described by Eq. (1) is a composite model, in the following sense. We may introduce independent complex-valued latent components ykfn Nc(0, wfkhkn) and write yfn = PK k=1 ykfn. Marginalizing the components from this simple Gaussian additive model leads to Eq. (1). In this perspective, the GCM implicitly assumes the data STFT Y to be a sum of elementary STFT components Yk = {ykfn}fn . In the GCM, the components can be reconstructed after estimation of W and H , using any statistical estimator. In particular, the minimum mean square estimator (MMSE), given by the posterior mean, reduces to so-called Wiener filtering:

ˆykfn = wfkhkn

[WH]fn yfn. (7)

The components may then be STFT-inversed to obtain temporal reconstructions that form the output of the overall signal decomposition approach.

Of course, the same principle applies to LRTFS. The synthesis coefficients αfn may equally be written as a sum of latent components, such that αfn = P k αkfn, with αkfn Nc(0, wfkhkn). Denoting by αk the column vector of dimension M with coefficients {αkfn}fn, Eq. (5) may be written as

k Φαk + e = X

k ck + e , (8)

where ck = Φαk. The component ck is the temporal expression of spectral pattern wk, the kth column of W. Given estimates of W and H, the components may be reconstructed in various way. The equivalent of the Wiener filtering approach used traditionally with the GCM would consist in computing ˆc MMSE k = ΦˆαMMSE k , with ˆαMMSE k = E{αk|x, W, H}. Though the expression of ˆαMMSE k is available in closed form it requires the inversion of a too large matrix, of dimensions T T (see also Section 3.2). We will instead use ˆck = Φˆαk with ˆαk = E{αk|ˆα, W, H}, where ˆα is the available estimate of α. In this case, the coefficients of ˆαk are given by

ˆαkfn = wfkhkn

[WH]fn ˆαfn. (9)

3 Estimation in LRTFS

We now consider two approaches to estimation of W, H and α in the LRTFS model defined by Eqs. (2)-(4). The first approach, described in the next section is maximum joint likelihood estimation (MJLE). It relies on the minimization of log p(x, α|W, H, λ). The second approach is maximum marginal likelihood estimation (MMLE), described in Section 3.2. It relies on the minimization of log p(x|W, H, λ), i.e., involves the marginalization of α from the joint likelihood, following the principle of SBL. Though we present MMLE for the sake of completeness, our current implementation does not scale with the dimensions involved in the audio signal processing applications presented in Section 5, and large-scale algorithms for MMLE are left as future work.

3.1 Maximum joint likelihood estimation (MJLE)

Objective. MJLE relies on the optimization of

CJL(α, W, H, λ) def = log p(x, α|W, H, λ) (10)

λ x Φα 2 2 + DIS(|α|2|v) + log(|α|2) + M log π , (11)

where we recall that v is the vectorized version of WH and where DIS(A|B) = P

ij d IS(aij|bij) is the IS divergence between nonnegative matrices (or vectors, as a special case), with d IS(x|y) = (x/y) log(x/y) 1. The first term in Eq. (11) measures the discrepancy between the raw signal and its approximation. The second term ensures that the synthesis coefficients are approximately low-rank. Unexpectedly, a third term that favors sparse solutions of α, thanks to the log function, naturally appears from the derivation of the joint likelihood. The objective function (11) is not convex and the EM algorithm described next may only ensure convergence to a local solution.

EM algorithm. In order to minimize CJL, we employ an EM algorithm based on the architecture proposed by Figueiredo & Nowak [6]. It consists of rewriting Eq. (5) as

β e1 , (12) x = Φz + e2 , (13)

where z acts as a hidden variable, e1 Nc(0, I), e2 Nc(0, λI βΦΦ ), with the operator denoting Hermitian transpose. Provided that β λ/δΦ, where δΦ is the largest eigenvalue of ΦΦ , the likelihood function p(x|α, λ) under Eqs. (12)-(13) is the same as under Eq. (5). Denoting the set of parameters by θJL = {α, W, H, λ}, the EM algorithm relies on the iterative minimization of

Q(θJL| θJL) = Z

z log p(x, α, z|W, H, λ)p(z|x, θJL)dz , (14)

where θJL acts as the current parameter value. Loosely speaking, the EM algorithm relies on the idea that if z was known, then the estimation of α and of the other parameters would boil down to the mere white noise denoising problem described by Eq. (12). As z is not known, the posterior mean value w.r.t z of the joint likelihood is considered instead.

The complete likelihood in Eq. (14) may be decomposed as

log p(x, α, z|W, H, λ) = log p(x|z, λ) + log p(z|α) + log p(α|WH). (15)

The hidden variable posterior simplifies to p(z|x, θJL) = p(z|x, λ). From there, using standard manipulations with Gaussian distributions, the (i + 1)th iteration of the resulting algorithm writes as follows.

E-step: z(i) = E{z|x, λ(i)} = α(i) + β λ(i) Φ (x Φα(i)) (16)

M-step: (f, n), α(i+1) fn = v(i) fn v(i) fn + β z(i) fn (17)

(W(i+1), H(i+1)) = arg min W,H 0

fn DIS |α(i+1) fn |2|[WH]fn (18)

T x Φα(i+1) 2 F (19)

In Eq. (17), v(i) fn is a shorthand for [W(i)H(i)]fn . Eq. (17) is simply the application of Wiener filtering to Eq. (12) with z = z(i). Eq. (18) amounts to solving a NMF with the IS divergence; it may be solved using majorization-minimization, resulting in the standard multiplicative update rules given in [3]. A local solution might only be obtained with this approach, but this is still decreasing the negative log-likelihood at every iteration. The update rule for λ is not the one that exactly derives from the EM procedure (this one has a more complicated expression), but it still decreases the negative log-likelihood at every iteration as explained in [6].

Note that the overall algorithm is rather computationally friendly as no matrix inversion is required. The Φα and Φ x operations in Eq. (16) correspond to analysis and synthesis operations that can be realized efficiently using optimized packages, such as the Large Time-Frequency Analysis Toolbox (LTFAT) [7].

3.2 Maximum marginal likelihood estimation (MMLE)

Objective. The second estimation method relies on the optimization of

CML(W, H, λ) def = log p(x|W, H, λ) (20)

α p(x|α, λ)p(α|WH)dα (21)

It corresponds to the type-II maximum likelihood procedure employed in [4, 5]. By treating α as a nuisance parameter, the number of parameters involved in the data likelihood is significantly reduced, yielding more robust estimation with fewer local minima in the objective function [5].

EM algorithm. In order to minimize CML, we may use the EM architecture described in [4, 5] that quite naturally uses α has the hidden data. Denoting the set of parameters by θML = {W, H, λ}, the EM algorithm relies on the iterative minimization of

Q(θML| θML) = Z

α log p(x, α|W, H, λ)p(α|x, θML)dα, (22)

where θML acts as the current parameter value. As the derivations closely follow [4, 5], we skip details for brevity. Using rather standard results about Gaussian distributions the (i + 1)th iteration of the algorithm writes as follows.

E-step : Σ(i) = (Φ Φ/λ(i) + diag(v(i 1)) 1) 1 (23)

α(i) = Σ(i)Φ x/λ(i) (24)

v(i) = E{|α|2|x, v(i), λ(i)} = diag(Σ(i)) + |α(i)|2 (25)

M-step : (W(i+1), H(i+1)) = arg min W,H 0

fn DIS v(i) fn|[WH]fn (26)

x Φα(i) 2 2 + λ(i) XM

m=1(1 Σ(i) mm/v(i) m ) (27)

The complexity of this algorithm can be problematic as it involves the computation of the inverse of a matrix of size M in the expression of Σ(i). M is typically at least twice larger than T, the signal length. Using the Woodbury matrix identity, the expression of Σ(i) can be reduced to the inversion of a matrix of size T, but this is still too large for most signal processing applications (e.g., 3 min of music sampled at CD quality makes T in the order of 106). As such, we will discard MMLE in the experiments of Section 5 but the methodology presented in this section can be relevant to other problems with smaller dimensions.

4 Multi-resolution LRTFS

Besides the advantage of modeling the raw signal itself, and not its STFT, another major strength of LRTFS is that it offers the possibility of multi-resolution modeling. The latter consists of representing a signal as a sum of t-f atoms with different temporal (and thus frequency) resolutions. This is for example relevant in audio where transients, such as the attacks of musical notes, are much shorter than sustained parts such as the tonal components (the steady, harmonic part of musical notes). Another example is speech where different classes of phonemes can have different resolutions. At even higher level, stationarity of female speech holds at shorter resolution than male speech. Because traditional spectral factorizations approaches work on the transformed data, the time resolution is set once for all at feature computation and cannot be adapted during decomposition.

In contrast, LRTFS can accommodate multiple t-f bases in the following way. Assume for simplicity that x is to be expanded on the union of two frames Φa and Φb, with common column size T

and with t-f grids of sizes Fa Na and Fb Nb, respectively. Φa may be for example a Gabor frame with short time resolution and Φb a Gabor frame with larger resolution such a setting has been considered in many audio applications, e.g., [8, 9], together with sparse synthesis coefficients models. The multi-resolution LRTFS model becomes

x = Φaαa + Φbαb + e (28)

(f, n) {1, . . . , Fa} {1, . . . , Na}, αa,fn Nc([Wa Ha]fn) , (29) (f, n) {1, . . . , Fb} {1, . . . , Nb}, αb,fn Nc([Wb Hb]fn) , (30)

and where {αa,fn}fn and {αb,fn}fn are the coefficients of αa and αb, respectively.

By stacking the bases and synthesis coefficients into Φ = [Φa Φb] and α = [αT a αT b ]T

and introducing a latent variable z = [z T a z T b ]T , the negative joint log-likelihood log p(x, α|Wa, Ha, Wb, Hb, λ) in the multi-resolution LRTFS model can be optimized using the EM algorithm described in Section 3.1. The resulting algorithm at iteration (i + 1) writes as follows.

E-step: for ℓ= {a, b}, z(i) ℓ = α(i) ℓ + β

λΦ ℓ(x Φaα(i) a Φbα(i) b ) (31)

M-step: for ℓ= {a, b}, (f, n) {1, . . . , Fℓ} {1, . . . , Nℓ}, α(i+1) ℓ,fn = v(i) ℓ,fn v(i) ℓ,fn + β z(i) fn (32)

for ℓ= {a, b}, (W(i+1) ℓ , H(i+1) ℓ ) = arg min Wℓ,Hℓ 0

fn DIS |α(i+1) ℓ,fn |2|[WℓHℓ]fn (33)

λ(i+1) = x Φaα(i+1) a Φbα(i+1) b 2 2/T (34)

The complexity of the algorithm remains fully compatible with signal processing applications. Of course, the proposed setting can be extended to more than two bases.

5 Experiments

We illustrate the effectiveness of our approach on two experiments. The first one, purely illustrative, decomposes a jazz excerpt into two layers (tonal and transient), plus a residual layer, according to the hybrid/morphological model presented in [8, 10]. The second one is a speech enhancement problem, based on a semi-supervised source separation approach in the spirit of [11]. Even though we provided update rules for λ for the sake of completeness, this parameter was not estimated in our experiments, but instead treated as an hyperparameter, like in [5, 6]. Indeed, the estimation of λ with all the other parameters free was found to perform poorly in practice, a phenomenon observed with SBL as well.

5.1 Hybrid decomposition of music

We consider a 6 s jazz excerpt sampled at 44.1 k Hz corrupted with additive white Gaussian noise with 20 d B input Signal to Noise Ratio (SNR). The hybrid model aims to decompose the signal as

x = xtonal + xtransient + e = Φtonalαtonal + Φtransientαtransient + e , (35)

using the multi-resolution LRTFS method described in Section 4. As already mentionned, a classical design consists of working with Gabor frames. We use a 2048 samples-long ( 46 ms) Hann window for the tonal layer, and a 128 samples-long ( 3 ms) Hann window for the transient layer, both with a 50% time overlap. The number of latent components in the two layers is set to K = 3.

We experimented several values for the hyperparameter λ and selected the results leading to best output SNR (about 26 d B). The estimated components are shown at Fig. 1. When listening to the signal components (available in the supplementary material), one can identify the hit-hat in the first and second components of the transient layer, and the bass and piano attacks in the third component. In the tonal layer, one can identify the bass and some piano in the first component, some piano in the second component, and some hit-hat ring in the third component.

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

0 1 2 3 4 5 0

Figure 1: Top: spectrogram of the original signal (left), estimated transient coefficients log |αtransient| (center), estimated tonal coefficients log |αtonal| (right). Middle: the 3 latent components (of rank 1) from the transient layer. Bottom: the 3 latent components (of rank 1) from the tonal layer.

5.2 Speech enhancement

The second experiment considers a semi-supervised speech enhancement example (treated as a single-channel source separation problem). The goal is to recover a speech signal corrupted by a texture sound, namely applauses. The synthesis model considered is given by

x = Φtonal αspeech tonal + αnoise tonal + Φtransient αspeech transient + αnoise transient + e, (36)

with αspeech tonal Nc 0, Wtrain tonal Hspeech tonal , αnoise tonal Nc 0, Wnoise tonal Hnoise tonal , (37)

and αspeech transient Nc 0, Wtrain transient Hspeech transient , αnoise transient Nc 0, Wnoise transient Hnoise transient . (38)

Wtrain tonal and Wtrain transient are fixed pre-trained dictionaries of dimension K = 500, obtained from 30 min of training speech containing male and female speakers. The training data, with sampling rate 16k Hz, is extracted from the TIMIT database [12]. The noise dictionaries Wnoise tonal and Wnoise transient are learnt from the noisy data, using K = 2. The two t-f bases are Gabor frames with Hann window of length 512 samples ( 32 ms) for the tonal layer and 32 samples ( 2 ms) for the transient layer, both with 50% overlap. The hyperparameter λ is gradually decreased to a negligible value during iterations (resulting in a negligible residual e), a form of warm-restart strategy [13].

We considered 10 test signals composed of 10 different speech excerpts (from the TIMIT dataset as well, among excerpts not used for training) mixed in the middle of a 7 s-long applause sample. For every test signal, the estimated speech signal is computed as

ˆx = Φtonal ˆαspeech tonal + Φtransient ˆαspeech transient (39)

Noisy signal: long window STFT analysis

0 1 2 3 4 5 6 7 0

Noisy signal: short window STFT analysis

0 1 2 3 4 5 6 7 0

Denoised signal: Tonal Layer

0 1 2 3 4 5 6 7 0

Denoised signal: Transient Layer

0 1 2 3 4 5 6 7 0

Figure 2: Time-frequency representations of the noisy data (top) and of the estimated tonal and transient layers from the speech (bottom).

and a SNR improvement is computed as the difference between the output and input SNRs. With our approach, the average SNR improvement other the 10 test signals was 6.6 d B. Fig. 2 displays the spectrograms of one noisy test signal with short and long windows, and the clean speech synthesis coefficients estimated in the two layers. As a baseline, we applied IS-NMF in a similar setting using one Gabor transform with a window of intermediate length (256 samples, 16 ms). The average SNR improvement was 6 d B in that case. We also applied the standard OMLSA speech enhancement method [14] (using the implementation available from the author with default parameters) and the average SNR improvement was 4.6 d B with this approach. Other experiments with other noise types (such as helicopter and train sounds) gave similar trends of results. Sound examples are provided in the supplementary material.

6 Conclusion

We have presented a new model that bridges the gap between t-f synthesis and traditional NMF approaches. The proposed algorithm for maximum joint likelihood estimation of the synthesis coefficients and their low-rank variance can be viewed as an iterative shrinkage algorithm with an additional Itakura-Saito NMF penalty term. In [15], Elad explains in the context of sparse representations that soft thresholding of analysis coefficients corresponds to the first iteration of the forwardbackward algorithm for LASSO/basis pursuit denoising. Similarly, Itakura-Saito NMF followed by Wiener filtering correspond to the first iteration of the proposed EM algorithm for MJLE.

As opposed to traditional NMF, LRTFS accommodates multi-resolution representations very naturally, with no extra difficulty at the estimation level. The model can be extended in a straightforward manner to various additional penalties on the matrices W or H (such as smoothness or sparsity). Future work will include the design of a scalable algorithm for MMLE, using for example message passing [16], and a comparison of MJLE and MMLE for LRTFS. Moreover, our generative model can be considered for more general inverse problems such as multichannel audio source separation [17]. More extensive experimental studies are planned in this direction.

Acknowledgments

The authors are grateful to the organizers of the Modern Methods of Time-Frequency Analysis Semester held at the Erwin Schr oedinger Institute in Vienna in December 2012, for arranging a very stimulating event where the presented work was initiated.

[1] P. Smaragdis, C. F evotte, G. Mysore, N. Mohammadiha, and M. Hoffman. Static and dynamic source separation using nonnegative factorizations: A unified view. IEEE Signal Processing Magazine, 31(3):66 75, May 2014.

[2] D. D. Lee and H. S. Seung. Learning the parts of objects with nonnegative matrix factorization. Nature, 401:788 791, 1999.

[3] C. F evotte, N. Bertin, and J.-L. Durrieu. Nonnegative matrix factorization with the Itakura Saito divergence. With application to music analysis. Neural Computation, 21(3):793 830, Mar. 2009.

[4] M. E. Tipping. Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1:211 244, 2001.

[5] D. P. Wipf and B. D. Rao. Sparse bayesian learning for basis selection. IEEE Transactions on Signal Processing, 52(8):2153 2164, Aug. 2004.

[6] M. Figueiredo and R. Nowak. An EM algorithm for wavelet-based image restoration. IEEE Transactions on Image Processing, 12(8):906 916, Aug. 2003.

[7] Z. Pr uˇsa, P. Søndergaard, P. Balazs, and N. Holighaus. LTFAT: A Matlab/Octave toolbox for sound processing. In Proc. 10th International Symposium on Computer Music Multidisciplinary Research (CMMR), pages 299 314, Marseille, France, Oct. 2013.

[8] L. Daudet and B. Torr esani. Hybrid representations for audiophonic signal encoding. Signal Processing, 82(11):1595 1617, 2002.

[9] M. Kowalski and B. Torr esani. Sparsity and persistence: mixed norms provide simple signal models with dependent coefficients. Signal, Image and Video Processing, 3(3):251 264, 2009.

[10] M. Elad, J.-L. Starck, D. L. Donoho, and P. Querre. Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA). Journal on Applied and Computational Harmonic Analysis, 19:340 358, Nov. 2005.

[11] P. Smaragdis, B. Raj, and M. V. Shashanka. Supervised and semi-supervised separation of sounds from single-channel mixtures. In Proc. 7th International Conference on Independent Component Analysis and Signal Separation (ICA), London, UK, Sep. 2007.

[12] TIMIT: acoustic-phonetic continuous speech corpus. Linguistic Data Consortium, 1993.

[13] A. Hale, W. Yin, and Y. Zhang. Fixed-point continuation for ℓ1-minimization: Methodology and convergence. SIAM Journal on Optimisation, 19(3):1107 1130, 2008. [14] I. Cohen. Noise spectrum estimation in adverse environments: Improved minima controlled recursive averaging. IEEE Transactions on Speech and Audio Processing, 11(5):466 475, 2003.

[15] M. Elad. Why simple shrinkage is still relevant for redundant representations? IEEE Transactions on Information Theory, 52(12):5559 5569, 2006.

[16] M. W. Seeger. Bayesian inference and optimal design for the sparse linear model. The Journal of Machine Learning Research, 9:759 813, 2008.

[17] A. Ozerov and C. F evotte. Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation. IEEE Transactions on Audio, Speech and Language Processing, 18(3):550 563, Mar. 2010.