TITLE OF INVENTION

METHOD FOR TRANSFORMING INPUT SIGNAL

TECHNICAL FIELD

[0001] This invention relates generally to signal processing, and more particularly to transforming an input signal to an output signal using a dynamic model, where the signal is an audio (speech) signal.

BACKGROUND ART

[0002] A common framework for modeling dynamics in non-stationary signals is a hidden Markov model (HMM) with temporal dynamics. The HMM is the de facto standard for speech recognition. A discrete-time HMM models a sequence of N observed (acquired) random variables

def def

{½ } ^{= χ} ί Ν ~ { ^χ 1> ^χ 2>" ·> ^χ Ν}> i- ^e -' ^s ig ^na l samples, by conditioning probability distributions on the sequence of unobserved random state variables {h _n } . Two constraints are typically defined on the HMM.

[0003] First, the state variables have first-order Markov dynamics. This means that p(h _n | h\. _n _ ) = p(h _n \ h _n _i) , where the p(h _n \ h _n _\) are known as transition probabilities. The transition probabilities are usually constrained to be time-invariant.

[0004] Second, each sample x _n , given the corresponding state h _n , is independent of all other hidden states h _n > , n ≠ n , so that » where the p(x _n \ h _n ) are known as observation probabilities. In many speech applications, the states h _n are discrete, and observations x _n are F -dimensional vector-valued continuous acoustic features,

def def

xn ⁼ { ^■ -/,(«)} ⁼ { ^χ η ·> ^χ 2η ·>· · -> ^x Fn}>

where the parentheses indicate that n is not iterated. Typical frequency features are short-time log power spectra, where /indicates a frequency bin.

def

[0005] Defining initial probabilities ^(/¾ | /¾) = p(h ), the joint distribution of the random variables of the HMM is

N

/>({·½}>¾}) = \ ^h n )P( ^h n I Vl)- (!)

[0006] Linear Dynamical Systems

[0007] A related model is a linear dynamical system used in Kalman filters. The linear dynamical system is characterized by states and observations that are continuous, vector-valued, and jointly Gaussian distributed

h _n = Ah _n _ ₁ + e _n , (2) ^v n = ^Bh n ^{+ v} n > ( ³ )

K K

where h _n €Ξ R (or h _n E C ) is the state at time n, K the dimension of the state space, A is a state transition matrix, ε _η is additive Gaussian transition noise, F F

V„ER (or V _¾ GC ) is the observation at time n, F is the dimension of the observation (or feature) space, B is an observation matrix, v _n is additive

Gaussian noise, and R is real.

[0008] Non-Negative Matrix Factorization

[0009] In the context of audio signal processing, the signal is typically processed using a sliding window and a feature vector representation that is often a magnitude or power spectrum of the audio signal. The features are nonnegative. In order to discover repeating patterns in the signal in an unsupervised way, nonnegative matrix factorization (NMF) is extensively used.

[0010] For a nonnegative matrix V of dimensions F x N , a rank-reduced approximation is

V « WH ,

where W and H are nonnegative matrices of dimensions F x K and K N , respectively. The a roximation is typically obtained from a minimization

where d (x \ y) is a positive function scalar cost function with a unique minimum at x = y . [0011] Itakura-Saito Nonnegative Matrix Factorization (IS-NMF)

[0012] For the audio signal, where the matrix Vis the power spectrogram of a complex-valued short-time Fourier transform (STFT) matrix X, conventional methods have used the Itakura-Saito distance, which measures the difference between the actual and approximated spectrum, as the cost function, because the cost function implies a latent model of superimposed zero-mean Gaussian components that is relevant for audio signals. More precisely, let Xj^ be the complex -valued STFT coefficient at frame n and frequency/ , and x^ = _η , k where cfkn ^{' N} c ( > ^w fr ^h kn )'

[0013] Then,

= D _IS (\ X \ ² \ WH) ₊ cst, (5) where =\ x ^ \ .

[0014] The model can also be expressed as k [0015] It is equivalent to assume that | x |^ is exponentially distributed with parameter β^ _η and uniform phase

Exponential^ fk ^h kn ) > ( ⁶ )

ZJC^ : Unifor i-π ,+π) . (7)

[0016] Smooth IS-NMF

