System and method for training a class-specific hidden Markov model using a modified Baum-Welch algorithm

[0002] The invention described herein may be manufactured and used by or for the Government of the United States of America for governmental purposes without the payment of any royalties thereon or therefore.

[0003] (1) Field of the Invention

[0004] This invention relates to systems and methods for modeling physical phenomena, and more particularly to a system and method for modeling physical phenomena, such as speech, using a class-specific implementation of the Baum-Welch algorithm for estimating the parameters of a class-specific hidden Markov model (HMM).

[0005] (2) Description of the Prior Art

[0006] By way of example of the state of the art, reference is made to the following papers, which are incorporated herein by reference. Not all of these references may be deemed to be relevant prior art.

[0007] P. M. Baggenstoss, “Class-specific features in classification,” IEEE Trans. Signal Processing, December 1999.

[0008] S. Kay, “Sufficiency, classification, and the class-specific feature theorem,” to be published IEEE Trans Information Theory.

[0009] B. H. Juang, “Maximum likelihood estimation for mixture multivariate stochastic observations of Markov chains,” AT&T Technical Journal, vol. 64, no. 6, pp. 1235-1249, 1985.

[0010] L. R. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition,” Proceedings of the IEEE, vol. 77, pp. 257-286, February 1989.

[0011] L. E. Baum, “A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains,” Ann. Math. Stat., vol. 41, pp. 164-171, 1970.

[0012] E. H. Lehmann, Theory of Point Estimation, New York: Wiley, 1983.

[0013] S. Kay, Modern Spectral Estimation; Theory and Applications. Prentice Hall, 1988.

[0014] E. J. Hannan, Multiple Time Series, Wiley, 1970.

[0015] M. H. Quenouille, “The joint distribution of serial correlation coefficients,” Ann. Math. Stat., vol. 20, pp. 561-571, 1949.

[0016] Many systems, e.g., communication, data processing and other information systems, can be described or characterized in terms of a series of transitions through a set of states. Hidden Markov models (HMMs) have found applications in modeling physical phenomena characterized by a finite number of states. Often these states represent distinct physical phenomena. In speech, for example, the human voice is characterized by distinct physical phenomena or modes, e.g., voiced speech, fricatives, stops, and nasal sounds. When applied to speech processing applications, the speech modes or components are first modeled by HMMs using an algorithm to estimate parameters for the HMMs (referred to as the training phase). The trained HMMs can then be used to determine which speech components are present in a speech signal (referred to as the recognition phase).

