Roland Badeau
email: roland.badeau@enst.fr

Rémy Boyer
email: remy.boyer@enst.fr

Bertrand David
email: bertrand.david@enst.fr

EDS PARAMETRIC MODELING AND TRACKING OF AUDIO SIGNALS

Despite the success of parametric modeling in various fields of digital signal processing, the Fourier analysis remains a prominent tool for many audio applications. This paper aims at demonstrating the usefulness of the Exponentially Damped Sinusoidal (EDS) model both for analysis/synthesis and tracking purposes.

INTRODUCTION

The advantages of this model are on the one hand to overcome the Fourier resolution limit related to windowing and on the other hand to enhance the classical sinusoidal model used in speech and audio coding. The main drawbacks of the EDS-based methods consist in the complexity of the algorithms and in the assumption of non time-varying parameters. In this paper, applications of some recent enhancements of the algorithms to audio signal processing are presented, in order to both reduce the complexity and track the parameter variations with time.

Modeling context. Recently, many efforts have been made to achieve a powerful representation of an audio signal as speech or music, for a compression purpose [START_REF] Iso-Mpeg | Call for proposals for new tools for audio coding[END_REF][START_REF] Purnhagen | HILN-the MPEG-4 parametric audio coding tools[END_REF]. More specifically, in parametric audio coding, it is worthwhile to have compact (sparse) representations of the signal: the model order M (i.e. the number of elementary components) must be far less than N , the length of the analysis window in samples. One way to obtain a more compact representation is to increase the parameter N and keep unchanged the model order. Unfortunately, for large N , the audio signal can no longer be considered as a quasi-stationary signal. In this case, the basic sinusoidal model [START_REF] Mcaulay | Speech analysis and synthesis based on a sinusoidal representation[END_REF], which tends to represent the audio signal as a sum of constant-amplitude components, becomes ineffective. Consequently, the EDS model was introduced in the audio modeling context [START_REF] Nieuwenhuijse | Robust exponential modeling of audio signal[END_REF][START_REF] Boyer | Non-stationary signal parametric modeling techniques with an application to low bitrate audio coding[END_REF]. In this work, we compare the sinusoidal and EDS models with the same total number of model parameters. Note that keeping a satisfactory algorithmic complexity implies setting a maximal bound to the parameter N .

Tracking context. The EDS model relies on the assumption of non varying signal parameters within the observation window. A more realistic modeling of musical signals should include slow variations of the parameters. Tracking these time variations would have interesting applications, such as:

• evaluating the degree of stationarity of the audio signal,

• detecting model breaks, which characterize transient sounds,

• developing more realistic synthesis techniques.

A reference method in frequency tracking is the Sintrack algorithm introduced by P. Duvaut [6]. This method relies on a fast linear prediction technique, which makes it useful for real-time estimation and tracking of damped sinusoids in noise. However, its lack of robustness results in repeated re-initializations which increase the computational cost. Concurrently, subspace-based high resolution methods, despite their higher computational complexity, prove to be much more reliable than linear prediction. Therefore, adaptive subspace estimation may offer interesting outlooks for frequency tracking.

Contents. This paper is organized as follows.

Section 2 introduces the EDS model and presents subspacebased high resolution methods for the estimation and tracking of the model parameters. Some synthesis techniques are proposed both in a static and an adaptive context, with an application to pitch modification.

Section 3 shows the application of these methods to coding, tracking and re-synthesis of audio signals.

Finally, section 4 summarizes the main conclusions of this paper.

THEORETICAL BACKGROUND

The EDS model defines the discrete signal as

x(t) = M m=1
am exp(dmt) cos (2πfmt + φm), t ∈ {0, . . . N -1}

