DSP编程

开发平台：
Visual C++

DHMM_LHS.cpp：源码内容
							//	DHMM_LHS.cpp:
//		Implementation of the DHMM_LHS Module.
//		That is the transform of previous DHMM Code by LiHuSheng.
//
//	Created 2001/08, By DongMing, MDSR.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "DHMM_LHS.h"
#include "DHMM_GL.h"
#include <math.h>
//#define KMEAN_PRECISION		(1.0E-3)
//#define Epsilon_B			(1.0E-6)
//#define EPSILON_CALC		(1.0E-6)
//#define MAX_CALC_NUM		300
//////////////////////////////////////////////////////////////////////
//	Private Function Head
void DHMM_VQ(
			 int DMNU_Train_num,					// 训练集人数
			 double **DMNU_Codebook,				// VQ码本
			 int Codeword_num,				// 码字个数
			 double K_means_precision,		// K平均精度
			 int *DMNU_Train_Speaker,			// 训练集说话人标号
			 int Dimension,					// 特征维数
			 int DMNU_sex,						// 性别标志
			 int DMNU_num,						// 测试集的组号(num=0~Test_time)
			 int DMAD_n_Code_Word_Num,
			 DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Word,
			 DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Book
			 );
void DHMM_train(
				int Status_num,				// 状态数
				int Train_num,				// 训练集人数
				int *DMNU_Train_Speaker,			// 训练集说话人标号
				double **DMNU_Codebook,			// VQ码本
				double *Pi,					// 初始概率
				double **A,					// A矩阵
				double **B,					// B矩阵
				int Output_num,				// VQ码本大小（码字个数）
				int DMNU_Cepstrum_order,			// 特征维数
				int DMNU_digit,					// 词条标号
				FILE *DMNU_Fp1,					// 模型文件指针
				FILE *DMNU_Fp2,					// ？？？
				int DMNU_sex,					// sex在该函数中没有引用
				int DMNU_Mode,					// 程序中出现的Mode参数的取值包括：
											// WHOLE_WORD
											// WORD_HEAD1
											// WORD_HEAD2
											// WORD_TAIL
				WORD_SAMPLE * DMAD_pu_Word_Sample,
				DHMM_MODEL * DMAD_pu_DHMM_Model
				);
double DHMM_Calculate(
					  double ***Alpha,
					  double ***Belta,
					  double *Pi,
					  double **A,
					  double **B,
					  int **Out,			// 用于训练的各句子的观察序列（输出序列）
					  int Sample_num,		// 训练集人数
					  int Output_num,		// VQ码本大小（码字数，即观察值的个数）
					  int *T,				// 用于训练的各句子的帧数
					  int n,				// 状态数
					  double **P			// ？？？
					  );
double DHMM_Priority(double *Pi, double **A, double **B, int *Out, int T, int n);
double Viterbi(
			   int *VQ_Out,				// VQ后的一个识别样本
			   int Frame_num,				// 帧数
			   double **B,					// B矩阵
			   double **A,					// A矩阵
			   int Status_num,				// 状态数
			   int *Status_sequence		// 状态序列
			   );
