Roland Badeau

Mark D Plumbley

Multichannel high resolution NMF for modelling convolutive mixtures of non-stationary signals in the time-frequency domain

I. INTRODUCTION

N ONNEGATIVE matrix factorisation was originally intro- duced as a rank-reduction technique, which approximates a non-negative matrix V ∈ R F ×T as a product V ≈ W H of two non-negative matrices W ∈ R F ×S and H ∈ R S×T with S < min(F, T ) [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF]. In audio signal processing, it is often used for decomposing a magnitude or power TF representation, such as a Fourier or a constant-Q transform (CQT) spectrogram. The columns of W are then interpreted as a dictionary of spectral templates, whose temporal activations are represented in the rows of H. Several applications to audio have been addressed, such as multi-pitch estimation [START_REF] Raczyński | Multipitch analysis with harmonic nonnegative matrix approximation[END_REF]- [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF], automatic music transcription [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF], [START_REF] Bertin | Enforcing harmonicity and smoothness in Bayesian non-negative matrix factorization applied to polyphonic music transcription[END_REF], musical instrument recognition [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multiway Data Analysis and Blind Source Separation[END_REF], and source separation [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[END_REF]- [START_REF] Liutkus | Informed source separation using latent components[END_REF].

In the literature, several probabilistic models involving latent components have been proposed to provide a probabilistic framework to NMF. 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 nonnegative 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] Schmidt | Non-negative matrix factorization with Gaussian process priors[END_REF]- [START_REF] Virtanen | Bayesian extensions to nonnegative 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 modelling a magnitude or power TF representation, and simply ignore the phase information. In an application of source separation, the source estimates are then obtained by means of Wiener-like filtering [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[END_REF]- [START_REF] Liutkus | Informed source separation using latent components[END_REF], which consists in applying a mask to the magnitude TF representation of the mixture, while keeping the phase field unchanged. It can be easily shown that this approach cannot properly separate sinusoidal signals lying in the same frequency band, which means that the frequency resolution is limited by that of the TF transform. In other respects, the separated TF representation is generally not consistent, which means that it does not correspond to the TF transform of a temporal signal, resulting in artefacts such as musical noise. Therefore enhanced algorithms are needed to reconstruct a consistent TF representation [START_REF] Roux | Consistent Wiener filtering for audio source separation[END_REF]. In the same way, in an application of model-based audio synthesis, where there is no available phase field to assign to the sources, reconstructing consistent phases requires employing ad-hoc methods [START_REF] Griffin | Signal reconstruction from short-time Fourier transform magnitude[END_REF], [START_REF] Le Roux | Fast signal reconstruction from magnitude STFT spectrogram based on spectrogram consistency[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 clearly not relevant in the case of sinusoidal or impulse signals for instance, and it is not consistent with the existence of spectral or temporal dynamics. Indeed, in the case of wide sense stationary (WSS) processes, spectral dynamics (described by the power spectral density) is closely related to temporal correlation (described by the autocovariance function). Reciprocally, in the case of uncorrelated processes (all samples are uncorrelated with different variances), temporal dynamics induces spectral correlation. In other respects, further dependencies in the TF domain may be induced by the TF transform, due to spectral and temporal overlap between TF bins.

In order to overcome the assumption of independent TF bins, Markov models have been introduced for taking the local dependencies between contiguous TF bins of a magnitude or power TF representation into account [START_REF] Ozerov | Factorial scaled hidden Markov model for polyphonic audio representation and source separation[END_REF]- [START_REF] Mysore | Non-negative hidden Markov modeling of audio with application to source separation[END_REF]. However, these models still ignore the phase information. Conversely, the complex NMF model [START_REF] Kameoka | Complex NMF: A new sparse representation for acoustic signals[END_REF], [START_REF] Le Roux | Complex NMF under spectrogram consistency constraints[END_REF], which was explicitly designed to represent phases alongside magnitudes in a TF representation, is based on a deterministic framework that does not represent statistical correlations.

In order to model both phases and correlations within frequency bands in a principled way, 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]High resolution NMF for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF] a new model called high resolution (HR) NMF. We showed that this model offers an improved frequency resolution, able to separate sinusoids within the same frequency band, and an improved synthesis capability, able to restore missing TF observations. It can be used with 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 implemented in the TF domain by following the exact approach proposed in [START_REF] Badeau | Probabilistic time-frequency sourcefilter 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 fast variational EM algorithm. This paper further develops a previous work presented in [START_REF]Multichannel HR-NMF for modelling convolutive mixtures of non-stationary signals in the time-frequency domain[END_REF], by providing a theoretical ground for the TF implementation of convolution.

The paper is structured as follows. The filter bank used to compute the TF representation is presented in Section II. We then show in Section III how convolutions in the original time domain can be accurately implemented in the TF domain. The multichannel HR-NMF model is introduced in section IV, and the variational EM algorithm is derived in section V. An application to a stereophonic piano signal is presented in section VI. Finally, conclusions are drawn in section VII.

II. TIME-FREQUENCY ANALYSIS

In the literature, the STFT [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF] is often preferred to other existing TF transforms, because under some smoothness assumptions it allows the approximation of linear filtering by multiplying each column of the STFT by the frequency response of the filter. Instead, we propose to use the more general and flexible framework of perfect reconstruction (PR) filter banks [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF]. Indeed, we will show in Section III that PR actually permits us to accurately implement convolutions in the TF domain.

We thus consider a filter bank [START_REF] Vaidyanathan | Multirate Systems and Filter Banks[END_REF], which transforms an input signal x(n) ∈ l ∞ (F) in the original time domain n ∈ Z (where F = R or C and l ∞ (F) denotes the space of bounded sequences on F) into a 2D-array

x(f, t) ∈ l ∞ (F) ∀f ∈ [0 . . . F -1] in the TF domain (f, t) ∈ [0 . . . F -1] × Z. More precisely, x(f, t) is defined as x(f, t) = (h f * x)(Dt)
, where D is the decimation factor, * denotes standard convolution, and h f (n) is an analysis filter of support [0 . . . N -1] with N = LD and L ∈ N. The synthesis filters hf (n) of same support [0 . . . N -1] are designed so as to guarantee PR. This means that the output, defined as

x ′ (n) = F -1 f =0 t∈Z hf (n -Dt)x(f, t), satisfies x ′ (n) = x(n -N ), which corresponds to an overall delay of N samples. Let H f (ν) = n∈Z h f (n)e -2iπνn
(with an upper case letter) denote the discrete time Fourier transform (DTFT) of overlap provided that |t 1t 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 1f 2 | ≥ K, due to high rejection in the stopband.

h f (n) over ν ∈ R. Considering that the time supports of h f (Dt 1 -n) and h f (Dt 2 -n) do not ↓ F ↓ F h0 hF -1 . . . . . . ↑ F ↑ F h0 hF -1 . . . . . . TTF Time-domain transformation TTD x(n) y(n-N) x(f, t) y(f, t) (a

III. TF IMPLEMENTATION OF CONVOLUTION

In this section, we consider a stable filter of impulse response g(n) ∈ l 1 (F) (where l 1 (F) denotes the space of sequences on F whose series is absolutely convergent) and two signals x(n) ∈ l ∞ (F) and y(n) ∈ l ∞ (F), such that y(n) = (g * x)(n). Our purpose is to directly express the TF representation y(f, t) of y(n) as a function of x(f, t), i.e. to find a TF transformation T TF in Figure 1(a) such that if the input of the filter bank is x(n), then the output is y(n-N ) (y is delayed by N samples in order to take the overall delay of the filter bank into account). The following developments further investigate and generalize the study presented in [START_REF] Badeau | Probabilistic time-frequency sourcefilter decomposition of non-stationary signals[END_REF], which focused on the particular case of critically sampled PR cosine modulated filterbanks. The general case of stable linear filters is first addressed in section III-A, then the particular case of stable recursive filters is addressed in section III-B.

A. Stable linear filters

The PR property of the filter bank implies that the relationship between y(f, t) and x(f, t) is given by the transformation T TF described in the larger frame in Figure 1(b), 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 ). The resulting mathematical expression is given in Proposition 1.

Proposition 1. Let g(n) ∈ l 1 (F) be the impulse response of a stable linear filter, and

x(n) ∈ l ∞ (F) and y(n) ∈ l ∞ (F) two signals such that y(n) = (g * x)(n).
Let y(f, t) and x(f, t) be the TF representations of these signals as defined in Section II. Then

y(f, t) = ϕ∈Z τ ∈Z c g (f, ϕ, τ ) x(f -ϕ, t -τ ) (1) f t τ ϕ c g (f, ϕ, τ ) x(f, t) * y(f, t) Fig. 2. TF implementation of convolution where ∀f ∈ [0 . . . F -1], ∀ϕ ∈ Z, ∀τ ∈ Z, c g (f, ϕ, τ ) = (h f * hf-ϕ * g)(D(τ + L)), ( 2
)
with the convention ∀f / ∈ [0 . . . F -1], h f = 0.
Proof: Equations ( 1) and ( 2) are obtained by successively substituting equations y

(n) = (g * x)(n) and x(n) = F -1 f =0 t∈Z hf (n-D(t-L))x(f, t) into y(f, t) = (h f * y)(Dt).
Remark 1. As mentioned in section II, if |ϕ| ≥ K, then subbands f and fϕ do not overlap, thus c g (f, ϕ, τ ) can be neglected.

Equation [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] shows that a convolution in the original time domain is equivalent to a 2D-convolution in the TF domain, which is stationary w.r.t. time, and non-stationary w.r.t. frequency, as illustrated in Figure 2.

B. Stable recursive filters

In this section, we introduce a parametric family of TF filters based on a state space representation, and we show a relationship between these TF filters and equation (1).

Definition 1. Stable recursive filtering in TF domain is defined by the following state space representation

: ∀f ∈ [0 . . . F -1], t ∈ Z, z(f, t) = x(f, t) - Qa τ =1 a g (f, τ )z(f, t -τ ) y(f, t) = P b ϕ=-P b τ ∈Z b g (f, ϕ, τ ) z(f -ϕ, t -τ ) (3) where Q a ∈ N, P b ∈ N, and ∀f ∈ [0 . . . F -1], x(f, t) ∈ l ∞ (F) is the sequence of input variables, z(f, t) ∈ l ∞ (F) is the sequence of state variables, and y(f, t) ∈ l ∞ (F) is the sequence of output variables. The autoregressive term a g (f, τ ) ∈ F is a causal sequence of support [0 . . . Q a ] w.r.t.
τ (with a g (f, 0) = 1), having only simple poles lying inside the unit circle. The moving average term b g (f, ϕ, τ ) ∈ F is a sequence of finite support w.r.t. τ , and ∀f ∈

[0 . . . F -1], ∀ϕ ∈ [-P b . . . P b ], b g (f, ϕ, τ ) = 0 provided that f -ϕ / ∈ [0 . . . F -1]. Proposition 2. If g(n) ∈ l 1 (F)
is the impulse response of a causal and stable recursive filter, then the TF input/output system defined in Proposition 1 admits the state space representation [START_REF] Smaragdis | Relative pitch tracking of multiple arbitrary sounds[END_REF], where

P b = K -1 and ∀f ∈ [0 . . . F -1], ∀ϕ ∈ [-P b , P b ], b g (f, ϕ, τ ) is a sequence of support [-L + 1 . . . -L + 1 + Q b ] w.r.t. τ , where Q b ≥ 2L + Q a -1.
Proposition 2 is proved in Appendix A.

Proposition 3. In Definition 1, equation (3) can be rewritten in the form of equation

(1), where ∀f ∈ [0 . . . F -1], ∀τ ∈ Z, c g (f, ϕ, τ ) = 0 if |ϕ| > P b , and ∀f ∈ [0 . . . F -1], ∀ϕ ∈ [-P b . . . P b ], filter c g (f, ϕ, τ
) is defined as the only stable (bounded-input, bounded-output) solution of the following recursion:

∀τ ∈ Z, Qa u=0 a g (f -ϕ, u)c g (f, ϕ, τ -u) = b g (f, ϕ, τ ). (4) Proposition 3 is proved in Appendix A. Remark 2.
In Definition 1, a g (f, τ ) and b g (f, ϕ, τ ) are overparametrised compared to g(n) in Proposition 1. Consequently, if a g (f, τ ) and b g (f, ϕ, τ ) are arbitrary, then it is possible that no filter g(n) exists such that equation (2) holds, which means that this state space representation does no longer correspond to a convolution in the original time domain. In this case, we will say that the TF transformation defined in equation ( 3) is inconsistent.

IV. MULTICHANNEL HR-NMF

In this section we present the multichannel HR-NMF model, initially introduced in [START_REF]Multichannel HR-NMF for modelling convolutive mixtures of non-stationary signals in the time-frequency domain[END_REF]. Before defining HR-NMF in the TF domain in section IV-B, we first provide an intuitive interpretation of this model in the time domain.

A. HR-NMF in the time domain

The HR-NMF model of a multichannel signal y m (n) ∈ F is defined for all channels m ∈ [0 . . . M -1] and times n ∈ Z, as the sum of S source images y ms (n) ∈ F plus a Gaussian noise n m (n) ∈ F: y m (n) = n m (n)+ S-1 s=0 y ms (n). Moreover, each source image y ms (f, t) for any s ∈ [0 . . . S -1] is defined as y ms (n) = (g ms * x s )(n), where g ms is the impulse response of a causal and stable recursive filter, and x s (n) is a Gaussian process 1 . Additionally, processes x s and n m for all s and m are mutually independent. In order to make this model identifiable, we will further assume that the spectrum of x s (n) is flat, because the variability of source s w.r.t. frequency can be modelled within filters g ms for all m. Thus filter g ms represents both the transfer from source s to sensor m and the spectrum of source s. In section IV-B, recursive filters g ms will be directly implemented in the TF domain via equations ( 7) and ( 8), following Definition 12 .

B. HR-NMF in the TF domain

The multichannel HR-NMF model of TF data

y m (f, t) ∈ F 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 source images y ms (f, t) ∈ F plus a 2D-white noise n m (f, t) ∼ N F (0, σ 2 y ), (5)
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 :

y m (f, t) = n m (f, t) + S-1 s=0 y ms (f, t). (6)
Each source image y ms (f, t) for s ∈ [0 . . . S -1] is defined as

y ms (f, t) = P b ϕ=-P b Q b τ =0 b ms (f, ϕ, τ ) z s (f -ϕ, t -τ ) (7)
where

P b , Q b ∈ N, b ms (f, ϕ, τ ) = 0 if f -ϕ / ∈ [0 . . . F -1]
, and the latent components z s (f, t) ∈ F are defined as follows:

• ∀t ∈ [0 . . . T -1], x s (f, t) ∼ N F (0, σ 2 xs (t)) and z s (f, t) = x s (f, t) - Qa τ =1 a s (f, τ )z s (f, t -τ ) (8)
where Q a ∈ N and a s (f, τ ) defines a stable autoregressive filter,

• ∀t ∈ [-Q z . . . -1] where Q z = max(Q b , Q a ), z s (f, t) ∼ N (µ s (f, t), 1/ρ s (f, t)). (9)
Moreover, the random variables n m (f 1 , t 1 ) and

x s (f 2 , t 2 ) for all s, m, f 1 , f 2 , t 1 , t 2 are assumed mutually independent. Additionally, ∀m ∈ [0 . . . M -1], ∀f ∈ [0 . . . F -1], ∀t ∈ [-Q z . . . -1], y m (f, t) is unobserved, and ∀s ∈ [0 . . . S -1],
the prior mean µ s (f, t) ∈ F and the prior precision (inverse variance) ρ s (f, t) > 0 of the latent variable z s (f, t) are considered to be fixed parameters.

The set θ of parameters to be estimated consists of:

• the autoregressive parameters a s (f, τ . This model encompasses the following special cases:

) ∈ F for s ∈ [0 . . . S -1], f ∈ [0 . . . F -1], τ ∈ [1 . . . Q a ] (we further define a s (f, 0) = 1), • the moving average parameters b ms (f, ϕ, τ ) ∈ F for m ∈ [0 . . . M -1], s ∈ [0 . . . S -1], f ∈ [0 . . . F -1], ϕ ∈
• If M = 1, σ 2 y = 0 and P b = Q b = Q a = 0, then equa- tion (6) reduces to y 0 (f, t) = S-1 s=0 b 0s (f, 0, 0)x s (f, t), thus y 0 (f, t) ∼ N F (0, V f t ), where matrix V of co- efficients V f t is defined by the NMF V = W H with W f s = |b 0s (f, 0, 0)| 2 and H st = σ 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 f t = |y 0 (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 = 1 and P b = Q b = 0, then y 0 (f, t) follows the monochannel 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 non-stationary 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

a s (f, τ ) for all s ∈ [0 . . . S -1], f ∈ [0 . . . F -1] and τ ∈ [1 . . . Q a ],
and the NMF V = W H. , y m (f, t) can be written in the form y m (f, t) = Qa τ =1 α mτ z τ (f ) t , where z 1 (f ) . . . z Qa (f ) are the roots of the polynomial z Qa + Qa τ =1 a 0 (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]. Because it generalizes both IS-NMF and ESM models to multichannel data, the model defined in equation ( 6) is called multichannel HR-NMF.

• If S = 1, σ 2 y = 0, P b = 0, σ 2 x0 (t) = 1 ∀t ∈ [0 . . . T -1], and µ s (f, t) = 0 and ρ s (f, t) = 1 ∀t ∈ [-Q z . . . -1], then ∀m ∈ [0 . . . M -1], ∀f ∈ [0 . . . F -1], y m (f, t) is an autoregressive moving average (ARMA) process [25, Section 3.6]. • If S = 1, σ 2 y = 0, P b = 0, Q a > 0, Q b = Q a -1, ∀t ∈ [-Q z . . . -1], µ 0 (f, t) = 0, ρ 0 (f, t) ≫ 1,

V. VARIATIONAL EM ALGORITHM

In early works that focused on monochannel HR-NMF [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 non-stationary signals in the time-frequency domain[END_REF], in order to estimate the model parameters we proposed to resort to an expectation-maximization (EM) algorithm based on a Kalman filter/smoother. The approach proved to be appropriate for modelling audio signals in applications such as source separation and audio inpainting. However, its computational cost was high, dominated by the Kalman filter/smoother, and prohibitive when dealing with high-dimensional signals.

In order to make the estimation of HR-NMF faster, we then proposed two different strategies. The first approach aimed to improve the convergence rate, by replacing the M-step of the EM algorithm by multiplicative update rules [START_REF] Badeau | Multiplicative updates for modeling mixtures of non-stationary signals in the time-frequency domain[END_REF]. However we observed that the resulting algorithm presented some numerical stability issues. The second approach aimed to reduce the computational cost, by using a variational EM algorithm, where we introduced two different variational approximations [START_REF] Badeau | Variational Bayesian EM algorithm for modeling mixtures of non-stationary signals in the time-frequency domain (HR-NMF)[END_REF]. We observed that the mean field approximation led to both improved performance and maximal decrease of computational complexity.

In this section, we thus generalize the variational EM algorithm based on mean field approximation to the multichannel HR-NMF model introduced in section IV-B, as proposed in [START_REF]Multichannel HR-NMF for modelling convolutive mixtures of non-stationary signals in the time-frequency domain[END_REF]. 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], novelties also include a reduced computational complexity and a parallel implementation.

A. 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, determined by a set θ of parameters. Let F be a set of probability density functions (PDFs) 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 set of parameters θ, the variational free energy is defined as L(q; θ) = ln p(y, z; θ) q(z)

q . ( 10
)
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 ) (11)
• Maximization (E)-step (update θ):

θ i = argmax θ L(q ⋆ ; θ). (12)
In the case of multichannel HR-NMF, θ has been specified in section IV-B. We further define

δ m (f, t) = 1 if y m (f, t) is observed, otherwise δ m (f, t) = 0, in particular δ m (f, t) = 0 ∀(f, t) / ∈ [0 . . . F -1] × [0 . . . T -1].
The complete set of variables 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 z s (f, t) for s ∈ [0 . . . S -1], f ∈ [0 . . . F -1], and t ∈ [-Q z . . . T -1].
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 (z s (f, t)). (13)
With this particular factorization of q(z), the solution of ( 11) is such that each PDF q sf t is Gaussian: z s (f, t) ∼ N F (z s (f, t), γ zs (f, t)).

B. Variational free energy

Let α = 1 if F = C, and α = 2 if F = R. Let D y = M-1 m=0 F -1 f =0 T -1 t=0
δ m (f, t) be the number of observations, and

I(f, t) = ½ {0≤f <F, 0≤t<T } , e ym (f, t) = δ m (f, t) y m (f, t) - S-1 s=0 y ms (f, t) , x s (f, t) = I(f, t) Qa τ =0 a s (f, τ )z s (f, t -τ ) .
Then using equations ( 5) to ( 9), the joint log-probability distribution L = log(p(y, z; θ)) of the complete set of variables satisfies

-αL = -α (ln(p(y|z; θ)) + ln(p(z; θ))) = (D y + SF (T + Q z )) ln(απ) +D y ln(σ 2 y ) + 1 σ 2 y M-1 m=0 F -1 f =0 T -1 t=0 |e ym (f, t)| 2 + S-1 s=0 F -1 f =0 -1 t=-Qz ln( 1 ρs(f,t) ) + S-1 s=0 F -1 f =0 -1 t=-Qz ρ s (f, t)|z s (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) |x s (f, t)| 2 .
Thus the variational free energy defined in [START_REF] Liutkus | Informed source separation using latent components[END_REF] satisfies

-αL = D y ln(απ) -SF (T + Q z ) +D y 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)γ zs (f, t)) +ρ s (f, t) γ zs (f, t) + |z s (f, t) -µ s (f, t)| 2 + S-1 s=0 F -1 f =0 T -1 t=0 ln σ 2 xs (t) γz s (f,t) + γx s (f,t)+|xs(f,t)| 2 σ 2 xs (t) (14) where ∀f ∈ [0 . . . F -1], ∀t ∈ [0 . . . T -1], γ ey m (f, t) = δ m (f, t) S-1 s=0 γ yms (f, t), γ yms (f, t) = P b ϕ=-P b Q b τ =0 |b ms (f, ϕ, τ )| 2 γ zs (f -ϕ, t -τ ), e ym (f, t) = δ m (f, t) y m (f, t) - S-1 s=0 y ms (f, t) , y ms (f, t) = P b ϕ=-P b Q b τ =0 b ms (f, ϕ, τ ) z s (f -ϕ, t -τ ), γ xs (f, t) = I(f, t) Qa τ =0 |a s (f, τ )| 2 γ zs (f, t -τ ) , x s (f, t) = I(f, t) Qa τ =0 a s (f, τ )z s (f, t -τ ) .