[0017] In smooth variants of IS-NMF, an inverse-gamma or gamma random walk is assumed for independent rows of H. More precisely, the following model has been considered: where kn is nonnegative multiplicative innovation random variable with mode 1, such as f _te : G(o,a - l) _or

£ _kn IG(a,a + \)

where by convention gamma and inverse-gamma are σ(χ \ α, β) =—— x ^a~l exp- βχ , and (8)

T(a) [0018] Models Combining HMMs and NMF

[0019] If HMMs and NMF are combined, then the restriction that only one discrete state can be active at a time is mhereted from the HMMs. This means that multiple model are required for multiple source, leading to potential issues to computational tractability.

[0020] U.S. 7,047,047 describes denoising a speech signal using an estimate of a noise-reduced feature vector and a model of an acoustic environment. The model is based on a non-linear function that describes a relationship between the input feature vector, a clean feature vector, a noise feature vector and a phase relationship indicative of mixing of the clean feature vector and the noise feature vector.

[0021] U.S. 8,015,003 describes denoising a mixed signal, e.g., speech and noise, using a NMF constrained by a denoising model. The denoising model includes training basis matrices of a training acoustic signal and a training noise signal, and statistics of weights of the training basis matrices. A product of the weights of the basis matrix of the acoustic signal and the training basis matrices of the training acoustic signal and the training noise signal is used to reconstruct the acoustic signal.

[0022] In general, the prior art methods that focus on slow-changing noise, are inadequate for fast-changing nonstationary noise, such as experienced by using a mobile telephone in a noisy environment. [0023] Although HMMs can handle speech dynamics, HMMs often lead to combinatorial issues due to the discrete state space, which is computationally complex, especially for mixed signals from several sources. In conventional HMM approaches it is also not straightforward to handle gain adaptation.

[0024] NMF solves both the computational and gain adaptation issues. However, NMF does not handle dynamic signals. Smooth IS-NMF attempts to handle dynamics. However, the independence assumption of the rows of H is not realistic, as the activation of a spectral pattern at frame n is likely to be correlated with the activation of other patterns at a previous frame n-\.

[0025] It is an object of the invention to solve inherent problems associated with signal and data processing using HMMs and NMF frameworks.

SUMMARY OF INVENTION

[0026] It is an object of the invention to transform an input signal to an output signal when the input signal is a non-stationary signal, and more

specifically a mixture of signals. Therefore, the embodiments of the invention provide a non-negative linear dynamical system model for processing the input signal, particularly a speech signal that is mixed with noise. In the context of speech separation and speech denoising, our model adapts to signal dynamics on-line, and achieves better performance than conventional methods.

[0027] Conventional models for signal dynamics frequently use hidden Markov models (HMMs) or non-negative matrix factorization (NMF). [0028] HMMs lead to combinatorial problems due to the discrete state space, are computationally complex, especially for mixed signals from several sources. In conventional HMM approaches it is also not straightforward to handle gain adaptation.

[0029] NMF solves both the computational complexity and gain adaptation problems. However, NMF does not take advantage of past observations of a signal to model future observations of that signal. For signals with predictable dynamics, this is likely to be suboptimal.

[0030] Our model has advantages of both the HMMs and the NMF. The model is characterized by a continuous non-negative state space. Gain adaptation is automatically handled during inference. The complexity of the inference is linear in the number of signal sources, and dynamics are modeled via a linear transition matrix.

[0031] Specifically the input signal, in the form of a sequence of feature vectors, is transformed to the output signal by first storing parameters of a model of the input signal in a memory.

[0032] Using the vectors and the parameters, a sequence of vectors of hidden variables is inferred. There is at least one vector h _n of hidden variables hi _n for each feature vector X _n , and each hidden variable is nonnegative.

[0033] The output signal is generated using the feature vectors, the vectors of hidden variables, and the parameters. Each feature vector X _n is dependent on at least one of the hidden variables h _n for the same n. The hidden variables are related according to fy/i where j and / are summation

indices. The parameters include non-negative weights i, and ε _{t n} are independent non-negative random variables.

BRIEF DESCRIPTION OF DRAWINGS

[0034] Fig. 1 is a flow diagram for transforming an input signal to an output signal;