(1) where x(t) is the discrete signal observed in the window t ∈ {0, . . . N -1}, M is the order of the model, am ∈ R * + are the amplitudes, dm are the real valued damping factors, fm ∈ [- 1 2 , 1 2 [ are the frequencies and φm ∈ [-π, π[ denote initial phases. Equation (1) can equivalently be rewritten with the complex amplitudes αm = 1 2 am exp(iφm) and the complex poles zm = exp(dm + i2πfm) as in equation (2):

x(t) = M m=1 αmz t m + α * m z * m t . (2)
In section 2.1, EDS-based analysis/synthesis methods are presented in a block processing context (with constant model parameters). In section 2.2, it will be shown how these methods can be adapted to track slow variations of these parameters.

Block signal processing

The estimation of the model parameters is achieved in two steps: first the frequencies and damping factors are computed using a high resolution (HR) method, from which the amplitudes and initial phases are deduced by minimizing a least squares (LS) criterium. The estimated parameters are then used to re-synthesize the signal.

Initialization: U S = I2M 0 (L-2M )×2M
For n = 1, 2, . . . until convergence iterate:

A(n) = H U S (n -1) fast matrix product A(n) = U S (n)R(n)
skinny QR factorization

Subspace-based signal analysis

Define the L × L real Hankel data matrix H (with N = 2L -1) as

H =      x(0) x(1) . . . x(L -1) x(1)
x(2) . . . x(L) . . . . . . . . . . . .

x(L -1) x(L) . . . x(N -1)      . (3)
Suppose that 2M ≤ L. Then this matrix can be decomposed as H = E A E T , where

A = Diag(α1, . . . , αM , α * 1 , . . . , α * M )
and E is the L × 2M Vandermonde matrix

E =      1 . . . 1 1 . . . 1 z1 . . . zM z * 1 . . . z * M . . . . . . . . . . . . . . . . . . z L-1 1 . . . z L-1 M z * 1 L-1 . . . z * M L-1      . (4)
H has a 2M -dimensional range space, spanned by the fullrank matrix E. This range space fully characterizes the signal poles, even in presence of an additive white noise [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], and thus is referred to as the signal subspace. An orthonormal basis U S of this space can be obtained from the eigenvalue decomposition (EVD) of H . Indeed, since H is a rank-deficient symmetric real matrix, there exist a L × 2M orthonormal real matrix U S and a 2M × 2M diagonal real matrix Λ such that H = U S Λ U S T . The columns of U S thus span the signal subspace. In the presence of an additive white noise, the columns of U S are defined as the 2M -dominant eigenvectors of H (i.e. the eigenvectors associated to the 2M eigenvalues which have the highest magnitudes).

These dominant eigenvectors can be computed using the classical EVD algorithm called orthogonal iteration [START_REF] Golub | Matrix computations[END_REF] 1 (cf. table 1), which involves an auxiliary matrix A. The Hankel structure of the matrix H can be taken into account to make the algorithm faster by computing the first-step matrix product using Fast Fourier Transforms, which requires only O(LM log(L)) operations [START_REF] Golub | Matrix computations[END_REF] 2 . Then the second step can be achieved in O(LM 2 ) operations [START_REF] Golub | Matrix computations[END_REF] 3 . Since in practice this algorithm converges in a few iterations, the overall process requires O(LM (M + log(L))) operations.

1 Chapter 8, section 2.4. 2 Chapter 4, section 7.7. 3 Chapter 5, section 2.

Estimation of the frequencies and damping factors

The poles {zm, z * m } 1≤m≤M can be calculated by exploiting the rotational invariance property of the signal subspace. More precisely, define E ↓ (respectively E ↑ ) the matrix extracted from E by deleting the last (respectively the first) row. These matrices satisfy the equation

E ↑ = E ↓ D (5)
where

D = diag(z1, z2, . . . , zM , z * 1 , z * 2 , . . . , z * M ).
Since the matrices U S and E span the same subspace, there exist an invertible matrix C such that

U S = E C -1 . (6)
As for E, let U S ↓ (respectively U S ↑ ) be the matrix extracted form U S by deleting the last (respectively the first) row. Then equations ( 5) and ( 6) yield

U S ↑ = U S ↓ Φ (7)
where Φ = C D C -1 . The Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) method [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF] consists in:

• computing the matrix Φ = U S ↓ † U S ↑ (where the symbol † denotes the Moore-Penrose pseudo-inverse; this computation requires O(LM 2 ) operations),

• extracting the estimated poles ẑm as the eigenvalues of Φ (which can be achieved in O(M 3 ) operations).

Finally, for m = 1, . . . , M , the m th estimated frequency and damping factor can be deduced using fm = angle(ẑm) 2π and dm = ln |ẑm|.

Estimation of the amplitudes and initial phases

The complex amplitudes {αm} 1≤m≤M can be determined by minimizing the LS criterion minα x -Eα 2 2 , where • x = [x(0), . . . , x(L -1)] T are the signal samples,

• α = [α1, . . . , αM , α * 1 , . . . , α * M ] T are complex amplitudes. The solution to this criterion is

α = E † x. (8)
Hence, for m = 1, . . . , M , the m th estimated real amplitude and initial phase are âm = 2| αm| and φm = angle( αm).

Note that the full computation of E † can be avoided since equation [START_REF] Duvaut | Traitement du signal[END_REF] shows that

E † = C -1 U S T where C = U S T E. Thus, α = C -1 U S T x can be computed in O(LM 2 ) opera- tions.

Re-synthesis

Once the model parameters have been estimated, the signal can be reconstructed using equation [START_REF] Purnhagen | HILN-the MPEG-4 parametric audio coding tools[END_REF]. Thus, the estimated signal sample at time t is

x(t) = M m=1 (xm(t) + x * m (t)) (9)
where xm(t) = αm ẑt m is the m th complex damped sinusoid. Note that equation ( 9) can be implemented in O(LM ) operations.

In a block processing context, some interpolation techniques are required in order to force the continuity of the parameters between consecutive blocks [START_REF] Serra | Spectral modeling synthesis : a sound system based on a deterministic plus stochastic decomposition[END_REF].

Pitch-scale modification

An immediate application of the EDS model is a frequency-scale modification of the signal, which just consists in multiplying the estimated frequencies fm by a same factor β. Thus, the frequency of the m th complex damped sinusoid in the modified signal is

f s m = β fm so that the corresponding pole is ẑs m = exp( dm + i 2π f s m ) = ẑm exp(i 2π(β -1) fm) (10)
Therefore, equation ( 9) becomes

xs (t) = M m=1 (x s m (t) + xs m * (t)) (11)
where xs m (t) = αm (ẑ s m ) t . Note that this pitch modification method is no more computationally demanding than the exact re-synthesis.

Adaptive signal processing

The section transposes the HR methods presented above in an adaptive context. It will be shown that tracking the slow variations of the model parameters leads to a very simple re-synthesis method.

Model parameters tracking

The Sintrack method for frequency estimation and tracking [START_REF] Duvaut | Traitement du signal[END_REF] consists in a two-steps estimation: the Matrix Pencil HR method [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF] is first applied to obtain the initial parameters, and the tracking is then achieved using an adaptive Least Mean Square (LMS) algorithm, the frequencies and damping factors being extracted from the roots of a backward prediction polynomial [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF]. When the prediction error exceeds a certain threshold, the algorithm switches back to the initialization step.

Although this method has proved to be successful on musical signals [START_REF] David | Sintrack analysis for tracking components of musical signals[END_REF], the lack of robustness of the LMS algorithm results in an intensive use of the Matrix Pencil method, which is very timeconsuming.

To avoid this increase of complexity, the prediction polynomial tracking can be replaced by a signal subspace tracking, since subspace-based HR methods are known to give more reliable estimates of the signal poles than linear prediction. Subspace tracking has been intensively studied in the fields of adaptive filtering, source localization or parameter estimation. A first class of tracking algorithms is based on the projection approximation hypothesis [START_REF] Abed-Meraim | Fast orthonormal PAST algorithm[END_REF]; an other one relies on EVD or SVD tracking techniques, derived from classical EVD or SVD algorithms.

Table 2: Sequential iteration EVD algorithm

Initialization: U S = I2M 0 (L-2M )×2M
For each time step t iterate:

A(t) = H (t) U S (t -1) fast matrix product A(t) = U S (t)R(t) skinny QR factorization
For example, the orthogonal iteration algorithm of Table 1 can be adapted to track the dominant eigenvectors of a sliding-window matrix

H (t) =      x(t -(L -1)) . . . x(t) x(t -(L -2)) . . . x(t + 1) . . . . . . . . . x(t) . . . x(t + L -1)      (12)
just by replacing the iteration index n in table 1 by the discrete time index t [14] (cf. table 2). Thus, only one iteration is completed at each time step. Once the signal subspace basis U S is computed, the standard ESPRIT method can be applied. However, for the sake of computational efficiency, adaptive implementations of ESPRIT have been developed [START_REF] Strobach | Fast recursive subspace adaptive ESPRIT algorithms[END_REF], which require O(LM 2 ) or O(LM ) operations at each time step.

Finally, the estimation of the amplitudes and initial phases can be achieved as in section 2.1.3. Equation ( 8) now becomes

α(t) = E(t) † x(t) (13)
where E(t) is the Vandermonde matrix of the estimated poles at time t, αm(t) and ẑm(t) denote the estimated m th complex amplitude and pole at time t, and x(t) = [x(t), . . . , x(t + L -1)] T . Since this estimation involves the matrix E defined in equation (4) for a time window [0 . . . L -1], it must be noted that αm(t) now is the complex amplitude of the m th damped sinusoid at time t.

Re-synthesis

In an adaptive context, since the complex amplitudes of the damped sinusoids are estimated at each time step, equation ( 9) stands with xm(t) = αm(t). Therefore, the re-synthesis of the signal at each time step just consists in summing the complex amplitudes, which only requires O(M ) operations.

Pitch scale modification

Let ϕm(t) be the phase shift between the m th estimated damped sinusoid and the m th synthesized damped sinusoid at time t, so that equation [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF] stands with xs m (t) = xm(t) exp(i ϕm(t)) = αm(t) exp(i ϕm(t)). [START_REF] Strobach | Square hankel SVD subspace tracking algorithms[END_REF] Since these sinusoids satisfy the following recurrences

• xm(t) = xm(t -1) ẑm(t),

• xs m (t) = xs m (t -1) ẑs m (t), equations ( 10) and [START_REF] Strobach | Square hankel SVD subspace tracking algorithms[END_REF] show that ϕm(t) can be recursively updated using the following scheme:

ϕm(t) = ϕm(t -1) + 2π(β -1) fm(t). (15)
Then, xs m (t) can be computed using equation ( 14), from which the synthesized sample xs (t) can be deduced using equation [START_REF] Kumaresan | Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise[END_REF]. Note that this pitch modification method has the same complexity as the exact re-synthesis.

EXPERIMENTAL RESULTS

This section illustrates first the enhancement of the coding quality using the EDS rather than a simple sinusoidal model, then the tracking and re-synthesis of musical signals. The study deals with two piano tones, C5 and G5, sampled at 11025 Hz.

It is to be noticed that in real audio signal applications, the data matrix H is never rank-deficient, because of the presence of noise. Moreover, the rank-truncation order 2M is unknown, and must be chosen carefully. Indeed, over-estimating M is harmless, but under-estimating M often generates biases in the estimates of the frequencies and damping factors. Then, L must be chosen much greater than M , in order to enforce the robustness of the HR method. On the other hand, the higher L is, the more this method is computationally demanding. Therefore, audio signals with a high number of sinusoids (typically low-pitched sounds) may first be decomposed into several sub-band signals (via filtering/decimating, as proposed in [START_REF] Laroche | The use of the Matrix Pencil method for the spectrum analysis of musical signals[END_REF]), before applying the HR method. In the examples proposed below, this pre-processing is useless, since the chosen piano tones have few sinusoidal components. modeled signals. Figure 1 shows a strong pre-echo (energy before the onset) with the sinusoidal model. Moreover, the global variation of the attack is wrongly estimated. Thanks to the exponentially time-varying amplitudes, the EDS model provides a better modeling since it creates a short pre-echo and offers a good reproduction of the attack. After several structural considerations, a frequency aspect is introduced in the analysis by using the polyphase 32-bands pseudo-QMF filter-bank of MPEG1-audio [START_REF] Banderburg | ISO-MPEG-1 audio: a generic standard for coding of high-quality digital audio[END_REF], which provides a uniform partition of the frequency axis. After that, a power and a SNR measure are computed in each sub-band, noted SNRT F .

EDS vs sinusoids

Figure 3 shows the better SNRT F values of the EDS model. Note that several sub-bands are not reliable due to their weak power (see figure 4).

Adaptive signal processing

This section illustrates first the frequency tracking of a musical signal composed of two piano tones, then the synthesis method of section 2.1.4. The parameters tracking method presented in section 2.2.1 has been tested on a piano signal: the C5 tone of figure 1-a is played at time t = 0s, then the G5 is played at time t = 0.36s, while the C5 is maintained (the time waveform and the spectrogram of this signal are plotted in figure 5). Figure 6 shows the result of the tracking. The model order and the window length where M = 16 and L = 160, and the sinusoids energies

Frequency tracking

(Em = α 2 m 1-|zm| 2L 1-|zm| 2 )
are represented on a logarithmic scale using gray levels for the plot. Since the number of sinusoids is over-estimated, it can be seen that spurious poles are detected in the low frequency band (below 1000Hz), which actually corresponds to the highest level of noise in the original signal.

Re-synthesis and pitch scale modification

The synthesis method proposed in section 2.1.4 gave excellent results on the piano tones: the synthesized sounds were perceptually very similar to the original ones.

The hearing sensation is particularly well reproduced at the attack of both sounds. This may be related to the spurious poles detection mentioned above. Their number and energy are greater at the attack, which allows a good representation of the mechanical noise. It is well known that this impact noise occurring during the action of the hammer on the strings is of great importance for the naturalness of the sound. This method could be directly implemented without any further modification. On the opposite, the pitch-scale modification requires additional work. Indeed, the recursion on the phase shift between the estimated and the pitch-shifted signal in equation ( 15) relies on several implicit assumptions, such as:

• the number of frequencies is constant through time,

• each pole characterizes one single time-varying frequency, which is present in the whole signal,

• ẑm(t) matches ẑm(t -1) (ie. the m th frequency trajectory is known).

In real audio signals, however, the frequencies may appear or disappear, so that their number changes throw time. Moreover, spurious frequencies are sometimes detected, and should be eliminated. Consequently, tracking the poles trajectories is a difficult problem.

In the literature, several strategies were proposed to track sinusoids in the presence of noise in a block processing context [START_REF] Mcaulay | Speech analysis and synthesis based on a sinusoidal representation[END_REF][START_REF] Serra | Spectral modeling synthesis : a sound system based on a deterministic plus stochastic decomposition[END_REF][START_REF] Levine | Audio representations for data compression and compressed domain processing[END_REF]. These methods were designed in association with frequency estimators based on the Short Time Fourier Transform (STFT), but they can easily be adapted to the EDS model and the HR methods.

Finally, the pitch-scale modification technique proposed in section 2.2.3 in combination with these classical frequency matching strategies proved to be successful on the piano tones.

Note that once the poles trajectories are estimated, the discrimination should be made between

• the harmonics (related to the pitch of the sound),

• the remaining poles, which model the signal noise.

A realistic pitch scale modification should change the frequencies of the first class and leave the second class unchanged. Of course, the classification of the poles would require additional work.

CONCLUSIONS

The EDS model is a useful tool for audio signals modeling. It leads to a better representation of signal frames than the undamped sinusoidal model for a coding purpose. The use of a HR algorithm achieves an accurate estimation which can be efficiently updated by tracking the signal subspace through time. Moreover, tracking the model parameters offers very interesting outlooks for signal re-synthesis and modification.

Figure 1 :

 1 Figure 1: Time-shape waveforms : original signal, sinusoidal modeling with MSIN = 16 and EDS modeling with MEDS = 12.

Figure 2 :

 2 Figure 2: Fourier spectra : (a) original signal, (b) sinusoidal modeling with MSIN = 16, (c) EDS modeling with MEDS = 12.

Figure 3 :

 3 Figure 3: SNRT F in dB.

Figure 4 :

 4 Figure 4: Power in sub-bands in dB.

Figure 5 :

 5 Figure 5: Time waveform and spectrogram of the piano signal for the frequency tracking test.

Figure 6 :

 6 Figure 6: Frequency tracking of the piano signal.

Table 1 :

 1 Orthogonal iteration EVD algorithm


Roland Badeau
email: roland.badeau@enst.fr

Bertrand David

WEIGHTED MAXIMUM LIKELIHOOD AUTOREGRESSIVE AND MOVING AVERAGE SPECTRUM MODELING

Keywords: Autoregressive moving average processes, Spectral domain analysis, Maximum likelihood estimation

We propose new algorithms for estimating autoregressive (AR), moving average (MA), and ARMA models in the spectral domain. These algorithms are derived from a maximum likelihood approach, where spectral weights are introduced in order to selectively enhance the accuracy on a predefined set of frequencies, while ignoring the other ones. This is of particular interest for modeling the spectral envelope of harmonic signals, whose spectrum only contains a discrete set of relevant coefficients. In the context of speech processing, our simulation results show that the proposed method provides a more accurate ARMA modeling of nasal vowels than the Durbin method.

INTRODUCTION

Parametric spectrum modeling is a prominent tool in time series analysis [START_REF] Kay | Modern Spectral Estimation, Theory and Application[END_REF]. The popular autoregressive moving average models are particularly useful for predicting the future values of a time series from its past samples. This is why the estimation task is naturally performed in the time domain. However, performing this task in the spectral domain may be relevant to process signals whose spectrum only partly satisfies the parametric model. For instance, a given region of the spectrum can be corrupted by an extraneous signal. In the case of harmonic spectra, the parametric model is appropriate for representing the envelope of the observed spectrum, which only contains a discrete set of relevant coefficients. The problem of partial spectrum modeling was first addressed by A. El-Jaroudi and J. Makhoul [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF], who introduced the Discrete All-Pole (DAP) method for estimating AR models. The DAP method included a frequencydependent weighting of the whole frequency range, so that the spectral accuracy could be enhanced in certain frequency regions compared to others. A similar approach was later presented in [START_REF] Weruaga | All-Pole Estimation in Spectral Domain[END_REF]. In [START_REF] Villavicencio | Improving LPC spectral envelope extraction of voiced speech by true-envelope estimation[END_REF], the True Envelope Linear Prediction Coding (TELPC) method was proposed for estimating the envelope of harmonic speech signals. Finally, the problem of ARMA modeling in a frequency subband was addressed in [START_REF] Broersen | Time series analysis in a frequency subband[END_REF]. In this paper, we generalize the DAP approach to the estimation of MA and ARMA models. Our estimator involves the iterative maximization of a weighted maximum likelihood function, which is performed by using a new algorithm, whose convergence is proved in the case of AR modeling.

The research leading to this paper was supported by the French GIP ANR under contract ANR-06-JCJC-0027-01, Décompositions en Éléments Sonores et Applications Musicales -DESAM, and by the Groupe des Écoles des Télécommunications (GET), under contract TAM-TAM. We also acknowledge Gaël Richard, Professor at ENST, who kindly lent his voice for the signal analyzed in section 4, and Jean-Louis Durrieu, Ph.D. student at ENST, for its feedback about this work.

PRINCIPLE

Given an observed power spectrum Y (ν), where the normalized frequency ν belongs to a subset of [0, 1), we propose a Weighted Maximum Likelihood (WML) approach for estimating a general ARMA model (which includes AR and MA models). The Power Spectral Density (PSD) of an ARMA process of order (P, Q) is of the form

S(ν) = σ 2 ∑ Q q=0 b q e -i2πνq 2 ∑ P p=0 a p e -i2πν p 2 = σ 2 B(ν) A(ν) , (1)
where σ > 0, a 0 = 1, and b 0 = 1. As usually assumed in the literature, we suppose that ∀ν, Y (ν) follows an exponential distribution of parameter S(ν): its probability density function

p ν satisfies log(p ν (Y (ν))) = log 1 S(ν) - Y (ν) S(ν) .
Assuming that the values Y (ν) for all frequencies are independent, we define the weighted loglikelihood function of the whole observed spectrum as a function of σ 2 and of the vectors a = [a 0 , . . . , a P ] T and b = [b 0 , . . . , b Q ] T :

L(σ 2 , a, b) = 1 0 1 + log Y (ν) S(ν) - Y (ν) S(ν) dµ(ν), ( 2
)
where we added the constant1 0 (1 + log (Y (ν)))dµ(ν), so that this expression is similar to that of the Itakura-Saito distance [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF]. Spectral weights are introduced via measure µ, which can be any positive measure such that 1 0 dµ(ν) = 1 (for instance, µ can typically be a discrete measure whose support uniformly samples the interval [0, 1)). Substituting equation (1) into equation [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF], and maximizing the log-likelihood with respect to parameter σ , yields the estimate

(σ ⋆ ) 2 = 1 0 Y (ν)A(ν) B(ν) dµ(ν). ( 3
)
Then substituting equations ( 1) and (3) into equation (2) yields

L(a, b) = 1 0 log Y (ν)A(ν) B(ν) dµ(ν)-log 1 0 Y (ν)A(ν) B(ν) dµ(ν) .
(4) The concavity of function log shows that L is non-positive, and attains the value 0 if and only if

Y (ν)A(ν) B(ν)
is constant in the support of measure µ. Thus maximizing the log-likelihood is equivalent to maximizing the flatness of the ratio Y (ν)A(ν) B(ν) . Since L is unchanged by inverting and conjugating any roots of the polynomials A(z) = ∑ P p=0 a p z p and B(z) = ∑ Q q=0 b q z q , we can look for vectors a and b which belong to the set S of vectors of first coefficient 1, and such that the polynomial roots are strictly outside the unit circle 1 .

MAXIMIZATION OF THE LOG-LIKELIHOOD

Autoregressive models

In this section, we investigate the case of AR models, for which function L can be straightforwardly optimized 2 . For all vectors a ∈ S , we parameterize the continuous spectrum in the form

A(ν) = P ∑ p=0 a p e -i2πν p 2 = a T T (ν)a, (5)
where ∀ν ∈ [0, 1), T (ν) is the Toeplitz matrix whose coefficient of index (p, q) is cos(2πν(pq)).

R(a ⋆ , µ) = R(a ⋆ , µ YA ⋆ ) (6)
which involves the (P + 1)

× (P + 1) Toeplitz matrices R(a, µ) = 1 0 T (ν) A(ν) dµ(ν) (7) R(a, µ YA ) = 1 0 Y (ν) T (ν) dµ(ν) 1 0 Y (ν)A(ν) dµ(ν) . ( 8
)
The proof of proposition 3.1 can be found in Appendix A, and relies on the gradient of L: 1 2 ∇L(a) = (R(a, µ) -R(a, µ YA )) a. In proposition 3.1, notation R(a, µ YA ) is to be understood as follows: for any measure µ on the interval [0, 1), and any positive function W of period 1, we define the measure dµ W (ν) =

W (ν) 1 0 W (ξ ) dµ(ξ ) dµ(ν).
In the case W = YA, the notation R(a, µ YA ) is compatible with R(a, µ). Similar notations will be used throughout the paper.

In practice, L always admits a global maximum a ⋆ in S , except in the singular case where the polynomial A ⋆ (z) has at least one root on the unit circle (see [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF] for more details). Below, we propose an iterative algorithm which solves equation [START_REF] Durbin | The Fitting of Time Series Models[END_REF] in S . Since ∇L(a ⋆ ) = 0, it can be noted that a ⋆ is a fixed point in S of the function

φ (a) = R(a, µ YA ) -1 R(a, µ)a. (9)
Then a first order expansion yields

L(φ (a)) -L(a) = (φ (a) -a) T ∇L(a) + o( φ (a) -a ) = 1 2 ∇L(a) T R(a, µ YA ) -1 ∇L(a) + o( φ (a) -a ) ≥ 0 in a neighborhood of a ⋆ .
(10) This suggests an ascent method for maximizing function L, which consists in recursively applying function φ to an initial point in S . In practice however, the polynomial defined by vector φ (a) may have some roots inside the unit circle. Thus φ (a) should be remapped into S by replacing those roots by their inverse conjugate. This remapping, which is denoted P S below, does not modify the value of function L. Moreover, we show in Appendix B that forcing vector a in S dramatically improves the convergence rate of the recursion. Finally, our algorithm consists of the following steps: Finite Impulse Response (FIR) filters, the vectors in S correspond to minimum phase filters. 2 Note that the proposed WML estimator is equivalent to the DAP method in the case of AR modeling. However, our new algorithm for computing the optimal solution has an enhanced convergence rate, as proved in Appendix B.

• Initialization : a = [1, 0 . . . 0] T • Repeat until convergence: a ← P S R(a, µ YA ) -1 R(a, µ)a
The convergence of this algorithm is analyzed in the following proposition, which is proved in Appendix B.

Proposition 3.2 (Convergence of the algorithm for AR estimation).

1. A global maximum a ⋆ of function L in S is a locally stable equilibrium point of the discrete dynamical system formed by the recursion a ← P S (φ (a)).

2. If µ is the Lebesgue measure, the recursion converges in one iteration only.

In practice, we observed that this algorithm globally converges to the global maximum (whatever the initial point is). Moreover, if µ is a measure close enough to the Lebesgue measure 3 , the algorithm converges in a few iterations, as shown in Appendix B.

Moving average models

R(b ⋆ , µ) = R(b ⋆ , µ Y /B ⋆ ) ( 11
)
which involves the (Q + 1) × (Q + 1) Toeplitz matrices R(b, µ) = 1 0 T (ν) B(ν) dµ(ν) (12) R(b, µ Y /B ) = 1 0 Y (ν) T (ν) B(ν) 2 dµ(ν) 1 0 Y (ν) B(ν) dµ(ν) . (13)
Proposition 3.3 can be proved in the same way as proposition 3.1. Note however that this proposition does not guarantee that the solutions of equation ( 11) are global maxima. In order to maximize function L, we propose the following iterative algorithm:

• Initialization : b = [1, 0 . . . 0] T • Repeat until convergence: b ← P S R(b, µ) -1 R(b, µ Y /B ) b
As in section 3.1, it can be verified that this recursion constitutes an ascent method for maximizing function L. In practice, we observed that this algorithm globally converges to the global maximum of function L; however our proof of convergence, which could not be included in this paper, relies on additional assumptions regarding the observed spectrum. Again, if µ is a measure close enough to the Lebesgue measure, the convergence rate is enhanced.

Autoregressive moving average models

We propose the following iterative algorithm for maximizing function L, which consists in interlacing the two updates introduced above:

• Initialization : a = b = [1, 0 . . . 0] T
• Repeat until convergence: 3 For instance, we say that a discrete measure, which uniformly samples the interval [0, 1) with constant weights, is "close to" the Lebesgue measure, in the sense that the integral of any continuous function with respect to this discrete measure approximates the integral of the same function with respect to the Lebesgue measure, when the number of samples is high enough.

-Repeat: a ← P S R(a, µ YA/B ) -1 R(a, µ)a

-Repeat: b ← P S R(b, µ) -1 R(b, µ YA/B ) b
where R(a, µ) and R(b, µ) are defined in equations ( 7) and (12), and

R(a, µ YA/B )) = 1 0 Y (ν) T (ν) B(ν) dµ(ν) 1 0 Y (ν)A(ν) B(ν) dµ(ν) R(b, µ YA/B ) = 1 0 Y (ν)A(ν) T (ν) B(ν) 2 dµ(ν) 1 0 Y (ν)A(ν) B(ν) dµ(ν)
.

The convergence was not proved, but observed in practical cases 4 .

SIMULATION RESULTS

We address the problem of estimating the spectral envelope of a voiced speech signal. In speech processing, AR modeling is a prominent tool, since it constitutes a simplified model of the speech production system. However, it is not well adapted to fit the spectra of nasal sounds, which contain zeros due to the coupling between the vocal and nasal cavities. Below we focus on a French nasal vowel, which is denoted /Ẽ/ in the International Phonetic Alphabet. This vowel was pronounced at a fundamental frequency F 0 = 120 Hz, and recorded at a sampling frequency of 8000 Hz. The spectrum of this original sound was then obtained by successively extracting a signal of length 60 ms, applying a Hann window, and computing a Digital Fourier Transform (DFT), zero-padded to N = 1024 samples. Then an ARMA model was estimated by applying the proposed WML method with weights equal to 1 for the DFT samples corresponding to multiples of F 0 , and 0 elsewhere 5 . A good fit was obtained with P = 8 and Q = 7. Finally, the vowel was resynthesized at F 0 = 100 Hz and 300 Hz respectively, by filtering a train of periodic pulses, and the spectra of those two synthetic signals were computed in the same way as that of the original one. The first experiment illustrated in Figure 1-a shows the inadequacy of AR models for modeling nasal vowels. The spectrum of the synthetic signal at F 0 = 100 Hz, represented by the solid gray line, is modeled by an AR model of order P = 15, represented by the solid black line 6 . The AR model was estimated with the DAP method [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF]. The formants are well represented by this AR model; however the model does not fit the zeros present at 1250, 2250, and 3900 Hz. The second experiment illustrated in Figure 1-b confirms that conversely, an ARMA model of orders P = 8 and Q = 7 can efficiently represent the spectral envelope. The dashed black line represents the result obtained with the classical Durbin method [START_REF] Durbin | The Fitting of Time Series Models[END_REF], and the solid black line represents the result obtained with the WML method. Both methods provide an accurate estimate of the original ARMA model. In the last experiment illustrated in Figure 1-c, both methods are applied to the synthetic signal at F 0 = 300 Hz, with P = 8 and Q = 7. The WML method still correctly estimates the positions of the zeros, even though one of them falls between two harmonics, whereas Durbin's method only provides a smooth approximation of the spectral envelope.

CONCLUSIONS

We proposed new algorithms for estimating AR, MA, and ARMA models in the spectral domain, when only a subset of the original spectrum is observed. This technique is derived from a Weighted Maximum Likelihood (WML) approach, and can be seen as a generalization of the DAP method, which was initially designed for AR models only [START_REF] El-Jaroudi | Discrete all-pole modeling[END_REF]. In the simple case of AR modeling, we proved that the proposed algorithm converges to the optimal solution, and that the convergence rate is enhanced by remapping the poles at each iteration. In the case of MA modeling, the proposed algorithm still converges to the optimal solution, although the proof could not be included in the paper. The case of ARMA models seems more complex, and we noted that the practical convergence requires a more careful implementation. In all cases, the convergence toward the optimal solution requires that the number of distinct parameters remains lower than the number of non-zero spectral weights. Finally, the algorithm for estimating ARMA models was successfully applied in the context of speech analysis. Our results first confirmed that ARMA models are more appropriate for representing nasal vowels than AR models, and the WML method proved to provide a much more accurate estimation of the envelope of high-pitched harmonic spectra than the classical Durbin method.

APPENDIX A. CHARACTERIZATION OF THE AR ESTIMATOR

Proof of proposition 3.1. Any critical point a ⋆ of the C 1 function L is such that the gradient ∇L(a ⋆ ) is zero. By substituting equation [START_REF] Broersen | Time series analysis in a frequency subband[END_REF] into equation [START_REF] Villavicencio | Improving LPC spectral envelope extraction of voiced speech by true-envelope estimation[END_REF] where B(ν) ≡ 1, and by differentiating the resulting expression, we obtain 1 2 ∇L(a) = (R(a, µ) -R(a, µ YA )) a, where ma- trices R(a, µ) and R(a, µ YA ) were defined in equations ( 7) and (8). Note that R(a, µ) -R(a, µ YA ) is a symmetric Toeplitz matrix. Let {r k (a)} k∈{-P...P} be the series formed by the coefficients in the first column and the first row of this matrix; this series is symmetric.

Since 1 2 ∇L(a ⋆ ) = 0, matrix R(a ⋆ , µ) -R(a ⋆ , µ YA ) is singular. How- ever classical results regarding singular Toeplitz matrices [START_REF] Badeau | High resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF] show that the series {r k (a ⋆ )} k∈{-P...P} is a linear combination of complex exponentials (possibly modulated by polynomials), whose poles are the roots of the polynomial A ⋆ (z). If a ⋆ ∈ S , all these roots are strictly outside the unit circle, and the symmetry of this series yields only one solution: r k (a ⋆ ) = 0 ∀k ∈ {-P . . . P}, which proves equation [START_REF] Durbin | The Fitting of Time Series Models[END_REF]. Now let show that any solution of equation ( 6) is a global maximum of function L. For any vector a, substituting equations ( 6) to (8) in to equation (4) yields

L(a) -L(a ⋆ ) = 1 0 log A(ν) A ⋆ (ν) dµ(ν) -log 1 0 A(ν)
A ⋆ (ν) dµ(ν) ≤ 0 because of the concavity of function log .

B. CONVERGENCE ANALYSIS FOR THE AR MODEL

Proof of assertion 1. in Prop. 3.2. Any equilibrium point of the dynamical system a ← P S (φ (a)) in S is the global maximum of function L. Now let prove that this equilibrium point is locally stable. According to the general theory of discrete dynamical systems, the local stability of an equilibrium point a ⋆ of a recursion of the form a ← f (a) is guaranteed if the Jacobian of f at point a ⋆ has a spectral radius strictly smaller than one. Here, we will prove that the Jacobian J φ (a ⋆ ) of function φ has all its eigenvalues of magnitude strictly lower than 1, except that associated to the eigenvector a ⋆ , which is equal to 1. This shows that the direction of the dynamical system locally converges to the direction of a ⋆ . Since in the neighborhood of a ⋆ ∈ S , the remapping P S consists in normalizing a vector so that its first coefficient becomes 1, we will then conclude that the dynamical system locally converges to a ⋆ . First, differentiating equation (9) shows that the Jacobian of φ satisfies J φ (a) = R(a, µ YA ) -1 R(a, µ) -P(a, µ) + 2u(a, µ)u(a, µ YA ) T (14) where

u(a, ν) = T (ν)a A(ν) (15) P(a, µ) = 2 1 0 u(a, ν)u(a, ν) T dµ(ν) (16) u(a, µ) = R(a, µ) a = 1 0 u(a, ν)dµ(ν) (17) u(a, µ YA ) = R(a, µ YA ) a. ( 18
)
Then noting that T (ν) can be written in the form T (ν) = ℜ e(ν)e(ν) H where e(ν) = [1, e i2πν , . . . , e i2πνP ] T , it can be verified that

P(a, µ) = R(a, µ) + H(a, µ), (19)
where the Toeplitz matrix R(a, µ) was defined in equation [START_REF] Badeau | High resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF], and H(a, µ) is a Hankel matrix:

