hmmutils.c

matlab 例程

/************************
功能：HMM文件的读写等操作
***********************/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"nrutil.h"
#include"hmm.h"
#define PI 3.14
/************************
函数名称：ReadHMM
功能：读取HMM结构
参数：fp:文件指针
      phmm:HMM结构指针，保存HMM结构
返回值：无
************************/
void ReadHMM(FILE *fp,struct HMM *phmm)
{
	int i,j,k;
	fscanf(fp,"M=%d\n",&(phmm->M));
    fscanf(fp,"N=%d\n",&(phmm->N));
	fscanf(fp,"A:\n");
	phmm->A=(double **)dmatrix(1,phmm->N,1,phmm->N);
	for(i=1;i<=phmm->N;i++){
		for(j=1;j<=phmm->N;j++){
			fscanf(fp,"%lf",&(phmm->A[i][j]));
		}
		fscanf(fp,"\n");
	}
	fscanf(fp,"B:\n");
	phmm->B=(double **)dmatrix(1,phmm->N,1,phmm->M);
    for(j=1;j<=phmm->N;j++){
		for(k=1;k<=phmm->M;k++){
			fscanf(fp,"%lf",&(phmm->B[j][k]));
		}
		fscanf(fp,"\n");
	}
    fscanf(fp,"pi:\n");
	phmm->pi=(double *)dvector(1,phmm->N);
	for(i=1;i<=phmm->N;i++)
		fscanf(fp,"%lf",&(phmm->pi[i]));
}
/**********************************
函数：FreeHMM
功能：释放HMM结构
参数：phmm:HMM结构指针
返回值：无
*********************************/
void FreeHMM(struct HMM *phmm)
{
	free_dmatrix(phmm->A,1,phmm->N,1,phmm->N);
	free_dmatrix(phmm->B,1,phmm->N,1,phmm->M);
	free_dvector(phmm->pi,1,phmm->N);
}
/************************************
函数：InitHMM
功能：初始化HMM结构
参数：phmm:HMM结构指针
      N:状态数
	  M：可观察值的个数
	  seed:随机数种子
返回值：无
*********************/
void InitHMM(struct HMM *phmm,int N,int M)
{
	int i,k;
	

    phmm->N=N;
    phmm->M=M;
    
    phmm->A=(double **)dmatrix(1,phmm->N,1,phmm->N);

	phmm->A[1][1]=0.95;
    phmm->A[1][2]=0.05;
    phmm->A[2][1]=0.05;
    phmm->A[2][2]=0.95;
	
	phmm->B=(long double **)dmatrix(1,phmm->N,1,phmm->M);
    phmm->u=(double *)dvector(1,phmm->N);
    phmm->sigma=(double *)dvector(1,phmm->N);
	for(i=1;i<=phmm->N;i++)
	{
     phmm->u[i]=0;
	}
    phmm->sigma[2]=50.0;
    phmm->sigma[1]=5.0;
    for(k=0;k<phmm->M;k++)  
	{
		phmm->B[1][k]=(long double)exp((double) (-((k-phmm->u[1])*(k-phmm->u[1]))/(2*phmm->sigma[1]*phmm->sigma[1])));
        phmm->B[1][k] = phmm->B[1][k]/sqrt((double)(PI*2.0*phmm->sigma[1]*phmm->sigma[1]));
	}
	for(k=0;k<phmm->M;k++)  
	{
		phmm->B[2][k]=(long double)exp((double) (-((k-phmm->u[2])*(k-phmm->u[2]))/(2*phmm->sigma[2]*phmm->sigma[2])));
        phmm->B[2][k] = phmm->B[2][k]/sqrt((double)(PI*2.0*phmm->sigma[2]*phmm->sigma[2]));
	}
		
	phmm->pi=(double *)dvector(1,phmm->N);
	
	for(i=1;i<=phmm->N;i++)
	{
		phmm->pi[i]=0.5;
		
	}
     	
}