[0035] Fig. 2 is a flow diagram of a method for determining parameters of a dynamic model according to embodiment of the invention; and

[0036] Fig. 3 is a flow diagram of a method for enhancing a speech signal using the dynamic model according to embodiments of the invention.

DESCRIPTION OF EMBODIMENTS

[0037] Introduction

[0038] The embodiments of our provide a model for transforming and processing dynamic (non-stationary) signal and data that has advantages of HMMs and NMF based models. [0039] The model is characterized by a continuous non-negative state space. Gain adaptation is automatically handled on-line during inference. Dynamics of the signal are modeled using a linear transition matrix A. The model is a non-negative linear dynamical system with multiplicative non-negative innovation random variables 8 _n . The signal can be a non-stationary linear signal, such as an audio or speech signal, or a multi-dimensional signal. The signal can be expressed in the digital domain as data. The innovation random variable is described in greater detail below.

[0040] The embodiments also provide applications for using the model.

Specifically, the model can be used to process an audio signal acquired from several sources, e.g., the signal is a mixture of speech and noise (or other acoustic interference) and the model is used to enhance the signal by, e.g., reducing noise. When we say "mixed," we mean that the speech and noise are acquired by a single sensor (microphone).

[0041] However, it is understood that the model can also be used for other non-stationary signals and data that have characteristicsthat vary over time, such as economic or financial data, network data and signals, or signals, medical signals, or other signals acquired from natural phenomena. The parameters include non-negative weights <¾/, and ε ι _>η are independent non-negative random variables, the distributions of which also have parameters. The indices i, j, I, and n are described below.

[0042] General Method [0043] As shown in Fig. 1, parameters 101 of a model of an input signal 102 are stored in a memory 103.

[0044] The input signal is received as a feature vectors x _n 104 of salient characteristics of the signal. The features are of course application and signal specific. For example, if the signal is an audio signal, the features can be log power spectra. It is understood that the different type of of features that can be used is essentially unlimited for many types of different signals and data that can be processed by the method according to the invention.

[0045] The method infers 110 a sequence of vectors of hidden variables 111. The inference is based on the feature vector 104, the parameters, a hidden variable relationship 130, and a relationship 140 of observations to hidden variables. There is at least one vector h _n of hidden variables for each feature vector X _n . Each hidden variable is nonnegative.

[0046] An output signal 122 corresponding to the input signal is generated 120 to form the feature vectors, the vectors of hidden variables, and the

parameters.

[0047] General Method Details

[0048] In our method, each feature vector X _n is dependent on at least one of the hidden variables hi _>n for the same n. The hidden variables are related according to a hidden variable relationship _n = ^ j / £/ _n h j _η _ 130, where / and /

JJ are summation indices. The stored parameters include non-negative weights Ci , and ε n are independent non-negative random variables. This formulation enables the model to represent statistical dependency over time in a structured way, so that the hidden variables for the current frame, n, are dependent on those of the previous frame, n-1 with a distribution that is determined by the combination of

Ci , and the parameters of the distribution of the weights ε ι _>η . The weight ε ι _>η , for example, may be Gamma random variables with shape parameter cc and inverse scale parameter β .

[0049] In one embodiment, C ₍ = d(i, l )ai j , where d j are

non-negative scalars, so that δ is a

Kronecker delta. In this case, if the weights ε ι _>η , are Gamma random variables with shape parameter a and inverse scale parameter β , then the conditional distribution of h _n given he

hidden states vector, is pi i b ^a _i - where Gamma(x \ a,b) = x ^a e ^x is the gamma distribution for random

Γ»

00

—t z— 1

variable * with shape a, inverse scale b, and T(z) = Je t dt is the gamma

0 function. This embodiment is designed to conform to the simplicity of the basic structure of a conventional linear dynamical system, but differs from prior art by the non-negative structure of the model, and the multiplicative innovation random variables.