R(a, µ) = ℜ 1 0 e(ν)e(ν) H |a T e(ν)| 2 dµ(ν) (20) H(a, µ) = ℜ 1 0 e(ν)e(ν) T (a T e(ν)) 2 dµ(ν) . (21)
Then substituting equation (19) into equation ( 14), we obtain

J φ (a) = -R(a, µ YA ) -1 H(a, µ) + 2φ (a)u(a, µ YA ) T .
Applying this equality to the equilibrium point a ⋆ yields

J φ (a ⋆ ) = a ⋆ u(a ⋆ , µ) T -ε(a ⋆ ), (22)
where ε(a

⋆ ) = R(a ⋆ , µ) -1 H(a ⋆ , µ) -u(a ⋆ , µ) u(a ⋆ , µ) T . (23)
Then it can be verified that ε(a ⋆ )a ⋆ = 0, and u(a ⋆ , µ) T ε(a ⋆ ) = 0 T . Thus zero is an eigenvalue of ε(a ⋆ ), associated to the right and left eigenvectors a ⋆ and u(a ⋆ , µ), which additionally satisfy u(a ⋆ , µ) T a ⋆ = 1. In other respects, noting that matrix R(a ⋆ , µ) ± H(a ⋆ , µ) is positive semi-definite, it can be proved that all the other eigenvalues of ε(a ⋆ ) have a magnitude strictly lower than 1. This finally proves the local stability of the dynamical system.