//////////////////////////////////////////////////////////////////////
//	API functions
int DHMM_VQ_Train_Code_Book_LHS(DYNA_2DIM_DOUBLE_ARRAY d2dda_Code_Word, int n_Code_Word_Num, int n_Code_Word_Dim,
								DYNA_2DIM_DOUBLE_ARRAY d2dda_Initial_Code_Book, DYNA_2DIM_DOUBLE_ARRAY d2dda_Code_Book, int n_Code_Book_Size)
{
	DHMM_VQ(0, NULL, n_Code_Book_Size, KMEAN_PRECISION, 0, n_Code_Word_Dim, 0, 0,
		n_Code_Word_Num, d2dda_Code_Word, d2dda_Code_Book);
	return 0;
}
int DHMM_Model_Train_DHMM_Model_LHS(WORD_SAMPLE * pu_Word_Sample, int n_Word_Sample_Num,
									DHMM_MODEL * pu_DHMM_Model)
{
	DHMM_train(pu_DHMM_Model->n_State_Num, n_Word_Sample_Num, NULL, NULL,
		pu_DHMM_Model->pdPi, pu_DHMM_Model->d2dda_A, pu_DHMM_Model->d2dda_B, pu_DHMM_Model->n_Code_Book_Size, 0,
		0, NULL, NULL, 0, 0, pu_Word_Sample, pu_DHMM_Model);
	return 0;
}
int DHMM_Recog_Viterbi_LHS(DHMM_MODEL * pu_DHMM_Model,
						   WORD_SAMPLE * pu_Word_Sample,
						   double * pd_Max_Likelihood, int * pn_Status_Sequence)
{
	(*pd_Max_Likelihood) = Viterbi(pu_Word_Sample->pn_VQed_Feature_Sequence, pu_Word_Sample->n_Feature_Sequence_Len,
		pu_DHMM_Model->d2dda_B, pu_DHMM_Model->d2dda_A, pu_DHMM_Model->n_State_Num,
		pn_Status_Sequence);
	return 0;
}
//////////////////////////////////////////////////////////////////////
//	Original Code
//////////////////////////////////////////////////////////////////////
//	DMNU = DongMing doesn't use this when transforming.
//	DMAD = DongMing addes this when transforming.
static void DHMM_VQ(
			 int DMNU_Train_num,					// 训练集人数
			 double **DMNU_Codebook,				// VQ码本
			 int Codeword_num,				// 码字个数
			 double K_means_precision,		// K平均精度
			 int *DMNU_Train_Speaker,			// 训练集说话人标号
			 int Dimension,					// 特征维数
			 int DMNU_sex,						// 性别标志
			 int DMNU_num,						// 测试集的组号(num=0~Test_time)
			 int DMAD_n_Code_Word_Num,
			 DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Word,
			 DYNA_2DIM_DOUBLE_ARRAY DMAD_d2dda_Code_Book
			 )
{
	int    i, j, p, /* n, */ m /* q,t,h */;
//	int    Seg[3];					// 除了在Read_Feature中初始化，Seg在该函数中并没有被引用
	int    Nearer;
	int    *Counter;				// 长度：Codeword_num
	int    Total_num;				// 训练样本的帧数总和
//	int    Acoustic_Mode;
//	char   Feature_file[50];		// 一个训练样本(utterance)的特征数据文件名
	double **Codeword;				// size: Codeword_num*Dimension
	double **Sum;					// size: Codeword_num*Dimension
	double **s;						// size: Total_num*Dimension
//	double **Feature;				// size: Frame_num*Dimension （一个训练样本（utterance）的特征数据）
	double Min_Dist;
	double Dist;
	double New_D;
	double Old_D;
//	int    Frame_num;				// 一个训练样本（utterance）的帧数
//	char   Codefile_name[2][40] = {{"\work\dspv\comp"},{"\work\dspv\comp"}};
//	char   Savefile[50];
//	FILE * Fp;
	// 计算训练样本的帧数总和
/*
	Total_num = 0;
	for(n=0;n<Train_num;n++)
		for(i=0;i<MODEL_NUM;i++)
		{
			sprintf(Feature_file,"\work\dspv\mm%02d\mm%d.mfc",Train_Speaker[n],i);
			// 每个说话人一个目录，每个utterance一个文件
			Total_num += Get_Length(Feature_file);
		}
*/
	Total_num = DMAD_n_Code_Word_Num;
/*
	if((Codeword = new double *[Codeword_num]) == NULL)
	{
		printf("Not enough memory for Codeword !n");
		exit(-1);
	}
*/
	if((Sum = new double *[Codeword_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Sum !n");
		ASSERT(0);
	}
	for(i=0; i<Codeword_num; i++)
	{
/*
		if((Codeword[i] = new double [Dimension]) == NULL)
		{
			printf("Not enough memory for Codeword[%d] !n", i);
			exit(-1);
		}
*/
		if((Sum[i] = new double [Dimension]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Sum[%d] !n", i);
			ASSERT(0);
		}
	}
	Codeword = DMAD_d2dda_Code_Book;
		
	if((Counter  = new int [Codeword_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Counter !n");
		ASSERT(0);
	}
/*	
	if((s  = new double * [Total_num]) == NULL)
	{
		printf("Not enough memory for s !n");
		exit(-1);
	}
	for(i=0; i<Total_num; i++)
		if((s[i]  = new double [Dimension]) == NULL)
		{
			printf("Not enough memory for s !n");
			exit(-1);
		}
*/
	// 读取所有的训练样本，逐帧写入s中，形成训练向量集
/*
	q = 0;
	for(i=0;i<Train_num;i++)
		for(n=0;n<MODEL_NUM;n++)
		{
			sprintf(Feature_file,"\work\dspv\mm%02d\mm%d.mfc",Train_Speaker[i],n);
			Read_Feature(Feature_file,Dimension,WHOLE_WORD,Frame_num,Feature,Seg,Acoustic_Mode);
			for(t=0;t<Frame_num;t++)
				for(h=0;h<Dimension;h++)
					s[q+t][h] = Feature[t][h];
				q += Frame_num;
				Free_Feature(Feature,Frame_num);
		}
*/
	s = DMAD_d2dda_Code_Word;
	// k-means Clusteration
	
	for(p=0; p<Dimension; p++) Codeword[0][p] = 0;
	for(i=0; i<Total_num; i++)
		for(p=0; p<Dimension; p++)
			Codeword[0][p] += s[i][p];
	for(p=0; p<Dimension; p++) Codeword[0][p] /= Total_num;
	
//	printf("  Clustering...n");
	for(m=2; m<=Codeword_num; m*=2)
	{
		DEBUG_PRINTF("tM = %d.n", m);
		for(i=m/2-1; i>=0; i--)
			for(j=0; j<Dimension; j++)
			{
				Codeword[i*2+1][j] = Codeword[i][j] * 0.9;
				Codeword[i*2][j]   = Codeword[i][j] * 1.1;
			}
			
		Old_D = 1.0e20;
		New_D = 1.0e15;
		while(((Old_D - New_D) / New_D) > K_means_precision)
		{
			for(j=0; j<m; j++)
			{
				for(p=0; p<Dimension; p++) Sum[j][p] = 0;
				Counter[j] = 0;
			}
			
			Old_D = New_D;
			New_D = 0;
			
			for(i=0; i<Total_num; i++)
			{
				Dist = 0;
				for(p=0; p<Dimension; p++)
					Dist += (Codeword[0][p] - s[i][p]) * (Codeword[0][p] - s[i][p]);
				Min_Dist = Dist;
				Nearer = 0;
				for(j=1; j<m; j++)
				{
					Dist = 0;
					for(p=0; p<Dimension; p++)
						Dist += (Codeword[j][p] - s[i][p]) * (Codeword[j][p] - s[i][p]);
					if(Dist < Min_Dist)
					{
						Min_Dist = Dist;
						Nearer = j;
					}
				}	//	To find out the nearest center to s[i][].
				for(p=0; p<Dimension; p++) Sum[Nearer][p] += s[i][p];	//	To compute center.
				Counter[Nearer] ++;	//	Class Number.
				New_D += Min_Dist;
			}
			
			for(j=0; j<m; j++)
			{
				if(Counter[j] > 0)
					for(p=0; p<Dimension; p++)
						Codeword[j][p] = Sum[j][p] / Counter[j];
			}
		}
	}
	
	// 训练结束，保存码本到Codebook参数中
/*
	for(i=0;i<Codeword_num;i++)
		for(j=0;j<Dimension;j++)
			Codebook[i][j] = Codeword[i][j];
		
	// 写入到码本文件中
	sprintf(Savefile,"%s%d.cod",Codefile_name[sex],num);
	Fp = fopen(Savefile, "wb");
	if(Fp == NULL)
	{
		printf("File open error with %s !n", Codefile_name);
		exit(-1);
	}
	fwrite(&Codeword_num, sizeof(int), 1, Fp);
	fwrite(&Dimension, sizeof(int), 1, Fp);
	for(i=0; i<Codeword_num; i++)
		fwrite(Codeword[i], sizeof(double), Dimension, Fp);
	fclose(Fp);
*/
	// 释放内存
	for(i=0; i<Codeword_num; i++)
	{
//		delete Codeword[i];
		delete Sum[i];
	}
//	delete Codeword;
	delete Sum;
	delete Counter;
	
//	for(i=0; i<Total_num; i++) delete s[i];
//	delete s;
}
static void DHMM_train(
				int Status_num,				// 状态数
				int Train_num,				// 训练集人数
				int *DMNU_Train_Speaker,			// 训练集说话人标号
				double **DMNU_Codebook,			// VQ码本
				double *Pi,					// 初始概率
				double **A,					// A矩阵
				double **B,					// B矩阵
				int Output_num,				// VQ码本大小（码字个数）
				int DMNU_Cepstrum_order,			// 特征维数
				int DMNU_digit,					// 词条标号
				FILE *DMNU_Fp1,					// 模型文件指针
				FILE *DMNU_Fp2,					// ？？？
				int DMNU_sex,					// sex在该函数中没有引用
				int DMNU_Mode,					// 程序中出现的Mode参数的取值包括：
											// WHOLE_WORD
											// WORD_HEAD1
											// WORD_HEAD2
											// WORD_TAIL
				WORD_SAMPLE * DMAD_pu_Word_Sample,
				DHMM_MODEL * DMAD_pu_DHMM_Model
				)
{
	int     i, j, k,m;
//	int     Acoustic_Mode;
//	int     Seg[3];
	int     **Output;			// for each word item model being trained, the output sequence (observation sequence) of each utterance of this model. （每个元素保存该帧的码字标号，每一行是一个utterance，且各行不等长）
								// 1st dimension: index of the training utterance; 
								// 2nd dimension: index of frame
	int     *Frame_num;			// 各训练句子的帧数
	
	double  ** P;				// Status_num*Output_num矩阵（干什么用的？？？）
	double  **Old_A;
	double  **Old_B;
	double  ***Alpha;
	double  ***Belta;
//	double  **Feature;			// feature data of one utterance
	double  w, u;
//	char    File_Name[40];
	
//	time_t  t;
	
//	srand((unsigned) time(&t));
	
	// 分配内存
	if((Old_A = new double *[Status_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Old_A !n");
		ASSERT(0);
	}
	for(i=0; i<Status_num; i++)
	{
		if((Old_A[i] = new double [Status_num]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Old_A[%d] !n", i);
			ASSERT(0);
		}
	}
	
	if((Old_B = new double *[Status_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Old_B !n");
		ASSERT(0);
	}
	for(i=0; i<Status_num; i++)
	{
		if((Old_B[i] = new double [Output_num]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Old_B[%d] !n", i);
			ASSERT(0);
		}
	}
	if((P = new double*[Status_num]) == NULL)
	{
		DEBUG_PRINTF("Not eoungh memory for P[]n");
		ASSERT(0);
	}
	for(i=0;i<Status_num;i++)
		if((P[i] = new double[Output_num]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for P[%d]n",i);
			ASSERT(0);
		}
		
		
	// Initialize Alpha & Belta
	if((Alpha = new double **[Train_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Alpha !n");
		ASSERT(0);
	}
	
	if((Belta = new double **[Train_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Belta !n");
		ASSERT(0);
	}
	
	if((Output = new int *[Train_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Output !n");
		ASSERT(0);
	}
	
	if((Frame_num = new int [Train_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Frame_num !n");
		ASSERT(0);
	}
	
	// 初始化初始概率Pi:状态0概率为1，其余状态概率为0
	Pi[0] = 1;
	for(i=1; i<Status_num; i++) Pi[i] = 0;
	
	for(i=0; i<Status_num; i++)
	{
		// 初始化A矩阵
		for(j=0; j<Status_num; j++) A[i][j] = 0;
		if(i < Status_num-1)
		{
			A[i][i]   = 0.5;
			A[i][i+1] = 0.5;
		}
		else
		{
			A[i][i]   = 1.0;
		}
		// 归一化A矩阵
		w = 0;
		for(j=0; j<Status_num; j++) w = w + A[i][j];
		for(j=0; j<Status_num; j++) A[i][j] = A[i][j] / w;
		
		w = 0;
		// 初始化B矩阵
		for(j=0; j<Output_num; j++)
		{
			B[i][j] = 1.0 / Output_num * (0.9 + ((rand()%1000) / 5000.0));
			w = w + B[i][j];
		}
		// 归一化B矩阵
		for(j=0; j<Output_num; j++) B[i][j] = B[i][j] / w;
	}
	
	k = 0;
	for(m=0;m<Train_num;m++)
	{
		// 初始化第m个训练集说话人的特征数据文件名
//		sprintf(File_Name,"\work\dspv\mm%02d\mm%d.mfc",Train_Speaker[m],digit);
		// 读取该说话人的该词条的训练数据
//		Read_Feature(File_Name,Cepstrum_order,Mode,Frame_num[k],Feature,Seg,Acoustic_Mode);
		Frame_num[k] = DMAD_pu_Word_Sample[k].n_Feature_Sequence_Len;
		
		// 根据帧数为Alpha,Belta和Output分配内存
		// allocate memory for Alpha
		if((Alpha[k] = new double *[Frame_num[k]]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Alpha[%d] !n", k);
			ASSERT(0);
		}
		for(i=0; i<Frame_num[k]; i++)
			if((Alpha[k][i] = new double [Status_num]) == NULL)
			{
				DEBUG_PRINTF("Not enough memory for Alpha[%d][%d] !n", k, i);
				ASSERT(0);
			}
			
		// allocate memory for Belta
		if((Belta[k] = new double *[Frame_num[k]]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Belta[%d] !n", k);
			ASSERT(0);
		}
		for(i=0; i<Frame_num[k]; i++)
			if((Belta[k][i] = new double [Status_num]) == NULL)
			{
				DEBUG_PRINTF("Not enough memory for Belta[%d][%d] !n", k, i);
				ASSERT(0);
			}
			
		// allocate memory for Output
/*
		if((Output[k] = new int [Frame_num[k]]) == NULL)
		{
			printf("Not enough memory for Output[%d] !n", k);
			exit(-1);
		}
		
		// VQ
		for(i=0; i<Frame_num[k]; i++)
		{
			Output[k][i] = DHMM_Match(Feature[i], Codebook, Output_num, Cepstrum_order);
		}
		
		// 释放Feature的内存
		Free_Feature(Feature,Frame_num[k]);
*/		
		Output[k] = DMAD_pu_Word_Sample[k].pn_VQed_Feature_Sequence;
		k ++;
	}
	
	w = 1;
	m = 0;
	k = MAX_CALC_NUM;				// 最多训练300遍
	while(m < k)
	{
		u = DHMM_Calculate(
			Alpha,
			Belta,
			Pi,
			A,
			B,
			Output,
			Train_num,
			Output_num,
			Frame_num,
			Status_num,
			P
			);
		m ++;
		if((u > w) && (fabs((u-w)/w) < EPSILON_CALC)) break;
		else w = u;
	}
	DEBUG_PRINTF("tM = %d.n", m);
/*	
	fwrite(Pi, sizeof(double), Status_num, Fp1);
	for(i=0; i<Status_num;i++)
		fwrite(A[i], sizeof(double), Status_num, Fp1);
	for(i=0; i<Status_num; i++)
		fwrite(B[i], sizeof(double), Output_num, Fp1);
	if(Mode & WHOLE_WORD)
	{
		// P的各元素是B的各元素加1，保存这个干什么？？？
		for(i=0; i<Status_num; i++)
			fwrite(P[i], sizeof(double), Output_num, Fp2);
	}
*/	
	for(k=0; k<Train_num; k++)
	{
		for(i=0; i<Frame_num[k]; i++)
			delete Alpha[k][i];
		delete Alpha[k];
		
		for(i=0; i<Frame_num[k]; i++)
			delete Belta[k][i];
		delete Belta[k];
		
//		delete Output[k];
	}
	
	for(i=0; i<Status_num; i++) delete Old_A[i];
	delete Old_A;
	for(i=0; i<Status_num; i++) delete Old_B[i];
	delete Old_B;
	for(i=0;i<Status_num;i++)
		delete P[i];
	delete P;
	
	
	delete Alpha;
	delete Belta;
	
	delete Output;
	delete Frame_num;
	
}
static double DHMM_Calculate(
					  double ***Alpha,
					  double ***Belta,
					  double *Pi,
					  double **A,
					  double **B,
					  int **Out,			// 用于训练的各句子的观察序列（输出序列）
					  int Sample_num,		// 训练集人数
					  int Output_num,		// VQ码本大小（码字数，即观察值的个数）
					  int *T,				// 用于训练的各句子的帧数
					  int n,				// 状态数
					  double **P			// ？？？
					  )
{
	int    i, j, k, m;
	double w;
	double *R;
	double **New_A;
	double **New_B;
	double **Ct;					// to prevent underflow of Alpha and Belta, Ct is the scaling coefficient
	double Residue;
	
	// 分配内存
	// allocate memory for New_A
	if((New_A = new double *[n]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for New_A !n");
		ASSERT(0);
	}
	for(i=0; i<n; i++)
	{
		if((New_A[i] = new double [n]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for New_A[%d] !n", i);
			ASSERT(0);
		}
	}
	
	// allocate memory for New_B
	if((New_B = new double *[n]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for New_B !n");
		ASSERT(0);
	}
	for(i=0; i<n; i++)
	{
		if((New_B[i] = new double [Output_num]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for New_B[%d] !n", i);
			ASSERT(0);
		}
	}
	
	if((R = new double [n]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for R !n");
		ASSERT(0);
	}
	
	if((Ct = new double *[Sample_num]) == NULL)
	{
		DEBUG_PRINTF("Not enough memory for Ct !n");
		ASSERT(0);
	}
	for(i=0; i<Sample_num; i++)
		if((Ct[i] = new double [T[i]]) == NULL)
		{
			DEBUG_PRINTF("Not enough memory for Ct[%d] !n", i);
			ASSERT(0);
		}
		
	// Calculate Alpha & Belta.
	
	for(m=0; m<Sample_num; m++)
	{
		w = 0;
		for(i=0; i<n; i++)	//	n, State Number.
		{
			Alpha[m][0][i] = Pi[i] * B[i][ Out[m][0] ];
			w = w + Alpha[m][0][i];
		}
		Ct[m][0] = 1.0 / w;								// scaling coefficient
		for(i=0; i<n; i++)
			Alpha[m][0][i] = Alpha[m][0][i] * Ct[m][0];		// scaling
		
		// calculate Alpha
		for(k=1; k<T[m]; k++)
		{
			Ct[m][k] = 0;
			for(i=0; i<n; i++)
			{
				w = 0;
				for(j=0; j<n; j++)
					w = w + Alpha[m][k-1][j] * A[j][i];
				w = w * B[i][ Out[m][k] ];
				Alpha[m][k][i] = w;
				Ct[m][k] = Ct[m][k] + w;
			}
			Ct[m][k] = 1.0 / Ct[m][k];
			for(i=0; i<n; i++)
				Alpha[m][k][i] = Alpha[m][k][i] * Ct[m][k];
		}
		//在顾良的程序中，通过End_state_mode来区分最后一个状态是否为吸收态(ENDSTATE_ABSORB)，现在不再区分了，而是认为不是ENDSTATE_ABSORB.（参见顾良程序R_DHMM.CPP）
		for(i=0; i<n; i++) Belta[m][ T[m]-1 ][i] = Ct[m][ T[m]-1 ];
		
		// calculate Belta
		for(k=T[m]-2; k>=0; k--)
			for(i=0; i<n; i++)
			{
				w = 0;
				for(j=0; j<n; j++)
					w = w + A[i][j] * B[j][ Out[m][k+1] ] * Belta[m][k+1][j];
				Belta[m][k][i] = w * Ct[m][k];
			}
	}
	
	
	// Re-estimate B[i][j], i = 0 to n-1, j = 0 to m-1.
	for(i=0; i<n; i++)
		for(j=0; j<Output_num; j++)
			New_B[i][j] = 0;
		
	for(m=0; m<Sample_num; m++)
	{
		w = 0;
		for(i=0; i<n; i++)
			w += Alpha[m][T[m]-1][i];	//	w: 观察序列概率，即P(O)
		
		for(k=0; k<T[m]; k++)
			for(i=0; i<n; i++)	//	n, State Number.
			{
				New_B[i][Out[m][k]] += Alpha[m][k][i] * Belta[m][k][i] /
					Ct[m][k] / w;
			}
	}
	// 算这个P干什么用？？？
	for(i=0;i<n;i++)
		for(k=0;k<Output_num;k++)
			P[i][k] = 1 + New_B[i][k];
		
	for(i=0; i<n; i++)
	{
		w = 0;
		for(j=0; j<Output_num; j++) w += New_B[i][j];
		for(j=0; j<Output_num; j++) New_B[i][j] /= w;		// 归一化
		
		for(j=0; j<Output_num; j++)
			if(New_B[i][j] < Epsilon_B)						// Epsilon_B=1.0e-6
				New_B[i][j] = Epsilon_B;
			
		w = 0;
		for(j=0; j<Output_num; j++) w += New_B[i][j];
		for(j=0; j<Output_num; j++) New_B[i][j] /= w;		// 归一化
	}
	
	for(i=0; i<n; i++)									// 更新B矩阵
		for(j=0; j<Output_num; j++)
		{
			B[i][j] = New_B[i][j];
		}
		
		
	Residue = 0;
	for(m=0; m<Sample_num; m++)
		Residue += DHMM_Priority(Pi, A, B, Out[m], T[m], n);
	
	for(i=0; i<n; i++) delete New_A[i];
	delete New_A;
	
	for(i=0; i<n; i++) delete New_B[i];
	delete New_B;
	
	delete R;
	for(i=0; i<Sample_num; i++) delete Ct[i];
	delete Ct;
	
	return Residue;
}
//	Calculate P(O|lamuda)
static double DHMM_Priority(double *Pi, double **A, double **B, int *Out, int T, int n)
{
	int   i, j, k;
	double  w, p;
	double  *Alpha1;
	double  *Alpha2;
	
	if((Alpha1 = new double [n]) == NULL) { DEBUG_PRINTF("Not enough memory for !n"); ASSERT(0); }
	if((Alpha2 = new double [n]) == NULL) { DEBUG_PRINTF("Not enough memory for !n"); ASSERT(0); }
	
	p = 0;
	w = 0;
	for(i=0; i<n; i++)
	{
		Alpha1[i] = Pi[i] * B[i][ Out[0] ];
		w += Alpha1[i];
	}
	for(i=0; i<n; i++) Alpha1[i] /= w;
	p += log(w);
	
	for(k=1; k<T; k++)
	{
		for(i=0; i<n; i++)
		{
			w = 0;
			for(j=0; j<n; j++)
				w = w + Alpha1[j] * A[j][i];
			w = w * B[i][ Out[k] ];
			Alpha2[i] = w;
		}
		
		if(k < T-1)
		{
			w = 0;
			for(i=0; i<n; i++) w += Alpha2[i];
			for(i=0; i<n; i++) Alpha1[i] = Alpha2[i] / w;
			p += log(w);
		}
	}
	
	w = 0;
	for(i=0; i<n; i++) w += Alpha2[i];
	p += log(w);
	
	delete Alpha1;
	delete Alpha2;
	
	return p;
}
static double Viterbi(
			   int *VQ_Out,				// VQ后的一个识别样本
			   int Frame_num,				// 帧数
			   double **B,					// B矩阵
			   double **A,					// A矩阵
			   int Status_num,				// 状态数
			   int *Status_sequence		// 状态序列
			   )
{
	int **Predecessor;		// Status_num*Frame_num矩阵：Predecessor[j][i]中保存第i帧状态为j的情况下，第i-1帧的状态
	int  i, j, k;
	
	double  w1, w2,p;
	double  **Current_Pr;		// Status_num*Frame_num矩阵：即Viterbi搜索中的delta概率
	
	if((Predecessor = new int*[Status_num]) == NULL)
	{
		printf("Not enough memory for Predecessorn");
		exit(0);
	}
	for(i=0;i<Status_num;i++)
		if((Predecessor[i] = new int[Frame_num]) == NULL)
		{
			printf("Not enough memory for Predecessor[%d]n",i);
			exit(0);
		}
		
	if((Current_Pr = new double*[Status_num]) == NULL)
	{
		printf("Not enough memory for Current_Prn");
		exit(0);
	}
	for(i=0;i<Status_num;i++)
		if((Current_Pr[i] = new double[Frame_num]) == NULL)
		{
			printf("Not enough memory for Current_Pr[%d]n",i);
			exit(0);
		}
		
	/* Viterbi decoding */
	/* the first frame */
	// 初始化delta[0][0]：第一帧状态0的概率为1，其余状态的概率为0。（这里概率取了对数）
	Current_Pr[0][0] = 0;
	for(i=1;i<Status_num;i++)
		Current_Pr[i][0] = -1e50;
	
	/* forward */
	for(i=1;i<Frame_num;i++)
	{
		Predecessor[0][i] = 0;			// 状态0的前一状态只能是状态0
		Current_Pr[0][i] = log10(A[0][0] * B[0][VQ_Out[i]]) + Current_Pr[0][i-1];	// 计算delta[i][0]
		for(j=1;j<Status_num;j++)		// 当前状态是j
		{
			w1 = log10(A[j-1][j] * B[j][VQ_Out[i]]) + Current_Pr[j-1][i-1];	// 前一状态是j-1
			w2 = log10(A[j][j] * B[j][VQ_Out[i]]) + Current_Pr[j][i-1];		// 前一状态是j
			// 对于前一状态的两种可能：j-1和j，取delta[i][j]更大的（j为当前状态）
			if(w1 >= w2)
			{
				Predecessor[j][i] = j - 1;
				Current_Pr[j][i] = w1;
			}
			else
			{
				Predecessor[j][i] = j;
				Current_Pr[j][i] = w2;
			}
		}
	}
	
	/* Maximum Likelihood */
	// 计算最大输出概率
	p = -1e50;
	k = -1;
	for(i=0;i<Status_num;i++)
	{
		if(Current_Pr[i][Frame_num-1] > p)
		{
			p = Current_Pr[i][Frame_num-1];
			k = i;
		}
	}
	
	/*Backward for Status-Sequence */
	Status_sequence[Frame_num-1] = Status_num - 1; //k
	for(i=Frame_num-2;i>=0;i--)
		Status_sequence[i] = Predecessor[Status_sequence[i+1]][i+1];
	
	/* free the memory */
	for(i=0;i<Status_num;i++)
	{
		delete Current_Pr[i];
		delete Predecessor[i];
	}
	delete Current_Pr;
	delete Predecessor;
	
	return p;
}

						