[0050] In another embodiment, j = δ {m(i, j), I j , where are non-negative scalars, δ is the Kronecker delta,

ise ' ^ai m(i, j) is a one-to-one mapping from each combination of i and j to an index corresponding to / , (e.g., m(i, j) = (i - l)K + j, where K is a number of elements in the hidden variable h _n ) so that h{ _n =∑ ^a i,j ^£ m(i,j),nhj,n-l · T is embodiment j

enables flexibility in modeling the signal, because each transition can be inferred independently.

[0051] Another embodiment that is important to modeling multiple sources comprises partitioning hidden variables hi _n into S groups, where each group corresponds to one independent source in a mixture. Likewise, the non-negative random variables ε \ _n are partitioned according to the same S groups. This can be accomplished by a special case of the parameters where i = 0 when i _{z ?w} , and hj _n are not in the same group or when h _{j n} and ε\ _η are not associated with the same group. When the hidden variables are ordered accordingly, this gives a block structure, where each block corresponds to the model for one of the signal sources. [0052] In our embodiments, the hidden variables are related 140 to feature variables via a non-negative feature V f , of the signal indexed by feature /and

∑ (v) , (v)

C f i I ^{n £} l n ' j

C ^(V) (V)

where f ,i,l is a non-negative scalar, and ^ are independent non-negative random variables, and /, and / are indices of different components.

[0053] In a more constrained embodiment

wf _} i are non-negative scalars, where δ is the Kronecker delta, and s are the Gamma distributed random variables, so that the observation model based, at least in part, on

P( ^v f ,n \ h _n ) = Gamma(vf _n \ α ^ ^ν β^ /∑wf _t ih _iin ) ,

i

where Vf _>n is non-negative feature of the signal at frame n and frequency/, a ( ^v ) and are positive scalars, and Wf _ti are non-negative scalars.

[0054] In applications where the features xf, _n are complex spectrogram values of the input signal, a frame n and frequency/, the observation model can , which is the power in frame n, and frequency/. Thus, an

observation model can be formed based on

where V -1 is the unit imaginary number, and Θ j _n = Zjcy _n is a phase for a frame n and frequency/.

[0055] In another embodiment, we select the parameter & — 1 , so that the gamma distribution reduces to an exponential distribution as a special case. In this case, if the phases ej _n are distributed uniformly, then we obtain the observation model

P( ^x f ,n \ ^h n ) = ^N c (0,∑ ^w f ,i ^h i,n )

i

where N _c is a complex Gaussian distribution. This observation model corresponds to the Itakura-Saito nonnegative matrix factorization described above, and is combined in our embodiments with the non-negative dynamical system model.

[0056] Another embodiment uses an observation model for cascade of transformations of the same t e:

(U) (v) _£ (u) (v) where /' and ^c j are non-negative scalars,and Γ, ^η and V n are independent non-negative random variables, and i, are indices.

[0057] The method for inferring the hidden variables depends on the model parameterization for each embodiment.

[0058] Model Parameters

[0059] As shown in Fig. 2, from the input signal 102, we obtain the model parameters 101 as follows. The input signal can be considered a training signal, although it should be understood that the method can be adaptive to the signal, and "learn" the parameters on-line. The input signal can also be in the form of a digital signal or data.

[0060] For example, the training signal is a speech signal, or a mixed signal from multiple acoustic sources, perhaps including non-stationary noise, or other acoustic interference. The signal is processed as frames of signal samples. The sampling rate and number of samples in each frame is application specific. It is noted that the updating 230 described below for processing the current frame n is dependent on a previous frame n-1. For each frame we determine 210 a feature vector x _n representation. For an audio input signal, frequency features such as log power spectra could be used.

[0061] Parameters of the model are initialized 220. The parameters can include basis functions W, a transition matrix A, activation matrix H, and a fixed shape parameter and an inverse scale parameter β of a continuous gamma distribution parameter, and various combinations of these parameters depending on the particular application. For example in some applications, updating H and β are optional. In a variational Bayes (VB) method, H is not used. Instead an estimate of the posterior distribution of H is used and updated. If a maximum a-posteriori (MAP) estimation, then updating β is optional.

[0062] During each iteration of the method, the activation matrix, the basis function, the transition matrix, and the gamma parameter are updated 231-134. It should again be noted that the set of parameters to be updated is also application specific.

[0063] A termination condition 260, e.g., convergence or a maximum number of iterations, is tested after the updating 230. If true, store the parameters in a memory, otherwise if false, repeat at step 230.

[0064] The above steps of the general method and the parameter