Proof of assertion 2. in Prop. Similar considerations yield H(a, µ) = u(a, µ)u(a, µ) T . Consequently, matrix ε(a ⋆ ) defined in equation ( 23) is zero. Thus the algorithm converges in one iteration only, which was expected, since u(a, µ) does not depend on a ∈ S .

If µ is a measure close to the Lebesgue measure, the spectral radius of ε(a ⋆ ) remains much lower than 1. Thus the discrete dynamical system locally converges to a ⋆ at a high convergence rate.

Proposition 3 . 1 (

 31 Characterization of the WML AR estimator). If the C 1 function L admits a critical point a ⋆ in S , then this point is the global maximum of L, characterized by the equation

Proposition 3 . 3 (

 33 Characterization of the WML MA estimator). If L admits a critical point b ⋆ in S , then this point satisfies the equation

Fig. 1 .

 1 Fig. 1. Estimation of the spectral envelope of a nasal vowel

3 . 2 .•

 32 If µ(ν) is the Lebesgue measure, Substituting equation (20) into equation (17) yieldsu(a, µ) = ℜ Γ is the unit circle and z = e 2iπν . The k-th coefficient is u k (a, µ) = ℜ 1 2iπ Γ f k (z)dz , where f k (z) = z k-1 P ∑ p=0 a p z p . If a ∈S and k > 0, the holomorphic function f k has no singularity inside Γ, thus u k (a, µ) = 0. If a ∈ S , f 0 has one singularity inside Γ (z = 0), thus the residue theorem yields u 0 (a, µ) = 1.• Equation (21) shows that H(a, µ) can be written in the form H(a, µ) = ℜ 1 2iπ Γ e(z)e(z) T (a T e(z))

For simplicity, the notation S does not account for the vector dimension, which can be P or Q when appropriate. When viewed as the coefficients of

Note that within one iteration of this algorithm, we perform several updates of a and b, because we observed that if only one update was performed for each vector, the algorithm could alternate between two sub-optimal solutions. Computing a few updates within each iteration fixed the

problem.[START_REF] Broersen | Time series analysis in a frequency subband[END_REF] Note that if F 0 is greater than the spectral width of the Hann window, the values of the smoothed periodogram at multiples of F 0 can still be considered independent, as assumed in

section 2.[START_REF] Durbin | The Fitting of Time Series Models[END_REF] We chose P = 15, so that the total number of parameters used for modeling the observed spectrum is the same as that of the original ARMA model.


Roland Badeau

Angélique Drémeau

VARIATIONAL BAYESIAN EM ALGORITHM FOR MODELING MIXTURES OF NON-STATIONARY SIGNALS IN THE TIME-FREQUENCY DOMAIN (HR-NMF)

Keywords: Nonnegative Matrix Factorization, High Resolution methods, Expectation-Maximization algorithm, Variational inference

Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)

INTRODUCTION

Nonnegative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] is a powerful tool for decomposing mixtures of non-stationary signals in the timefrequency (TF) domain. However, unlike the high resolution (HR) methods [START_REF] Badeau | High resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF] dedicated to mixtures of complex exponentials, its spectral resolution is limited by that of the underlying TF representation. Following previous works which aimed at providing a probabilistic framework for NMF [START_REF] Mikkel | Non-negative matrix factorization with Gaussian process priors[END_REF][START_REF] Smaragdis | Blind Speech Separation, chapter Probabilistic decompositions of spectra for sound separation[END_REF][START_REF] Virtanen | Bayesian extensions to nonnegative 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], we introduced in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] a unified probabilistic model called HR-NMF, that permits to overcome this limit by taking both phases and local correlations in each frequency band into account. It can be used with both complex-valued and realvalued TF representations (like the short-time Fourier transform or the modified discrete cosine transform). Moreover, we showed that HR-NMF generalizes some very popular models: 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, and the exponential sinusoidal model (ESM), commonly used in HR spectral analysis of time series [START_REF] Badeau | High resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF]. In [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF], HR-NMF was estimated with the expectation-maximization (EM) algorithm, which involves time-demanding Kalman filtering and smoothing. In this paper, we introduce two faster algorithms based on variational inference, and compare the performance of the three algorithms.

This paper is organized as follows. In Section 2, we present the HR-NMF model, as introduced in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF]. We recall the basics of the This work is supported by the French National Research Agency (ANR) as a part of the DReaM project (ANR-09-CORD-006-03) and partly supported by the Quaero Program, funded by OSEO.

variational Bayesian EM algorithm in Section 3, before particularizing it to the HR-NMF model in Section 4. Section 5 is devoted to experimental results, and conclusions are drawn in Section 6.

The following notation will be used throughout the paper:

• M * : conjugate of matrix (or vector) M ,

• M : transpose of matrix (or vector) M ,

• M H : conjugate transpose of matrix (or vector) M , • [M ; N ]: vertical concatenation of M and N ,

• c =: equality up to an additive constant, • h * m: discrete convolution of times series h and m,

