dhmm_lhs.cpp

语音识别配套的VQ及DHMM模型训练程序(C语言)

		{
			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 Predecessor\n");
		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_Pr\n");
		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;
}