C. Variational EM algorithm for multichannel HR-NMF

According to the mean field approximation, the maximizations in equations ( 11) and ( 12) are performed for each scalar parameter in turn [START_REF] Mackay | Information Theory, Inference, and Learning Algorithms[END_REF]. The dominant complexity of each iteration of the resulting variational EM algorithm is 4M F ST ∆f ∆t, where ∆f = 1 + 2P b and ∆t = 1 + Q z . 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.

1) E-step:

For all s ∈ [0 . . . S -1], f ∈ [0 . . . F -1], t / ∈ [-Q z , -1], let ρ s (f, t) = 0.
Considering the mean field approximation [START_REF] Virtanen | Bayesian extensions to nonnegative matrix factorisation for audio signal modelling[END_REF], the E-step defined in equation ( 11) leads to the updates described in Table I (where * denotes complex conjugation). Note that z s (f, t) has to be updated after γ zs (f, t).

2) M-step: The M-step defined in (12) leads to the updates described in Table II. The updates of the four parameters can be processed in parallel.

VI. SIMULATION RESULTS

In this section, we present a basic proof of concept of the multichannel HR-NMF model introduced in section IV-B. The following experiments deal with a single source (S = 1) formed of a real piano sound sampled at 11025 Hz. A 1.25msshort stereophonic signal (M = 2) has been synthesized by filtering the monophonic recording of a C3 piano note with two room impulse responses simulated using the Matlab code