• N F (µ, R): real (if F = R) or circular complex (if F = C)
multivariate normal distribution of mean µ and covariance matrix R.

HR-NMF TIME-FREQUENCY MIXTURE MODEL

The HR-NMF mixture model of TF data x(f, t) ∈ F (where F = R or C) is defined for all discrete frequencies 1 ≤ f ≤ F and times 1 ≤ t ≤ T as the sum of K latent components c k (f, t) ∈ F plus a white noise n(f, t) ∼ N F (0, σ 2 ):

x(f, t) = n(f, t) + K k=1 c k (f, t) (1)
where

• c k (f, t) = P (k,f ) p=1 a(p, k, f ) c k (f, t -p) + b k (f, t) is ob- tained by autoregressive filtering of a non-stationary signal b k (f, t) ∈ F (where a(p, k, f ) ∈ F and P (k, f ) ∈ N is such that a(P (k, f ), k, f) = 0), • b k (f, t) ∼ N F (0, v k (f, t)) where v k (f, t) is defined as v k (f, t) = w(k, f ) h(k, t), (2)
with w(k, f ) ≥ 0 and h(k, t) ≥ 0, • processes n and b1 . . . bK are mutually independent.

Moreover, ∀(k, f ) ∈ {1 . . . K} × {1 . . . F }, the random vec- tors c k (f, 0) = [c k (f, 0); . . . ; c k (f, -P (k, f ) + 1)
] are assumed to be independent and distributed according to the prior distribution

c k (f, 0) ∼ N F (µ k (f ), Q k (f ) -1 )
, where the mean µ k (f ) and the precision matrix Q k (f ) are fixed parameters 1 . Lastly, we assume that ∀f ∈ {1 . . . F }, ∀t ≤ 0, x(f, t) is unobserved.

Let c denote the set {c k (f, t)} (k,f,t) , x denote the set {x(f, t)} (f,t) and θ the set of model parameters σ

2 , {a(p, k, f )} (p,k,f ) , {w(k, f )} (k,f ) p(x; θ). (4)
The solution of ( 3)-( 4) can be found by means of an EM algorithm. We proposed in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] an efficient implementation, using a Kalman filter/smoother in the E-step. However, the computational cost remains high, dominated by the complexity of the Kalman filter/smoother, and may be prohibitive when dealing with large dimensions. We propose here an alternative, based on the variational Bayesian EM (VB-EM) algorithm, which uses a mean-field approximation of the posterior p(c|x; θ ) to reach a good compromise between quality and complexity of the MAP estimation (3).

VARIATIONAL BAYESIAN EM ALGORITHM

Variational inference [START_REF] Beal | Variational Algorithms for Approximate Bayesian Inference[END_REF][START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF] is now a classical approach for estimating a probabilistic model involving both observed variables x and latent variables c, parametrized by θ. Let F be a set of probability density functions (PDF) over the latent variables c. For any PDF q ∈ F and any function f (c), we note f q = f (c)q(c)dc.

Then for any PDF q ∈ F and any parameter θ, the log-likelihood L(θ) = ln(p(x; θ)) can be decomposed as

L(θ) = DKL(q||p(c|x; θ)) + L(q; θ) (5)
where

DKL(q||p(c|x; θ)) = ln q(c) p(c|x; θ) q (6)
is the Kullback-Leibler divergence between q and p(c|x; θ), and

L(q; θ) = ln p(c, x; θ) q(c) q (7)
is called the variational free energy. Moreover, L(q; θ) can be further decomposed as L(q; θ) = E(q; θ) + H(q), where

E(q; θ) = ln (p(c, x; θ)) q , (8)
and H(q) =ln (q(c)) q is the entropy of distribution q. Since DKL(q||p(c|x; θ)) ≥ 0, L(q; θ) is a lower bound of L(θ). The variational Bayesian EM algorithm is a recursive algorithm for estimating θ. It consists of the two following steps at each iteration i:

• E-step (update q): q = argmin q∈F DKL(q||p(c|x; θi-1)) = argmax q∈F L(q; θi-1) (9) • M-step (update θ): θi = argmax θ L(q ; θ) = argmax θ E(q ; θ). (10)
F defines a set of constraints leading to a particular approximation of the posterior distribution p(c|x; θi-1). We note that:

• In the standard EM algorithm, q is not constrained, thus q = p(c|x; θi-1) and DKL(q ||p(c|x; θi-1)) = 0. Therefore L(θi) ≥ L(q ; θi) ≥ L(q ; θi-1) = L(θi-1), which proves that the log-likelihood is non-decreasing. • In the general case, L(θ) is no longer guaranteed to be non-decreasing, but its lower bound L(q; θ) is still nondecreasing.

VARIATIONAL BAYESIAN EM FOR HR-NMF

Considering the HR-NMF model defined in [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF],

∀(k, f ) ∈ {1 . . . K}× {1 . . . F }, let c kf denote the set {c k (f, t)} t∈{-P (k,f )+1...T } . More- over, let α = 1 if F = C, and α = 2 if F = R. Then α ln(p(c, x)) = α K k=1 F f =1 ln(p(c kf )) +α F f =1 T t=1 δ(f, t) ln(p(x(f, t)|c1(f, t) . . . cK (f, t))) = -KF T + K k=1 F f =1 P (k, f ) ln(απ) - K k=1 F f =1 (c k (f, 0) -µ k (f )) H Q k (f )(c k (f, 0) -µ k (f )) + K k=1 F f =1 ln(det(Q k (f ))) + T t=1 ln(ρ k (f, t)) - K k=1 F f =1 T t=1 ρ k (f, t) c k (f, t) - P (k,f ) p=1 a(p, k, f )c k (f, t -p) 2 - F f =1 T t=1 δ(f, t) ln(απσ 2 ) + 1 σ 2 x(f, t) - K k=1 c k (f, t) 2 (11) where • δ(f, t) = 1 if x(f, t) is observed, and δ(f, t) = 0 else (in particular δ(f, t) = 0 ∀t < 1 and ∀t > T ), • ρ k (f, t) = 1 v k (f,t) if t ∈ {1 . . . T },
and ρ k (f, t) = 0 else. In the following subsections, we will first recall the EM-based algorithm presented in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] as a particular case of the variational procedure ( 9)-(10) (Sections 4.1 and 4.2) and then propose two different alternatives to this costly approach, based on two mean-field approximations, i.e. two different definitions of F (Sections 4.3 and 4.4). These three algorithms only differ in the E-step, but they share the same implementation of the M-step.

M-step

The M-step defined in equation [START_REF] Wainwright | Graphical Models, Exponential Families, and Variational Inference[END_REF] consists in maximizing E(q ; θ) w.r.t. the model parameters θ. First, equations ( 8) and (11) yield

αE(q ; θ) c = - F f =1 T t=1 δ(f, t) ln(απσ 2 ) + e(f, t)/σ 2 - K k=1 F f =1 T t=1 ln(w(k, f )h(k, t)) + a(k, f ) H S(k, f, t)a(k, f ) w(k, f )h(k, t) , (12)
where

• e(f, t) = δ(f, t) x(f, t) - K k=1 c k (f, t) 2 q , • S(k, f, t) = c k (f, t) * c k (f, t) q , • c k (f, t) = [c k (f, t); . . . ; c k (f, t -P (k, f ))], • a(k, f ) = [1; -a(1, k, f); . . . ; -a(P (k, f ), k, f)].
We note that e(f, t) and S(k, f, t) can be computed as

• e(f, t) = δ(f, t) x(f, t) - K k=1 m k (f, t) 2 + K k=1 Γ k (f, t) , • S(k, f, t) = R k (f, t) * + m k (f, t) * m k (f, t) ,
where we have defined:

• m k (f, t) = c k (f, t) q , • Γ k (f, t) = |c k (f, t) -m k (f, t)| 2 q , • m k (f, t) = c k (f, t) q , • R k (f, t) = (c k (f, t) -m k (f, t)) (c k (f, t) -m k (f, t)) H
q . The maximization of E(q ; θ) in equation ( 12), w.r.t. σ 2 , a(p, k, f ), w(k, f ), and h(k, t), can then be performed as in the M-step presented in [START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF], using the current estimations of m k (f, t) and R k (f, t) derived from the E-steps as presented in the next sections.

E-step in the exact EM algorithm

As mentioned in Section 3, in the exact EM algorithm q is not constrained, thus the solution of ( 9) is given by q = p(c|x; θ), and the variational free energy L(q , θi) is equal to the log-likelihood L(θi). In [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF], we showed that the posterior distribution p(c|x; θ) is Gaussian, and that its first and second order moments, as well as the value of L(θi), can be computed by means of Kalman filtering and smoothing. The resulting E-step can symbolically be written as:

for 1 ≤ f ≤ F do {m k (f, t), R k (f, t)} 1≤k≤K 1≤t≤T = Kalman ({x(f, t)} 1≤t≤T )
end for Its computational complexity was shown to be O(F T K 3 (1 + P ) 3 ), where P = max k,f P (k, f ).

E-step with structured mean field approximation

If K > 1, we assume that F , introduced in Section 3, is the set of PDFs which can be factorized in the form

q(c) = K k=1 F f =1 q kf (c kf ). (13)
Using this particular factorization for q(c), the solution of (9) satisfies (see [START_REF] Beal | Variational Algorithms for Approximate Bayesian Inference[END_REF]):

∀(k, f ) ∈ {1 . . . K} × {1 . . . F }, ln(q kf (c kf )) c = ln(p(c, x)) (l,g) =(k,f ) q lg . (14)
Then, reformulating equation (11) and using equation ( 14), we get

α ln(p(c, x)) c = α ln(p(c kf )) - T t=1 δ(f, t) σ 2 |c k (f, t) -ĉk (f, t)| 2 , ( 15
) with ĉk (f, t) = x(f, t) - l =k
m l (f, t). We observe that q kf is the posterior distribution of a HR-NMF model of order K = 1, where the posterior means of all components other than k have been subtracted to the observed data x(f, t). Hence q kf is Gaussian, and its first and second moments can be computed by applying the Kalman filter/smoother presented in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] to ĉk (f, t) instead of x(f, t). The resulting E-step can symbolically be written as:

for 1 ≤ f ≤ F do for 1 ≤ k ≤ K do ∀1 ≤ t ≤ T , ĉk (f, t) = x(f, t) - l =k m l (f, t) {m k (f, t), R k (f, t)} 1≤t≤T = Kalman({ĉ k (f, t)} 1≤t≤T ) end for end for
The complexity of this procedure is O(F T K(1 + P ) 3 ) instead of O(F T K 3 (1 + P ) 3 ) for the "classical" E-step.

In order to evaluate this algorithm, we are also interested in computing the variational free energy L. After some straightforward calculations, we note that the entropy H(q kf ) satisfies αH(q kf ) = (T + P (k, f )) (ln(απ) + 1)

+ T t=1 ln(det(R k (f, t))) - T -1 t=1 ln(det(R k (f, t))), (16)
where

R k (f, t) is the P (k, f ) × P (k, f ) top-left submatrix of R k (f, t).
Thus equations ( 7), ( 11) and ( 16) yield