/**********************************
函数拷贝：CopyHMM
功能：拷贝HMM结构
参数：phmm1:HMM结构指针
      phmm2:HMM结构指针
返回：无
**********************************/
void CopyHMM(struct HMM *phmm1,struct HMM *phmm2)
{
  int i,j,k;
  phmm2->M=phmm1->M;
  phmm2->N=phmm2->N;
  phmm2->A=(double **)dmatrix(1,phmm2->N,1,phmm2->N);
  for(i=1;i<=phmm2->N;i++)
	  for(j=1;j<=phmm2->N;j++)
		  phmm2->A[i][j]=phmm1->A[i][j];
  phmm2->B=(double **)dmatrix(1,phmm2->N,1,phmm2->M);
  for(j=1;j<=phmm2->N;j++)
	  for(k=1;k<=phmm2->M;k++)
		  phmm2->B[j][k]=phmm1->B[j][k];
  phmm2->pi=(double *)dvector(1,phmm2->N);
  for(i=1;i<=phmm2->N;i++)
	  phmm2->pi[i]=phmm1->pi[i];
}
/*********************************
函数：PrintHMM
功能：保存HMM结构
参数：fp:文件指针
      phmm:HMM结构指针
返回值：无
**************************/
void PrintHMM(FILE *fp,struct HMM *phmm)
{
	int i,j,k;
	fprintf(fp,"M=%d\n",phmm->M);
    fprintf(fp,"N=%d\n",phmm->N);
	fprintf(fp,"A:\n");
	for(i=1;i<=phmm->N;i++)
	{
		for(j=1;j<=phmm->N;j++)
		{
			fprintf(fp,"%f",phmm->A[i][j]);
		}
        fprintf(fp,"\n");
	}
	fprintf(fp,"B:\n");
	for(j=1;j<=phmm->N;j++)
	{
		for(k=1;k<=phmm->M;k++)
		{
			fprintf(fp,"%f",phmm->B[j][k]);
		}
        fprintf(fp,"\n");
	}
     	fprintf(fp,"pi:\n");
	for(i=1;i<=phmm->N;i++)
	{
		
        fprintf(fp,"%f",phmm->pi[i]);
	}
    fprintf(fp,"\n\n");
}
/***************************
函数：forward.c
功能：给定观察值序列和HMM模型，利用前向算法求概率
参数：phmm:指向HMM的指针
      T：观察值序列长度
	  O：观察值序列
	  alpha:运算中的临时数组
	  pprob:返回值，要求的概率
****************************/
void Forward(struct HMM *phmm,int T,int **O,double **alpha,double *pprob,int k)
{
	int i,j;
	int t;

	double sum;
	/*初始化*/
	for(i=1;i<=phmm->N;i++)
		alpha[1][i]=phmm->pi[i]*phmm->B[i][O[1][k]];
	/*递归*/
	for(t=1;t<T;t++){
		for(j=1;j<=phmm->N;j++){
			sum=0.0;
			for(i=1;i<=phmm->N;i++)
				sum+=alpha[t][i]*(phmm->A[i][j]);
			    
			    alpha[t+1][j]=sum*(phmm->B[j][O[t+1][k]]);
		}
	}
	/*终止*/
	*pprob=0.0;
	for(i=1;i<=phmm->N;i++)
		*pprob+=alpha[T][i];
}
/**********************
功能：前向算法估计参数
参数：phmm:指向HMM的指针
      T：观察值序列长度
	  O：观察值序列
	  alpha:运算中的临时数组
	  scale:比例因子
	  pprob:返回值，要求的概率
	  ****************************/
