Roland Badeau

Mark D Plumbley

MULTICHANNEL HR-NMF FOR MODELLING CONVOLUTIVE MIXTURES OF NON-STATIONARY SIGNALS IN THE TIME-FREQUENCY DOMAIN

Keywords: Non-stationary signal modelling, Timefrequency analysis, Separation of convolutive mixtures, Multichannel signal analysis, Variational EM algorithm

In the literature, several probabilistic models involving latent components have been proposed for modelling time-frequency (TF) representations of audio signals (such as spectrograms), notably in the nonnegative matrix factorization (NMF) literature. Among them, the recent high resolution (HR)-NMF model is able to take both phases and local correlations in each frequency band into account, and its potential has been illustrated in applications such as source separation and audio inpainting. In this paper, HR-NMF is extended to multichannel signals and to convolutive mixtures. A fast variational expectation-maximization (EM) algorithm is proposed to estimate the enhanced model. This algorithm is applied to a stereophonic piano signal, and proves capable of accurately modelling reverberation and restoring missing observations.

INTRODUCTION

Following previous works which aimed to provide a probabilistic framework to NMF [START_REF] Schmidt | Non-negative matrix factorization with Gaussian process priors[END_REF][START_REF] Smaragdis | Blind Speech Separation[END_REF][START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF][START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], the HR-NMF model introduced in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF]High resolution NMF for modeling mixtures of nonstationary signals in the time-frequency domain[END_REF] offers an improved frequency resolution (it can separate sinusoids within the same frequency band) and an improved synthesis capability (it can restore missing TF observations). It is suitable for both complex-valued and real-valued TF representations, such as the short-time Fourier transform (STFT) and the modified discrete cosine transform (MDCT). It also generalizes some popular models, such as the Itakura-Saito NMF model (IS-NMF) [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], autoregressive (AR) processes [START_REF] Hayes | Statistical Digital Signal Processing And Modeling[END_REF], and the exponential sinusoidal model (ESM), commonly used in HR spectral analysis of time series [START_REF] Hayes | Statistical Digital Signal Processing And Modeling[END_REF].

In this paper, HR-NMF is extended to multichannel signals and to convolutive mixtures. Contrary to the multichannel NMF [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] where convolution was approximated, convolution is here accurately modelled in the TF domain by following the exact implementation proposed in [START_REF] Badeau | Probabilistic time-frequency source-filter decomposition of non-stationary signals[END_REF]. Consequently, correlations within and between frequency bands are both taken into account. In order to estimate this multichannel HR-NMF model, we propose a variational EM algorithm, whose novelties compared to [START_REF] Badeau | Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF] include a reduced computational complexity and a parallel implementation.

The paper is structured as follows. The multichannel HR-NMF model is introduced in section 2, and the variational EM algorithm is derived in section 3. An application to a stereophonic piano signal is presented in section 4. Finally, conclusions are drawn in section 5.

MULTICHANNEL HR-NMF MODEL

The multichannel HR-NMF model of TF data ym(f, t) ∈ F (where

F = R or C) is defined for all channels m ∈ [0 ...M -1], discrete frequencies f ∈ [0 ...F -1], and times t ∈ [0 ...T -1], as the sum of S components yms(f, t) ∈ F plus a white noise nm(f, t) ∼ N F (0,σ 2 y ), where N F (0,σ 2 y ) denotes a real (if F = R) or circular complex (if F = C) normal distribution of mean 0 and variance σ 2 y : ym(f, t)=nm(f, t)+ S-1 s=0 yms(f, t). (1)
Each component yms(f, t) for any s ∈ [0 ...S -1] is defined as

yms(f, t)= P b ϕ=-P b Q b τ =0 bms(f, ϕ, τ) zs(f -ϕ, t -τ )
where bms(f, ϕ, τ)=0if fϕ/ ∈ [0 ...F -1], and the latent component zs(f, t) is defined as follows:

•∀ t ∈ [0 ...T -1], xs(f, t) ∼N F (0,σ 2 xs (t)) and zs(f, t)=xs(f, t) - Qa τ =1 as(f, τ)zs(f, t -τ )
where as(f, τ) defines a causal and stable autoregressive filter,

•∀ t ∈ [-Qz ... -1] where Qz =m a x ( Q b ,Qa), zs(f, t) ∼ N (µs(f, t), 1/ρs(f, t)).
Moreover, all scalar random variables nm(f, t) and xs(f, t)

are assumed independent. Besides, ∀m ∈ [0 ...M -1], ∀f ∈ [0 ...F -1], ∀t ∈ [-Qz ... -1], ym(f, t
) is unobserved, and ∀s ∈ [0 ...S -1], the prior mean µs(f, t) and the prior precision (inverse variance) ρs(f, t) of the latent variable zs(f, t) are considered as known (fixed) parameters.

The set θ of parameters to be estimated consists of:

• the autoregressive parameters as(f, τ) ∈ F for s ∈ [0 ...S-1], τ ∈ [1 ...Qa] and f ∈ [0 ...F -1] (we further define as(f, 0) = 1), • the moving average parameters bms(f, ϕ, τ)

∈ F for m ∈ [0 ...M -1], s ∈ [0 ...S -1], f ∈ [0 ...F -1], ϕ ∈ [-P b ...P b ], and τ ∈ [0 ...Q b ],
• the variance parameters σ 2 y > 0 and σ 2 xs (t) > 0 for s ∈ [0 ...S -1] and t ∈ [0 ..

.T -1].

This model encompasses the following particular cases:

• If M =1 , σ 2 y =0and P b = Q b = Qa =0 , equation (1)
becomes y0(f, t)= S-1 s=0 b0s(f, 0, 0)xs(f, t), thus y0(f, t) ∼ N F (0, V ft ), where V is defined by the NMF V = WHwith W fs = |b0s(f, 0, 0)| 2 and Hst = σ 2 xs (t). The maximum likelihood estimation of W and H is then equivalent to the minimization of the Itakura-Saito (IS) divergence between matrix V and spectrogram V (where V ft = |y0(f, t)| 2 ), hence this model is referred to as IS-NMF [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF].

• If M =1and P b = Q b =0 , y0(f, t) follows the HR-NMF
model [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF]High resolution NMF for modeling mixtures of nonstationary signals in the time-frequency domain[END_REF][START_REF] Badeau | Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF] involving variance σ 2 y , autoregressive parameters as(f, τ), and same NMF V = WH.

• If S =1 , σ 2 y =0 , P b =0 , ∀t ∈ [0 ...T -1], σ 2 x 0 (t)=1 , and ∀t ∈ [-Qz ... -1], µs(f, t)=0and ρs(f, t)) = 1, then ∀m ∈ [0 ...M -1], ∀f ∈ [0 ...F -1], ym(f, t
) is an autoregressive moving average (ARMA) process [START_REF] Hayes | Statistical Digital Signal Processing And Modeling[END_REF].

• If S =1 , σ 2 y =0 , P b =0 , Qa > 0, Q b = Qa -1, ∀t ∈ [-Qz ... -1], µ0(f, t)=0 , ρ0(f, t) → +∞, and σ 2 0 (t)= {t=0} , then ∀m ∈ [0 ...M -1], ∀f ∈ [0 ...F -1], ym(f, t) can be written in the form ym(f, t)= Qa τ =1 αmτ z t
τ , where z1 ...zQ a are the roots of the polynomial z Qa -Qa τ =1 as(f, τ)z Qa-τ . This corresponds to the Exponential Sinusoidal Model (ESM) commonly used in HR spectral analysis of time series [START_REF] Hayes | Statistical Digital Signal Processing And Modeling[END_REF].

VARIATIONAL EM ALGORITHM

In order to estimate the multichannel HR-NMF model introduced in section 2, we derive below a variational EM algorithm.

Review of variational EM algorithm

Variational inference [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF] is now a classical approach for estimating a probabilistic model involving both observed variables y and latent variables z, parametrized by θ. Let F be a set of probability density functions (PDF) over the latent variables z. For any PDF q ∈F and any function f (z), we note f q = f (z)q(z)dz. Then for any parameter θ, the variational free energy is defined as

L(q; θ)= ln p(y, z; θ) q(z) q (2)
The variational EM algorithm is a recursive algorithm for estimating θ. It consists of the two following steps at each iteration i:

• Expectation (E)-step (update q): q ⋆ = argmax q∈F L(q; θi-1) (3)
• Maximization (E)-step (update θ):

θi = argmax θ L(q ⋆ ; θ). (4)
In the case of multichannel HR-NMF, θ has been defined in section 2, and we define δm(f,

t)=1if ym(f, t) is observed, other- wise δm(f, t)=0(in particular δm(f, t)=0∀(f, t) / ∈ [0 ...F - 1] × [0 ...T -1]
), so that the complete data consists of:

• the set y of observed variables y m(f, t) for m ∈ [0 ...M-1]
and for all f and t such that δm(f, t)=1,

• the set z of latent variables zs(f, t) for s ∈ [0 ...S -1], f ∈ [0 ...F -1], t ∈ [-Qz ...T -1].
Besides, we use a mean field approximation [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]: F is defined as the set of PDFs which can be factorized in the form

q(z)= S-1 s=0 F -1 f =0 T -1 t=-Qz q sf t (zs(f, t)). (5)
With this particular factorization of q(z), the solution of ( 3) is such that each PDF q sf t is Gaussian: zs(f, t) ∼N F (zs(f, t),γz s (f, t)).

Variational free energy

Let α =1if F = C, and α =2if F = R. The joint probability distribution of the complete data can be written as

-α ln(p(y, z; θ)) = (D + SF(T + Qz)) ln(απ) +D ln(σ 2 y )+ 1 σ 2 y M -1 m=0 F -1 f =0 T -1 t=0 |ey m (f, t)| 2 + S-1 s=0 F -1 f =0 -1 t=-Qz -ln(ρs(f, t)) + ρs(f, t)|zs(f, t) -µs(f, t)| 2 + S-1 s=0 F -1 f =0 T -1 t=0 ln(σ 2 xs (t)) + 1 σ 2 xs (t) |ex s (f, t)| 2
where we have introduced

D = M -1 m=0 F -1 f =0 T -1 t=0 δm(f, t) I(f, t)= 0≤f<F 0≤t<T
ey m (f, t)=δm(f, t) ym(f, t) -