αL(q; θ) = KF T + K k=1 F f =1 P (k, f ) -trace(Q k (f )R k (f, 0)) - K k=1 F f =1 (m k (f, 0) -µ k (f )) H Q k (f )(m k (f, 0) -µ k (f )) + K k=1 F f =1 ln(det(Q k (f ))) + T t=1 ln(ρ k (f, t)) - K k=1 F f =1 T t=1 ρ k (f, t) a(k, f ) H S(k, f, t)a(k, f ) - F f =1 T t=1 δ(f, t) ln(απσ 2 ) + e(f, t)/σ 2 + K k=1 F f =1 T t=1 ln(det(R k (f, t))) - T -1 t=1 ln(det(R k (f, t))) (17) where m k (f, 0) is the P (k, f ) × 1 top subvector of m k (f, 0).

E-step with mean field approximation

If P > 0, we further assume that F is the set of PDFs which can be factorized in the form

q(c) = K k=1 F f =1 T t=-(P (k,f )-1)
q kf t (c k (f, t)). ( 18
)
With this particular factorization of q(c), the solution of ( 9) satisfies (see [START_REF] Beal | Variational Algorithms for Approximate Bayesian Inference[END_REF]): ∀(k, f, t)∈{1 . . . K}×{1 . . . F }×{-P (k, f )+1 . . . T },

ln(q kf t (c k (f, t))) c = ln(p(c, x)) (l,g,u) =(k,f,t) q lgu . ( 19
)
Let us define the filter of impulse response h kf , such that h kf (0)= 1, h kf (p) = -a(p, k, f ) ∀p ∈ {1 . . . P (k, f )}, and h kf (p) = 0 everywhere else, and the filter hkf (p) = h kf (-p) * . After some straightforward calculations, equations ( 11) and ( 19)

yield ∀(k, f, t) ∈ {1 . . . K}×{1 . . . F }×{-P (k, f )+1 . . . T }, q kf t (c k (f, t)) ∼ N F (m k (f, t), Γ k (f, t)) , where 2 Γ k (f, t) = δ(f, t) σ 2 + q k (f, t) + | hkf | 2 * ρ k (f, t) -1 (with q k (f, t) = Q k (f ) (1-t,1-t) if -P (k, f ) + 1 ≤ t ≤ 0 and q k (f, t) = 0 else), and
m k (f, t) = m k (f, t) + Γ k (f, t) -q k (f, t) H (m k (f, 0) -µ k (f )) + δ(f,t) σ 2 (x(f, t) - K l=1 m l (f, t)) -hkf * (ρ k (f, t)(h kf * m k (f, t)))
2 | hkf | 2 denotes the filter whose coefficients are the square magnitude of the corresponding coefficients of hkf . (where 3 . The computational complexity of the E-step is thus further reduced from O(KF T (1 + P ) 3 ) to O(KF T (1 + P )), which is linear w.r.t. all dimensions.

q k (f, t) is the (1 -t) th column of Q k (f ) if -P (k, f ) + 1 ≤ t ≤ 0 and q k (f, t) = [0; . . . ; 0] else)
Finally, note that the variational free energy L(q; θ) can be calculated as in equation ( 17), where R k (f, t) becomes a diagonal matrix of diagonal coefficients Γ k (f, t) . . . Γ k (f, t -P (k, f )).

SIMULATION RESULTS

The VB-EM algorithm aims to maximize the free energy. As we mentioned in Section 3, the log-likelihood is thus no longer guaranteed to increase, while remaining an indicator of the estimation quality. It can then be interesting to evaluate the influence of the approximations ( 13) and (18) on the maximization of the log-likelihood. To this end, we consider a fully observed TF data x(f, t) generated according to model (1) with T =20, F =3, P (k, f )=3 ∀(k, f ) and K =2 (and random parameters θ), and compare the performance of the three algorithms described respectively in Subsections 4.2, 4.3 and 4.4 with regard to the maximization of the log-likelihood. Figure 1 presents the value of the log-likelihood at each iteration of the three algorithms. Interestingly, we can observe that although focusing on the maximization of the free energy, the VB-EM algorithm permits here to increase the log-likelihood, whatever the considered approximation (mean-field or structured mean-field). In addition, as intuitively expected, the most constrained factorization (18) leads to a lesser increase of the log-likelihood. In practice however, this expected quality loss is not tangible. As an example of the good behavior of the VB-EM approach, we focus here on a simple case of source separation, where the observation is the whole STFT x(f, t) (of dimensions F = 400 and T = 44) of a 1.05 s-long piano sound sampled at 11025 Hz, containing three notes, C3, C4 and C5, starting respectively at 0 ms, 260 ms and 525 ms, and lasting until the end of the sound. Within this scenario, we aim at separating K =3 components c k (f, t) of order P (k, f ) = 2 in the frequency band f which corresponds to the first harmonic of C5, the second harmonic of C4 and the fourth harmonic of C3 (around 540 Hz). These three sinusoidal components (whose real parts are represented as red solid lines in Figure 2) have very close frequencies, making them hardly separable. We compare then three different approaches, namely, the HR-NMF model estimated by means of the EM algorithm, the HR-NMF model estimated by means of the VB-EM algorithm using the mean-field approximation (18) and the IS-NMF model [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]. 3 In this equation, although the term m k (t) appears several times in the right-hand side, it can be easily verified that its contributions add up to zero. Two important observations can be made about Figure 2. As previously noticed in [START_REF] Badeau | Gaussian modeling of mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF][START_REF] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF], IS-NMF (in black dash-dotted lines), which involves Wiener filtering, is not able to properly separate the components when they overlap. As a comparison, the components estimated by HR-NMF (blue dashed lines and magenda dotted lines) better fit the ground truth. We see on this example that the EMbased and VB-EM-based approaches lead to very similar results: the separated components are often merged. More precisely, we measured an averaged mean squared error of 0.0161 for IS-NMF, 0.0016 for the VBEM-based HR-NMF and 0.0006 for the EM-based HR-NMF on the whole set of frequencies and components reconstructed within this experiment. The slight quality loss due to the mean-field approximation is largely compensated by a significant computation time saving: with a 2.20GHz CPU processor and 8Go RAM, the CPU time required to run the E-step in the exact EM approach with a Matlab implementation is 19.5s, while 1.9s is enough for the E-step with mean-field approximation.

CONCLUSIONS

This paper introduces two novel methods as alternatives to estimate 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] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF]. These methods are based on the variational Bayesian EM algorithm and two different mean-field approximations. Their low complexities allow using the HR-NMF model in high-dimensional problems without altering the good quality of the estimation. We illustrated these good properties with a simple example of source separation. In future work, we will investigate other kinds of structured and unstructured mean field approximations, as well as a fully Bayesian approach involving uninformative or informative priors for the various model parameters. We will also apply variational inference to the future extensions of the HR-NMF model (e.g. involving convolutive and multichannel mixtures).

Fig. 1 .

 1 Fig. 1. Maximization of the log-likelihood and the variational free energy over the iterations.

Fig. 2 .

 2 Fig. 2. Separation of three sinusoidal components.

In practice we choose µ k (f ) = [0; . . . ; 0] and Q k (f ) -1 = ξI, where I is the identity matrix and ξ is small relative to 1, in order to both enforce the causality of the latent components and avoid singular matrices.


Roland Badeau

Mark D Plumbley

PROBABILISTIC TIME-FREQUENCY SOURCE-FILTER DECOMPOSITION OF NON-STATIONARY SIGNALS

Keywords: Probabilistic modelling, Non-stationary processes, Time-frequency analysis, Source-filter models, Nonnegative matrix factorisation

Probabilistic modelling of non-stationary signals in the timefrequency (TF) domain has been an active research topic recently. Various models have been proposed, notably in the nonnegative matrix factorization (NMF) literature. In this paper, we propose a new TF probabilistic model that can represent a variety of stationary and non-stationary signals, such as autoregressive moving average (ARMA) processes, uncorrelated noise, damped sinusoids, and transient signals. This model also generalizes and improves both the Itakura-Saito (IS)-NMF and high resolution (HR)-NMF models.

INTRODUCTION

In the literature, several probabilistic models involving latent components have been proposed for modelling TF representations of audio signals such as spectrograms. Such models include NMF with additive Gaussian noise [START_REF] Schmidt | Non-negative matrix factorization with Gaussian process priors[END_REF], probabilistic latent component analysis (PLCA) [START_REF] Smaragdis | Probabilistic decompositions of spectra for sound separation[END_REF], NMF as a sum of Poisson components [START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF], and NMF as a sum of Gaussian components [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]. Although they have already proved successful in a number of audio applications such as source separation [START_REF] Smaragdis | Probabilistic decompositions of spectra for sound separation[END_REF][START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF] and multipitch estimation [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], most of these models still lack of consistency in some respects. Firstly, they focus on a magnitude or power TF representation, and simply ignore the phase information. As a consequence, reconstructing a consistent TF representation, including the phase field, often proves to be difficult [START_REF] Le Roux | Fast signal reconstruction from magnitude STFT spectrogram based on spectrogram consistency[END_REF][START_REF] Roux | Consistent Wiener filtering for audio source separation[END_REF]. Secondly, these models generally focus on the spectral and temporal dynamics, and assume that all time-frequency bins are independent. This assumption is not consistent with the existence of signal dynamics, because spectral dynamics always induces temporal dependencies, and temporal dynamics always induces spectral dependencies. In addition, further dependencies in the TF domain may be induced by the TF transform, due to spectral and temporal overlap between TF bins. In this paper, we address this problem by introducing a new probabilistic model called probabilistic time-frequency source-filter decomposition (PTFSFD), which aims to take both phases within TF bins and correlations between TF bins into account. It consists of a sum of source-filter models, where each component is obtained (in the original time domain) by successively applying a multiplication and a convolution to a white noise.

The paper is structured as follows. The filter bank used to compute the TF representation is introduced in Section 2. We then show in Sections 3 and 4 how convolutions and multiplications in the original time domain can be accurately implemented in the TF domain. The PTFSFD model is introduced in Section 5, and some examples are provided in Section 6. Finally, conclusions are drawn in Section 7.

DEFINITION OF THE FILTER BANK

In order to properly define the PTFSFD model in the TF domain, we need to accurately implement convolutions and multiplications. In the literature, the Short Time Fourier Transform (STFT) [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF] is often considered as a convenient TF transform, because under some smoothness assumptions it allows the approximation of convolution by multiplying each column of an STFT by the same spectrum, and of multiplication by multiplying each row of an STFT by the same waveform. However, STFT is an oversampled filter bank, and it transforms a real random process into a complex TF random field, whose distribution is generally approximated as circularly symmetric 1 . Instead, we propose to use a critically sampled perfect reconstruction (PR) cosine-modulated filter bank [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF], which has several desired properties. Firstly, uniform sampling in both time and frequency will allow an accurate implementation of convolutions and multiplications in the TF domain (cf. Sections 3 and 4). Secondly, real TF distributions will permit us to avoid the need for any circularlysymmetric assumption. The combination of PR, paraconjugate2 analysis/synthesis filters and critical sampling will al-

↓ F ↓ F h 0 h F -1 . . . . . . ↑ F ↑ F h0 hF -1 . . . . . . T TF

Time-domain transformation

Fig. 1. Analysis/synthesis filter bank low us to ensure the filter bank does not introduce any spurious correlation artefacts, concept that will be referred to as preservation of whiteness (PW), and explained in Section 5.1.