void ForwardWithScale(struct HMM *phmm,int T,int **O,double **alpha,double *scale,double *pprob,int k)
/*pprob是对数概率*/
{int i,j;
int t;

double sum;
/*初始化*/
scale[1]=0.0;
for(i=1;i<=phmm->N;i++){
	alpha[1][i]=phmm->pi[i]*phmm->B[i][O[1][k]];
	scale[1]+=alpha[1][i];
}
for(i=1;i<=phmm->N;i++)
  alpha[1][i]=(alpha[1][i])/(scale[1]);
/*递归*/
for(t=1;t<T;t++){
	scale[t+1]=0.0;
		for(j=1;j<=phmm->N;j++){
			sum=0.0;
			for(i=1;i<=phmm->N;i++)
				sum+=alpha[t][i]*(phmm->A[i][j]);
			    
			    alpha[t+1][j]=sum*(phmm->B[j][O[t+1][k]]);
		        scale[t+1]+=alpha[t+1][j];
		}
		for(j=1;j<=phmm->N;j++)
			alpha[t+1][j]/=scale[t+1];
	}
	/*终止*/
	*pprob=0.0;
	for(t=1;t<=T;t++)
		*pprob+=log(scale[t]);
}
/***************************
函数：Backward.c
功能：给定观察值序列和HMM模型，利用后向算法求概率
参数：phmm:指向HMM的指针
      T：观察值序列长度
	  O：观察值序列
	  beta:运算中的临时数组
	  pprob:返回值，要求的概率
****************************/
void Backward(struct HMM *phmm,int T,int **O,double **beta,double *pprob,int k)
{
	int i,j;
	int t;

	double sum;
	/*初始化*/
	for(i=1;i<=phmm->N;i++)
		beta[T][i]=1.0;
	/*递归*/
	for(t=T-1;t>=1;t--){
		for(i=1;i<=phmm->N;i++){
			sum=0.0;
			for(j=1;j<=phmm->N;j++)
				sum+=phmm->A[i][j]*(phmm->B[j][O[t+1][k]])*beta[t+1][j];
			    beta[t][i]=sum;
		}
	}
	/*终止*/
	*pprob=0.0;
	for(i=1;i<=phmm->N;i++)
		*pprob+=beta[1][i];
}
/***************************
函数：BackwardWithScale.c
功能：给定观察值序列和HMM模型，利用后向算法求概率
参数：phmm:指向HMM的指针
      T：观察值序列长度
	  O：观察值序列
	  beta:运算中的临时数组
	  scale:比例因子数组
	  pprob:返回值，要求的概率
****************************/
void BackwardWithScale(struct HMM *phmm,int T,int **O,double **beta,double *scale,double *pprob,int k)
{
	int i,j;
	int t;

	double sum;
	/*初始化*/
	for(i=1;i<=phmm->N;i++)
		beta[T][i]=1.0/scale[T];
	/*递归*/
	for(t=T-1;t>=1;t--){
		for(i=1;i<=phmm->N;i++){
			sum=0.0;
			for(j=1;j<=phmm->N;j++)
				sum+=phmm->A[i][j]*(phmm->B[j][O[t+1][k]])*beta[t+1][j];
			    beta[t][i]=sum/scale[t];
		}
	}
	/*终止*/
	*pprob=0.0;
	for(i=1;i<=phmm->N;i++)
		*pprob+=beta[1][i];
}
/*********************************
函数：Viterbi.c
功能：给定HMM和观察值序列，求最可能的状态
参数：phmm:HMM结构指针
      T:观察值的个数
	  O：观察值序列
	  deta,psi:中间变量
	  q:求的最佳状态序列
	  pprob:概率
*********************************/
#define VITHUGE 100000000000.0

void Viterbi(struct HMM *phmm,int T,int **O,double **delta,int **psi,int *q,double *pprob,int k)
{
	int i,j;
	int t;
	int maxvalind;
	double maxval,val;


	/*初始化*/
	for(i=1;i<=phmm->N;i++)
	{
		delta[1][i]=phmm->pi[i]*(phmm->B[i][O[1][k]]);
		psi[1][i]=0;
	}
	/*递归*/
	for(t=2;t<=T;t++)
	{
		for(j=1;j<=phmm->N;j++)
		{
			maxval=0.0;
			maxvalind=1;
			for(i=1;i<=phmm->N;i++)
			{
				val=delta[t-1][i]*(phmm->A[i][j]);
				if(val>maxval){
					maxval=val;
					maxvalind=i;
				}
			}
			delta[t][j]=maxval*(phmm->B[j][O[t][k]]);
			psi[t][j]=maxvalind;
		}
	}
/*终止*/
	*pprob=0.0;
	q[T]=1;
	for(i=1;i<=phmm->N;i++)
	{
		if(delta[T][i]>*pprob)
		{
			*pprob=delta[T][i];
			q[T]=i;
		}
	}