S-1 s=0 yms(f, t) ex s (f, t)=I(f, t) Qa τ =0
as(f, τ)zs(f, tτ )

Then the variational free energy defined in (2) becomes -αL(q; θ)=D ln(απ) -SF(T + Qz)

+D ln(σ 2 y )+ M -1 m=0 F -1 f =0 T -1 t=0 γe ym (f,t)+|ey m (f,t)| 2 σ 2 y + S-1 s=0 F -1 f =0 -1 t=-Qz -ln(ρs(f, t)γz s (f, t)) +ρs(f, t) γz s (f, t)+|zs(f, t) -µs(f, t)| 2 + S-1 s=0 F -1 f =0 T -1 t=0 ln(σ 2 xs (t)) -ln(γz s (f, t)) + γx s (f,t)+|ex s (f,t)| 2 σ 2 xs (t) where ∀f ∈ [0 ...F -1], ∀t ∈ [0 ...T -1], γe ym (f, t)=δm(f, t) S-1 s=0 P b ϕ=-P b Q b τ =0 |bms(f, ϕ, τ)| 2 γz s (f -ϕ, t -τ ) ey m (f, t)=δm(f, t) ym(f, t) - S-1 s=0 y ms (f, t) y ms (f, t)= P b ϕ=-P b Q b τ =0 bms(f, ϕ, τ) zs(f -ϕ, t -τ ) γx s (f, t)=I(f, t) Qa τ =0 |as(f, τ)| 2 γz s (f, t -τ ) ex s (f, t)=I(f, t) Qa τ =0
as(f, τ)zs(f, tτ )

Variational EM algorithm for multichannel HR-NMF

In this section, the maximizations in equations ( 3) and ( 4) are performed for each scalar parameter in turn. The resulting dominant complexity of each iteration of this variational EM algorithm is 4MFST∆F ∆T , where ∆f =1+2P b , and ∆t =1+Qz. However we highlight a possible parallel implementation, by making a difference between parfor loops which can be implemented in parallel, and for loops which have to be implemented sequentially.

E-step

For all s ∈ [0 ...S -1], f ∈ [0 ...F -1], t/ ∈ [-Qz, -1], let ρs(f, t)=0.
Considering the mean field approximation ( 5), the Estep defined in (3) leads to the updates described in Table 1 (where * denotes complex conjugation). Note that the updates of γz s (f, t) and zs(f, t) have to be processed sequentially.

M-step

The M-step defined in (4) leads to the updates described in Table 2.

The updates of the four parameters can be processed in parallel.

parfor s ∈ [0 ...S-1], f ∈ [0 ...F-1], t ∈ [-Qz ...T-1] do γz s (f, t) -1 = ρs(f, t)+ Qa τ =0 |as(f, τ)| 2 I(f,t+τ ) σ 2 xs (t+τ ) + M -1 m=0 P b ϕ=-P b Q b τ =0 |bms(f + ϕ, ϕ, τ )| 2 δm(f +ϕ,t+τ )

SIMULATION RESULTS

In this section, we present a basic proof of concept of the proposed multichannel HR-NMF model. The following experiments deal with a unique source (S =1 ) formed of a real piano sound sampled at 11025 Hz. A stereophonic signal (M =2 ) of length 1.25 ms has been synthesized by filtering the monophonic recording of a C3 piano note with two simulated room impulse responses. The TF representation y m(f, t) of this signal has then been computed by applying a critically sampled perfect reconstruction cosine modulated filterbank (F = R) with F =2 0 1frequency bands, involving filters of length 8F =1 6 0 8samples. The resulting TF representation, of dimensions F × T with T =7 7 , is displayed in Figure 1. In this configuration, the interchannel time difference is greater than the length of time frames, which suggests that the length Q b of filters bms(f, ϕ, τ) should be chosen greater than zero. In the following experiments, we have set µs(f, t)=0and ρs(f, t)=1 0 5 (which means that the observed sound is preceded by silence). The variational EM algorithm is initialized with zs(f, t)=0 , γz s (f, t)=σ 2 y = σ 2 xs (t)=1 , as(f, τ)= {τ =0} , and bms(f, ϕ, τ)= {ϕ=0,τ =0} . In order to illustrate the capability of the multichannel HR-NMF model to synthesize realistic audio data, we address the case of missing observations: all TF points within the red frame in Figure 1 In each experiment, 100 iterations of the algorithm are performed, and the restored signal is returned as y ms (f, t).

σ 2 y = 1 D M -1 m=0 F -1 f =0 T -1 t=0 γe ym (f, t)+|ey m (f, t)| 2 parfor s ∈ [0 ...S -1], t ∈ [0 ...T -1] do σ 2 xs (t)= 1 F F -1 f =0 γx s (f, t)+|ex s (f, t)| 2 end parfor for τ ∈ [1 ...Qa] do parfor s ∈ [0 ...S -1], f ∈ [0 ...F -1] do as(f, τ)= T -1 t=0 1 σ 2 xs (t) (zs(f,t-τ) * (as(f,τ)zs(f,t-τ)-ex s (f,t))) T -1 t=0 1 σ 2 xs (t) (γz s (f,t-τ )+|zs(f,t-τ )| 2 ) end parfor end for for s ∈ [0 ...S -1], ϕ ∈ [-P b ...P b ], τ ∈ [0 ...Q b ] do parfor m∈[0 ...M-1], f ∈[max(0,ϕ) ...F-1+min(0,ϕ)] do bms(f, ϕ, τ)= T -1 t=0 zs(f -ϕ,t-τ) * (δm(f,t)bms(f,ϕ,τ)zs(f-ϕ,t-τ)+ey m (f,t)) T -1 t=0 δm(f,t)(γz s (f -ϕ,t-τ)+|zs(f-ϕ,t-τ)| 2 ) end parfor end for Table 2: M-step
In the first experiment, a multichannel HR-NMF with Qa = Q b = P b =0is estimated. As highlighted in section 2, this is equivalent to modeling the two channels by two IS-NMF models [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] having distinct spectral atoms and sharing the same temporal activation, or by a multichannel NMF of rank 1 [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. The resulting TF representation y ms (f, t) is displayed in Figure 2. It can be noticed that wherever ym(f, t) is observed (δm(f, t)=1 ), y ms (f, t) does not accurately fit ym(f, t), because the length Q b of filters bms(f, ϕ, τ) has been underestimated. In other respects, the missing observations (δm(f, t)=0) could not be restored (y ms (f, t) is zero inside the red frame), because the correlations between contiguous TF coefficients in ym(f, t) have not been taken into account. In the second experiment, a multichannel HR-NMF model with Qa =2 , Q b =3 , and P b =1is estimated. The resulting TF representation y ms (f, t) is displayed in Figure 3. It can be noticed that wherever ym(f, t) is observed, y ms (f, t) better fits ym(f, t). Besides, the missing observations have been correctly estimated.

CONCLUSIONS

In this paper, the HR-NMF model [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF]High resolution NMF for modeling mixtures of nonstationary signals in the time-frequency domain[END_REF] has been extended to multichannel signals and to convolutive mixtures. The new multichannel HR-NMF model accurately represents convolution in the TF domain [START_REF] Badeau | Probabilistic time-frequency source-filter decomposition of non-stationary signals[END_REF], and also takes the correlations over frequencies into account. In order to estimate this model from real audio data, a variational EM algorithm has been proposed, whose novelties compared to [START_REF] Badeau | Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF] include a reduced computational complexity and a parallel implementation. This algorithm has been successfully applied to a stereophonic piano signal, and has been capable of accurately modelling reverberation and restoring missing observations.

In future work, the sparsity of the model parameters in the TF domain could be enforced by using sparse Baysian learning [START_REF] Wipf | Sparse bayesian learning for basis selection[END_REF]. Some other desirable properties such as harmonicity and temporal or spectral smoothness could also be enforced by introducing some prior distributions of the parameters. Similarly to the high spectral resolution, a high temporal resolution could be achieved by extending the model as proposed in [START_REF] Badeau | Probabilistic time-frequency source-filter decomposition of non-stationary signals[END_REF]. Other Bayesian estimation techniques such as Markov chain Monte Carlo (MCMC) methods and message passing algorithms [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF] might prove more effective than the variational EM algorithm. Lastly, the proposed approach could be used in a variety of applications, such as source separation, source coding, audio inpainting, and automatic music transcription.

  s∈[0 ...S-1], f0∈[0 ...∆F -1], t0∈[-Qz ...-Qz+∆T -1] do parforf -f 0 ∆f ∈ [0 ...⌊ F -1-f 0 ∆ f ⌋], t-t 0 ∆t ∈ [0 ...⌊ T -1-t 0 ∆t ⌋] do zs(f, t)=zs(f, t)γz s (f, t) ρs(f, t)(zs(f, t) -µs(f, t)) + Qa τ =0as(f, τ) * ex s (f,t+τ ) + ϕ, ϕ, τ ) * ey m (f +ϕ,t+τ )

Figure 1 :

 1 Figure 1: Input stereo signal ym(f, t).

  are unobserved (δm(f, t)=0 ∀t ∈ [26 ...50], and δm(f, t)=1for all other t in [0 ...T -1]).

Figure 2 :

 2 Figure 2: Stereo signal y ms (f, t) estimated with filters of length 1.

Figure 3 :

 3 Figure 3: Stereo signal y ms (f, t) estimated with longer filters.

Table 1 :

 1 E-step

October 20-23, 2013, New Paltz, NY

* This work was undertaken while Roland Badeau was visiting the Centre for Digital Music, partly funded by EPSRC Platform Grant EP/K009559/1, and by the French National Research Agency (ANR) as a part of the DReaM project (ANR-09-CORD-006-03).

† Mark D. Plumbley is funded by EPSRC Leadership Fellowship EP/G007144/1.


Romain Hennequin

Roland Badeau

Bertrand David

SCALE-INVARIANT PROBABILISTIC LATENT COMPONENT ANALYSIS

Keywords: Music signal processing, non-negative matrix factorization, probabilistic latent component analysis, shiftinvariant decomposition

In this paper, we present a new method for decomposing musical spectrograms. This method is similar to shift-invariant Probabilistic Latent Component Analysis, but, when the latter works with constant Q spectrograms (i.e. with a logarithmic frequency resolution), our technique is designed to decompose standard short time Fourier transform spectrograms (i.e. with a linear frequency resolution). This makes it possible to easily reconstruct the latent signals (which can be useful for source separation).