We thus consider in this paper a critically sampled cosinemodulated filter bank [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF], which transforms an input signal

x(n) ∈ R in the original time domain n ∈ Z into a 2D-array x(f, t) ∈ R in the TF domain (f, t) ∈ [0 . . . F -1] × Z. More precisely, x(f, t) is defined as x(f, t) = (h f * x)(F t),
where * denotes the standard convolution of non-periodic series over Z,

h f (n) = h(n) cos π F (f + 1 2 )(n + φ) is a real- valued analysis filter, F is an odd number so that φ = F +1
2 is a whole number, and h(n) is a prototype window of support [0 . . . N -1] with N = LF and L ∈ N, whose frequency passband is [-1 2F , 1 2F ] (modulo 1). The synthesis filters are defined as hf (n) = h f (N -n). This analysis/synthesis filter bank, represented in Figure 1 (T TF is defined as the TFdomain identity within Section 2), is designed so as to guarantee PR, which requires that L be even [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF]. This means that the output, defined as

x ′ (n) = f t hf (n -F t)x(f, t), satisfies x ′ (n) = x(n -N ).
The same PR property also guarantees that if x(f, t) is the input of the synthesis/analysis filter bank represented in Figure 2 (T TD is defined as the time-domain identity within Section 2), then the output is

x ′ (f, t) = x(f, t -L). Let H(ν) = n∈Z h(n)e -2iπνn
(with an upper case letter) denote the discrete time Fourier transform (DTFT) of h(n) over ν ∈ R. Considering that the time supports of h(F t 1 -n) and h(F t 2 -n) do not overlap provided that |t 1 -t 2 | ≥ L, we similarly define a whole number K, such that the overlap between the frequency supports of H f1 (ν) and H f2 (ν) can be neglected provided that |f 1 -f 2 | ≥ K, due to high rejection in the stopband. We note that h f (n) and x(f, t) are actually defined for all frequencies f ∈ Z, with the following properties:

• ∀n, t ∈ Z, h f (n) and x(f, t) are 2F -periodic w.r.t. f , • ∀n, t ∈ Z, h f (n) and x(f, t) are symmetric w.r.t. F -1 2 : ∀f ∈ Z, x(f, t) = x(2F -1 -f, t).

TF IMPLEMENTATION OF A CONVOLUTION

In this section, we consider a filter of impulse response

g(n) ∈ R and two signals x(n) ∈ R and y(n) ∈ R, such that y(n) = (g * x)(n).
Our purpose is to directly express y(f, t) as a function of x(f, t), i.e. to find a TF transformation T TF in Figure 1 such that the output of the filter bank

↓ F ↓ F h 0 h F -1 . . . . . . ↑ F ↑ F h0 hF -1 . . . . . .

TF-domain transformation

T TD Fig. 2. Synthesis/analysis filter bank is y(n -N ) when its input is x(n). achieve this, we use the PR property of this critically sampled filter bank, which implies that the unique solution T TF is that described in the larger frame in Figure 2, where the input is x(f, t), the output is y(f, t), and transformation T TD is defined as the time-domain convolution by g(n + N ). We then obtain

y(f, t) = F -1 ϕ=-F τ ∈Z c g (f, ϕ, τ ) x(f -ϕ, t -τ ), (1)
where

∀ϕ ∈ [-F . . . F -1], we have ∀f ∈ [0 . . . F -1], c g (f, ϕ, τ ) = (h f * hf-ϕ * g)(F (τ + L)), (2)
and ∀f ∈ [F . . . 2F -1], c g (f, ϕ, τ ) = c g (2F -1-f -ϕ, τ ),
where we have used the convention ∀ϕ /

∈ [0 . . . F -1], h ϕ = 0. This definition of c g (f, ϕ, τ ) is extended to all f ∈ Z and ϕ ∈ Z by 2F -periodicity. Note that ∀ϕ ∈ [-F . . . -K] ∪ [K . . . F -1]
, subbands f and f -ϕ do not overlap, thus c g (f, ϕ, τ ) can be neglected. Equation [START_REF] Schmidt | Non-negative matrix factorization with Gaussian process priors[END_REF] shows that a convolution in the original time domain is equivalent to a 2Dconvolution in the TF domain, which is stationary w.r.t. time, and circular of period 2F but non-stationary w.r.t. frequency.

In the following definition, we propose a parametric model of c g (f, ϕ, τ ). In Appendix A we prove that this model can represent any causal and stable ARMA filter g(n). [START_REF] Schmidt | Non-negative matrix factorization with Gaussian process priors[END_REF], where filters {c g (f, ϕ, τ )} 0≤f <F -F ≤ϕ≤F -1 are defined as the only stable solutions of the following equation: ∀f, ϕ, τ ,

Definition 1. The ARMA model of convolution in the TF domain is defined by

Qa u=0 a g (f -ϕ, u)c g (f, ϕ, τ -u) = b g (f, ϕ, τ ). (3)
In equation [START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF], the autoregressive term

a g (f -ϕ, u) is a causal sequence of support [0 . . . Q a ] (where Q a ∈ N) w.r.t.
u, having only simple poles lying inside the unit circle, the moving average term b

g (f, ϕ, τ ) is a sequence of support [-L + 1 . . . -L + 1 + Q b ] (where Q b ≥ 2L + Q a -1) w.r.t. τ , and ∀f ∈ [0 . . . F -1], ∀ϕ ∈ [-F . . . F -1], if f -ϕ / ∈ [0 . . . F -1] or |ϕ| > P b with P b = K -1, then b g (f, ϕ, τ )=0.
Remark 1. If a g (f, u) and b g (f, ϕ, τ ) are arbitrary, then there may be no filter g(n) such that equation (2) holds, which means that the operation defined in Definition 1 does no longer correspond to a convolution in the original time domain. In this case, we will say that the model is inconsistent.

TF IMPLEMENTATION OF A PRODUCT

In this section, we consider a sequence σ(n) ∈ R and two signals w

(n) ∈ R and x(n) ∈ R, such that x(n) = σ(n)w(n).
As in the previous section, our purpose is to directly express x(f, t) as a function of w(f, t), i.e. to find a TF transformation T TF in Figure 1 such that the output of the filter bank is x(n -N ) when its input is w(n). To achieve this, we also use the PR property of the filter bank, which implies that the unique solution is that represented inside the larger frame in Figure 2, where the input is w(f, t), the output is x(f, t), and transformation T TD is defined as the time-domain convolution by δ(n + N ) (δ is the Kronecker symbol: δ(u) = 1 if u = 0, δ(u) = 0 otherwise) followed by the product by σ(n).

Mathematically, we obtain x(f, t) = n∈Z σ(n)h f (F t - n) F -1 ϕ=0 τ ∈Z h ϕ (F τ -n)w(ϕ, τ ), which is equivalent to x(f, t) = (-1) f t F -1 ϕ=-F τ ∈Z p σ (t, ϕ, τ ) ×(-1) (f -ϕ)(t-τ ) w(f -ϕ, t -τ ) (4)
where

p σ (t, ϕ, τ ) = (-1) ϕt 2 m∈Z σ(F t -m)h(m) ×h(m -F τ ) cos π F ϕ(m + φ) + π 2 τ . (5)
Here, it is important to note that ∀t ∈ Z,

• ∀τ ∈ Z, p σ (t, ϕ, τ ) is 2F -periodic w.r.t. ϕ, • if τ is even, then p σ (t, ϕ, τ ) is symmetric w.r.t. ϕ, • if τ is odd, then p σ (t, ϕ, τ ) is antisymmetric w.r.t. ϕ, ∀τ ∈ Z, p σ (t, ϕ, -τ ) = (-1) τ p σ (t + τ, ϕ, τ ), • ∀τ / ∈ [-L + 1 . . . L -1], p σ (t, ϕ, τ ) = 0.
Equation [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] shows that a product in the original time domain is equivalent to a 2D-convolution in the TF domain, which is circular and stationary w.r.

α σ (t -τ, ω)p σ (t, ϕ -ω, τ ) = u∈Z β σ (t, ϕ -2F u, τ ). (6

)

In equation [START_REF] Roux | Consistent Wiener filtering for audio source separation[END_REF], the autoregressive term α σ (t -τ, ω) is a non-periodic symmetric sequence of bounded support [-P α . . . P α ] (where P α < F -K) w.r.t. ω, the moving average term β σ (t, ϕ, τ ) is a non-periodic sequence of bounded support [-P β . . . P β ] (where K + P α ≤ P β < F ) w.r.t. ϕ, and

• if τ is even, then β σ (t, ϕ, τ ) is symmetric w.r.t. ϕ, • if τ is odd, then β σ (t, ϕ, τ ) antisymmetric w.r.t. • ∀τ / ∈ [-Q β . . . Q β ] with Q β = L -1, β σ (t, ϕ, τ ) = 0.
Remark 2. If α σ (t, ω) and β σ (t, ϕ, τ ) are arbitrary, then there may be no sequence σ(n) such that equation ( 5) holds, which means that the operation defined in Definition 2 does not correspond to a multiplication in the original time domain. In this case, we will say that the model is inconsistent.

TF PROBABILISTIC MODELLING

The PTFSFD model introduced in Section 5.2 relies on the fundamental PW property presented in Proposition 1.

Preservation of whiteness

Proposition 1 (Preservation of whiteness). The filter bank defined in Section 2 guarantees that a stochastic process w(n) is white noise of variance σ 2 if and only if its TF counterpart w(f, t) is 2D-white noise of same variance σ 2 .

Proof.

If w(n) is white noise of variance σ 2 , then ∀f 1 , f 2 ∈ [0 . . . F -1], ∀t 1 , t 2 ∈ Z, E(w(f 1 , t 1 )w(f 2 , t 2 )) = σ 2 × n∈Z h f1 (F t 1 -n)h f2 (F t 2 -n) = σ 2 δ(f 1 -f 2 )δ(t 1 - t 2 )
, where E denotes mathematical expectation. The second equality comes from the PR property of the synthesis/analysis filter bank (cf. Figure 2, where T TD is defined as the timedomain identity). Thus w(f, t) is 2D-white noise of variance

σ 2 . Reciprocally, if w(f, t) is 2D-white noise of vari- ance σ 2 , then E(w(n 1 )w(n 2 )) = σ 2 F -1 f =0 t∈Z h f (F t - n 1 )h f (F t -n 2 ) = σ 2 δ(n 1 -n 2 )
, where the second equality comes from the PR property of the analysis/synthesis filter bank (cf. Figure 1, where T TF is defined as the TF-domain identity). Thus w(n) is white noise of variance σ 2 .

Definition of the PTFSFD model

Consider a mixture signal y(n) = R-1 r=0 y r (n), where each component y r (n) follows a source-filter model: y r = g r * x r , where g r (n) an impulse response of the form defined in Appendix A, x r = σ r × w r , where σ r (n) is a sequence of the form defined in Appendix B, and w r (n) is white noise of variance 1. Then in the TF domain, y(f, t) satisfies the PTFSFD model presented in Definition 3.

Definition 3. The PTFSFD model expresses a TF representation