[0017] For the classical hidden Markov model (HMM), all observations are assumed to be realizations of a random statistical model that depends on the Markov state. Although the statistical models, i.e. the observation probability density functions (PDF's), are different for each state, they are defined on the same observation space. The dimension of this observation space needs to be high to adequately observe the information content of the data for all states. The high dimension requires a large number of observations (or training samples) and leads to poor performance with limited amounts of training data.

[0018] In speech, for example, a different set of parameters control the uttered sound during each of the speech modes. Furthermore, a distinct type of signal processing or feature extraction is best suited to estimate the corresponding parameters of each mode. But, since one cannot know a priori which mode is in effect at a given instant of time and cannot change the observation space accordingly, it is necessary to operate in a unified observation space.

[0019] This requires a feature set that carries enough information for the estimation of all modes. This in turn leads to dimensionality issues since there is only a finite amount of data with which to train the observation PDF estimates. In effect, the observation PDF's of each state are represented using a feature set with higher dimension than would be necessary if the other states did not exist. The amount of data required to estimate a PDF is exponentially dependent on feature dimension. Given limitations of computer storage and available data, feature dimensions above a certain point are virtually impossible to accurately characterize. As a result, one may be forced to use a subset of the intended feature set to reduce dimension or else suffer the effects of insufficient training data.

[0020] Consider a hidden Markov model (HMM) for a process with N states numbered S₁through S_N. Let the raw data be denoted X[t], for time steps t=1,2, . . . , T. The parameters of the HMM, denoted λ, comprise the state transition matrix A={a_ij}, the state prior probabilities u_j, and the state observation densities b_j(X), where i and j range from 1 to N. These parameters can be estimated from training data using the Baum-Welch algorithm, as disclosed in the papers by Rabiner and Juang. But, because X[t] is often of high dimension, it may be necessary to reduce the raw data to a set of features z[t]=T(X[t]). We then define a new HMM with the same A and u_jbut with observations z[t], t=1,2, . . . T and the state densities b_j(z) (we allow the argument of the density functions to imply the identity of the function, thus b_j(X) and b_j(z) are distinct).

[0021] This is the approach used in speech processing today where z[t] are usually a set of cepstral coefficients. If z[t] is of low dimension, it is practical to apply probability density function (PDF) estimation methods such as Gaussian Mixtures to estimate the state observation densities. Such PDF estimation methods tend to give poor results above dimensions of about 5 to 10 unless the features are exceptionally, i.e., “well-behaved” are close to independent or multivariate Gaussian. In human speech, it is doubtful that 5 to 10 features can capture all the relevant information in the data. Traditionally, the choices have been (1) use a smaller and insufficient features set, (2) use more features and suffer PDF estimation errors, or (3) apply methods of dimensionality reduction. Such methods include linear subspace analysis, projection pursuit, or simply assuming the features are independent (a factorable PDF). All these methods involve assumptions that do not hold in general.

[0022] The class-specific method was recently developed as a method of dimensionality reduction in classification, as disclosed in U.S. patent application Ser. No. 09/431,716 entitled “Class Specific Classifier.” Unlike other methods of dimension reduction, it is based on sufficient statistics and results in no theoretical loss of performance due to approximation. Because of the exponential relationship between training data size and dimension, even a mere factor of 2 reduction in dimension can result in a significant difference.

[0023] Accordingly, one object of the present invention is to reduce the number of data samples needed for training HMMs.

[0024] Another object of the present invention is to extend the idea of dimensionality reduction in classification to the problem of HMM modeling when each state of the HMM may have its own minimal sufficient statistic.

[0025] A further object of the present invention is to modify the Baum-Welch algorithm used to estimate parameters of class-specific HMMs.

[0026] The foregoing objects are attained by the method and system of the present invention. The present invention features a method of training a class-specific hidden Markov model (HMM) used for modeling physical phenomena characterized by a finite number of states. The method comprises the steps of receiving training data forming an observation sequence; estimating parameters of the class-specific HMM from the training data using a modified Baum-Welch algorithm, wherein the modified Baum-Welch algorithm uses likelihood ratios with respect to a common state and based on the sufficient statistics for each state; and storing the parameters of the class-specific HMM for use in processing signals representing the physical phenomena.

[0027] The step of estimating parameters of the class-specific HMM preferably includes the step of conducting a plurality of iterations of the Baum-Welch algorithm, wherein each iteration of the Baum-Welch algorithm includes a class-specific forward procedure for calculating forward probabilities, a class-specific backward procedure for calculating backward probabilities, and HMM reestimation formulas for updating the parameters of the class-specific HMM based upon the forward probabilities and the backward probabilities.

[0028] The present invention also features a system for training a class-specific hidden Markov model (HMM) used for modeling physical phenomena characterized by a finite number of states. The system includes means for receiving training data represented as an observation sequence; means for estimating parameters of the class-specific HMM from the training data using a modified Baum-Welch algorithm, wherein the modified Baum-Welch algorithm uses likelihood ratios that compare each of the states to a common state; and means for storing the parameters of the class-specific HMM for use in processing signals representing the physical phenomena.

[0029] According to one example, the HMMs are used in speech processing and the physical phenomena includes modes of speech.

[0030] These and other features and advantages of the present invention will be better understood in view of the following description of the invention taken together with the drawings wherein:

[0031] FIG. 1 is a schematic block diagram of a conventional system for training HMMs used in speech processing;

[0032] FIG. 2 is a flow chart illustrating the method of training HMMs in accordance with the prior art Baum-Welch algorithm;

[0033] FIG. 3 is a flow chart illustrating a method of training HMMs for use in signal processing, according to the present invention; and

[0034] FIGS. 4 and 5 are diagrams illustrating state classification error probabilities as a function of number of training records for the method of the present invention compared to conventional methods.

[0035] The method and system, according to the present invention, is used to train class-specific hidden Markov models (HMM's) for use in modeling physical phenomena characterized by a finite number of states. In the exemplary embodiment, the system and method is used in a speech processing system to train class-specific HMMs modeling different modes or components of speech. The method of the present invention can also be applied to other types of systems for modeling other types of physical phenomena.

[0036] A conventional speech processing system 10, FIG. 1, generally includes a preprocessor 12, a storage medium 14, a training processor 16, and a recognition processor 18. The preprocessor 12 receives a “raw” input speech signal 20 and transforms the speech signal 20 into acoustic speech units, in the form of digital data, which are used to generate the HMMs. The preprocessor 12 samples the speech signal, extracts the features from the time series, and feeds them to a desired storage space (e.g., on the storage medium 14). The preprocessor can include any conventional preprocessing device known in the art since the method of the present invention is not dependent upon the structure of the preprocessor 12.

[0037] The training processor 16 is used to generate the estimated parameters of the HMMs and can include any suitable computer known in the art programmed to perform the training method of the present invention. The training processor 16 stores the HMMs and parameters for use by the speech recognition processor 18. The speech recognition processor 18 can also include any suitable computer known in the art programmed to perform speech recognition using the HMMs generated and stored by the training processor. The preprocessor 12,. a storage medium 14, a training processor 16, and a recognition processor 18 can be implemented as separate devices or as separate modules within a single computer.

[0038] FIG. 2 shows a method of training HMMs according to the prior art using the Baum-Welch algorithm. Training data is received in step 22. A feature set is computed in step 24. This computation extracts features which will allow one state of the HMM to be distinguished from another state. Step 26 calculates forward probabilities and step 28 calculates backward probabilities. HMM reestimation formulas are applied to the probability data function in step 30. Guassian mixture reestimation formulas are applied in step 32. In step 34, the algorithm determines if the probability data function has improved over the previous iteration. Steps 26-34 are performed again if the probability data function is still improving. Once the PDF has converged, the trained HMM with estimated parameters can be stored as in step 36.

[0039] The method of training the class-specific HMMs, according to the present invention, is illustrated in FIG. 3. First, the training data is received, step 110. In the exemplary embodiment, the training data is the digital data stored in the storage medium 14 by the preprocessor 12. Next, the parameters of the class-specific HMM are estimated from the training data using a modified Baum-Welch algorithm, steps 115-160, as will be described in greater detail below. In the exemplary embodiment, this algorithm is performed by the training processor 16. The trained class-specific HMM with the estimated parameters is then stored, step 170, for example, in the storage medium 14 to be used by the speech recognition processor 18.

[0040] The modified Baum-Welch algorithm is preferably an iterative algorithm that calculates new estimates of the HMM parameters given an observation sequence and a previous estimate of the HMM parameters. The Baum-Welch algorithm is generally composed of two parts: the forward/backward procedure and the reestimation of parameters. In step 115, the modified Baum-Welch algorithm of the present invention first extracts features which will allow distinction of one state or class from the common class. The current invention then applies a state-specific forward procedure to calculate forward probabilities, step 120, and a class-specific backward procedure to calculate backward probabilities, step 130. The modified Baum-Welch algorithm then applies HMM reestimation formulas based upon the forward and backward probabilities to update the parameters, step 140. The method can also apply Gaussian mixture reestimation formulas, step 150. The steps 120-150 of the Baum-Welch algorithm are repeated until convergence.

[0041] As can be seen in FIG. 3, features and PDFs are computed for each state. The number of states can be any number from 1 to N. After features are calculated for a particular state, the Baum-Welch algorithm is used to refine PDF for the state.

[0042] The modified Baum-Welch algorithm allows each state to be observed using a different set of features rather than relying on a common feature set. Each feature set (or observation space) is chosen to be a sufficient statistic for discrimination of the given state from a common “white-noise” state. Comparison of likelihood values is possible through the use of likelihood ratios with respect to the common statistical model. The new algorithm is the same in theory as the algorithm based on a common feature set, but without the necessity of estimating high-dimensional probability density functions (PDFs).

[0043] In speech, for example, the word “Act” has the voiced sound “A” and unvoiced sounds “c” and “t”. The voiced and unvoiced sounds are produced by distinctly different mechanisms. The class-specific method considers each sound separately in a set of binary tests to a common (noise-only) condition. Let us ignore for the moment the relationship between these sounds in time and consider only the static classification problem. We consider the binary test of “A” vs. white noise, the test of “c” vs. white noise, and the test of “t” vs. white noise. Each of these simpler tests may have a distinct sufficient or approximately sufficient class-specific statistic. The sufficient statistic for these simpler tests are almost certainly lower in dimension than a sufficient statistic for the combined (3-way classification) problem. It can be shown that the union of the three class-specific statistics is a sufficient statistic for the combined problem based on a test that in theory has the same probability of error as the optimal test based on the raw data. This union of features is a sufficient statistic for the HMM as well. A class-specific HMM can be constructed in a manner that involves only PDF estimates involving the separate class-specific feature vectors, not the union.

[0044] To establish some notation, the classic HMM is reviewed using the notation from the Juang paper (with the exception that the PDF's are represented as p(·) and b_j(·) instead of f(·) and f_j(·)). Consider a set of state occurrences θΔ{θ₀, θ₁, . . . θ_T} which are the result of a Markov chain with state priors u_j, j=1,2, . . . , N and state transition matrix A. Rather than observing the states θ directly, we observe the stochastic outputs SΔ{s₀, s₁, . . . , s_T} which are realizations from a set of state PDF's b_j(s), j=1,2, . . . , N such that: $\begin{array}{cc}p\left(S❘\theta \right)=\prod _t{b}_{{\theta }_i}\left(s_t\right)& \left(1\right)\end{array}$

[0045] Note that s is used to denote an observation at an unspecified time, while the bold notation S denotes the collection of all observations from t=1, . . . , T. The complete set of parameters λ defining the HMM are $\begin{array}{cc}\lambda \underset{\_}{\underset{\_}{\Delta }}\left\{\begin{array}{ccc}u\underset{\_}{\underset{\_}{\Delta }}{u}_j,& j\in \left[1\dots N\right],& \\ A\underset{\_}{\underset{\_}{\Delta }}{a}_{\mathrm{ij}},& j\in \left[1\dots N\right]& j\in \left[1\dots N\right],\\ b_j& j\in \left[1\dots N\right],& \end{array}\right\}& \left(2\right)\end{array}$

[0046] The class-specific HMM is defined as follows. Let the observations s_tbe each passed into a bank of N feature transformations 115. For each observation s_t, we have features, z_j,t=T(s_t), j=1, . . . , N. These features are sufficient statistics, as described below, and are of lower dimension than s_t. The goal is to derive an estimation algorithm for the Markov model based on these lower-dimensional statistics. The parameter vector λ for the class-specific HMM is denoted

[0047] λ^z{u,A,{b_j^zjε[1 . . . N]}}  (3)

[0048] where b_j^z(·) is the PDF of z_j=T_j(s) when state j is true. This does not require the distribution of z_jwhen state k is true (k≠j)).

[0049] To apply the class-specific method, we need to define a condition that is common to all states. Let the PDF of s in each state j be dependent on an underlying unknown parameter set λ_jwith distribution f_j(λ_j), thus

[0050] b_j(s)=∫_λb_j(s|λ_j)f_j(λ_j)dλ_j, jε[1 . . . N].  (4)

[0051] In most, if not all, signal processing applications, the observations s. are corrupted by noise. Let the parameter set for each state include an amplitude parameter a_j

[0052] λ_j=[a_jλ′_j]  (5)

[0053] such that there is a common noise-only condition H₀defined by

[0054] H₀: a₁=a₂= . . . a_N=0.  (6)

[0055] Let the distribution of s under H₀be denoted b₀(s):

[0056] b₀(s)Δb₁(s|[0λ′₁]=b₂(s|[0λ′₁]= . . . =b_N(s|[0λ′_N]).  (7)

[0057] We assume s_tare independent under H₀. Thus, $\begin{array}{cc}p\left(S❘H_0\right)=\prod _{t=1}^T{b}_0\left(s_t\right).& \left(8\right)\end{array}$

[0058] While this structure does not seem to fit many problems of interest, any problem can be modified to include an amplitude parameter even if the amplitude parameter is never zero. Furthermore, the noise-only condition H₀can have an arbitrarily small assumed variance because H₀is only a theoretical tool and does not need to approximate any realistic situation. The choice of H₀, however, affects the choice of sufficient statistics.

[0059] For each state j, let there exist a sufficient statistic for z=T_j(s) such that the Neyman-Fisher factorization holds:

[0060] b_j(s|λ_j)₃₂g_j(z_j;λ_j)h_j(s).  (9)

[0061] It is well known, and follows from the above factorization, that the likelihood ratio is invariant if written in terms of a sufficient statistic: $\begin{array}{cc}\frac{b_j\left(s❘{\lambda }_j\right)}{b_0\left(s\right)}=\frac{b_j^z\left(z_j;{\lambda }_j\right)}{b_{0,j}^z\left(z_j\right)}& \left(10\right)\end{array}$

[0062] where b_j^z(z_j;λ_j) and b_0,j^z(z_j) are the PDFs of z_junder state j and H₀. Marginalizing both sides over λ_j, $\begin{array}{cc}\frac{{\int }_{{\lambda }_j}{b}_j\left(s❘{\lambda }_j\right)f_j\left({\lambda }_j\right)d{\lambda }_j}{b_0\left(s\right)}=\frac{{\int }_{{\lambda }_j}{b}_j^z\left(z_j;{\lambda }_j\right)f_j\left({\lambda }_j\right)d{\lambda }_j}{b_{0,j}^z\left(z_j\right)}& \left(11\right)\end{array}$

[0063] we obtain the relationship $\begin{array}{cc}\frac{b_j\left(s\right)}{b_0\left(s\right)}=\frac{b_j^z\left(z_j\right)}{b_{0,j}^z\left(z_j\right)}& \left(12\right)\end{array}$

[0064] Thus, we say that z_jis a sufficient statistic for the binary hypothesis test for state j vs. H₀. The distributions of the statistics under H₀:b₀(s) and b_0,j^z(z_j) are assumed to be known. Often we let H₀be white Gaussian noise for which it may be possible to derive the densities analytically.

[0065] We assume the following Gaussian mixture representation for each b_j^z(z_j): $\begin{array}{cc}{b}_j^z\left(z_j\right)=\sum _{k=1}^M{c}_{\mathrm{jk}}N\left(z_j,{\mu }_{\mathrm{jk}}^z,U_{\mathrm{jk}}^z\right),1\le j\le N.& \left(13\right)\end{array}$

[0066] where K is the mode number

[0067] N(z_j,μ, U)Δ(2π)^−P/2|U|^−1/2exp{−½(s−μ)′U⁻¹(s−μ)},  (14)

[0068] and P_jis the dimension of z_j. Note that there is no problem assuming that M is independent of j as long as M is sufficiently large for the worst case.

[0069] Central to the Baum-Welch algorithm is the ability to calculate the probability that a state is true at a particular time given an observation sequence. These probabilities depend on the forward and backward probabilities, which in turn depend on the state likelihood functions. It is possible to re-formulate the problem to compute these probabilities using low-dimensional PDF's obtained by replacing likelihood functions with likelihood ratios based on the sufficient statistics for each state. When compared to the standard feature-based HMM using a common feature set, the parameters may be estimated using fewer training samples or for a fixed number of training samples, more accuracy and robustness can be obtained.

[0070] The class-specific Baum-Welch algorithm is now described in detail. The HMM can be used as a classifier by calculating p(S|λ). This is accomplished efficiently by the forward procedure. A far more difficult problem is the estimation of λ from training data: $\begin{array}{cc}\underset{\lambda \in \Lambda }{\max }p\left(S❘\lambda \right).& \left(15\right)\end{array}$

[0071] A well-behaved iterative algorithm based on the Expectation Maximization (EM) method is available as defined by Juang. Formulas for updating the parameters λ at each iteration are called the reestimation formulas. The derivation by Juang is well known for the case when b_j(s) are Gaussian mixtures. The derivation of Juang is modified to estimate the parameters of Gaussian mixture representations for b_j^z(z_j).

[0072] We begin with the global density for the HMM (Juang eq. 3) $\begin{array}{cc}p\left(S❘\lambda \right)=\sum _{\theta \in \Theta }{u}_{{\theta }_0}\prod _{t=1}^T{a}_{{\theta }_{t-1^{{\theta }_t}}}{b}_{{\theta }_t}\left(s_t\right).& \left(16\right)\end{array}$

[0073] We may re-write the above as $\begin{array}{cc}p\left(S❘\lambda \right)=\frac{p\left(S❘\lambda \right)}{p\left(S❘H_0\right)}p\left({SH}_0\right)=\left[\sum _{\theta \in \Theta }{u}_{{\theta }_0}\prod _{t=1}^T{a}_{{\theta }_{t-1^{{\theta }_t}}}\frac{b_{{\theta }_t}\left(s_t\right)}{b_0\left(s_t\right)}\right]p\left(S❘H_0\right)& \left(17\right)\end{array}$

[0074] using equation (8). By the sufficiency relationship of equation (12), we have $\begin{array}{cc}p\left(S❘\lambda \right)=\left[\sum _{\theta \in \Theta }{u}_{{\theta }_0}\prod _{t=1}^T{a}_{{\theta }_{t-1^{{\theta }_t}}}\frac{b_j^z\left(T_j\left(s_t\right)\right)}{b_{j,0}^z\left(T_j\left(s_t\right)\right)}\right]p\left(S❘H_0\right).& \left(18\right)\end{array}$

[0075] Thus, $\begin{array}{cc}p\left(S❘\lambda \right)=\left[\sum _{\theta \in \Theta }{u}_{{\theta }_0}\prod _{t=1}^T{a}_{{\theta }_{t-1^{{\theta }_t}}}\frac{\sum _{k=1}^M{c}_{\mathrm{jk}}N\left(T_j\left(s_t\right),{\mu }_{\mathrm{jk}}^z,U_{\mathrm{jk}}^z\right)}{b_{j,0}^z\left(T_j\left(s_t\right)\right)}\right]p\left(S❘H_0\right).& \left(19\right)\end{array}$

[0076] We define $\begin{array}{cc}{b}_{\mathrm{jk}}^z\left(s,{\lambda }^z\right)\underset{\_}{\underset{\_}{\Delta }}\frac{N\left(T_j\left(s\right),{\mu }_{\mathrm{jk}}^z,U_{\mathrm{jk}}^z\right)}{b_{j,0}^z\left(T_j\left(s_t\right)\right)}p\left(S❘H_0\right)& \left(20\right)\end{array}$

[0077] (note that b_jk^z(·) and b_0,j^z(·) where j=1,2, . . . , N and k=1,2, . . . , M are unrelated). We then have $\begin{array}{cc}p\left(S❘{\lambda }^z\right)=\sum _{\theta \in \Theta }{u}_{{\theta }_0}\prod _{t=1}^T\left[a_{{\theta }_{t-1^{{\theta }_t}}}\sum _{k=1}^M{c}_{{\theta }_{t}k}{b}_{{\theta }_{t}k}^z\left(s_t,{\lambda }^z\right)\right]& \left(21\right)\end{array}$

[0078] which is the same as (Juang eq. 8). We then may follow the proof of Juang, provided the necessary requirements of b_jk^z(sλ^z) are met.

[0079] Notice that b_jk^z(s,λ^z) depends on λ^zonly through the multivariate Gaussian density. Thus, b_jk^z(s,λ^z) meets the necessary log-concavity and elliptical symmetry requirements necessary for the reestimation formulas that follow. We can proceed in the proof of Juang, until it is necessary to differentiate b_jk^z(s_t,λ^z) with respect to λ^z. At that point, the additional terms in equation (20), not dependent on λ^z, do not affect the solution of the maximization; we only need to let z_j,ttake the place of s_t. For recursive computation of the forward and backward probabilities, we could use $\begin{array}{cc}{b}_j\left(s_t\right)=b_0\left(s_t\right)\frac{b_j^z\left(T_j\left(s_t\right)\right)}{b_{j,0}^z\left(T_j\left(s_t\right)\right)},& \left(22\right)\end{array}$

[0080] which follows directly from equation (12). However, this requires knowledge of s_t. If we adopt the class-specific approach, we would like to forget s_taltogether. If we use instead $\begin{array}{cc}{b}_j^c\left(s_t\right)=\frac{b_j^z\left(T_j\left(s_t\right)\right)}{b_{0,j}^z\left(T_j\left(s_t\right)\right)},& \left(23\right)\end{array}$

[0081] the resulting forward and backward probabilities, which we rename α_t^c(i) and β_t^c(j) are related to α_t(i) and β_t(j) through: $\begin{array}{cc}{\alpha }_t^c\left(j\right)=\frac{{\alpha }_t\left(j\right)}{\prod _{u=1}^t{b}_0\left(s_u\right)}\mathrm{and}& \left(24\right)\\ {\beta }_t^c\left(j\right)=\frac{{\beta }_t\left(j\right)}{\prod _{u=t+1}^t{b}_0\left(s_u\right)},& \left(25\right)\end{array}$

[0082] which yields at termination $\begin{array}{cc}\sum _{j=1}^N{\alpha }_T^c\left(j\right)=\frac{p\left(S❘{\lambda }^z\right)}{p\left(S❘H_0\right)}.& \left(26\right)\end{array}$

[0083] Since these quantities are used with normalization, it may be verified that the added factors have no affect on the reestimation formulas.

[0084] Accordingly, the modified Baum-Welch algorithm is as follows:

[0085] The class-specific forward procedure

[0086] 1. Initialization: $\begin{array}{cc}{\alpha }_1^c\left(i\right)=u_i\frac{b_i^z\left(z_{i,1}\right)}{b_{0,i}^z\left(z_{i,1}\right)},1\le i\le N& \left(27\right)\end{array}$

[0087] 2. Induction: $\begin{array}{cc}{\alpha }_{t=1}^c\left(j\right)=\left[\sum _{i=1}^N{\alpha }_t^c\left(i\right)a_{\mathrm{ij}}\right]\frac{b_j^z\left(z_{j,t+1}\right)}{b_{0,j}^z\left(z_{j,t+1}\right)},1\le t\le T-1,1\le j\le N& \left(28\right)\end{array}$

[0088] 3. Termination: $\begin{array}{cc}\frac{p\left(S❘{\lambda }^z\right)}{p\left(S❘H_0\right)}=\sum _{i=1}^N{\alpha }_T^c\left(i\right)& \left(29\right)\end{array}$

[0089] The class-specific backward procedure

[0090] 1. Initialization:

[0091] β_T^c(i)=1  (30)

[0092] 2. Induction: $\begin{array}{cc}{\beta }_t^c\left(i\right)=\sum _{j=1}^N{a}_{\mathrm{ij}}\frac{b_j^z\left(z_{j,t+1}\right)}{b_{0,j}^z\left(z_{j,t+1}\right)}{\beta }_{t+1}^c\left(j\right),t=T-1,T-2,\dots 11\le i\le N& \left(31\right)\end{array}$

[0093] HMM Reestimation formulas Define γ_t(j) as p(θ_t=j|S). We have $\begin{array}{cc}{\gamma }_t\left(j\right)=\frac{{\alpha }_t^c\left(j\right){\beta }_t^c\left(j\right)}{\sum _{i=1}^N{\alpha }_t^c\left(i\right){\beta }_t^c\left(j\right)}& \left(32\right)\end{array}$

[0094] Let $\begin{array}{cc}{\xi }_t\left(i,j\right)=\frac{{\alpha }_t^c\left(i\right)a_{\mathrm{ij}}\frac{b_j^z\left(z_{j,t+1}\right)}{b_{0,j}^z\left(z_{j,t+1}\right)}{\beta }_{t+1}^c\left(j\right)}{\sum _{i=1}^N\sum _{m=1}^N{\alpha }_t^c\left(i\right)a_{\mathrm{im}}\frac{b_m^z\left(z_{m,t+1}\right)}{b_{0,m}^z\left(z_{m,t+1}\right)}{\beta }_{t+1}^c\left(m\right)}& \left(33\right)\end{array}$

[0095] The updated state priors are

[0096] û_i=γ_l(i).  (34)

[0097] The updated state transition matrix is $\begin{array}{cc}{\hat{a}}_{\mathrm{ij}}=\frac{\sum _{t=1}^{T-1}{\xi }_t\left(i,j\right)}{\sum _{t=1}^{T-1}{\gamma }_t\left(i\right)}.& \left(35\right)\end{array}$

[0098] Keep in mind that γ_t(j), ξ_t(i,j), û_i, and û_ijwill be identical to those estimated by the conventional approach if equation (12) is true. There is no correspondence, however, between the Gaussian mixture parameters of the two methods except that the densities they approximate obey equation (12).

[0099] Gaussian mixture reestimation formulas

[0100] Let $\begin{array}{cc}{\gamma }_t^c\left(j,m\right)={\gamma }_t\left(j\right)\left[\frac{c_{\mathrm{jm}}N\left(z_{j,t},{\mu }_{\mathrm{jm}}^z,U_{\mathrm{jm}}^z\right)}{b_j^z\left(z_{j,t}\right)}\right]& \left(36\right)\\ {\hat{c}}_{\mathrm{jm}}^z=\frac{\sum _{t=1}^T{\gamma }_t^c\left(j,m\right)}{\sum _{t=1}^T\sum _{l=1}^M{\gamma }_t^c\left(j,l\right),}& \left(37\right)\\ {\hat{\mu }}_{\mathrm{jm}}^z=\frac{\sum _{t=1}^T{\gamma }_t^c\left(j,m\right)z_{j,t}}{\sum _{t=1}^T{\gamma }_t^c\left(j,m\right)}& \left(38\right)\\ \mathrm{and}& \\ {\hat{U}}_{\mathrm{jm}}^z=\frac{\sum _{t=1}^T{\gamma }_t^c\left(j,m\right)\left(z_{j,t}-{\mu }_{\mathrm{jm}}^z\right){\left(z_{j,t}-{\mu }_{\mathrm{jm}}^z\right)}^{\prime }}{\sum _{t=1}^T{\gamma }_t^c\left(j,m\right)}& \left(39\right)\end{array}$

[0101] Since the class-specific approach is based on sufficient statistics, one must find the sufficient statistics. In real-world problems, the sufficient statistics are not well defined because no mathematical formulas are known for the densities. If we look at conventional classifiers and HMM's, we often make the “leap” from raw data to a feature set without being overly concerned if the feature set is sufficient or not. The same is true for the class-specific method. We should be no more concerned about exact sufficiency than we are for conventional feature-based approaches. Sufficiency is used for the theoretical development only—the theory then serves as a guide to the practice.

[0102] One way to arrive at approximate sufficiency is to develop features from which the original time-series (from a given state) can be reconstructed to a suitable fidelity. If we were to build a classifier to distinguish human speech from other natural sounds, we could say that linear prediction coefficients are approximately sufficient for human speech because the speech could be intelligibly reconstructed from the features. The degree of sufficiency needed in any given problem depends on the application. The LPC coefficients would not be sufficient for a lie detector which needs more subtle information.

[0103] The relationship of the class-specific Baum-Welch algorithm to the conventional algorithm is described below. The conventional Baum-Welch algorithm computes $\begin{array}{cc}\underset{\lambda }{\max }p\left(S❘\lambda \right),& \left(40\right)\end{array}$

[0104] while the class-specific counterpart (see eq. 16) computes $\begin{array}{cc}\underset{{\lambda }^z}{\max }g\left(Z❘{\lambda }^z\right)p\left(S❘H_0\right)& \left(41\right)\end{array}$

[0105] where

[0106] Z₁,z₁. . . , z_N}  (42)

[0107] and

[0108] z_jj(s₁),T_j(s₂), . . . , T_j(s_T)_}.  (43)

[0109] In addition to the relationship between the parameter estimates of the two algorithms, the methods are also related as classifiers. If we had the true parameter values of λ, λ^z, so that equation (12) held, it is clear from equation (16) that

[0110] p(S|λ)=g(Z|λ^z)p(S|H₀).  (44)

[0111] Thus, given M classes and M sets of parameters λ_m, λ_m^z, m=1, . . . , M, the classifiers $\begin{array}{cc}\arg \underset{m=1}{\overset{M}{\max }}p\left(S❘{\lambda }_m\right)& \left(45\right)\\ \mathrm{and}& \\ \arg \underset{m=1}{\overset{M}{\max }}g\left(Z❘{\lambda }_m^z\right)p\left(S❘H_0\right)& \left(46\right)\end{array}$

[0112] give the same answer. Thus, S and Z are equivalent statistics for classification purposes or Z is sufficient for the problem. If we combine the class-specific features into a common feature vector Z by forming the union of the class-specific features,

[0113] Z=U{z₁, . . . ,z_n},  (47)

[0114] then we could use Z in place of s in a conventional HMM. Any common feature vector from which we can derive the class-specific features is sufficient because the class-specific features are sufficient. Thus, a side benefit of the class-specific method is the identification of a sufficient common feature set Z. If we consider the static M-ary classification problem, it is a well known result that the ratios $\begin{array}{cc}\frac{p\left(X❘H_1\right)}{p\left(X❘H_0\right)},\frac{p\left(X❘H_2\right)}{p\left(X❘H_0\right)},\dots ,\frac{p\left(X❘H_M\right)}{p\left(X❘H_0\right)}& \left(48\right)\end{array}$

[0115] are sufficient for the M-ary classification problem. The sufficiency property is used to reduce each of these ratios separately by an appropriate choice of sufficient statistics: $\begin{array}{cc}\frac{p\left(z_1❘H_1\right)}{p\left(z_1❘H_0\right)},\frac{p\left(z_2❘H_2\right)}{p\left(z_2❘H_0\right)},\dots ,\frac{p\left(z_M❘H_M\right)}{p\left(z_M❘H_0\right)}& \left(49\right)\end{array}$

[0116] Thus, clearly z₁, z₂, . . . , z_mtogether are sufficient for the M-ary problem.

[0117] To test the concept of the modified Baum-Welch algorithm (or Specific HMM), a Markov model was created. The following transition matrix and state priors were used: $\begin{array}{cc}A=\left[\begin{array}{cccccc}0.7& 0.3& 0& 0& 0& 0\\ 0& 0.7& 0.3& 0& 0& 0\\ 0& 0& 0.7& 0.3& 0& 0\\ 0& 0& 0& 0.7& 0.3& 0\\ 0& 0& 0& 0& 0.7& 0.3\\ 0.1& 0& 0& 0& 0& 0.9\end{array}\right]u=\left[\begin{array}{c}0.167\\ 0.167\\ 0.167\\ 0.167\\ 0.167\\ 0.167\end{array}\right]& \left(50\right)\end{array}$

[0118] For each state, there is a corresponding signal model. For each time step, a state is chosen randomly according to the Markov model. Then, a 256-sample segment of a time series is created from the statistical model corresponding to the chosen state. Features are then computed from each segment and used as an observation vector of an HMM. The following signal models were used:

[0119] An impulse of duration 2 samples occurring on samples 1 and 2 of the time-series with additive Gaussian noise of variance 1.

[0120] An impulse of duration 2 samples occurring on samples 2 and 3 of the time-series with additive Gaussian noise of variance 1.

[0121] Gaussian noise of variance greater than 1.

[0122] Sinewave of frequency 0.100 radians/sec in additive Gaussian noise of variance 1.

[0123] Sinewave of frequency 0.101 radians/sec in additive Gaussian noise of variance 1.

[0124] Second-order autoregressive noise of variance 1.

[0125] The criteria for selecting these synthetic signals were (a) the signals are easy to describe and produce, (b) a sufficient or approximately sufficient statistic was known for distinguishing each signal from H₀, (c) these statistics had a known density under H₀, (d) the signals and statistics were diverse and statistically dependent (otherwise a factorization of the PDF is possible). The common signal class, H₀, is iid Gaussian noise of zero mean and unit variance, denoted N (0,1). For each model, a sufficient statistic is chosen for distinguishing the signal from H₀. These are denoted z₁through z₆.

[0126] Each signal is described in detail below including the sufficient statistic for distinguishing the signal from H₀, and formulas for the PDF of each sufficient statistic under H₀. Signal classes 1 and 2: impulses in additive noise.

[0127] Signal class 1 follows the model

[0128] x_t=2(δ[t]+δ[t−1])+n_t,  (51)

[0129] where δ[t]=1 if t=0 and equals zero otherwise. The independent additive noise n_tis distributed N(0,1). A sufficient statistic for distinguishing this signal from H₀is

[0130] z₁=x₀+x₁.  (52)

[0131] The distribution of z₁under H₀is clearly normal with mean zero and variance 2:

[0132] log p(z₁|H₀)=−0.5 log (4π)−0.25z₁².  (53)

[0133] Signal class 2 follows the model

[0134] x_t=2(δ[t−1]+δ[t−2])+n_t.  (54)

[0135] A sufficient statistic for distinguishing this signal from H₀is

[0136] z₂=x₁+x₂.  (55)

[0137] The distribution of z₂under H₀is identical to the distribution of z₁under H₀.

[0138] Signal class 3: random noise with variance greater than 1.

[0139] Signal type 3 is iid Gaussian noise of zero mean and variance 1.7. A sufficient statistic for distinguishing this signal from H₀is $\begin{array}{cc}{z}_3=\log \left\{\sum _{t=1}^N{x}_t^2\right\}.& \left(56\right)\end{array}$

[0140] The distribution of z₃under H₀is Chi-squared with N degrees of freedom:

[0141] log p(z₃|H₀)=−log Λ(N/2)−(N/) log (2)+(N/2)z₃−exp(z₃)/2.  (57)

[0142] Signal classes 4 and 5: Sinewaves in additive noise.

[0143] Signal class 4 and 5 follow the model

[0144] x_t=n_t+asin (ω_kt+φ),  (58)

[0145] where φ is a uniformly distributed random phase and a=0.4 and ω_kdepends on the signal class k: ω₄=0.1, ω₅=0.101 radians per sample. The additive noise n_tis distributed N (0, 1). A sufficient statistic for distinguishing these signals from H₀is $\begin{array}{cc}{z}_k=\log \left\{{\left|\sum _{t=1}^N{x}_t{e}^{-{\mathrm{j\omega }}_{k}t}\right|}^2\right\}k=4,5& \left(59\right)\end{array}$

[0146] The distribution of z₄and z₅under H₀is the distribution of the log of a scaled χ²(2) random variable. More precisely, $\begin{array}{cc}\log p\left(z_k❘H_0\right)=-\log N-\frac{\exp \left\{z_k\right\}}{N}+z_{k}k=4,5.& \left(60\right)\end{array}$

[0147] Signal class 6: Second-order autoregressive noise of variance 1.

[0148] Signal type 6 is autoregressive noise of variance 1. We have

[0149] x_t=[−a₁x_t−1−a₂x_t−2+n_t]a  (61)

[0150] where n_tis distributed N(0,1). The autoregressive coefficients are fixed at values of a₁=−0.75, a₂=0.78. The theoretical output autocorrelation may be solved for analytically, and gives r₀=3.1048 α². Thus, to create a first-order AR process of variance 1, we require α=0.5675. Because the variance of the resulting process is 1, it is distinguished from H₀primarily by its first 2 autocorrelation lags. Thus, an approximate sufficient statistics to distinguish this signal type from H₀are the first 2 autocorrelation function estimates. The theoretical distribution of the normalized circular autocorrelation estimates {{tilde over (r)}_{1, {tilde over (r)}2}, . . . , {tilde over (r)}_p}, where P is the model order,

[0151] {tilde over (r)}_k={circumflex over (r)}_k/{circumflex over (r)}₀, k=1,2, . . . ,P,  (62)

[0152] and $\begin{array}{cc}{\hat{r}}_k=\frac{1}{N}\sum _{i=1}^N{x}_i{x}_{\left(i-k\right)\mathrm{mod}N},& \left(63\right)\end{array}$

[0153] are known. There is a 1:1 transformation that can be used to transform {{tilde over (r)}_k} into statistics that are better-behaved statistically (the distribution appears to closer to Gaussian and may be approximated with a Gaussian Mixture of low model order). Such a transformation consists of first converting {{tilde over (r)}_k}, k=1, . . . , {dot over (P)} into reflection coefficients {{circumflex over (K)}₁, {circumflex over (K)}₂, . . . , {circumflex over (K)}_P}, then applying the transformation $\begin{array}{cc}{K}_k^{\prime }=\log \left(\frac{1-K_k}{1+K_k}\right).& \left(64\right)\end{array}$

[0154] In summary, statistic z₆is produced as follows (with P=2):

[0155] 1. Compute {{circumflex over (r)}₀, {circumflex over (r)}₁, {circumflex over (r)}₂, . . . {circumflex over (r)}_P}.

[0156] 2. Normalize. Let {tilde over (r)}_k={circumflex over (r)}_k/{circumflex over (r)}₀, k=1,2, . . . , P.

[0157] 3. Convert {{tilde over (r)}_k} to reflection coefficients {K_k}. For P=1 this is trivial: K₁=−{tilde over (r)}₁.

[0158] 4. Let z₆={K′₁, K′₂, . . . , K′_P} where $K_k^{\prime }=\log \left(\frac{1-K_k}{1+K_k}\right).$

[0159] Under H₀, we may find the exact distribution of {{tilde over (r)}_k, k=1,2, . . . , P}, however this becomes numerically unstable for large N. A simpler asymptotic form is available. The general form for the asymptotic (large N) distribution of {K′₁, . . . , K′_P} is the following. Let $\begin{array}{cc}{K}_k=\left(\frac{1-\exp K_k^{\prime }}{1+\exp K_k^{\prime }}\right).& \left(65\right)\end{array}$

[0160] The log of the Jacobian of this transformation is $\begin{array}{cc}\log J=-\sum _{k=1}^P\log \left(\frac{2}{\left|K_k^2-1\right|}\right),& \left(66\right)\end{array}$

[0161] and the distribution of {K₁, K₂, . . . , K_P} under H₀is approximated for large N by the formula as follows:

[0162] Let $\begin{array}{cc}{c}_0=\log \Gamma \left(\frac{N}{2}\right)-\frac{\log \pi }{2}-\log \Gamma \left(\frac{N-1}{2}\right)& \left(67\right)\\ c_0=\log \Gamma \left(N+1\right)-N\log 2-\log \Gamma \left(\frac{N-1}{2}\right)-\log \Gamma \left(\frac{N+3}{2}\right)& \left(68\right)\end{array}$

[0163] Let n_ebe the largest integer less than or equal to P/2, and n₀be the smallest integer greater than or equal to P/2. Then $\begin{array}{cc}\log p\left(K_1,\dots ,K_p\right)=\sum _{i=1}^{n_e}\left\{c_e+2\log \left(1+K_{2i}\right)+\frac{N-3}{2}\log \left(1-K_{2i}^2\right)\right\}+=\sum _{i=1}^{n_o}\left\{c_o+\log \left(1+K_{2i-1}\right)+\frac{N-3}{2}\log \left(1-K_{2i-1}^2\right)\right\}.& \left(69\right)\end{array}$

[0164] Then,

[0165] log p(z₆|H₀)=log J+log p(K₁, . . . ,K_P).  (70)

[0166] To illustrate the strength of the Class-Specific algorithm, it was compared with two classical implementations. The classical implementations used a 7-dimensional feature vector consisting of all features used by the Class-Specific (CS) approach. Two variations of the classical approach were created by assuming the Gaussian Mixtures had (1) a general covariance matrix (“classical”) or CL approach) or (2) a diagonal structure which assumes all the features are statistically independent (“independence assumption” or IA approach). We will refer to these three methods by the abbreviations CS, CL, and IA.

[0167] The multiple record implementation of the HMM described in Section V.(B.) of the Rabiner paper (page 273) was used. All training and testing data was organized into records with a constant length of 99 data segments per record (T_r=99). Each segment consisted of 256 time samples. Features were computed on each segment (a segment is associated with a time-step or observation of the HMM).

[0168] In the experiment, algorithm performance was determined as a function of R. In the experiments, we used R=1,2,5,10,20,40,80, 160, and 320 records. To measure algorithm performance of all implementations, the appropriate version of the Baum-Welch algorithm was first run to convergence on the available training data. Next, the Viterbi algorithm was used to determine the most likely state sequence on a separate pool of 640 records of testing data. The Viterbi algorithm (described in the Rabiner paper on pp. 263-264) was appropriately modified for the CS method by substituting likelihood ratios in place of likelihoods, as described above. The number of state errors divided by the number of opportunities (640 times 99) was used as a performance metric. For more accurate determination of algorithm performance, sixteen independent trials were made at each value of R.

[0169] The details of the algorithm as applied in the tests are now discussed in greater detail. Algorithm scaling as detailed in the Rabiner paper (pp. 272-273), is needed in any practical HMM implementation, especially in the class-specific HMM. To prevent singular covariance matrices, a small value, ρ_ji², i=1,2, . . . , P_jshould be added to the diagonal elements of the matrix U_jm^z. In practice, it prevents collapsing covariance estimates due to insufficient data. Values of ρ_ji²are as follows: z₁, z₂: 1×10⁻⁴, z₂through z₆: 1×10⁻⁶.

[0170] One difficulty arises in trying to determine the number of state errors. If the algorithm is initialized randomly, there is no predefined mapping between the states and the models. The CS approach was always observed to converge in such a way that the state and model indexes coincide. This is because the states are treated specially with respect to the sufficient statistic used. But for the CL and IA implementations of the HMM, there would be no way to predict the mapping if the initialization was random. Therefore, it would be necessary to search for the best mapping on each trial by minimizing the error count over all mappings. This approach was adopted initially, but because of “catastrophic” errors, as discussed below, an “assisted” initialization method was used instead. This solved the unknown mapping problem.

[0171] Two initialization methods are used, “unassisted” for the CS approach and “assisted” for the CL and IS approaches. The “unassisted” method was to initialize with many mixture components with means located at randomly selected data points, and covariances equal to the measured sample covariance of the ensemble of training data. The HMM probabilities u and a_ijare initialized to uniform values of {fraction (1/N)}. The “assisted” method and the reasons for it are discussed below.

[0172] During training, components are removed (pruned) if the corresponding c_jm^zis “small”, indicating an insufficient number of samples “used” by the component. The number of mixture components remaining after convergence thus increased with the amount of training. However, components were also merged together if it was determined that they were too “close”. The measure of closeness was based on percentage of mutual overlap. As a result of “merging”, the number of components reached a certain limit and stopped increasing. A full discussion of the pruning and merging rules is not in the scope of this paper and is of little academic interest; however, it may be of interest to state the average numbers of mixing components remaining after convergence. This is provided in Table 1 as follows.

[0173] The CL and IA approaches encountered severe problems with initialization. More often than not, a poor initialization resulted in the Baum-Welch algorithm finding the incorrect stationary point. This resulted in catastrophic errors. This appeared to be caused by an inability to distinguish two or more states. The number of catastrophic errors decreased as the number of training records increased.

[0174] The CS method did not encounter any catastrophic errors like this. To study only the PDF estimation issue apart from any initialization issues, it was decided to “assist” the CL and IA approaches by providing a good initialization point. This initialization point was found by training separately on separate labeled data from each of the states, then using these PDF's as the state PDF's with the known u and A. In short, the algorithm was initialized with the true parameter values. This “good” parameter set was used as the starting point for all trials of the CL and IA approaches. Doing this provides somewhat of a “lower bound” on error performance. The CS results were unassisted.

[0175] A plot of the average error probability for the CL and CS methods as a function of R is shown in FIG. 4. The three lines represent maximum, minimum, and median P_eout of 16 trials at each R. From this graph, it is clear that the CL can achieve the same P_e, however requires far more data. Now consider FIG. 5 in which we add results of the IA method and show only the median error curves. The IA approach has a higher limiting error rate, due to the built-in assumption of independence, which is not valid. It has, however, nearly identical convergence behavior as the CS method. This is expected since the dimension of the PDF estimates is similar for the two approaches.

[0176] Accordingly, the present invention provides a method and system for training class-specific HMMs using a modified Baum-Welch algorithm. The class-specific implementation of the Baum-Welch algorithm has been developed and verified to work in principle. A computer simulation used a 6-state HMM with six synthetic signal models. The class-specific feature dimensions were either 1 or 2. When the performance was measured in terms of state errors in the most likely state sequence, the class-specific method greatly outperformed the classic HMM approach which operated on a combined 7-dimensional feature set. The improved performance was due to the lower dimension of the feature spaces made possible by knowledge of the appropriate feature set for each state. It also outperformed a method that assumed the features were statistically independent. In this case, the improved performance comes from the fact that the class-specific approach does not compromise theoretical performance when reducing dimension.

[0177] In light of the above, it is therefore understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described.