INTRODUCTION

Non-negative decompositions are widely used for audio spectrograms processing: non-negative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] and PLCA [START_REF] Shashanka | Sparse overcomplete latent variable decomposition of counts data[END_REF] are both used to decompose spectrograms with applications such as source separation [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF][START_REF] Smaragdis | Supervised and semisupervised separation of sounds from single-channel mixtures[END_REF] and automatic music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF]. Shift-invariant models [START_REF] Schmidt | Nonnegative matrix factor 2-D deconvolution for blind single channel source separation[END_REF][START_REF] Smaragdis | Sparse and shift-invariant feature extraction from non-negative data[END_REF] decompose constant-Q spectrograms [START_REF] Brown | Calculation of a constant Q spectral transform[END_REF] with a single frequency template for each harmonic instrument: with a log-frequency resolution, a frequency shift corresponds to a transposition. Then each note of a single instrument can be modeled as a base template shifted to the right pitch.

Unfortunately, constant-Q transforms (CQT) are difficult to use for sound source separation: even if a near-perfect inverse CQT transform was recently proposed [START_REF] Schörkhuber | Constant-Q transform toolbox for music processing[END_REF], the variable resolution of the decomposition makes it difficult to apply time-frequency masking and source separation using CQT is still an open problem. Attempts were made to use shift-invariant decomposition of CQT for source separation [START_REF] Fitzgerald | Shifted non-negative matrix factorisation for sound source separation[END_REF] using a mapping between log-frequency and linear frequency resolution to avoid CQT inversion. This paper presents a new decomposition method inspired by shift-invariant decompositions, but which is designed to directly decompose STFT spectrograms. In a constant-Q spectrogram, a change of fundamental frequency approximately corresponds to a translation of the spectral template; in a standard STFT spectrogram, one can model such a change with an homothety of the spectral template. This approximation remains only valid for small modifications since:

• For real acoustic instruments, a note change usually leads to a different energy balance between the different harmonics.

Even if it is widely used in Shift-Invariant decomposition, the assumption that this balance does not change is only exact for a few electronic instruments. • The spectral width of harmonics (or partials) is affected by the homothety, in opposition with the real world case where it is mostly related to the analysis window type and length.

We call the new decomposition scale-invariant PLCA (SIPLCA). This scale-invariant model presents some new issues that were not encountered with shift-invariant models: an homothety on a set of integers does not yield integers. We propose a solution to this issue. In section 2, we remind the principle of standard PLCA and shift-invariant PLCA. We then present, in section 3, the new scaleinvariant model and derive an algorithm to estimate the parameters. Some examples are presented in section 4 and we propose an application of single notes repitching in a polyphonic signal. Finally, we draw conclusions in section 5.

PROBABILISTIC LATENT COMPONENT ANALYSIS

The model that we used is inspired by shift-invariant PLCA [START_REF] Smaragdis | Sparse and shift-invariant feature extraction from non-negative data[END_REF] which is a probabilistic drawing model. In PLCA decompositions [START_REF] Shashanka | Sparse overcomplete latent variable decomposition of counts data[END_REF], a non-negative spectrogram V f t is considered as an histogram obtained from a structured draw of a frequency random variable f ∈ {1, 2, . . . , F } and a time random variable t ∈ {1, 2, . . . , T } which follow the joint distribution P (f, t). One can design P (f, t) in different ways, which lead to different decompositions.

Standard PLCA

Standard PLCA (non shift-invariant) [START_REF] Shashanka | Sparse overcomplete latent variable decomposition of counts data[END_REF] leads to a decomposition similar to NMF. The draw is structured with a latent (hidden) random variable z which corresponds to a "component", assuming that f and t are independent given z:

P (f, t) = Z z=1 P (z)P (f, t|z) = Z z=1 P (z)P (f |z)P (t|z).
The histogram V f t is thus assumed to be obtained in the following way: first z is drawn following P (z) and then f and t are drawn following respectively P (f |z) and P (t|z).

In an NMF framework, P (f |z) corresponds to the spectral templates and P (t|z) corresponds to the activations of each component. P (z) is the relative weight of each component.

Shift-invariant PLCA

Shift-invariant PLCA introduces another latent random variable f ′ ∈ {1, 2, . . . , F ′ }, the base frequency, in reference of which a transposition variable τ ∈ Z is defined. f ′ and t are assumed independent given z, and τ and f ′ are also independent given z (but t and τ are not). f is obtained by a transposition of the base template: f = f ′ + τ . P (f, t) then takes the form:

P (f, t) = Z z=1 P (z) F ′ f ′ =1 PK (f ′ |z)PI (f -f ′ , t|z).
PK is called the kernel distribution: it corresponds to the base spectral templates which are shifted by the impulse distribution PI .

SCALE-INVARIANT PLCA

Model

In standard short-time Fourier transform spectrograms (with a linear frequency resolution), transposition is no longer a shift but approximately corresponds to an homothety of the spectral template: we thus model it with a multiplication by a scalar λ ∈ R + \{0}. Let X be a discrete random variable taking its values in {0, 1, 2, ...K} and λ a continuous positive random variable with a density function p. The density of u = λX is then:

p λX (u) = K k=1 P (X = k)p( u k ) k + δ(u)P (X = 0) (1)
where δ(u) is a Dirac delta function.

In our model, one assumes that the continuous frequency random variable fc ∈ R (the link between fc and the observed variable f will be made clearer later) is obtained by multiplying the base frequency f ′ ∈ {0, 1, . . . , F ′ } (which is independent of t given z) with the transposition factor λ ∈ R + \{0} (which depends on t but not on f ′ given z). Using (1) with u = fc, k = f ′ and K = F ′ , we get:

P (fc, t|z) = F ′ f ′ =1 PK (f ′ |z) f ′ PI fc f ′ , t|z + δ(fc)PK (0|z).
We use the notation PK for the kernel distribution and PI for the impulse distribution, as for shift-invariant PLCA. However they do not represent the same object. PK still corresponds to the base spectral templates but is now rescaled by the impulse distribution PI .

In this paper, we consider that PK (0|z) = 0 to avoid the singularity at the null frequency (there still can be energy in the frequency channel 0 introduced when scaling down the frequency template). We then get:

P (fc, t) = Z z=1 P (z) F ′ f ′ =1 PK (f ′ |z) f ′ PI fc f ′ , t|z .
The random variable fc is continuous, but the observed random variable f ∈ {0, 1, . . . , F } is discrete. Then we will suppose that f = round(fc) and consequently:

P (f, t) = f + 1 2 f -1 2 P (fc, t)dfc.
Consequently, ∀f ∈ {0, 1, . . . , F }, ∀t ∈ {1, . . . , T }:

P (f, t) = z,f ′ P (z)PK (f ′ |z) f ′ f + 1 2 f -1 2 PI fc f ′ , t|z dfc.
The parameters to be estimated are then:

θ = {P (z), PK (f ′ |z), PI (λ, t|z)}.
For practical purposes, one needs to discretize PI (which is a continuous density function with respect to λ) in some way, in order to estimate θ. We propose to perform this discretization by parameterizing PI assuming that λ → PI (λ, t|z) is piece-wise constant for all t and z. We select a family {λ k } k∈{1,...,K} (which does not depend on t and z). In this paper, we choose λ k = 2 k-k 0 12n st : this exponential discretization is chosen to fit a transposition scale in subdivisions of the tone (nst corresponds to the number of discretized values of λ for each semitone). We assume that PI is given by:

∀λ ∈ [λ k 2 -1 24n st , λ k 2 1 24n st ], PI (λ, t|z) = PI (λ k , t|z).
Moreover, we assume that PI is zero outside these intervals. The values PI (λ k , t|z) (for all k, t and z) then completely describe PI .

Then:

f + 1 2 f -1 2 PI ( fc f ′ , t|z)dfc = f ′ k f,f ′ max k=k f,f ′ min PI (λ k , t|z)δλ f,f ′ k where k f,f ′ min is chosen so that λ k f,f ′ min 2 -1 24n st < f -1 2 f ′ ≤ λ k f,f ′ min 2 1 24n st and k f,f ′ max is chosen so that λ k f,f ′ max 2 -1 24n st ≤ f + 1 2 f ′ < λ k f,f ′ max 2 1
24n st (with the following constraints:

1 ≤ k f,f ′ min ≤ K and 1 ≤ k f,f ′ max ≤ K): kmin = k0 - 1 2 + 12nst log 2 ( f -1 2 f ′ ) kmax = k0 + 1 2 + 12nst log 2 ( f + 1 2 f ′ )
where ⌈.⌉ is the ceiling function and ⌊.⌋ the floor function.

Thus δλ f,f ′ k is given by δλ f,f ′ k = min (λ k 2 1 24n st , f + 1 2 f ′ ) - max (λ k 2 -1 24n st , f -1 2 f ′ ). When δλ f,f ′ k
is not limited by constraints on f and f ′ , we will denote

δλ k = λ k 2 1 24n st -λ k 2 -1 24n st .
The parameters to be estimated are then:

θ = {P (z), PK (f ′ |z), PI (λ k , t|z)|z ∈ [[1, Z]], f ′ ∈ [[1, F ′ ]], k ∈ [[1, K]], t ∈ [[1, T ]]}.

Expectation-Maximization algorithm

We intend to estimate the value of the parameter θ that maximizes the log-likelihood of observing V f t :

L(( f , t)|θ) = i∈I log P (fi, ti), (2)
where f and t correspond to the draws of f and t (draws are indexed by i ∈ I = {1, ..., N } where N is the total number of draws). As the number of draws that leads to the value (f, t) is V f t , the loglikelihood can be rewritten:

L(( f , t)|θ) = F f =1 T t=1 V f t log P (f, t)
The estimation will be performed with the Expectation-Maximization (EM) algorithm with latent variables z and f ′ (it would be equivalent to consider z and λ as latent variables since f = λf ′ ). The completed log-likelihood is:

L(( f , t, z, f ′ )|θ) = F f =1 T t=1 V f t log P (f, t, z, f ′ ) = f,t V f t log P (z) + log PK (f ′ |z) + log f + 1 2 f -1 2 PI ( fc f ′ , t|z)dfc + c
where c is a constant that does not depend on θ.