y(f, t) (with 0 ≤ f ≤ F -1) as a sum of R uncorrelated components {y r (f, t)} 0≤r<R , where y r (f, t) = F -1 ϕ=-F τ ∈Z c gr (f, f -ϕ, t -τ )x r (ϕ, τ ), x r (f, t) = (-1) f t F -1 ϕ=-F τ ∈Z p σr (t, f -ϕ, t-τ ) ((-1) ϕτ w r (ϕ, τ )) ,
w r (f, t) is 2D-white noise of variance 1, and c gr (f, ϕ, τ ) and p σr (t, ϕ, τ ) are parametrised as in Definitions 1 and 2.

DEMONSTRATION EXAMPLES

The PTFSFD model introduced in Definition 3 encompasses the following particular cases:

1) ARMA process: If ∀t, p r (t, ϕ, τ )=δ 2F (ϕ)δ(τ ) (where δ 2F (ϕ) = u∈Z δ(ϕ -2F u)), then y r (n) can be any ARMA process (an example is given in Figures 3-1a and3-1b).

2) Uncorrelated noise: If ∀f, c gr (f, ϕ, τ ) = δ 2F (ϕ)δ(τ ), y r (n) belongs to a parametric family of uncorrelated noises (an example is represented in Figures 3-2a and3-2b).

3) Impulses: If ∀f, c gr (f, ϕ, τ ) = δ 2F (ϕ)δ(τ ) and p r (t, ϕ, τ ) is defined as in Example 1 in Appendix B with ε → 0, then y r (t) is a linear combination of impulses (an example is represented in Figures 3-3a and3-3b).

4) Damped sinusoids: If p r (t, ϕ, τ ) is defined as in Example 1 in Appendix B with P = 1 and ε → 0, and if c r (f, ϕ, τ ) is defined as in Appendix A with N 0 = 0, then y r (n) follows a real-valued exponential sinusoidal model (ESM) commonly used in HR spectral analysis of time series [START_REF] Hayes | Statistical Digital Signal Processing And Modeling[END_REF] (an example is represented in Figures 3-4a and3-4b).

The PTFSFD model also includes some NMF frameworks. Indeed, if W and H are F × R and R × T nonnegative matrices of entries W f r and H rt , if w r (f, t) is 2D-white Gaussian noise, and if p σr (t, ϕ, τ

) = √ H rt δ 2F (ϕ)δ(τ ) then: 5) if c gr (f, ϕ, τ )= W f r δ 2F (ϕ)δ(τ ), y(f, t) follows an IS-NMF model [4] of order R, defined by the NMF V =W H; 6) if c r (f, ϕ, τ ) = W f r δ 2F (ϕ)a f r (τ )
where a f r (τ ) is the impulse response of a causal and stable autoregressive filter, y(f, t) follows an 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] Badeau | High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] of order R, defined by the NMF V =W H and the autoregressive filters a f r .

Unlike cases (1)-( 4), IS-NMF and HR-NMF are inconsistent in the sense defined in Remarks 1 and 2 in Sections 3 and 4. Besides, IS-NMF cannot accurately generate damped sinusoids as in Figure 3-4a, and neither IS-NMF nor HR-NMF can accurately generate impulses as in Figure 3-3a.

CONCLUSIONS

In this paper, we proposed a new probabilistic model of nonstationary signals, which consists of a sum of source-filter models. Each component of the sum is obtained by successively applying a multiplication and a convolution to white noise. In order to properly express this model in the TF domain, we used a critically sampled cosine-modulated filter bank, which allows an accurate implementation in the TF domain of convolutions and multiplications in the original time domain, by means of 2D-filters. We then proposed an ARMA parametrisation of these 2D-filters, which can represent a variety of transformations, ranging from smooth to sharply selective in time or in frequency. We have demonstrated that the resulting probabilistic time-frequency source-filter decomposition (PTFSFD) model can represent a broad range of stationary and non-stationary signals, including ARMA processes and uncorrelated noise, and it has both high spectral 1) is an ARMA process, ( 2) is an uncorrelated noise, ( 3) is an impulse, ( 4) is a damped sinusoid 3and temporal resolutions. PTFSFD also generalizes and improves both IS-NMF and HR-NMF: it can accurately generate damped sinusoids, which IS-NMF cannot; and importantly it can also accurately generate impulses, which is not possible for either IS-NMF or HR-NMF.

Because audio signals are sparse in the time-frequency domain, we observed that the PTFSFD model involves a small number of non-zero parameters in practice. In future work, we will investigate enforcing this property, for instance by introducing an a priori distribution of the parameters inducing sparsity. The proposed approach could also be applied to other types of filter banks, e.g. STFT filter banks. We also plan to propose some algorithms for estimating this model from the samples of an input signal. In particular, the variational Bayesian methods already investigated in [START_REF] Badeau | Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF] for the HR-NMF model seem to be a promising approach, because of their low complexity and their capacity to efficiently estimate complex graphical models.

Acknowledgements. 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 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.

A. ARMA MODEL OF CONVOLUTION

We consider the TF implementation of convolution presented in Section 3, and we define g(n) as the impulse response of a causal and stable ARMA filter, having only simple poles. Then the partial fraction expansion of its transfer function shows that it can be written in the form g

(n) = g 0 (n) + P k=1 g k (n), where P ∈ N, g 0 (n) is a causal sequence of support [0 . . . N 0 -1] (with N 0 ∈ N), and ∀k ∈ [1 . . . P ], g k (n) = A k e δ k n cos(2πν k n + ψ k )1 n≥0 , where A k > 0, δ k < 0, ν k ∈ [0, 1 2 ], ψ k ∈ R. Then ∀f , equation (2) yields c g (f, ϕ, τ ) = P k=0 c g k (f, ϕ, τ ) with c g0 (f, ϕ, τ ) = (h f * hf-ϕ * g 0 )(F (τ + L)) and ∀k ∈ [1 . . . P ], c g k (f, ϕ, τ ) = e δ k F τ (A k (f, ϕ, τ ) cos(2πν k F τ ) + B k (f, ϕ, τ ) sin(2πν k F τ )), where we defined A k (f, ϕ, τ ) = A k N -1 n=-N +1 (h f * hf-ϕ )(n+N )e -δ k n cos(2πν k n-ψ k )1 n≤F τ and B k (f, ϕ, τ ) = A k N -1 n=-N +1 (h f * hf-ϕ )(n + N )e -δ k n sin(2πν k n -ψ k )1 n≤F τ . It can be easily proved that: • the support of c g0 (f, ϕ, τ ) is [-L+1 . . . L+⌈ N0-2 F ⌉] w.r.t. τ , • if τ ≤ -L, then c g0 (f, ϕ, τ ), A k (f, ϕ, τ ) and B k (f, ϕ, τ ) are zero, thus c g (f, ϕ, τ ) = 0, • if τ ≥ L, A k (f, ϕ, τ ) and B k (f, ϕ, τ ) do not depend on τ .
It results that ∀f, ϕ, c g (f, ϕ, τ -L + 1) is the impulse response of a causal and stable ARMA filter.

As a particular case, suppose that ∀k ∈ [1 . .

. P ], |δ k | ≪ 1. If τ ≥ L, then A k (f, ϕ, τ ) and B k (f, ϕ, τ
) can be neglected as soon as ν k does not lie in the supports of both H f (ν) and H f -ϕ (ν). Thus for each f and ϕ, there is a limited number P (f, ϕ) ≤ P (possibly 0) of c g k (f, ϕ, τ ) which contribute to c g (f, ϕ, τ ). In the general case, we can still consider without loss of generality that ∀f, ϕ, there is a limited number P (f, ϕ) ≤ P of c g k (f, ϕ, τ ) which contribute to c g (f, ϕ, τ ). We then define Q a 2max

f,ϕ P (f, ϕ)

and

Q b = 2L + Q a -1 + ⌈ N0-2 F ⌉.
Then it can be easily proved that c g (f, ϕ, τ ) matches Definition 1.

B. ARMA MODEL OF MULTIPLICATION

We consider the TF implementation of multiplication presented in Section 4. Suppose that σ 0 (n) is a band-limited sequence, so that there is 0 < ν 0 < 1 2 -K 2F such that the DTFT Σ 0 (ν) is zero ∀ν ∈ [ν 0 , 1 2 ] (we note that if ν 0 ≪ 1, then σ 0 (n) has slow temporal variations). Then the DTFT of the sequence σ 0 (F t -m)h(m)h(m -F τ ) is zero ∀ν ∈ [ν 0 + K 2F , 1 2 ]. Consequently, equation ( 5) yields p σ0 (t, ϕ, τ ) = 0 if K + ⌈2F ν 0 ⌉ ≤ |ϕ| ≤ F . Besides, let P ∈ N, and ∀k ∈ [1 . .

. P ], let σ k (n) = σ ′ k (n) 2(cos( π(n-φ) F )-cos(θ k ))
, where + ε 2 )sinc( n-n k 2F ) (where sinc(u) = sin(πu) πu ). In particular, ν 0 = 3 4F . Then when ε → 0, σ k (n) → A k δ(n -n k ). Asymptotically, we thus get σ(n) = P k=1 A k δ(n -n k ), i.e. a linear combination of impulses. Moreover, in this case, if either t or t -τ / ∈ ⌈ n k F ⌉ . . . ⌈ n k F ⌉ + L -1 , then p σ k (t, ϕ, τ ) = 0. Thus for each t and τ , there is a finite number P (t, τ ) ≤ P (possibly 0) of p σ k (t, ϕ, τ ) which contribute to p σ (t, ϕ, τ ).

In the general case, we can still consider without loss of generality that ∀t, τ , there is a limited number P (t, τ ) ≤ P of p σ k (t, ϕ, τ ) which contribute to p σ (t, ϕ, τ ). Moreover, we suppose that P α = max t,τ P (t, τ ) and ν 0 are such that P α < F -K and ν 0 < 1 2 -K+Pα 2F , and we define P β = K + P α -1+⌈2F ν 0 ⌉ < F . Then it can be easily proved that p σ (t, ϕ, τ ) matches Definition 2.

Fig. 3 .

 3 Fig. 3. Examples of signals generated by the PTFSFD model (F =15, L=8, K=4; y r (n) are in the left column and y r (f, t) in the right column): (1) is an ARMA process, (2) is an uncorrelated noise, (3) is an impulse, (4) is a damped sinusoid3

F θ k π ∈

 π R\Z and σ ′ k (n) is a band-limited sequence, such that the DTFT Σ ′ k (ν) is zero ∀ν ∈ [ν 0 , 1 2 ]. Then equation (5) yields p σ k (t, ϕ + 1, τ ) -2 cos(θ k )p σ k (t, ϕ, τ ) + p σ k (t, ϕ -1, τ ) = p σ ′ k (t, ϕ, τ ).Finally, we define σ(n) = P k=0 σ k (n). Example 1 (Linear combination of impulses). As an example, suppose that σ 0 (n) = 0, and ∀k ∈ [1 . . . P ], θ k = π(φ-n k ) F + ε (where n k ∈ Z, all n k are distinct, and ε > 0), and σ ′ k (n) = 2A k ε sin( π(2φ-n-n k ) 2F

A complex scalar random variable Z is circularly symmetric if ∀Ψ ∈ R, the random variables Z and Ze iΨ have the same probability distribution.

The paraconjugation of a discrete linear filter consists of the time reversal and complex conjugation of its impulse response.

The numerical values of the model parameters that we used to generate those signals are listed in our Matlab code, which is available upon request.