determination can be performed in a processor connected to a memory and input/output interfaces as know. Specialized microprocessors, and the like can also be used. It is understood that the signals processed by the method, e.g., speech or financial data, can be extremely complex. The method transforms the input signal into features which can be stored in the memory. The method also stores the model parameters and inferred hidden variables in the memory.

[0065] Model Parameters Details [0066] For simplicity of this description, we limit the notation to the embodiment where ^ = d(i, l)w† j , thew f are non-negative scalars, δ is a Kronecker delta, and £ j 0) _n are gamma distributed random variables, with

( ^v ) 1 Θ

parameter CC — I , and phases / ,« are distributed uniformly. In this case, our model is ^x fn · ^' ^ ^w flc^kn ) , (10) where Χ is the complex-valued STFT coefficient at frame n and frequency/, N _c is the complex Gaussian distribution, Wffc is the value of the k ^th basis function for the power spetrum at frequency /, h _n and

columns of the activation matrix H, respectively, A is the nonnegative K K transition matrix that models the correlations between the different patterns in successive frames n-1 and n, £ _n is a nonnegative innovation random variable, e.g., a vector of dimension K, and o denotes entry-wise multiplication. The smooth IS-NMF can be obtained as a particular case of our model by setting A = I _j , where l _K is the K x K identity matrix.

[0067] Advantages

[0068] A distinctive and advantageous property of our model is that more than one state dimension can be non-zero at a given time. This means that a signal simultaneously acquired from multiple sources by a single sensor can analyzed using a single model, unlike the prior art HMM which requires multiple models.

[0069] Gamma Model of Innovations

[0070] We use an independent gamma distribution for the innovation ^ kn , namely ρ(ε _ίη \ α, β) = β(α _ι , β _ι ).

[0071] It follows that h _n is conditionally gamma distributed, such that

P(K \ Ah _n _ ₁ ) = YPihin \ _ί , β _ί /[ΑΗ _η _ _{1 ί} ),

i

and in particular

[0072] For /¾, we use an independent scale-invariant noninformative

Jeffreys prior i.e., p(h ) = In Bayesian probability, the Jeffreys prior

k

is a non-informative (objective) prior distribution on a parameter space that is proportional to the square root of the determinant of Fisher information.

[0073] MAP Inference in the Gamma Innovation Model

[0074] The maximum a-posteriori (MAP) objective function is

C(W,H,A^) =

+ ( ^N - l)∑(log ^T ( i ) - ^a i ^lQ g ft )

-∑l g^(¾)

[0075] Scales

[0076] Scale-Ambiguity Between A and /?

[0077] A K x K nonnegative diagonal matrix with coefficients λ _} on its diagonal is Λ , thus,

CQV,H, AA, Afi) = C(W,H, A, fi),

which has a scale-ambiguity between A and β . When both A and ? are estimated, the scale-ambiguity can be corrected in a number of ways, for example by fixing β to arbitrary values or by normalizing the rows of A at every iteration 230 and rescaling β accordingly. For example, we can normalize the rows of the transition matrix A such that the rows sum to 1, or so that the maximum coefficient in every row is 1. In some embodiments, ¾ = αι , i.e., the model expectation of the innovation random variable is 1. [0078] Ill-Posedness of MAP

[0079] The scales of W and H are related by

Ο \- ¹ ,ΑΗ,Α,β) = C(W,H,\ ^~l AK, fi) + JV^ logA _j , where i

λ i is the i-th element of the diagonal of Λ .

[0080] Without further constraints, the minimization of the MAP objective leads to a degenerate solution such that || W ||→∞ and || H \\→ 0 . If we assume that all the diagonal elements of Λ are equal, such that Λ = λΐ , then C(W\ ^~l ,hH,A) = C(W,H,A) + iSV log A.

[0081] The MAP objective can be made arbitrarily small by decreasing the value of λ . Hence, the norm of W is controlled during optimization. This can be achieved by hard or soft constraints. The hard constraint is a regular constraint that must be satisfied, and the soft constraint is a cost functions expressing a preference.

[0082] Hard Constraint

[0083] We solve