The completed log-likelihood expectation is then:

Q(θ|θ (c) ) = f ′ ,z,f,t V f t P (z, f ′ |f, t, θ (c) ) log P (z) + log PK (f ′ |z) + log f + 1 2 f -1 2 PI ( fc f ′ , t|z)dfc + c (3)
where θ (c) is the current value of the parameter. We can get an expression for P (z, f ′ |f, t, θ (c) ) with respect to θ (c) using the Bayes theorem (E step):

P (z, f ′ |f, t, θ (c) ) = P (c) K (f ′ |z)P (c) (z) f + 1 2 f -1 2 P (c) I ( fc f ′ , t|z)dfc f ′ P (c) (f, t)
.

(4) The notation (.) (c) refers to values computed from the current parameter:

θ (c) = {P (c) (z), P (c) K (f ′ |z), P (c) I (λ, t|z)}.
The completed expectation (3) will be maximized (M step) with respect to θ (θ (c) being fixed). As θ is made up of probabilities that must sum to 1, the maximization is constrained. Thus we use the Lagrangian1 :

H(θ|θ (c) ) = Q(θ|θ (c) ) + z ρz   1 - f ′ PK (f ′ |z)   + µ 1 - z P (z) + z τz 1 - t λmax λ min PI (λ, t|z) (5)
where ρz, µ, and τz are Lagrange multipliers.

Update of P (z) and PK

(f ′ |z) ∂H(θ|θ (c) ) ∂P (z)
= 0 leads to the update rule of P (z):

P (z) ← f,f ′ ,t V f t P (z, f ′ |f, t, θ (c) ) z,f,f ′ ,t V f t P (z, f ′ |f, t, θ (c) ) . ( 6
)
In a similar way, we obtain the update rule of PK (f ′ |z):

PK (f ′ |z) ← f,t V f t P (z, f ′ |f, t, θ (c) ) f,f ′ ,t V f t P (z, f ′ |f, t, θ (c) ) . (7) 3.2.2. Update of PI (λ k , t|z)
Because of the expression of the Lagrangian H(θ|θ (c) ) with respect to PI (λ k , t|z), the update rule of PI (λ k , t|z) is more complex to derive. We consider the following "fixed-point" update rule (iterated several times) which hopefully will converge to a zero of

∂H ∂P I (λ k ,t|z) : PI (λ k , t|z) ← f,f ′ V f t P (z,f ′ |f,t,θ (c) )P I (λ k ,t|z)1 [k min ,kmax ] (k) τz δλ k k f,f ′ max k ′ =k f,f ′ min P I (λ k ′ ,t|z)δλ f,f ′ k ′ δλ f,f ′ k . (8)
In equation ( 8), the division by τz is a normalization. We have not managed to prove convergence of PI under several iterations of the update rule (8) to a zero of ∂H ∂P I . However we observed it in practice.

As Q(θ|θ (c) ) (defined in equation ( 3)) and the constraints are

C 1 , Q(θ|θ (c)
) is strictly concave and the constraints are affine, and the regularity conditions are satisfied, a stationary point of H(θ|θ (c) ) is necessarily the global maximum of Q(θ|θ (c) ) under the normalization constraints. Thus, Q(θ|θ (c) ) is effectively maximized at each iteration using the update rules ( 6), ( 7) and ( 8) and the EM algorithm makes the likelihood of the observed data increase at each iteration. More details about the calculation of the update rules can be found in [START_REF] Hennequin | Scale-invariant probabilistic latent component analysis[END_REF].

EXAMPLES AND APPLICATIONS

Toy example

In this section, we use our algorithm to decompose a simple spectrogram. This spectrogram is obtained by STFT (using a 1024 sample-long Hann window with 75% overlap and no zero-padding) of a short excerpt of synthesizer which plays the notes of a A major scale on two octaves (from A4 to A6) sampled at 11025Hz. The original spectrogram is pictured in figure 1(a). The decomposition provides the reconstructed spectrogram pictured in figure 1(b): the reconstructed spectrogram is very similar to the original one. The difference of maximum amplitudes between original and reconstructed spectrograms come from the normalization of P (f, t) (V f t is not normalized), but the dynamic remains the same in both spectrograms. High frequency harmonics of the reconstructed spectrogram are slightly larger than the original one.

Obtained kernel distribution PK is represented in figure 1(d): we can see that the factorized template is actually harmonic. For high values of the frequency index, amplitudes of the templates tend to be very small. This comes from the fact that our model considers that values of the spectrogram outside the observed frequencies are zeros whereas the model spectrogram P (f, t) can take positive values outside this range. This can be solved using the approach described in [START_REF] Smaragdis | Missing data imputation for time-frequency representations of audio signals[END_REF]. In practice, for signals of actual acoustic instruments, this is not a real issue, since the harmonic content of such signals is mostly smothered by noise from 5000Hz: thus a sampling rate of 22050Hz or more permits to reduce this issue. Impulse distribution PI is represented in figure 1(c): high probabilities clearly appear at actual notes relative positions. There are also some replicas of the notes at positions with high harmonic similitudes (octave, twelfth, double octave...). At onset times, PI takes high values for many values of the homothety factor λ: this is caused by the flat shape of the spectrum of onsets which is matched here with several rescaled harmonic templates.

Real audio data

In this section, we present the SIPLCA decomposition of the spectrogram of the 10 first seconds of the song Because by the Beatles with a single template (Z = 1). The decomposed signal consists of a polyphonic harpsichord introduction recorded in real conditions. The original stereo signal was transformed into a mono signal (summing both channels) and downsampled to 22050Hz. The spectrogram was calculated with a STFT using a 2048 sample-long Hann window with 75% overlap and no zero-padding.

The obtained impulse distribution PI is represented in figure 2: actually played notes are materialized by a rectangle in the figure. We can see that in all rectangles, PI takes high values.

The impulse distribution PI is thus very similar to the impulse distribution that can be obtained with shift-invariant decompositions. However, as our decomposition is computed on linear frequency resolution spectrograms, it has an important advantage: it is possible to generate time-frequency masks that can be directly used to separate different components with Wiener filtering. Thus, it is possible to isolate single notes in a polyphonic signal and to repitch them individually. Examples of processed sounds are available on the webpage [13].

CONCLUSION

In this paper, we proposed a new way of decomposing nonnegative music spectrograms: the decomposition is based on a few frequency templates that can be rescaled at each frame (which corresponds to a transposition). We presented examples of this decom- position on music spectrograms and showed how it can be used to modify individual notes in a polyphonic signal.

Future works will include a better representation of nonharmonic components of musical sounds, such as transients.

Figure 1 :

 1 Figure 1: Scale-invariant PLCA: decomposition of a diatonic scale.

Figure 2 :

 2 Figure 2: Impulse distribution of the introduction of Because.

The research leading to this paper was supported by the CONTINT program of the French National Research Agency (ANR), as a part of the DReaM project (ANR-09-CORD-006-03), and by the Quaero Programme, funded by OSEO, French State agency for innovation.

October 16-19, 2011, New Paltz, NY

We did not write non-negative constraints in the problem statement because our final solution guarantees that parameters remain non-negative, then these constraints are inactive.


Romain Hennequin
email: romain.hennequin@telecom-paristech.fr

Roland Badeau

Bertrand David

TIME-DEPENDENT PARAMETRIC AND HARMONIC TEMPLATES IN NON-NEGATIVE MATRIX FACTORIZATION

This paper presents a new method to decompose musical spectrograms derived from Non-negative Matrix Factorization (NMF). This method uses time-varying harmonic templates (atoms) which are parametric: these atoms correspond to musical notes. Templates are synthesized from the values of the parameters which are learnt in an NMF framework. This parameterization permits to accurately model some musical effects (such as vibrato) which are inaccurately modeled by NMF.

INTRODUCTION

The decomposition of audio signals in terms of elementary atoms has been a large field of research for years. Sparse decomposition techniques [START_REF] Mallat | Matching pursuit with time-frequency dictionaries[END_REF] use a redundant dictionary of vectors (called atoms) and try to decompose signals on a few of them (much less than the dimension of the space). When atoms are designed to better encompass some signal properties (for instance harmonic atoms for musical signals [START_REF] Gribonval | Harmonic decomposition of audio signals with matching pursuit[END_REF]), the elements become less generic but more meaningful with regard to the context, and then a supervised classification can be performed to cluster atoms corresponding to a real event in the signal [START_REF] Leveau | Instrument-specific harmonic atoms for midlevel music representation[END_REF].