parfor s ∈ [0 . . . S-1], f ∈ [0 . . . F -1], t ∈ [-Qz . . . T -1] do γz s (f, t) -1 = ρs(f, t) + Qa τ =0 I(f,t+τ )|as(f,τ )| 2 σ 2 xs (t+τ ) + M -1 m=0 P b ϕ=-P b Q b τ =0 δm (f +ϕ,t+τ )|bms(f +ϕ,ϕ,τ )| 2 σ 2 y end parfor for s ∈ [0 . . . S-1], f 0 ∈ [0 . . . ∆f -1], t 0 ∈ [-Qz . . .-Qz+∆t-1] do parfor f -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 τ =0 as(f,τ ) * xs(f,t+τ ) σ 2 xs (t+τ ) - M -1 m=0 P b ϕ=-P b Q b τ =0 bms (f +ϕ,ϕ,τ ) * ey m (f +ϕ,t+τ ) σ 2
y end parfor end for

TABLE I E-STEP OF THE VARIATIONAL EM ALGORITHM σ 2 y = 1 Dy 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) + |xs(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-τ )-xs(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 II M-STEP OF THE VARIATIONAL EM ALGORITHM

presented in [START_REF] Mcgovern | A model for room acoustics[END_REF] 3 . The TF representation y m (f, t) of this signal has then been computed by applying a critically sampled PR cosine modulated filter bank (F = R) with F = 201 frequency bands, involving filters of length 8F = 1608 samples. The resulting TF representation, of dimension F × T with T = 77, is displayed in Figure 3. In particular, it can be noticed that the two channels are not synchronous (the starting time in the left channel is ≈ 0.04s, whereas it is ≈ 0.02s in the right channel), which suggests that the order Q b of filters b ms (f, ϕ, τ ) should be chosen greater than zero.

In the following experiments, we have set µ s (f, t) = 0 and ρ s (f, t) = 10 5 . These values force z s (f, t) to be close to 3 Those impulse responses were simulated using 15625 virtual sources. The dimensions of the room were [20m, 19m, 21m], the coordinates of the two microphones were [19m, 18m 1.6m] and [15m, 11m, 10m], and those of the source were [5m, 2m, 1m]. The reflection coefficient of the walls was 0.3. zero ∀t ∈ [-Q z . . . -1] (since the prior mean and variance of z s (f, t) are µ s (f, t) = 0 and 1/ρ s (f, t) = 10 -5 ), which is relevant if the observed sound is preceded by silence. The variational EM algorithm is initialized with the neutral values z s (f, t) = 0,

γ zs (f, t) = σ 2 y = σ 2 xs (t) = 1, a s (f, τ ) = ½ {τ =0}
, and b ms (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. We suppose that all TF points within the frame in Figure 3 In each experiment, 100 iterations of the algorithm are performed, and the restored signal is returned as y ms (f, t).

In the first experiment, a multichannel HR-NMF with Q a = Q b = P b = 0 is estimated. Similarly to the example provided in section IV-B, this is equivalent to modelling the two channels by two rank-1 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 W and sharing the same temporal activation H, or by a rank 1 multichannel NMF [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 4. It can be noticed that wherever y m (f, t) is observed (δ m (f, t) = 1), y ms (f, t) does not accurately fit y m (f, t) (this is particularly visible in high frequencies), because the length Q b of filters b ms (f, ϕ, τ ) has been underestimated: the source to distortion ratio (SDR) in the observed area is 11.7dB. In other respects, the missing observations (δ m (f, t) = 0) could not be restored (y ms (f, t) is zero inside the frame, resulting in an SDR of 0dB in this area), because the correlations between contiguous TF coefficients in y m (f, t) have not been taken into account.

In the second experiment, a multichannel HR-NMF model with Q a = 2, Q b = 3, and P b = 1 is estimated. The resulting TF representation y ms (f, t) is displayed in Figure 5. It can be noticed that wherever y m (f, t) is observed, y ms (f, t) better fits y m (f, t): the SDR is 36.8dB in the observed area. Besides, the missing observations have been better estimated: the SDR is 4.8dB inside the frame. Actually, choosing P b > 0 was necessary to obtain this result, which means that the spectral overlap between frequency bands cannot be neglected.

VII. CONCLUSIONS

In this paper, we have shown that convolution can be accurately implemented in the TF domain, by applying 2Dfilters to a TF representation obtained as the output of a PR filter bank. In the particular case of recursive filters, we have also shown that filtering can be implemented by means of a state space representation in the TF domain. These results have then been used to extend the monochannel HR-NMF model initially proposed 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 non-stationary signals in the time-frequency domain[END_REF] to multichannel signals and convolutive mixtures. The resulting multichannel HR-NMF model can accurately represent the transfer from each source to each sensor, as well as the spectrum of each source. It 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, which has a reduced computational complexity and a parallel implementation 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]. This algorithm has been successfully applied to a stereophonic piano signal, and has been capable of modelling reverberation due to room impulse response, and restoring missing observations.

Because audio signals are sparse in the time-frequency domain, we observed that the multichannel HR-NMF 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 [START_REF] Wipf | Sparse Bayesian learning for basis selection[END_REF]. In order to deal with more realistic music signals, the estimation of HR-NMF should be performed in a more informed way, for instance by means of semi-supervised learning, or by using any kind of prior information about the sources. For instance, harmonicity and temporal or spectral smoothness could be enforced by introducing some prior distributions of the parameters, or by reparametrising the model. The model could also be extended in several ways, for instance by taking the correlations over latent components into account, or by using other types of TF transforms, e.g. wavelet transforms. 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.

ACKNOWLEDGMENT

This work was undertaken while Roland Badeau was visiting the Centre for Digital Music, partly funded by EPSRC Platform Grant EP/K009559/1. Mark D. Plumbley is funded by EPSRC Leadership Fellowship EP/G007144/1.

APPENDIX A TF IMPLEMENTATION OF STABLE RECURSIVE FILTERING

Proof of Proposition 2:

We consider the TF implementation of convolution given in Proposition 1, and we define g(n) as the impulse response of a causal and stable recursive filter, having only simple poles. Then the partial fraction expansion of its transfer function [START_REF] Cheng | Analysis of Linear Systems[END_REF] shows that it can be written in the form g(n) = g 0 (n) + where we have defined

A k (f, ϕ, τ ) = A k N -1 n=-N +1 (h f * hf-ϕ )(n + N ) ×e -δ k n cos(2πν k n -ψ k )½ n≤Dτ , B k (f, ϕ, τ ) = A k N -1 n=-N +1 (h f * hf-ϕ )(n + N ) ×e -δ k n sin(2πν k n -ψ k )½ n≤Dτ .
It can be easily proved that ∀f ∈ [0 . . . F -1], ∀ϕ ∈ Z,

• the support of c g0 (f, ϕ, τ ) is [-L + 1 . . . L + ⌈ N0-2 D ⌉],

• if τ ≤ -L, then c g0 (f, ϕ, τ ), A k (f, ϕ, τ ) and B k (f, ϕ, τ ) are zero, thus c g (f, ϕ, τ ) = 0,

• if τ ≥ L, then A k (f, ϕ, τ ) and B k (f, ϕ, τ ) don't depend on τ . Therefore ∀f ∈ [0 . . . F -1], ∀ϕ ∈ Z, c g (f, ϕ, τ -L + 1) is the impulse response of a causal and stable recursive filter, whose transfer function has a denominator of order 2Q and a numerator of order 2L + 2Q -1 + ⌈ N0-2 D ⌉]. As a particular case, suppose that ∀k ∈ [1 . . . Q], |δ 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 Q(f, ϕ) ≤ Q (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 ∈ [0 . . . F -1], ∀ϕ ∈ Z, there is a limited number Q(f, ϕ) ≤ Q of c g k (f, ϕ, τ ) which contribute to c g (f, ϕ, τ ). We then define Q a = 2max f,ϕ Q(f, ϕ) and

Q b = 2L + Q a -1 + ⌈ N0-2 D ⌉.
Then ∀f ∈ [0 . . . F -1], ∀ϕ ∈ Z, c g (f, ϕ, τ -L + 1) is the impulse response of a causal and stable recursive filter, whose transfer function has a denominator of order Q a and a numerator of order Q b . Considering Remark 1, we conclude that the input/output system described in equation ( 1) is equivalent to the state space representation (3), where P b = K -1.

Proof of Proposition 3:

We consider the state space representation in Definition 1, and we first assume that ∀f ∈ [0 . . . F -1], sequences x(f, t), y(f, t), and z(f, t) belong to l 1 (Z). Then the following DTFTs are well-defined: Y (f, ν) = t∈Z y(f, t)e -2iπνt , X(f, ν) = t∈Z x(f, t)e -2iπνt , B g (f, ϕ, ν) = τ ∈Z b g (f, ϕ, τ )e -2iπντ , A g (f, ν) =

where

C g (f, ϕ, ν) = B g (f, ϕ, ν) A g (f -ϕ, ν) (16)
is the frequency response of a recursive filter. Since this frequency response is twice continuously differentiable, then this filter is stable, which means that its impulse response c g (f, ϕ, τ ) =

  Multichannel high resolution NMF for modelling convolutive mixtures of non-stationary signals in the time-frequency domain Roland Badeau, Mark D. Plumbley Abstract-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 NMF (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. The new model can represent a variety of stationary and non-stationary signals, including autoregressive moving average (ARMA) processes and mixtures of damped sinusoids. 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. Index Terms-Non-stationary signal modelling,

)

  Fig. 1. Time-frequency vs. time domain transformations

  [-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]. We thus have θ = {σ 2 y , σ 2 xs , a s , b ms } s∈[0...S-1],m∈[0...M-1]

and σ 2 0

 2 (t) = ½ {t=0} (where ½ S denotes the indicator function of a set S), then ∀m ∈ [0 . . . M -1], ∀f ∈ [0 . . . F -1]

Fig. 3 .

 3 Fig. 3. Input stereo signal ym(f, t).

  are unobserved: δ m (f, t) = 0 ∀t ∈ [26 . . . 50] (which corresponds to the time range 0.47s-0.91s), and δ m (f, t) = 1 for all other t in [0 . . . T -1].

Fig. 4 .Fig. 5 .

 45 Fig. 4. Stereo signal y ms (f, t) estimated with filters of length 1.

Q

  k=1 g k (n), where Q ∈ N, g 0 (n) is a causal sequence of support [0 . . . N 0 -1] (with N 0 ∈ N), and ∀k ∈ [1 . . . Q], g k (n) = A k e δ k n cos(2πν k n + ψ k )½ n≥0 where A k > 0, δ k < 0, ν k ∈ [0, 1 2 ], ψ k ∈ R. Then ∀f ∈ [0 . . . F -1], equation (2) yields c g (f, ϕ, τ ) = Q k=0 c g k (f, ϕ, τ ) with c g0 (f, ϕ, τ ) = (h f * hf-ϕ * g 0 )(D(τ + L)) and ∀k ∈ [1 . . . Q], c g k (f, ϕ, τ ) = e δ k Dτ (A k (f, ϕ, τ ) cos(2πν k Dτ ) +B k (f, ϕ, τ ) sin(2πν k Dτ ))

Qaτ 1 Ag

 1 =0 a g (f, τ )e -2iπντ . Then applying the DTFT to equation (3) yields Z(f, ν) = (f,ν) X(f, ν) and Y (f, ν) = P b ϕ=-P b B g (f, ϕ, ν)Z(fϕ, ν). Therefore Y (f, ν) = P b ϕ=-P b C g (f, ϕ, ν)X(fϕ, ν),

The probability distributions of processes nm(n) and xs(n) will be defined in the TF domain in section IV-B.

More precisely, compared to the result of Proposition 2, processes zs(f, t) and xs(f, t) as defined in section IV-B are shifted L-1 samples backward, in order to write bms(f, ϕ, τ ) in a causal form. This does not alter the definition of HR-NMF, since equation (8) is unaltered by this time shift, and yms(f, t) is unchanged in equation[START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multiway Data Analysis and Blind Source Separation[END_REF].

1 0 C g (f, ϕ, ν)e +2iπντ dν belongs to l 1 (F). Equations ( 1) and ( 4) are then obtained by applying an inverse DTFT to (15) and (16). Finally, even if x(f, t), y(f, t), and z(f, t) belong to l ∞ (Z) but not to l 1 (Z), equations ( 1) and ( 3) are still well-defined, and the same filter c g (f, ϕ, τ ) ∈ l 1 (F) is still the only stable solution of equation ( 4).


A Jon Altuna

B Bernie Mulgrew

C Roland Badeau

D Vicente Atxa

A FAST ADAPTIVE METHOD FOR SUBSPACE BASED BLIND CHANNEL ESTIMATION

In this paper, a new fast adaptive blind channel estimation method is proposed using the subspace information from the correlation matrix. The algorithm is fully adaptive in the sense that both the subspace information and the optimization which leads to the channel estimation are computed adaptively. It is based on the recently proposed YAST subspace tracker which has been shown to outperform other methods both in terms of speed of convergence and computational complexity. A discussion on the convergence properties of the proposed algorithm is presented. We also propose a hybrid method which makes use of the YAST subspace tracker for initial fast convergence and the subspace information is then updated using the numerically stable OPAST subspace tracker.

INTRODUCTION

The increasing demand for high spectral efficiency in broadband wireless communications where multipath delay spread is a major problem has motivated the design of (semi-) blind channel estimation and equalization techniques which require little or no knowledge of the transmitted sequence. A number of (semi-) blind techniques have been proposed, which rely on some statistical properties or algebraic structure of the received signal [START_REF] Tong | Multichannel Blind Identification: From Subspace to Maximum Likelihood Methods[END_REF]. In particular, several subspace techniques have been proposed for blind and semi-blind channel identification after the pioneering work by Moulines et al. [START_REF] Cardoso | Subspace methods for the blind identification of multichannel FIR filters[END_REF] which has been applied to single-carrier [START_REF] Gorokhov | Semi-blind second order identification of convolutive channels[END_REF] and multicarrier systems [START_REF] Shengli Zhou | Subspacebased (semi-) blind channel estimation for block precoded space-time OFDM[END_REF]. These algorithms generally require the computation of an SVD on the received signal correlation matrix to determine the signal and noise subspace in a first step and then the channel impulse response is estimated as the solution to a quadratic form up to a scalar factor.

Recent developments in subspace tracking have shown that performance matching a detailed SVD can be obtained using subspace projection based algorithms [START_REF] Davila | Efficient, high performance, subspace tracking for time-domain data[END_REF] [START_REF] Badeau | Yet Another Subspace Tracker[END_REF]. These algorithms view the estimation of the signal subspace as an optimization problem and compute the signal subspace as the exact solution over a subspace of reduced dimensions. Essentially, an orthonormal matrix which spans the n-dimensional dominant subspace of the correlation matrix is adaptively estimated where the orthonormality of the estimated projection matrix is guaranteed at each iteration. In particular, the YAST algorithm proposed in [START_REF] Badeau | Yet Another Subspace Tracker[END_REF] clearly outperforms subspace trackers such as OPAST [START_REF] Abed-Meraim | Fast orthonormal PAST algorithm[END_REF] or plane-rotation based FST [START_REF] Rabideau | Fast, rankadaptive subspace tracking and applications[END_REF] with a global complexity O(mn) where m is the dimension of the observed data window. Recently, a promising subspace tracker known as FAPI [START_REF] Badeau | Fast approximated power iteration subspace tracking[END_REF] has also been presented, which outperforms OPAST both in terms of stability and speed of convergence, but YAST has a much faster initial convergence than FAPI.

This paper presents a fast fully-adaptive subspace-based blind channel estimation method. We first propose a method which uses the YAST subspace tracker to adaptively calculate the signal subspace from the correlation matrix and the quadratic form solution is also adaptively calculated using the iterative power method [START_REF] Doukopoulos | Adaptive power techniques for blind channel estimation in CDMA systems[END_REF]. However, we have observed that, although the YAST algorithm converges remarkably at a similar speed of a detailed SVD with an important reduction in computational complexity, the algorithm diverges when the iterations reach the dimension of the signal subspace. Nevertheless, we propose a hybrid approach which combines the YAST subspace tracker for initial fast convergence and stability is guaranteed by updating the subspace information using the numerically stable OPAST subspace tracker. The proposed method is numerically stable and achieves faster convergence than other methods proposed in the literature [START_REF] Doukopoulos | Adaptive power techniques for blind channel estimation in CDMA systems[END_REF] with a global complexity of O(mn).

The paper is organized as follows: in Section 2 the multichannel system model is presented and Section 3 presents the fully adaptive blind channel estimation method based on the YAST subspace tracker. Section 4 discusses convergence properties of the adaptive algorithm and proposes a hybrid method which guarantees stability. Finally, computer simulations are presented in Section 5 which highlight the ability of the proposed hybrid approach to track the impulse response of a frequency-selective propagation channel faster than other methods using different subspace trackers.

SYSTEM MODEL

The system is described as a multichannel system where the discrete time input sequence x(k) is received by K antennas or equivalently is sampled at a rate K-times faster than the source symbol rate. The received signal y (i) (k) at the ith subchannel is formulated as:

y (i) k = H (i) N x k + w (i) k (1)
where

x k = [x(k), x(k -1), • • • , x(k -N -M + 1)
] † is the (N + M ) × 1 vector which contains the input sequence x(k), N represents the dimension of the data window at each subchannel and M + 1 is the length of the subchannel impulse response. On the other hand, y

(i) k = [y (i) (k), y (i) (k -1), • • • , y (i) (k - N + 1)] † represents the N × 1 received signal vector and w (i) k = [w (i) (k), w (i) (k -1), • • • , w (i) (k -N + 1)] † is the noise vec- tor. The subchannel filtering matrix H (i) N is a N × (N + M ) Block-Toeplitz matrix associated with the ith subchannel impulse response h (i) k def = h (i) (k) as: H (i) N =       h (i) 0 • • • h (i) M 0 • • • • • • 0 0 h (i) 0 • • • h (i) M 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . 0 • • • • • • 0 h (i) 0 • • • h (i) M       (2)
Equivalently,

y k = HN x k + w k (3)
where

y k def = [y (0) † k , • • • , y (K-1) † k ] † is the received signal vector and HN = [H (0) † N , • • • , H (K-1) † N ] † is a KN × (N + M ) filtering matrix.

ADAPTIVE SUBSPACE-BASED BLIND CHANNEL ESTIMATION

The identification of the channel coefficients is based on the SVD of the KN ×KN covariance matrix Ry(k) of the received signal.

Provided that the filtering matrix HN meets the full column-rank condition, the observation window N is sufficiently large so that N ≥ M and the subschannels h

(i)
k do not share common zeros [START_REF] Cardoso | Subspace methods for the blind identification of multichannel FIR filters[END_REF], the singular vectors S(k) associated with the M +N highest singular values span the signal subspace and the singular vectors G(k) associated with the last KN -M -N singular values span the noise-subspace:

Ry(k) = HN Rx(k)H H N + Rw(k) = [S(k)G(k)] ΣS(k) 0 0 ΣG(k) S(k) H G(k) H (4)
where Rw(k) is the autocorrelation matrix of the stationary noise process.

Since the signal subspace is a linear space spanned by the filtering matrix HN , according to [START_REF] Cardoso | Subspace methods for the blind identification of multichannel FIR filters[END_REF], the subschannel impulse responses h

(i)
k can be estimated by maximizing the following quadratic form:

q ′ (H) = M +N -1 i=0 | Ŝi(k) H HN | 2 = H H M +N -1 i=0 Ŝi(k) Ŝi(k) H H = H H Q(k)H (5)
H = [H (0)T , . . . , H (K-1)T ] is a 1 × K(M + 1) channel vector and H (i) = [h (i) 0 , . . . , h (i) M ] T . Ŝi(k)
are the estimated signal subspace singular vectors and Ŝi(k) is a transformation matrix such that:

Ŝi(k) = ΥM+1[ Ŝi(k)], (6)
where ) ]. The solution to this maximization problem is the unit-norm singular vector associated with the maximum singular value of matrix Q(k) subject to the unit norm constraint H = 1.

H (i) N = ΥN [H (i
In particular, we propose an adaptive method for the estimation of the signal subspace based on maximizing the following criterion:

τ [ Ŝ(k)] = trace[ Ŝ(k) H Ry(k) Ŝ(k)] (7)
where Ŝ(k) is the estimated m × n orthonormal matrix which spans the signal subspace and m = KN and n = M + N . This maximization is subject to the recursive estimation of the received signal covariance matrix as:

Ry(k) = β Ry(k -1) + y k y H k (8)
The idea introduced in [START_REF] Davila | Efficient, high performance, subspace tracking for time-domain data[END_REF] and [START_REF] Badeau | Yet Another Subspace Tracker[END_REF] is to reduce the search space of the m × n orthonormal vectors Ŝ(k) which requires an overall complexity O(m 2 n) by maximizing the following alternative criterium instead:

τ [ Ŝ(k)] = trace[Rz(k)] = trace[U (k) H S(k) H Ry(k) S(k) U (k)] = trace[U (k) H Rz(k)U (k)] (9)
where

Ŝ(k) = S(k)U (k). The matrix S(k) = [ Ŝ(k -1), V (k)
] is an orthogonal matrix with p = 1 or p = 2 more columns appended to matrix Ŝ(k -1). On the other hand, V (k) = e k /norm(e k ) is the normalized prediction error, where:

e k = x k -Ŝ(k -1) Ŝ(k -1) H x k (10)
The maximization of eq. 9 reduces the search space from m × n to the (n + p) × n orthonormal matrix U (k) which spans the n-dimensional dominant subspace of Rz(k), but it would still require to compute an singular value decomposition with a complexity O(n 3 ). In order to avoid this, we take into account the fact that the orthonormal matrix U (k) also spans an invariant subspace of the inverse Rz(k)

-1 Rz(k) -1 = U (k) H Rz(k) -1 U (k) (11)
Thus, the YAST algorithm computes Rz(k) -1 and its p-dimensional dominant subspace which does not require any singular value decomposition. This optimization leads to a simple recursion for the adaptation rule of the signal subspace matrix [START_REF] Badeau | Yet Another Subspace Tracker[END_REF]:

Ŝ(k) = Ŝ(k -1) -ζ(k) f H k = Ŝ(k -1) -[V (k) + Ŝ(k -1) fk ( Ĩ + ρ(k)) -1 ] f H k (12)
where fk = φ(k) θ(k) H . For the sake of simplicity we study the adaptation rule for p = 1 rule. Then, φ(k) is built from the first m coefficients of the (m + 1) × 1 dominant singular vector of the matrix Rz(k)

-1 , and θ(k) and ρ(k) correspond to the sign and the absolute value of the last coefficient of this singular vector. The adaptation of the signal subspace matrix Ŝ(k) requires 4mn+ O(m + n 2 ) MAC (multiply-accumulate) operations.

On the other hand, the solution to the maximization of the quadratic form expressed in eq. 5 would require the computation of the dominant singular vector of the m1 × m1 matrix Q(k) with a complexity O(m 2 1 ), where m1 = K(M + 1). In this case, a further reduction in complexity can be achieved by tracking the most significant singular vector of matrix Q(k) using the power method [START_REF] Golub | Matrix computations[END_REF]. The power method is an iterative method for computing the dominant singular vector with a complexity in each iteration. The procedure is as follows:

Since the matrix Q(k) is a symmetric positive matrix, the dominant left singular vector Uq(k) of Q(k) can be directly obtained according to:

Uq(k) = Q(k)Uq(k -1) Uq(k) = Uq(k)/ Uq(k) (13)
Finally, the channel vector estimate can be obtained from the dominant left singular vector Uq(k) = [u

(0) q , • • • , u (K(M +1)) q ] us- ing H (i) = [u (i) q , u (K+i) q , • • • , u (KM +i) q ].
In summary, the adaptive blind channel estimation algorithm would require the following steps:

Table I Adaptive Blind Channel Estimation Recursion: for k = 1, 2, . . . Compute Ry(k) = β Ry(k -1) + y k y H k Update Ŝ(k) = Ŝ(k -1) -ζ(k) f H k (eq. 12) Compute Q(k) = M +N -1 i=0 Ŝi(k) Ŝi(k) H Update Uq(k) using Uq(k) = Q(k)Uq(k -1)
and Uq(k) = Uq(k)/ Uq(k) end

CONVERGENCE PROPERTIES AND HYBRID METHOD

In this section we consider the convergence properties of the YAST subspace tracker in the context of blind channel estimation. In particular, it is interesting to note that after k = M + N iterations, the signal subspace is somewhat estimated as the M + N -dimensional subspace of [ Ŝ(k -1), x k , x k+1 , . . . , x k+M +N -1 ] which maximizes the Rayleigh quotient. However, the family x k , x k+1 , . . . , x k+M +N -1 exactly spans the true signal subspace and therefore, after k = M + N iterations the estimated subspace matches the true signal subspace. This explanation only holds for noiseless signals.

On the other hand, we have observed that the YAST algorithm diverges immediately after convergence. This is caused by the loss of orthonormality of the updated matrix Ŝ(k) due to the fact that the prediction error e k = x k -Ŝ(k -1) Ŝ(k -1) H x k becomes neglegible at iteration k = M + N . Indeed, when the algorithm reaches convergence at k = M + N we obtain an orthonormal matrix Ŝ(k) which matches the exact SVD solution for the same iteration. In practice, it is observed that the normalized prediction error vector

V (k) = e k /norm(e k ) in S(k) = [ Ŝ(k -1), V (k)] is not exactly orthogonal to Ŝ(k -1).
In most cases, especially when the observation window N is small, the channel estimates are not able to converge to a reasonable value and we thus propose to use the YAST subspace tracker as an initialization with fast convergence for k = M +N iterations and after that, the numerically stable OPAST subspace tracker [START_REF] Abed-Meraim | Fast orthonormal PAST algorithm[END_REF] can be used for the adaptation of the signal subspace matrix.

The recursion used by the OPAST algorithm for the adaptation of the signal subspace is:

Ŝ(k) = Ŝ(k -1) + pk q H k (14) where pk = τ k Ŝ(k -1)q k + (1 + τ k q k
2 )p k . On the other hand, we define the following terms as:

p k = γ k e k , (15)
γ k = 1/(1 + x H k Ŝ(k -1)q k ), (16)
q k = 1/β(Z(k -1) Ŝ(k -1) H x k ) (17)
Z(k) = (1/β)Z(k -1) -γ k q k q H k (18)
and

τ k = 1 q k 2 1 1 + p k 2 q k 2 -1 . (19)
This adaptation rule requires 4mn+O(n 2 ) MAC (multiply-accumulate) operations. The pseudo-code of the proposed hybrid method is as follows

(k) = M +N -1 i=0 Ŝi(k) Ŝi(k) H Update Uq(k) using Uq(k) = Q(k)Uq(k -1) and Uq(k) = Uq(k)/ Uq(k) end
The computational complexity of the different adaptive subspace based channel estimation techniques is summarized in Table III:

Table III subspace estimation of estimation Uq(k) Detailed SVD O[(KN ) 3 ] [K(M + 1)] 2 YAST 4KN (M + N ) [K(M + 1)] 2 +O[KN + (M + N ) 2 ] OPAST 4KN (M + N ) [K(M + 1)] 2 +O[(M + N ) 2 ] Hybrid YAST/OPAST [K(M + 1)] 2

SIMULATION RESULTS

Computer simulations have been conducted in order to compare the performance of the proposed fully adaptive method with the blind subspace technique using exact SVD and an adaptive method using the OPAST subspace tracker Figure 1 shows the convergence of the NMSE averaged over 50 Monte-Carlo runs, using a fixed observation window size N = 20 and SNR=25 dB . It is observed that the adaptive method using the YAST algorithm the signal subspace estimation and the iterative power method for the quadratic form solution is able to converge to the solution provided by the exact SVD-based method (denoted BL Subs). However, it can be seen that the algorithm diverges immediately after convergence at k = M + N . On the other hand, we have also tested the performance of the adaptive technique using the OPAST subspace tracker and it is apparent that the convergence is slow compared to the YAST subspace tracker. We finally tested the proposed hybrid method which switches from the YAST adaptation rule to OPAST at k = M +N and it is shown that this combination enables global convergence of the algorithm while increasing the speed of convergence.

Another experiment was also conducted computing the NMSE over 100 Monte-Carlo runs using a fixed window size N = 20. It is observed in Figure 2 that the proposed hybrid adaptive technique achieves almost identical performance as the exact SVDbased technique at iterations k = 22 and k = 150, whereas YAST diverges after iteration k = 22 and the OPAST based adaptive method is not even able to converge at k = 150.

CONCLUSIONS

In this paper we propose a new method for adaptive subspacebased blind channel estimation using the YAST subspace tracker for initial signal subspace estimation and the solution of the quadratic form which leads to the blind channel estimate is also adaptively computed. It has been discussed and shown by means of computer simulations that although the YAST algorithm might diverge when the iterations reach the dimension of the signal subspace, it can still be used for initial fast convergence followed by other subspace trackers such as OPAST while keeping a global complexity of O(mn). The resulting hybrid method achieves a remarkable convergence speed and numerical stability.

Fig. 1 .

 1 Fig. 1. Convergence of channel estimation using an observation window size N = 20 and SNR=25 dB.

Fig. 2 .

 2 Fig. 2. NMSE versus SNR for different at iterations k = 22 (solid) and k = 150 (dashed).


Roland Badeau
email: roland.badeau@telecom-paristech.fr

GAUSSIAN MODELING OF MIXTURES OF NON-STATIONARY SIGNALS IN THE TIME-FREQUENCY DOMAIN (HR-NMF)

Keywords: Nonnegative Matrix Factorization, High Resolution methods, Expectation-Maximization algorithm, source separation, audio inpainting

Nonnegative Matrix Factorization (NMF) is a powerful tool for decomposing mixtures of non-stationary signals in the Time-Frequency (TF) domain. However, unlike the High Resolution (HR) methods dedicated to mixtures of exponentials, its spectral resolution is limited by that of the underlying TF representation. In this paper, we propose 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. This model is estimated with a recursive implementation of the EM algorithm, that is successfully applied to source separation and audio inpainting.

INTRODUCTION

Some very powerful tools have been recently introduced for modeling mixtures of non-stationary signal components in the TF domain. Among them, NMF [START_REF] Virtanen | Bayesian extensions to nonnegative matrix factorisation for audio signal modelling[END_REF] and Probabilistic Latent Component Analysis (PLCA) [START_REF] Smaragdis | Blind Speech Separation[END_REF] compute an approximate factorization of a magnitude or power TF representation, such as the spectrogram.

Since phases are generally discarded in these models, reconstructing the phase field requires employing ad-hoc methods [START_REF] Griffin | Signal reconstruction from short-time Fourier transform magnitude[END_REF]. To the best of our knowledge, apart the complex NMF which was designed in a deterministic framework [START_REF] Le Roux | Complex NMF under spectrogram consistency constraints[END_REF], the only probabilistic model that takes the phase field into account is the Itakura-Saito (IS)-NMF [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF]. Separating the signal components is then proven equivalent to Wiener filtering. The spectral resolution of IS-NMF is thus limited by that of the TF representation (sinusoids in the same frequency band cannot be properly separated).

In other respects, IS-NMF assumes that all TF coefficients are independent, which is not suitable for modeling sinusoids or transients for instance. Markov models have thus been proposed for taking the local dependencies between contiguous TF coefficients of a magnitude or power TF representation into account [START_REF] Dikmen | Gamma Markov random fields for audio source modeling[END_REF][START_REF] Mysore | Non-negative hidden Markov modeling of audio with application to source separation[END_REF].

In this paper, we introduce a unified model called HR-NMF, which natively takes both phases and local correlations in each frequency band into account. This approach avoids using a phase reconstruction algorithm, and we show that it overcomes the spectral resolution of the TF representation. It can be used with both complex-valued and real-valued TF representations (like the shorttime Fourier transform or cosine modulated filterbanks). This paper is organized as follows: HR-NMF is introduced in section 2, then our recursive implementation of the EM algorithm for estimating this model is presented in section 3. Section 4 is devoted to experimental results, and conclusions are drawn in section 5. The following notation will be used throughout the paper:

• x: scalar (normal letter),

• v = [v1; . . . ; vK ]: column vector (bold lower case letter),

• S: matrix (bold upper case letter),

• S (t,p) : (t, p)-th entry of matrix S (indexed upper case letter), • S * (resp. S H ): conjugate (resp. conjugate transpose) of S, • diag(.): (block)-diagonal matrix whose diagonal blocks are (.), • 1 (resp. 0): vector whose entries are all equal to 1 (resp. 0), • E[.] /x : conditional expectation of (.) given the observation x, • N (µ, R): real or circular complex normal distribution of mean µ and covariance matrix R, • for a given vector v of dimension D, and any subvector v of dimension D ≤ D (whose entries are a subset of those of v),

J v v denotes the D × D selection matrix such that v = J v v H v.

TIME-FREQUENCY MIXTURE MODEL

The mixture model x(f, t) is defined for all frequencies 1 ≤ f ≤ F and times 1 ≤ t ≤ T as the sum of K latent components c k (f, t) plus a white noise n(f, t) ∼ N (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) (and P (k, f ) ∈ N is such that a(P (k, f ), k, f ) = 0), • b k (f, t) ∼ N (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.

Since N denotes either the real or the circular complex normal distribution, model (1) is either real or complex-valued. Moreover, we assume that ∀k ∈ {1 . . . K}, ∀f ∈ {1 . . . F }, ∀t ≤ 0, c k (f, t) are independent and identically distributed random variables: c k (f, t) ∼ N (0, 1), and x(f, t) are unobserved. The parameters to be estimated are σ 2 , a(p, k, f ), w(k, f ), and h(k, t).

This time-frequency model generalizes some very popular models, widely used in various signal processing communities:

• If σ 2 = 0 and ∀k, f, P (k, f ) = 0, (1) becomes x(f, t) = K k=1 b k (f, t), thus x(f, t) ∼ N (0, V f t ), where V is defined by the NMF V = W H with W f k = w(k, f ) and H kt = h(k, t).
Maximum likelihood estimation of W and H is then equivalent to the minimization of the IS-divergence between the matrix model V and the spectrogram V (where

V f t = |x(f, t)| 2 )
, that is why 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].

• For given values of k and f , if ∀t, h(k, t) = 1, then c k (f, t) is an autoregressive process of order P (k, f ). • For given values of k and f , if P (k, f )≥1 and ∀t≥P (k, f )+1, h(k, t) = 0, then c k (f, t) can be written in the form c k (f, t) = P (k,f ) p=1
αp z t p , where z1 . . . z P (k,f ) are the roots of the polynomial z P (k,f ) -P (k,f ) p=1 a(p, k, f )z P (k,f )-p . This corresponds to the Exponential Sinusoidal Model (ESM) 1 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].

For these reasons, we refer to model (1) as HR-NMF.

EXPECTATION-MAXIMIZATION (EM) ALGORITHM

In order to estimate the model parameters, we apply the EM algorithm2 to the observed data x and the latent components c1 . . . cK . In order to handle the case of missing data, we define δ(f, t) = 1 if x(f, t) is observed, and δ(f, t) = 0 else.

Maximization Step (M-step)

The conditional expectation (given the observations) of the loglikelihood of the complete data is

Q = E /x [ln(p(c1 . . . cK , x))] = E /x [ln(p(x/c1 . . . cK ))] + K k=1 E /x [ln(p(c k ))] . We can thus write 3 Q c = Q0 + K k=1 Q k where Q0 = -F f =1 T t=1 δ(f, t) ln(σ 2 ) + e(f, t)/σ 2 , Q k = - 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) , e(f, t) = δ(f, t) E /x x(f, t) - K k=1 c k (f, t) 2 , (3)
a(k, f ) = [1; -a(1, k, f ); . . . ; -a(P (k, f ), k, f )], (4)
and ∀k, f , for all 0 ≤ p1, p2 ≤ P (k, f ),

S (p 1 ,p 2 ) (k, f, t) = E /x [c k (f, t -p1) * c k (f, t -p2)]. ( 5
)
Maximizing Q is thus equivalent to independently maximizing Q0 with respect to (w.r.t.) σ 2 and each Q k w.r.t. h(k, t), w(k, f ) and a(k, f ). Since the maximization of Q k does not admit a closed form solution, we propose to recursively maximize Q k w.r.t. (w(k, f ), a(k, f )) and w.r.t. h(k, t). We cannot provide here the full mathematical derivation of the M-step because of the page limit, but its pseudo-code is summarized in Table 1 (according to the notation introduced in section 1, J

a(k,f ) 1 = [1; 0; . . . ; 0]). Its complexity is O(F T KP 2 ), where P = max k,f P (k, f ).
Inputs after the E-step: δ(f, t), e(f, t), S(k, f, t), h(k, t)

σ 2 = F f =1 T t=1 e(f, t) F f =1 T t=1 δ(f, t) For k = 1 to K 4 ,
Repeat (as many times as wanted) 5 :

For f = 1 to F 4 , Σ(k, f ) = 1 T T t=1 S(k, f, t) h(k, t) α(k, f ) = Σ(k, f ) -1 J a(k,f ) 1 w(k, f ) = 1/(J a(k,f ) 1 H α(k, f )) a(k, f ) = w(k, f ) α(k, f ) End for f ; For t = 1 to T 4 , h(k, t) = 1 F F f =1 a(k, f ) H S(k, f, t)a(k, f ) w(k, f ) End for t; Normalization of the NMF: H k = max t (h(k, t)), w(k, f ) = H k w(k, f ), h(k, t) = h(k, t)/H k End repeat; End for k; Outputs: σ 2 , a(k, f ), w(k, f ), h(k, t).
Table 1: Pseudo-code of the M-step

Expectation Step (E-step)

The purpose of the E-step is to determine the a posteriori distribution of the latent components c k (f, t) given the observations x(f, t) (and more precisely, e(f, t) and S(k, f, t)). Since these random variables are mutually independent for different values of f , the Estep can process each f separately. As our recursive implementation of the E-step is inspired from Kalman filtering theory [START_REF] Bierman | Factorization methods for discrete sequential estimation[END_REF], we first introduce the Kalman representation of the HR-NMF model:

γ(f, t) = A(f ) γ(f, t -1) + b ′ (f, t), (6)
x(f, t) = u(f ) H γ(f, t) + n ′ (f, t), (7)
where

• K(f ) denotes the set {k ∈ {1 . . . K}/P (k, f ) ≥ 1}; • ∀k ∈ K(f ), c(k, f, t) = [c k (f, t); . . . ; c k (f, t-P (k, f )+1)]; • the state vector γ(f, t) contains c(k, f, t) for all k ∈ K(f ); • the state transition matrix is A(f ) = diag({A(k, f )} k∈K(f ) ), where 6 A(k, f ) = a(P (k, f ), k, f ) if P (k, f ) = 1, otherwise A(k, f )= a(1, k, f ) . . . a(P (k)-1, k, f ) a(P (k, f ), k, f ) diag(1) 0 • ∀f, t, vector c(f, t) contains c k (f, t) for all k ∈ K(f ), • the white process noise is b ′ (f, t) = J γ(f,t) c(f,t) b(f, t), where no- tation J γ(f,t) c(f,t) was introduced in section 1, and b(f, t) contains b k (f, t) ∀k ∈ K(f ); thus b(f, t) ∼ N (0, diag(v(f, t))), where v(f, t) contains v k (f, t) ∀k ∈ K(f ); • the observation matrix is u(f ) H , where u(f ) = J γ(f,t) c(f,t) 1; • the white observation noise is n ′ (f, t) = n(f, t) + k/P (k,f )=0 b k (f, t) ∼ N (0, σ 2 (f, t)), where σ 2 (f, t) = σ 2 + k/P (k,f )=0 v k (f, t). (8)
We cannot provide the full mathematical derivation of the Estep in this paper because of the page limit, but its pseudo-code is summarized in Table 2. Its overall computational complexity is O(F T K 3 P 3 ). In Table 2, we have used the following notation:

• ∀f, t, vector d(f, t) contains c k (f, t -P (k, f )) ∀k ∈ K(f ), • ∀f, t, vector c ′ (f, t) contains c k (f, t) for all k / ∈ K(f ), • ∀f, t, vector v ′ (f, t) contains v k (f, t) for all k / ∈ K(f ), • ∀f, t, ∀k∈K(f ), c(k, f, t)=[c k (f, t); . . . ; c k (f, t -P (k, f ))], • ∀f, t, γ(f, t) contains c(k, f, t) for all k ∈ K(f ),
• ∀f, t, and for any random vector v, v f,t is the conditional expectation of v given {x(f, 1) . .

. x(f, t)}. Besides, R f,t v is the conditional expectation of ( v f,t ) ( v f,t ) H given {x(f, 1) . . . x(f, t)}, where v f,t = v -v f,t . Similarly, for any vectors v1 and v2, R f,t v 1 ,v 2 is the conditional expectation of ( v f,t 1 ) ( v f,t
2 ) H given {x(f, 1) . . . x(f, t)}. Other letters denote temporary variables used in the computations.

APPLICATIONS

In order to illustrate the capabilities of HR-NMF, we consider two examples of straightforward applications. First, noticing that the E-step estimates c k (f, t), source separation will be addressed in section 4.1. Moreover, since c k (f, t) is estimated even at timefrequencies where the observation x(f, t) is missing, audio inpainting will be addressed in section 4.2. The following experiments deal with a real piano sound, composed of a C4 tone played alone at t = 0 ms, and a C3 tone played at t = 680 ms while the C4 tone is maintained. The sampling frequency is 8600 Hz, and x(f, t) is obtained by computing the STFT of the input signal with dimensions F = 400 and T = 60, using 90 ms-long Hann windows with 75% overlap (the corresponding spectrogram is plotted in Figure 1). HR-NMF is then computed by running 5 iterations of the EM algorithm with P (k, f ) = 1, after initialization with IS-NMF.

Source separation

In this first experiment, the whole STFT x(f, t) is observed, and we aim at separating K = 2 components c k (f, t) in the frequency band f which corresponds to the first harmonic of C4 and to the second harmonic of C3 (around 270 Hz). These two sinusoidal components (whose real parts are represented as red solid lines in Figure 2) have very close frequencies, which makes them hardly separable. As expected, IS-NMF, which involves Wiener filtering, is not able to properly separate the components when they overlap, from t = 680 ms to 1.36 s (the estimated components are represented as black dash-dotted lines). As a comparison, the components estimated by HR-NMF (blue dashed lines) better fit the ground truth.

Inputs after the M-step: 4 For t = 1 to T (forward pass): 5 Predict:

x(f, t), δ(f, t), σ 2 , A(f ), v k (f, t) Initialization: ∀f , γ f,0 (f, 0) = 0, R f,0 γ(f,0) = diag(1) ∀f, t, σ 2 (f, t) = σ 2 + k/P (k,f )=0 v k (f, t) For f = 1 to F ,
R f,t-1 γ(f,t),d(f,t) = A(f )R f,t-1 γ(f,t-1) J γ(f,t-1) d(f,t) R f,t-1 d(f,t) = J γ(f,t-1) d(f,t) H R f,t-1 γ(f,t-1) J γ(f,t-1) d(f,t) R f,t-1 b ′ (f,t) = J γ(f,t) c(f,t) diag(v(f, t)) J γ(f,t) c(f,t) H R f,t-1 γ(f,t) = A(f )R f,t-1 γ(f,t-1) A(f ) H + R f,t-1 b ′ (f,t) d f,t-1 (f, t) = J γ(f,t-1) d(f,t) H γ f,t-1 (f, t -1) γ f,t-1 (f, t) = A(f ) γ f,t-1 (f, t -1) Φ f,t-1 γ(f,t),d(f,t) = R f,t-1 γ(f,t) -1 R f,t-1 γ(f,t),d(f,t) Ψ f,t-1 d(f,t) = R f,t-1 d(f,t) -R f,t-1 d(f,t),γ(f,t) Φ f,t-1 γ(f,t),d(f,t) φ f,t-1 d(f,t) = d f,t-1 (f, t) -Φ f,t-1 γ(f,t),d(f,t) H γ f,t-1 (f, t) Update: µ(f, t) = R f,t-1 γ(f,t) u(f ) ε(f, t) = σ 2 (f, t) + u(f ) H µ(f, t) λ(f, t) = δ(f,t) ε(f,t) µ(f, t) R f,t γ(f,t) = R f,t-1 γ(f,t) -λ(f, t) µ(f, t) H ǫ f,t (f, t) = x(f, t) -u(f ) H γ f,t-1 (f, t) γ f,t (f, t) = γ f,t-1 (f, t) + λ(f, t) ǫ f,t (f, t) End for t;
For t = T down to 1 (backward pass): 5 Wiener filtering:

ǫ f,T (f, t) = x(f, t) -u(f ) H γ f,T (f, t) c ′f,T (f, t) = δ(f,t) σ 2 (f,t) v ′ (f, t) ǫ f,T (f, t) e ′ (f, t) = |ǫ f,T (f,t)| 2 +u(f ) H R f,T γ(f,t) u(f ) σ 2 (f,t) -1 e(f, t) = δ(f, t) σ 2 σ 2 σ 2 (f,t) e ′ (f, t) + 1 ∀k / ∈ K(f ), S(k, f, t) = v k (f, t) + δ(f, t) v k (f,t) 2 σ 2 (f,t) e ′ (f, t) Smoothing: R f,T γ(f,t),d(f,t) = R f,T γ(f,t) Φ f,t-1 γ(f,t),d(f,t) R f,T d(f,t) = Ψ f,t-1 d(f,t) + R f,T d(f,t),γ(f,t) Φ f,t-1 γ(f,t),d(f,t) R f,T γ(f,t) = J γ(f,t) γ(f,t) R f,T γ(f,t) J γ(f,t) γ(f,t) H + J γ(f,t) d(f,t) R f,T d(f,t) J γ(f,t) d(f,t) H +J γ(f,t) γ(f,t) R f,T γ(f,t),d(f,t) J γ(f,t) d(f,t) H +J γ(f,t) d(f,t) R f,T d(f,t),γ(f,t) J γ(f,t) γ(f,t) H R f,T γ(f,t-1) = J γ(f,t) γ(f,t-1) H R f,T γ(f,t) J γ(f,t) γ(f,t-1) d f,T (f, t) = φ f,t-1 d(f,t) + Φ f,t-1 γ(f,t),d(f,t) H γ f,T (f, t) γ f,T (f, t) = J γ(f,t) γ(f,t) γ f,T (f, t) + J γ(f,t) d(f,t) d f,T (f, t) γ f,T (f, t -1) = J γ(f,t) γ(f,t-1) H γ f,T (f, t) S(f, t) = (R f,T γ(f,t) + γ f,T (f, t) γ f,T (f, t) H ) * ∀k ∈ K(f ), S(k, f, t) = J γ(f,t) c(k,f,t) H S(f, t)J γ(f,t) c(k,f,t)
End for t; End for f ; Outputs: e(f, t), S(k, f, t), and c f,T k (f, t) if wanted. 3. Of course the C3 tone, which was unobserved, could not be recovered, but the C4 tone is correctly estimated. Moreover, the noise in the unobserved part has been removed. Actually we observed that a listening test does not permit to perceive any artifact in the signal synthesized from c1(f, t) by a standard overlap-add technique. As a comparison, note that IS-NMF is not capable of audio inpainting, because it does not take the correlations between contiguous TF coefficients into account (the missing coefficients are estimated as zeros).

CONCLUSIONS

In this paper, we introduced a new method for modeling mixtures of non-stationary signals in the time-frequency domain, which was successfully applied to source separation and audio inpainting. Compared to standard IS-NMF, the proposed approach natively takes both phases and local correlations in each frequency band into account. We showed that it achieves high resolution, which means that two sinusoids of different frequencies can be properly sepa-

  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.

Figure 1 :Figure 2 :

 12 Figure 1: Spectrogram of the input piano sound

Figure 3 :

 3 Figure 3: Recovery of the full C4 piano tone

Table 2 :

 2 Pseudo-code of the E-step

4.2. Audio inpainting

In this second experiment, the second part of the STFT (from t = 680 ms to 1.36 s) is unobserved, and in the first part (from

Actually HR-NMF also encompasses the more general Polynomial Amplitude Complex Exponentials (PACE) model introduced in[START_REF] Badeau | High resolution spectral analysis of mixtures of complex exponentials modulated by polynomials[END_REF].

[START_REF] Smaragdis | Blind Speech Separation[END_REF] The EM algorithm is an iterative method which alternates two steps called Expectation and Maximization. It is proved to increase the likelihood of the observed data at each iteration. In

section 3, we assume that w(k, f ), h(k, t), σ 2 and a(P (k, f ), k, f ) are non-zero. However, these parameters might become zero after some iterations, depending on the input data. Such singular cases should be addressed in a rigorous implementation of EM.3 c = denotes equality up to additive and multiplicative constants which do not depend on the model parameters to be estimated.

This loop can be processed in parallel.

This loop has to be processed sequentially.

Note that since ∀k ∈ K(f ), a(P (k, f ), k, f ) = 0, A(f ) is invertible.

Note that contrary to standard high resolution methods, the proposed approach is able to handle mixtures of amplitude-modulated sinusoids starting at different times; it also performs the clustering of these sinusoids into several components according to their temporal dynamics.