minC(W,H,A) si W≥ 0,H≥ 0,|| | w _k 1 using the change of variable W = WA , H = AH with Λ = diagfA^,. . .,λχ], and λ/ _ζ = PwfcPi, the norm-constraint can be relaxed by solving mm C(W,H,A) = D _IS (V \ WH) + S(AH) si W≥0,H≥ .

[0084] Soft Constraint (Penalization)

[0085] Another way we can control the norm of W is to add an appropriate penalty to the objective function, e.g.,

min CQV^^ + W W ^ si W≥0,H≥0

[0086] The soft constraint is typically simpler to implement than the hard constraint, but requires the tuning of λ .

[0087] Learning and Inference Procedures for MAP Estimation

[0088] We describe a majorization-minimization (MM) procedure. The MM is an iterative optimization procedure that can be applied to a convex objective function to determine maximums. That is, MM is a way to construct the objective function. MM determines a surrogate function that majorizes the objective function by driving the function to a local optimum. In our embodiments, the matrices H,A, and W are updated conditionnally on one and another. In the following, tildas (~) denote current parameter iterations.

[0089] Inequalities

[0090] For ¾ } such that = 1, we have

by Jensen's inequality. We can form an upper bound on log a by linearization, at any point φ ,

[0091] In particular, log ¾χ _έ < (log ¾¾ - 1) +— L- Ja _k x _k ,

* ^' k l ^a jXj k

[0092] Fit to Data

V fn

Pfk =

n v fn hkn n ^v n

Penalty terms

[0094] Update Rules

[0095] The MM framework includes majorizing the terms of the objective function with the previous inequalities, providing an upper bound of the objective function that is tight at the current parameters, and minimizing the upper bound instead of the original objective. This strategy applied to the minimization of the MAP objective with the soft constraint on the norm of W leads to the following updates 230 as shown in Fig. 2.

[0096] Update 231 Activation Matrix H

[0097] The columns of H are updated 231 sequentially. Left to right updates makes the update of h _n at iteration / dependent of ^ and r (1- -1) ' . The update of hfa involves rooting a polynomial of order 2, such that

where the values of a , b , c are given in the next table.

a

Si(n+l)

b 1 (Jeffreys) or 0 (uniform) c

~ ^h kn f ^j ¾, ₊ LA ^¾ i° ^,+1) ]

1 < n < N

a

&i(n+l) Skn . b

c

n = N

a fik

4kn +

Skn

b

c W

[0098] In particular, for the exponential innovation with expectation 1 (« = βΐ = 1), we obtain the following multiplicative updates: For n = 1,

i 8i(n+l) ^h kn

For 1 < n <

For n = N,

[0099] Update 232 Basis Function W

[0100] Update 233 Transition Matrix A

[0101] Variational EM Procedure for Maximum Likelihood Estimation

[0102] The activation partameter H is a latent variable to integrate from the joint likelihood. For generality, we assume the gamma distribution parameters β - { ¾ } ^t0 b ^e fr ^ee - The shape parameters α _¾ · are treated as fixed parameters. We minimize

Οθν,Α,β) =

[0103] This yields a better posed estimation problem because the set of parameters is of fixed-dimensionality w.r.t to the number of samples N.

Furthermore, the objective is now better posed in terms of scales. For any positive diagonal matrix Λ , we have

€( ,Α,β) =€( Α ^~1 ,ΑΑΑ ^~1 ,β)

so that the renormalization of solution W * only a renormalization of *

A . This is not true for the MAP approach.