Recently methods of data factorization were proposed to simultaneously extract atoms from the signal and provide a decomposition on these atoms giving more robustness to the diversity of signals: NMF [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] has been introduced both to reduce the dimensionality and to explain the whole data by a limited number of elementary parts, possibly more significant regarding the considered objects. For instance, thanks to the non-negativity constraint, NMF applied to musical spectrograms will hopefully decompose them into notes, percussive sounds, or rapid transients. This interesting behavior of NMF leads to a wide dissemination of the tool in the audio signal processing community, with a number of applications such as automatic music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF][START_REF] Paulus | Drum transcription with non-negative spectrogram factorization[END_REF][START_REF] Bertin | Blind signal decompositions for automatic transcription of polyphonic music: NMF and K-SVD on the benchmark[END_REF] and sound source separation [START_REF] Cichocki | New algorithms for nonnegative matrix factorization in applications to blind source separation[END_REF][START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorizationwith temporal continuity[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. Unfortunately, NMF is not well adapted for time-varying phenomena such as vibrato: NMF is a rank reduction technique relying on the frame-to-frame redundancy and slight variations of the fundamental frequency drastically increase the rank by eliminating this redundancy (frequency of high order partials extensively varies and successive frames are then quite dissimilar).

To address this issue, Smaragdis [START_REF] Smaragdis | Sparse and shift-invariant feature extraction from nonnegative data[END_REF] proposes a shift-invariant extension of NMF (in the theoretical background of Probabilis-tic Latent Component Analysis) to decompose constant-Q spectrograms in which transposition can be seen as a shift of the templates. However, the constant-Q spectrogram imposes a quantification of the frequency bins and thus fundamental frequency in a vibrato cannot be estimated with precision. Moreover, the constant-Q transform is not invertible and time-domain signal reconstruction from the decomposition is consequently a difficult task [START_REF] Fitzgerald | Towards an inverse constant q transform[END_REF].

When used on musical spectrograms, for instance in automatic transcription, "interesting" atoms should hopefully present a specific structure, since musical note spectra possess a harmonic structure. This harmonicity is generally desirable and templates are usually constrained to be harmonic in automatic transcription application [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF].

The main idea proposed in this paper is to analytically synthesize such harmonic atoms in a parametric way: this method provides a parametric representation of the harmonic atoms, which can depend on a fundamental frequency parameter, a chirp parameter, a decrease/increase parameter and so on.

The harmonicity constraint is quite strong: the shape of the atoms is restrictive, since they can only be harmonic. However the parameterization should allow for a unique (parametric) atom for a note, the fundamental frequency of which slightly varies (for instance during a vibrato), which was a real issue with standard NMF.

This method provides a representation which can be suitable for automatic transcription, further providing an accurate description of the slight frequency variations.

In section 2, the principle of NMF is briefly reminded and our model is presented as an extension of NMF. In section 3 an algorithm which provides the decomposition is presented. In section 4, the algorithm is used to decompose a musical excerpt showing its ability to deal with slight variations of fundamental frequency (vibrato) and it is then compared to NMF. Finally, conclusions are drawn in section 5.

PARAMETRIC TEMPLATE

Non-negative Matrix Factorization

Given an F × T non-negative matrix V and an integer R such that F R + RT ≤ F T , NMF approximates V by the product of an F × R non-negative matrix W and an R × T non-negative matrix H:

V ≈ WH, V f t ≈ Vft = R r=1 w f r hrt. ( 1
)
When V is the magnitude or power spectrogram of a musical signal, the templates that are redundant in multiple frame hopefully are most of the time the harmonic templates corresponding to musical tones. Thus, each column of W should correspond to a note and each row of H is the time activation associated to each note. However, this property is not assured and generally, further constraints are added [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF][START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF].

Model

In our model, templates are parameterized following a timedependent value θrt. Equation ( 1) is thus replaced by:

V f t ≈ Vft = R r=1 w θ rt f r hrt ( 2
)
where θrt is the parameter associated to the template r at time t.

It can be considered as the "state" of this template: template r will then be synthesized from the value of this parameter for each t.

The time-dependence of the parameter now allows for modeling time-varying phenomena such as vibrato.

In this paper, the parameter chosen is the instantaneous fundamental frequency (noted θrt = f rt 0 ) of the template: each template is a harmonic comb parameterized by its fundamental frequency.

The parametric template writes:

w f rt 0 f r = n h (f rt 0 ) k=1 a k g(f -kf rt 0 ). (3)
This template corresponds in the time domain to a windowed stationary periodic sound (i.e. a windowed sum of sine functions). The Fourier transform of a periodic signal of fundamental frequency f0 is a sum of Dirac distributions centered in kf0 (with k ∈ Z ⋆ ). Thus, when such a signal is windowed, its Fourier transform is the convolution of the previous sum with the Fourier transform of the window. As templates should be non-negative, we take the squared modulus of this Fourier Transform. To make the calculation simpler the modulus of the sum of harmonics is replaced by the sum of the squared modulus of each harmonic: g is the squared modulus of the Fourier transform of the window used in the STFT. The interference between two successive partials is thus neglected; this approximation can be made for sufficiently high fundamental frequencies (or equivalently, sufficiently long analysis windows). The choice of the squared modulus is natural to make function g differentiable, in order to permit standard minimization algorithm.

The expression of function g and its derivative is given in appendix A.1 for a Gauss window and in appendix A.2 for Hann and Hamming windows. Then equation ( 3) is obtained where a k are the amplitudes of each harmonic, and n h (f rt 0 ) the number of harmonics.

Amplitudes a k of each harmonic are supposed to be the same for every atom and will be learnt in an unsupervised way. It would have been possible to have a different set of partial amplitudes for each template, however, this choice extensively enhances octave (and twelfth, double octave...) ambiguities.

As for a standard NMF, a cost function which is a distance (or a divergence) between the observed spectrogram and the reconstructed spectrogram will be minimized:

C(Θ, H, A) = D(V| V) = f t d(V f t | Vft ) (4)
where Θ = (θrt) r∈ 1,R ,t∈ 1,T , H = (hrt) r∈ 1,R ,t∈ 1,T , and

A = (a k ) k∈ 1,K .
In this paper, a β-divergence is chosen for d: this is a frequently used class of divergences in NMF frameworks which encompasses a number of usual divergences (Euclidean distance for β = 2, Kullback-Liebler divergence for β = 1 and Itakura-Saito divergence for β = 0). The β-divergence is defined for β ∈ R\{0, 1} by:

d β (x, y) = 1 β(β -1) (x β + (β -1)y β -βxy β-1 ).
For β ∈ {0, 1}, the value of the β-divergence is the limit of the previous expression.

Then, the cost function is minimized according to hrt, f rt 0 and a k for r ∈ 1, R , t ∈ 1, T and k ∈ 1, K .

However, the cost function is highly non-convex with respect to (wrt) f rt 0 : it can be seen in figure 1 for fixed r and t. We can observe several local minima: two prevailing minima correspond to actually played notes, and several lies near the octaves, the suboctaves, the twelfths... of each note (This figure can be seen as the opposite of a spectral product-like function). Consequently, a global optimization seems to be doomed to failure: a template with a 800Hz initial fundamental frequency will certainly converge to 880Hz. To avoid this issue, rather than a unique generic template for the whole scale, a template by chromatic degree is introduced. Thus, in the proposed decomposition, there is a single harmonic template associated to each chromatic degree. The fundamental frequency of the r th template can vary over time around

f rt 0 ≈ f ref 0 2 r-1 12
, where f ref 0 is the fundamental frequency of the lowest template. Templates with a low activation are discarded afterwards.

Figure 1: Cost function wrt the fundamental frequency f rt 0 of the templates r: analyzed spectrum is a mix of two harmonic spectra with fundamental frequencies 440Hz and 523Hz

ALGORITHM

Minimization of C(Θ, H, A) can be done with a multiplicative descent algorithm similar to those generally used for NMF: for the fundamental frequencies, the choice of multiplicative update rules is motivated by the positiveness of the frequency parameter and their natural logarithmic distribution.

In a multiplicative descent algorithm, the update rules associated to one of the parameters λ (here, λ = f (r, t), λ = h(r, t) or λ = a k for some value of r, t and k) are usually obtained by expressing the partial derivative of the cost function wrt this parameter as a difference of two positive terms:

∂C ∂λ = P λ -M λ . (5)
The update rule is then:

λ ← λ M λ P λ . ( 6
)
This rule particularly ensures that λ remains non-negative and becomes constant if the partial derivative is zero.

The partial derivative of the cost function ( 4) wrt one of the parameter λ is:

∂C ∂λ = f t ∂d(V f t , Vft ) ∂λ = f t ∂ Vft ∂λ ∂d ∂y (V f t , Vft )
where ∂d ∂y stands for the partial derivative of d wrt to its second argument. With a β-divergence cost function, this partial derivative is:

∂d ∂y (V f t , Vft ) = V β-2 f t ( Vft -V f t ).
Thus:

∂C ∂λ = f t ∂ Vft ∂λ V β-2 f t ( Vft -V f t ).

Update of f0

In order to obtain the update rules for the parameter Θ, one needs the partial derivative of the cost function wrt the parameter θr 0 t 0 :

∂C ∂θr 0 t 0 = f t ∂ Vft ∂θr 0 t 0 V β-2 f t ( Vft -V f t ).
The partial derivative of the parametric spectrogram Vft (see equation ( 2)) wrt θr 0 t 0 is given by:

∂ Vft ∂θr 0 t 0 = δtt 0 hr 0 t 0 ∂w θ r 0 t 0 f r 0 ∂θr 0 t 0
where δ is the Kronecker delta.

When the parameter Θ is the fundamental frequency of each template at each time, the partial derivative of the template wrt this fundamental frequency is obtained from (3):

∂w f r 0 ∂f0 r 0 t 0 = - n h k=1 a k kg ′ (f -kf0 r 0 t 0 ).
The partial derivative of the cost function wrt f0 r 0 t 0 is then:

∂C ∂f r 0 t 0 0 = - f n h k=1 hr 0 t 0 a k kg ′ (f -kf0 r 0 t 0 )( V β-1 f t 0 -V β-2 f t 0 V f t 0 ).
When g has a single lobe (when the window used in STFT is a Gauss window), the following remark is not an issue. However, generally, g has several lobes, and g ′ changes sign in numerous points. Then, in order to make the calculation easier, only the support of the main lobe of g (noted Λ) is kept in the expression of the derivative:

g ′ (f ) ≈ g ′ a (f ) = g ′ (f )1Λ(t). (7)
Assuming that the main lobe has a unique local maximum (which is actually true for Hamming and Hann windows, see figure 2), function f → -g ′ a (f )/f is then non-negative (see figure 2) and one can write g ′ a (f -kf0 r 0 t 0 ) = -(f -kf0 r 0 t 0 )P (f -kf0 r 0 t 0 ), where P is a positive function.

(f ) = -g ′ (f ) f on Λ = [-2 T , 2
T ] for a T -long Hamming window. 2Using approximation [START_REF] Bertin | Blind signal decompositions for automatic transcription of polyphonic music: NMF and K-SVD on the benchmark[END_REF], the partial derivative of the cost function wrt f r 0 t 0 0 is expressible as a difference of two non-negative terms: ∂C ∂f r 0 t 0 0 ≈ Gr 0 t 0 -Fr 0 t 0 [START_REF] Cichocki | New algorithms for nonnegative matrix factorization in applications to blind source separation[END_REF] with:

Gr 0 t 0 = f,k hr 0 t 0 a k kP (f -kf0 r 0 t 0 ) V β-2 f t 0 (f Vft 0 +kf0 r 0 t 0 V f t 0 ),
Fr 0 t 0 = f,k hr 0 t 0 a k kP (f -kf0 r 0 t 0 ) V β-2 f t 0 (kf0 r 0 t 0 Vft 0 +f V f t 0 ).
This results in the update rule of f0:

f r 0 t 0 0 ← f r 0 t 0 0 Fr 0 t 0 Gr 0 t 0 . (9)
The spectrogram will then be decomposed using a single template per semitone. However the update rule of f0 does not ensure that for each template r, f rt 0 remains close to its original fundamental frequency and does not slip to the fundamental frequency band of another semitone. Thus, evolution of f rt 0 should be restricted to a fundamental frequency band around the corresponding semitone. In our algorithm, when the fundamental frequency leaves its allocated frequency band, we consider that the corresponding template should not be active at this time and set the corresponding activation to 0.

Update of H

Update rules of H are obtained in a way similar to standard NMF by computing the derivative of the cost function wrt hr 0 t 0 :

∂C ∂hr 0 t 0 = f t ∂ Vft ∂hr 0 t 0 V β-2 f t ( Vft -V f t ) (10)
= f w θ r 0 t 0 f r 0 V β-2 f t 0 ( Vft 0 -V f t 0 ) (11)
= Pr 0 t 0 -Mr 0 t 0 [START_REF] Fitzgerald | Towards an inverse constant q transform[END_REF] where both

Pr 0 t 0 = f w θ r 0 t 0 f r 0 V β-1 f t 0 and Mr 0 t 0 = f w θ r 0 t 0 f r 0 V β-2 f t 0 V f t 0 are positive terms.
Then, the update rule of hr 0 t 0 is:

hr 0 t 0 ← hr 0 t 0 Mr 0 t 0 Pr 0 t 0 . ( 13
)

Update of A

Update rules of A are obtained in a way similar to the previous section by computing the partial derivative of the cost function wrt a k :

∂C ∂a k = f t ∂ Vft ∂a k V β-2 f t ( Vft -V f t ). (14)
The partial derivative of Vft wrt a k is:

∂ Vft ∂a k = R r=1 g(f -kf rt 0 )hrt1 [1,n h (f rt 0 )] (k).
Noting r k the maximum value of r for which k ∈ [1, n h (f rt 0 )], the previous equation becomes:

∂ Vft ∂a k = r k r=1 g(f -kf rt 0 )hrt.
Thus, the partial derivative of the cost function wrt a k can be naturally expressed as a difference of two positive terms:

∂C ∂a k = P k -M k (15)
with:

P k = f t r k r=1 g(f -kf rt 0 )hrt V β-1 f t , M k = f t r k r=1 g(f -kf rt 0 )hrt V β-2 f t V f t .
The update rule of a k is then:

a k ← a k M k P k . ( 16
)

Standard NMF templates to model non-harmonic events

In a musical spectrogram, percussive events (generated by percussive instruments or by the onset of harmonic instruments) do not correspond to harmonic templates and are thus inaccurately taken into account in our model. To encompass this kind of event a standard NMF decomposition term can be added to the parametric spectrogram proposed in (2):

V f t ≈ Vft = R r=1 w θ rt f r hrt + R ′ r=1 w ′ f r h ′ rt . ( 17
)
Thus w ′ f r is not time-varying and should model percussive events. R ′ should be kept very low (in the examples of section 4, R ′ = 1) in order to avoid non-parametric templates representing harmonic events. Update rules of

W ′ = (w ′ f r ) f ∈ 1F r∈ 1R ′ and H ′ = (h ′ rt ) r∈ 1R ′ t∈ 1T
are standard NMF multiplicative updates for β-divergence which can be found in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]:

W ′ ← W ′ . ((W ′ H ′ ) .β-2 .V)H ′T (W ′ H ′ ) .β-1 H ′T (18) H ′ ← H ′ . W ′T ((W ′ H ′ ) .β-2 .V) W ′T (W ′ H ′ ) .β-1 (19)
where the dot and the fraction bar stands for element-wise operation (element-wise multiplication, element-wise exponent and element-wise division).

Constraints

As for standard NMF, penalty terms can be added in the cost function to favor certain properties of the decomposition. The update rules are obtained in the same way as presented previously. The partial derivative of the constraint term Cp wrt the parameter λ to be updated is expressed as a difference of two positive terms:

∂Cp ∂λ = P p λ -M p λ .
The update rule (6) of the parameter λ thus becomes:

λ ← λ M λ + M p λ P λ + P p λ
where P λ and M λ were defined in equation ( 5). Several constraints have been considered for the proposed decomposition:

• Sparsity constraint on the columns of H as proposed in [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorizationwith temporal continuity[END_REF]: only a few templates should be simultaneously active. • Uncorrelation constraint between activations of a template and its octave (and eventually twelfth, double octave...) as proposed in [START_REF] Zhang | A NMF algorithm for blind separation of uncorrelated signals[END_REF]: to avoid octave ambiguities, a template should not be active. • Smoothness constraint on the coefficients of the amplitude as proposed in [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorizationwith temporal continuity[END_REF]: the spectral shape of templates should be smooth.

It has been observed that smoothness and uncorrelation constraints noticeably improve the results.

The whole algorithm without penalty term is detailed in Algorithm 1. The sampling rate of the excerpt is fs = 11025Hz. We chose a 1024-sample long Hamming window with 75% overlap for the STFT. The decomposition was made with 72 templates (6 octaves) with fundamental frequencies every semitone from 55Hz (A0) to 3322Hz (G#6). Kullback-Liebler divergence (β-divergence with β = 1) was used. Two constraint terms were added: a uncorrelation constraint on activation of octave (and twelfths) and a smoothness constraint on A. All notes were played by the synthesizer with a slight vibrato.

Notes of the piece clearly appear with strong activation values. The maximum simultaneous polyphony is of 3 notes. Even when one note is played simultaneously with its octave, the activation values correspond to effectively played notes. At onset instant, numerous templates are active because no template is able to accurately represent an onset. To reduce this problem, one standard NMF template (non-parametric template) was added to the model in order to better represent the onsets (see section 3.4).

The reconstructed spectrogram is represented in figure 4. One can easily build a synthetic time/frequency representation of the activations that includes the varying fundamental frequency: for each template r, at each time t, a thin peak is generated in the time/frequency plane at (t, f rt 0 ) with amplitude equal to the activation hrt. Such a representation is given in figure 6 for the same excerpt in figure 5. This representation reveals the vibrato generated by the synthesizer.

Comparison with NMF

In this section we will compare our method with the method of decomposition proposed in [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF] which is based on NMF with templates that are imposed to be harmonic: in this method, templates are linear combinations of narrow band harmonic patterns. This method is used for transcription and decomposes Equivalent Rectangular Bandwidth (ERB) power spectrograms on 88 templates corresponding each one to a note in the chromatic scale. Smoothness constraints on the activations are added. The input data is not exactly the same (ERB spectrogram in [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF] and Standard STFT spectrogram in our method) but the original time-domain signal used is the same and the kind of decomposition provided is very similar in both cases (harmonic templates corresponding to notes The signal analyzed is a C diatonic scale played from C1 to B4 harmonized with the corresponding third in the scale (each note of the scale is simultaneously played with its third). The scale was played by a synthesizer.

The activation provided by our algorithm is given in figure 7: all notes played clearly appears. The amplitude shape of each note is clear, with a strong onset and a slight decay. One could notice the simultaneous activation of numerous templates at onset times, which is due to the percussive spectral shape of the spectrogram at these instants.

The activation provided by the algorithm of [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF] is given in figure 8 with the same dynamic as the activation of our algorithm: most notes are correctly located. However, there are several octave (or sub-octave, twelfth, sub-twelfth...) errors. Moreover the smoothness constraint does not make possible to see the amplitude shape of each note, the onset being spread and the decay being extended. The odd activation of the second template is probably due to the presence of non-harmonic events which are taken into ac-count by this template. Thus, our decomposition seems to enhance the representation of this signal. [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF] (amplitude scale is in dB).

CONCLUSION

In this paper, we presented a new way of decomposing musical spectrograms on a basis of parametric templates which correspond to musical notes. This decomposition provides a good representation of the different notes which are played in each column of the spectrogram. The decomposition being parametric is thus very flexible.

In future work, a supervised learning of atoms could be included in order to improve results and to decompose spectrograms containing several instruments. Moreover, the parameterization of templates could be extended in order to model inharmonicity (that occurs in piano notes), fast fundamental frequency variations (chirp) or fast decrease/increase of the amplitude of the harmonics, and thus to better fit to real audio spectrograms.

Figure 2 :

 2 Figure 2: Main lobe of g, derivative of g and positivity of P (f ) = -g ′ (f )

Algorithm 1 2 r- 1 12 f ref 0 for j = 1 rt 0 f ref 0 -

 121010 Figure5represents the activations hrt of the decomposition of the power spectrogram (represented in figure3) of an excerpt (four first bars) of the J.S. Bach's first prelude played by a synthesizer. The sampling rate of the excerpt is fs = 11025Hz. We chose a 1024-sample long Hamming window with 75% overlap for the STFT. The decomposition was made with 72 templates (6 octaves) with fundamental frequencies every semitone from 55Hz (A0) to 3322Hz (G#6). Kullback-Liebler divergence (β-divergence with β = 1) was used. Two constraint terms were added: a uncorrelation constraint on activation of octave (and twelfths) and a smoothness constraint on A. All notes were played by the synthesizer with a slight vibrato.Notes of the piece clearly appear with strong activation values. The maximum simultaneous polyphony is of 3 notes. Even when one note is played simultaneously with its octave, the activation values correspond to effectively played notes. At onset instant, numerous templates are active because no template is able to accurately represent an onset. To reduce this problem, one standard NMF template (non-parametric template) was added to the model in order to better represent the onsets (see section 3.4).The reconstructed spectrogram is represented in figure4.

Figure 3 :Figure 4 :

 34 Figure 3: Original spectrogram of the excerpt of J.S. Bach's first prelude

Figure 5 :

 5 Figure 5: Activations in the decomposition of the spectrogram of the excerpt of J.S Bach's first prelude. Color scale is in dB. Semitones correspond to relative MIDI note (relative to the lower pitch).

Figure 6 :

 6 Figure 6: Representation of the activations including fundamental frequencies (Two first bars of the excerpt of J.S. Bach's first prelude). Color scale is in dB.

Figure 7 :

 7 Figure 7: Representation of the activations provided by our algorithm (amplitude scale is in dB).

Figure 8 :

 8 Figure8: Representation of the activations provided by the algorithm of[START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian nonnegative matrix factorization applied to polyphonic music transcription[END_REF] (amplitude scale is in dB).

Proc. of the 13 th Int. Conference on Digital Audio Effects (DAFx-10), Graz, Austria ,September 6-10, 2010

In order to represent all curves on the same plot, the y-axis scale has been modified and is then arbitrary.

The research leading to this paper was supported by the Quaero Programme, funded by OSEO, French State agency for innovation and by the CONTINT program of the French National Research Agency (ANR), as a part of the DReaM project (ANR-09-CORD-006-03).

A. EXPRESSION OF g

In this section, we use the following definition of the Fourier transform of a continuous-time signal x:

x(t)e -i2πf t dt.

A.1. Gauss window

In this section, we will give the properties of the Gauss window, i.e. the window defined in the time-domain by:

where σ characterizes the width of the peak.

The Fourier transform of this window is:

Consequently, the expression of g for this type of window is:

The derivative of g is then:

For all frequencies g ′ (f ) f ≤ 0, which permits to easily write the partial derivative of the cost function wrt the fundamental frequency of a template at a given time as the difference of two positive terms (cf. equation ( 8)).

The Gauss window has good frequency properties (it has a single main lobe and it saturates the Heisenberg inequality) but are rarely used because of its infinite support and because it does not permit the perfect reconstruction of a signal from its spectrogram.

A.2. "Cosine" window

In this section, we will study the properties of windows defined in the time-domain by:

where T is the length of the window, 1 [0,T ] (t) is the indicator function of the interval [0, T ], and α + β = 1 (the maximum of the window is equal to 1). This class of windows encompasses:

• Hann window (sometimes referred as Hanning window), for α = β = 0.5.

• Hamming window, for α = 0.54 and β = 0.46.

The Fourier transform of this window is:

Thus, the expression of g for this type of window:

Remark: g can be C 1 -prolonged in 0 and ±T with g(0) = α 2 T 2 and g(±T ) = β 2 T 2 4 . The derivative of g is then:

For Hann and Hamming windows, the main lobe corresponds to frequencies f ∈ [-2 T , 2 T ]. For these frequencies, g ′ (f ) f ≤ 0, which permits to easily write the partial derivative of the cost function wrt the fundamental frequency of a template at a given time as a difference of two positive terms (cf. equation ( 8)), assuming that g is zero outside the main lobe.


Romain Hennequin

Roland Badeau

Bertrand David

NMF WITH TIME-FREQUENCY ACTIVATIONS TO MODEL NON STATIONARY AUDIO EVENTS

Keywords: music information retrieval, non-negative matrix factorization, unsupervised machine learning

Real world sounds often exhibit non-stationary spectral characteristics such as those produced by a harpsichord or a guitar. The classical Non-negative Matrix Factorization (NMF) needs a number of atoms to accurately decompose the spectrogram of such sounds. An extension of NMF is proposed hereafter which includes timefrequency activations based on ARMA modeling. This leads to an efficient single-atom decomposition for a single audio event. The new algorithm is tested on real audio data and shows promising results.

INTRODUCTION

Human cognition utilizes redundancies to understand visual and audio signals, and several techniques tend to mimic this behavior when decomposing and approximating signals (sounds or images for instance): Principal Component Analysis, Independent Component Analysis or NMF (Lee and Seung [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]) have been introduced both to reduce the dimensionality and to explain the whole data by a few meaningful elementary objects. Thanks to the non-negativity constraint, NMF is able to provide a significant picturing of the data: applied to musical spectrograms it will hopefully decompose them in elementary notes or impulses. The technique is widely used in audio signal processing, with a number of applications such as automatic music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF] and sound source separation [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity[END_REF].

However, the standard NMF is shown to be efficient when the elementary components (notes) of the analyzed sound are nearly stationary. In case of a noticeable spectral variability, the standard NMF will likely need several non-meaningful atoms to decompose a single event, which often leads to a necessary post-processing. To override this drawback, Smaragdis [START_REF] Smaragdis | Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs[END_REF] proposes a shift-invariant extension of NMF in which time/frequency templates are factorized from the original data: each atom then corresponds to a time-frequency musical event able to include spectral variations over time. This method gives good results, but does not permit any variation between different occurrences of the same event (atom), its duration and spectral content evolution being fixed.

Durrieu introduces in [START_REF] Durrieu | An iterative approach to monaural musical mixture de-soloing[END_REF] a source/filter model in a NMF framework in order to model the main melody of musical pieces. This model permits to efficiently take the strong spectral variations of the human voice into account. The filter model is constrained to be a linear combination of pass-band templates while the source is a harmonic template derived from physical modeling. In this paper an

The research leading to this paper was supported by the French GIP ANR under contract ANR-06-JCJC-0027-01, DESAM, and by the Quaero Program, funded by OSEO, French State agency for innovation. extension of NMF is proposed which presents some similarity with Durrieu's model but which includes AutoRegressive Moving Average (ARMA) models estimated from the data, and learns the sources (atoms) in a totally unsupervised way.

In section 2, we introduce the source/filter decomposition as an extension of NMF. In section 3, we derive an iterative algorithm similar to those used for NMF to compute this decomposition. In section 4, we present experiments of source/filter decomposition of the spectrogram of two sounds, and compare this decomposition to the standard NMF. Conclusions are drawn in section 5.

MODEL

NMF and extension

Given an F × T non-negative matrix V and an integer R such that F R+RT ≪ F T , NMF approximates V by the product of an F ×R non-negative matrix W and an R × T non-negative matrix H:

V ≈ WH i.e. V f t ≈ R r=1 w f r hrt (1)
This approximation is generally quantified by a cost function C(V, W, H) to be minimized with respect to (wrt) W and H, which is generally designed element-wise. In this article we will focus on a cost function built with the β-divergence d β (see [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] for its expression) which includes usual measures: Euclidian distance (β = 2), Kullback-Leibler divergence (β = 1) and Itakura-Saito divergence (β = 0).

When applied to power (squared magnitude) spectrograms, NMF factorizes data into a matrix (or basis) of frequency templates which are the R columns of W and a matrix H whose R rows are the temporal vectors of activations corresponding to each template. For a musical signal made of several notes played by the same instrument, it is hoped that the decomposition leads to spectral templates corresponding to single notes or percussive sounds. H will then display a representation similar to a "piano-roll" (cf. [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF]).

This factorization however does not permit to well represent a sound presenting a noticeable spectral evolution. For instance a single note of a plucked string instrument most of the time shows high frequency components which decrease faster than low frequency components. This characteristic is not well modeled with a single frequency template. Several templates are needed which results in a less meaningful decomposition.

To address this issue, we propose an extension of NMF where temporal activations become time/frequency activations. The factorization (1) becomes:

V f t ≈ Vft = R r=1 w f r hrt(f ) (2)
where the activation coefficients are now frequency dependent. To avoid an increase of the problem dimensionality the hrt(f ) coefficients are further parameterized by means of ARMA models (paragraph 2.2). Equation ( 2) can be interpreted with the help of the source/filter paradigm: the spectrum of each frame of the signal results from the combination of filtered templates (sources). hrt(f ) corresponds to the time varying filter associated to the source r. The decomposition thus benefits from the versatility of the source/filter model which is well suited for numerous sound objects.

AutoRegressive Moving Average (ARMA) Modeling

hrt(f ) is chosen following the general ARMA model:

h ARM A rt (f ) = σ 2 rt Q q=0 b q rt e -i2πν f q 2 P p=0 a p rt e -i2πν f p 2 where ν f = f -1 2F
is the normalized frequency associated to frequency index f ∈ {1, ..., F } (as audio signal are real valued, we only consider frequencies between 0 and the Nyquist frequency). b q rt are the coefficients of the MA part of the filter and a p rt those of the AR part. σ 2 rt is the global gain of the filter. For P = Q = 0, h ARM A rt (f ) no longer depends on f and the decomposition corresponds to a standard NMF with temporal activations σ 2 rt . Defining art = (a 0 rt , . . . , a P rt ) T and brt = (b 0 rt , . . . , b Q rt ) T , time/frequency activations can be rewritten as:

h ARM A rt (f ) = σ 2 rt b T rt T(ν f )brt a T rt U(ν f )art
where T(ν) is the (Q+1)×(Q+1) Toeplitz matrix with [T(ν)]pq = cos(2πν(p -q)) and U(ν) is similar to T(ν) but of dimension (P + 1) × (P + 1). MA only and AR only models are included by respectively taking P = 0 and Q = 0. It is worth noting that h ARM A rt (f ) is always non-negative while there exists neither nonnegativity constraint on b p rt nor on a q rt . The parameterized power spectrogram then becomes:

Vft = R r=1 w f r σ 2 rt b T rt T(ν f )brt a T rt U(ν f )art (3)

ALGORITHM

We choose a general β-divergence cost function:

C(W, A, B, Σ) = F f =1 T t=1 d β (V f t , Vft ) with [W] f r = w f r , [Σ]rt = σ 2 rt , [A]rtp = a p rt and [B]rtq = b q rt
and the expression of d β is given in [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF].

The partial derivative of the cost function wrt any variable y (y being any coefficient of W, Σ, A or B) is:

∂C(W, A, B, Σ) ∂y = F f =1 T t=1 V β-2 f t ( Vft -V f t ) ∂ Vft ∂y (4)
The expression of the gradient of C wrt a vector y of several coefficients of A or B is the same, replacing the partial derivative by a gradient ∇y in (4). This leads to update rules for a multiplicative gradient descent algorithm similar to those used in [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF][START_REF] Smaragdis | Non-negative matrix factor deconvolution; extraction of multiple sound sources from monophonic inputs[END_REF]. In such an iterative algorithm, the update rule associated to one of the parameters is obtained from the partial derivative of the cost function wrt this parameter, written as a difference of two positive terms. In the particular case of a scalar parameter, ∂C ∂y = Gy -Fy with Gy =

V β-1 f t ∂ Vft ∂y and Fy = V β-2 f t V f t ∂ Vft
∂y and the update rule for y is:

y ← y × Fy Gy (5)
This rule particularly ensures that y remains non-negative and becomes constant if the partial derivative is zero.

Update of frequency templates and filter gains

Following equation ( 5), we obtain the update rules of w f 0 r 0 and σ 2 r 0 t 0 :

w f 0 r 0 ← w f 0 r 0 T t=1 h ARM A r 0 t (f0) V β-2 f 0 t V f 0 t T t=1 h ARM A r 0 t (f0) V β-1 f 0 t (6) σ 2 r 0 t 0 ← σ 2 r 0 t 0 F f =1 w f r 0 h ARM A r 0 t 0 (f ) V β-2 f t 0 V f t 0 F f =1 w f r 0 h ARM A r 0 t 0 (f ) V β-1 f t 0 (7)

Update of filters

The update rules of the coefficients of the filters are derived in a similar way, but the updates are not element-wise, but rather "vectorwise": we derive an update rule for each brt and for each art.

Update of brt:

The gradient of the parameterized spectrogram Vft wrt br 0 t 0 is:

∇ b r 0 t 0 Vft = δt 0 t 2w f r 0 σ 2 r 0 t 0 a T r 0 t 0 U(ν f )ar 0 t 0 T(ν f )br 0 t 0
Then, by substituting this expression into equation ( 4) with y = br 0 t 0 , we obtain the gradient of the cost function wrt br 0 t 0 :

∇ b r 0 t 0 C = 2 F f =1 w f r 0 σ 2 r 0 t 0 V β-2 f t 0 ( Vft 0 -V f t 0 ) a T r 0 t 0 U(ν f )ar 0 t 0 T(ν f )br 0 t 0 = 2σ 2 r 0 t 0 Rr 0 t 0 -R ′ r 0 t 0 br 0 t 0
where:

Rr 0 t 0 = F f =1 w f r 0 V β-1 f t 0 a T r 0 t 0 U(ν f )ar 0 t 0 T(ν f ) R ′ r 0 t 0 = F f =1 w f r 0 V β-2 f t 0 V f t 0 a T r 0 t 0 U(ν f )ar 0 t 0 T(ν f )
Both matrices Rr 0 t 0 and R ′ r 0 t 0 are positive definite under mild assumptions. Then, we follow the approach given in [START_REF] Badeau | Weighted maximum likelihood autoregressive and moving average spectrum modeling[END_REF] and derive the following update rule for the MA part of the filter:

br 0 t 0 ← R -1 r 0 t 0 R ′ r 0 t 0 br 0 t 0 (8)
As Rr 0 t 0 and R ′ r 0 t 0 are both non singular, R -1 r 0 t 0 is well defined and br 0 t 0 is ensured to never be zero.

Update of art:

The update rules of art are derived in the same way as for brt. Thus, defining:

Sr 0 t 0 = F f =1 w f r 0 V β-1 f t 0 b T r 0 t 0 T(ν f )br 0 t 0 a T r 0 t 0 U(ν f )ar 0 t 0 2 U(ν f ) and S ′ r 0 t 0 = F f =1 w f r 0 V β-2 f t 0 V f t 0 b T r 0 t 0 T(ν f )br 0 t 0 a T r 0 t 0 U(ν f )ar 0 t 0 2 U(ν f )
we derive the following update rule for the AR part of the filter:

ar 0 t 0 ← S ′-1 r 0 t 0 Sr 0 t 0 ar 0 t 0 (9)

Description of the algorithm

The update rules ( 6), ( 7), ( 8) and ( 9) are applied successively to all the coefficients of W, all the coefficients of Σ, all the coefficients of B and all the coefficients of A. Between the updates of each of these matrices (and tensors), the parameterized spectrogram V is recomputed.

Identification: As for the standard NMF, the decomposition (3) which minimizes the cost function is not unique. To reduce identification problems, we impose constraints on W, Σ, B and A:

• for all r and t, we impose that brt and art (considered as polynomials) have all their roots inside the unit circle. • for all r, we impose wr = 1 for some norm . .

• for all r and t, we impose b 0 rt = 1 and a 0 rt = 1. Thus, at the end of each iteration of our algorithm, we transform br,t and ar,t by replacing roots outside the unit circle by the conjugate of their inverse and accordingly adapting the gain, normalize each column of W, divide br,t and ar,t by their first coefficient and update Σ in order not to change Vft by these modifications. All these transformations have no influence on the values of the parameterized spectrogram.

In the remainder of the article, we will refer to this algorithm by the expression "source/filter factorization".

Dimensionality

In the standard NMF with R atoms, the dimension of the parameterized spectrogram is dim W + dim H = R(F + T ). With our algorithm, the dimension of the parameters is:

dim W + dim Σ + dim A + dim B = RF + RT (P + Q + 1
). Thus, one should have RF + RT (P + Q + 1) ≪ F T , so P and Q must remain small.

One should notice that our decomposition permits to considerably reduce the number of atoms R needed to accurately fit the data when the parts of the sounds present strong spectral variations, which keeps the dimension of the parameters small.

EXAMPLES

In this section several experiments are presented to show that our algorithm is well adapted to decompose sounds having strong spectral variations. It is quite difficult to objectively compare the algorithm with other decomposition algorithms like NMF, as signal models are different. Thus, in this paper, we just present illustrations of the proposed decomposition on real audio data.

All the spectrograms used in these experiments are power spectrograms obtained from recorded signals by means of a short time Fourier transform (STFT).

In lack of any theoretical proof, we chose β = 0.5 since we empirically observed the monotonic decrease of the cost function and the convergence of the algorithm for 0 ≤ β ≤ 0.5 over a large set of tests. Moreover the results were more accurate with β = 0.5 than with β = 0 (Itakura-Saito divergence).

Algorithms (standard NMF and source/filter factorization) were initialized with random values (except for the filters which were initially flat) and were run until apparent convergence. The algorithms were reinitialized a hundred times (in order to come close to the "best" minimum possible).

Didgeridoo

Description of the excerpt

In this section our algorithm is applied to a short didgeridoo excerpt. The didgeridoo is an ethnic wind instrument of northern Australia. It makes a continuous modulated sound produced by the vibrations of the lips. The modulations result from the mouth and throat configuration with the help of which the player is able to control several resonances. Figure 1(a) represents the spectrogram of the excerpt: the sound produced is almost harmonic (with some noise) and a strong moving resonance appears in the spectrogram. We can thus consider that this signal is composed of a single event encompassing spectral variations, and try to decompose it with a single atom (R = 1). The sampling rate of the excerpt is fs = 11025Hz. We chose a 1024 samples long Hann window with 75% overlap for the STFT.

Experiment and results

The spectrogram of the excerpt is decomposed with a standard NMF algorithm for R = 1 atom and R = 5 atoms, and with source/filter factorization for R = 1 atom, with an order 3 AR modeling (Q = 0 and P = 3). Reconstructed spectrograms are respectively represented in figures 1(b), 1(c) and 1(d).

Albeit the didgeridoo is played alone in the analyzed spectrogram, the standard NMF needs many atoms to accurately decompose the power spectrogram. With 1 atom, NMF does not accurately represent the moving resonance (figure 1(b)). With 5 atoms, some spectral variations appears (figure 1(c)), but the resonance trajectory remains a bit unclear. Besides, the signal is not decomposed in a meaningful way (each atom is a part of the sound which has no perceptual meaning) and the dimensionality of the parameters is large (F R + RT = 3290).

In opposition to the standard NMF, source/filter factorization permits to accurately represent the spectral variability of the sound (figure 1(d)) with a single atom, keeping the dimensionality low (F R + T R(Q + 1) = 1093): the moving resonance of the original sound is well tracked, and the total error C is smaller than that of the standard NMF with R = 5. In this case, the decomposition is more efficient and relevant than the standard NMF.

Harpsichord

Description of the excerpt

In this section our algorithm is applied to a short harpsichord excerpt, composed of two different notes (C2 and E♭2): first, the C2 is played alone, then the E♭2, and at last, both notes are played simultaneously. The spectrogram of the extract is represented in figure 2(a). As for most of plucked string instruments, high frequency partials of a harpsichord tone decay faster than low frequency partials.

Experiment and results

The spectrogram of the excerpt was decomposed with a standard NMF algorithm for R = 2 atoms (1 atom per note) and R = 6 atoms, and with source/filter factorization for R = 2 atoms, with an ARMA modeling (Q = 1 and P = 1). Reconstructed spectrograms are respectively represented in figure 2(b), 2(c) and 2(d).

The standard NMF needs several atoms per note to accurately decompose the L-shaped power spectrograms: with only 2 atoms (1 per note played), the faster decay of high frequency content does not appear at all (figure 2(b)). With 6 atoms, the attenuation of high frequency partials appears (figure 2(c)), but each atom is a part of a note spectrum and has no real perceptual meaning.

The ARMA modeling included in our algorithm leads to a good description of the overall spectrogram shape. 2 atoms (1 per note) are enough to accurately fit the original short time spectrum: each atom is harmonic (figure 3(a)) and corresponds to one note while the decay of high frequency partials is clearly well described by the ARMA modeling (see time/frequency activations h ARM A rt (f ) in figure 3(b)). The dimensionality of the data provided by our algorithm (F R + T R(Q + P + 1) = 5704) remains lower than with a standard NMF with 6 atoms (F R + RT = 9804) and the global error C between the original and the reconstructed spectrogram is lower than that obtained with the standard NMF with R = 6.

Thus the decomposition provided by source/filter factorization seems to give a more meaningful representation of the given spectrogram than the one obtained with the standard NMF.

CONCLUSION AND FUTURE WORK

In this paper, we proposed a new iterative algorithm which is an extended version of the Non-negative Matrix Factorization based on a source/filter representation. We showed that this representation is particularly suitable to efficiently and meaningfully decompose non stationary sound objects including noticeable spectral variations.

In the future, this extended decomposition could be further developed to deal with small variations of the pitch (like vibrato), for instance by using a constant-Q spectrogram like in [START_REF] Schmidt | Nonnegative matrix factor 2-D deconvolution for blind single channel source separation[END_REF]. Moreover, quantitative evaluation of the method is under active consideration.

1 Fig. 1 . 2 Fig. 2 .

 1122 Fig. 1. Original power spectrogram of the extract of didgeridoo 1(a) and reconstructed spectrograms 1(b), 1(c) and 1(d)

Fig. 3 .

 3 Fig. 3. Source/filter decomposition (R = 2, Q = 1 and P = 1) of the power spectrogram of the harpsichord excerpt