[0104] For minimizing C(W, Α, β) , the EM procedure can be based on the complete dataset (V,H) , and on the iterative minimization of β(* I * ) = S _H ^l °gP( ^v > ^H I ^W P(H I V, 6 )dH ,

where Θ = {W, Α, β}. We do not use the posterior probability p(H \ V, Θ) . Instead, we use a variational EM procedure. For any probability density function q(H), the following inequality holds:

C(9)≤-(logp( \ WH)) _q - (logp(H \ A)) _q + (\ogq(H)) _q = B _q (0), where (') _q denotes the expectation under q (H) . Variational EM minimizes B _q (6) , instead of C(6) . At each iteration, the bound is first evaluated and tightened, given Wand A by minimizing B _q (0) over q, or more precisely, over the shape parameters of q , given a specific parametrized form, and then minimized with respect to (Θ) given q . Variational EM coincides with EM when q(H) = p(H \ θ) , in which case C(6) is decreased at every iteration. In other cases, variational EM conducts approximate inference. The validity depends on how well q(H) approximates the true posterior probability p(H \ θ) .

[0105] Derivation of the Bound

[0106] The expressions of log p(V \ WH) and log p(H \ A) show that the coefficients of H are coupled through ratios or logarithms of linear combinations . This makes expectations of

log p(V I WH) and log p{H \ A) very difficult to determine independently of the specific form of q(H) . [0107] Therefore, we majorize log p(V \ WH) and log p(H \ A) , to obtain a tractable bound. Using the above inequalities and assuming a factored form of the variational distribution, such that q(H) = J (/fo _j ) is an upper

kn

bound of C(W, Α,β), the function Β^( ,Α,β) = 0ogVyh-l)

+

+ f (l g¾>

i=l

+∑(l°g$(¾n)>

kn

Where φ _η are nonnegative coefficients s Vij _n are nonnegative coefficients such that

n _> ^/w are nonnegative coefficients, . denotes the set of all tuning parameters {φ _η ,v _t j _n , p _in ,ip }fi _n ij,

(·) denotes expectation w.r.t q , i.e., corresponds to (-)q . We remove subscript q to alleviate notations. W 201

[0108] The expression of the bound involves the expectation of hfa , and log h £ . These expectations are precisely the sufficient statistics of the generalized inverse-Gaussian (GiG), which is a practical convenience for q(H). We use

where

GIG(x I α,β, γ) =

and where K _a is a modified Bessel function of the second kind and X , β and γ are nonnegative scalars. Under the GIG distribution,

[0109] For any a , K _a+ i(x) = 2(a/x)K _a (χ) + Κ _α _ι(χ) , which leads to the alternative, implementation-efficient ex ression of

[0110] Optimization of the Bound [0111] We give the conditional updates of the various parameters of the bound. Update orders are described below.

[0112] Updates

Tuning parameters V

2 ^W fj ( ^h jn )>

∑ ^a ij ( ^h j(n-l) )

Variational distribution q

13] Parameters of interest

N n-2

aij = (21)

N

ai ∑ Pin ( ^h j(n-l))

n-2

[0114] Updating Order

[0115] We denote the set of of tuning parameter for frame n by ξ _η , i.e.,

[0116] As shown in Fig. 2, the following order of updates 230 leads to an efficient implementation.

[0117] At iteration (/) do

For n - Ι,. , .,Ν , W

[0118] Update 231 the activation parameters [q(h _n )]^ as a function of

[0119] Update

Update 232 the basis

ξ(21-1)

sition matrix A and the activation parameters, [q(H)

[0120] Under this updating order, the VB-EM procedure is:

Update q(H) .

(∑ )

Update ,Α,β

Determine the bound

[0121] Speech Denoising with the Dynamic Model

[0122] As shown in Fig. 3 for one embodiment, we use our method and model for speech enhancement, e.g., denoising. We construct our model parameters 101 for speech 306 by estimating bases W and the transition matrix^ on some speech (audio) training data 305 as described above. We denote the trained bases and transition matrix as and . where (s) is speech.

[0123] Similarly, we construct a noise model 307 with bases and transition matrix , and combining the two models 306-307 into the single model 300 by concatenating and into W = and and into A, where A is a block-diagonal matrix with A^ and

A ^} on the diagonal.

[0124] We can also train for noise on some noise training data, or we can fix the speech part of the model, and train for the noise part on the test data, thus making the noise part a general model that collects parts of the signal that cannot be modeled by the speech model. The simplest version of the later model uses a single basis for the noise, and uses an identity matrix as the transition matrix

[0125] After the model 300 is constructed, we can use the model to enhance an input audio signal x 301. We determine 310 a time-frequency feature

representation. We estimate 320 the parameters of the model 300 that vary, i.e., the activation matrix H ^{K }} for the speech and H ^} for the noise (n), and the bases and transition matrix for the noise.

[0126] Thus, we obtain a single model that combines speech, W^H^

(n) (ri)

and noise W ^,l >H ^> , which we then use to reconstruct 330 the complex STFT of the enhanced speech x 340 using

[0127] The time-domain signal can be reconstructed using a conventional overlap-add method, which evaluates a discrete convolution of a very long input signal with a finite impulse response filter

[0128] Extensions

[0129] Other complex models can also generated based on the above embodiments. [0130] Dirichlet Innovations

[0131] Instead of considering the innovation random variables ε _η to be gamma distributed, the innovation can be Dirichlet distributed, which is similar to a normalization of the activation parameter h _n .

[0132] HMM-Like Behavior

[0133] We can constrain h _n to be 1 -sparse during inference. [0134] Structured Variational Inference

[0135] Conventional variational inference assumes that the variational posteriror probabilities q(h _n ) are independent of each other, which, given a strong dependency relation between h _n and h _n _y, is likely to be very wrong. We can model the posteriror probability in terms of q(h _n \ h _n _\) . One possibility for such a q distribution uses a GIG distribution with parameters dependent on Ah _n _

[0136] Gamma Distribution of Innovation

[0137] The complex Gaussian model on the complex STFT coefficients in Eqn. (6) is equivalent to assuming that the power is exponentially distributed with parameter WH . We can extend the model by assuming that the power is gamma distributed, thus leading to a donut-shaped distribution for the complex coefficients.

[0138] Full Covariance of Innovation Random Variables

[0139] In linear dynamical systems, the innovation random variables can have a full-covariance. For positive random variables, one way to include the correlations is to transform an independent random vector with a non-negative matrix. This leads to the model,

K = (Ah _n _ _l ) o (Bf _n ) , where f _n is a nonnegative random vector of size / x 1 and B is a nonnegative matrix of dimension K x J . When B = IfcxK » ^mis simplifies to fn = ^£ n . This can be accom lished in the more general form of the model ctorized form: are the elements of B.

[0140] Transition Innovations

[0141] It can also be useful to model the transition between each of the components of h _n and h _n _\ using separate innovation random variables. This is analogous to the use of Dirichlet prior probabilities in discrete Markov models. One method would admit h _n - {A ° E _n )h _n _i, where E _n is a nonnegative innovations matrix of dimension K K . This can be accomplished in the more general form of the model hf _n = C( j / £ _n h _n _\ by setting the parameters u

i , where are the elements of A and m(i, j) is a one-to-one mapping from each combination of and j to an index corresponding to

/ . Then, the i j ^' -th element of E _n is ^m(i,j ),n .

[0142] Considering Other Innovation Types Besides Gamma

[0143] A log-normal, Poisson distribution leads to yet different types of dynamical systems.

[0144] Considering Other Divergences

[0145] We so far only considered the Itakura-Saito divergence. We can also use the KL-divergence, and different divergences for h _n \ _n _\ and for v | h .

[0146] Online Procedure

[0147] For real-time applications, only the signal up to the current time is used, e.g., an application where only the activation matrix H are estimated, or another application where all parameters are optimized. In the later application, we can perform a "warm" start with pretrained bases W and transition matrix A .

[0148] Multi-Channel Version [0149] Because our model relies on a generative model involving the complex STFT coefficients, the model can be extended to a multi-channel application. Optimization in this setting involves EM updates between mixing system and a source NMF procedure.

[0150] Effect of the Invention

[0151] The embodiments of the invention provide a non-negative linear dynamical system model for processing non-stationary signals, particularly speech signals mixed with noise. In the context of speech separation and speech denoising, our model adapts to signal dynamics on-line, and achieves better performance than conventional methods.

[0152] Conventional models for signal dynamics frequently use hidden Markov models (HMMs) or non-negative matrix factorization (NMF). HMMs lead to combinatorial problems due to the discrete state space, are computationally complex, especially for mixed signals from several sources, and make it difficult to handle gain adaptation. NMF solves both the computational complexity and gain adaptation problems. However, NMF does not take advantage of past observations of a signal to model future observations of that signal. For signals with predictable dynamics, this is likely to be suboptimal.

[0153] Our model has advantages of both the HMMs and the NMF. The model is characterized by a continuous non-negative state space. Gain adaptation is automatically handled during inference. The complexity of the inference is linear in the number of sources, and dynamics are modeled via a linear transition matrix.