sequence.cpp

HMM模型的C语言实现

/*
**	File:	sequence.cpp
**	功能：	序列的生成等函数
*/


#include"StdAfx.h"
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<string.h>
#include"nrutil.h"
#include"hmm.h"

/************************************
**函数名称：GenSequenceArray
**功能：生产随机状态
**参数：phmm,HMM结构指针
**	seed:随机数种子
**	T：观察序列长度
**	O:观察值序列
**	q:产生的序列***/

void GenSequenceArray(HMM*phmm,int seed,int T,int*O,int*q)//O  q 为动态分配数组
{
	int t=1;

	hmmsetseed(seed);

	q[1]=GenInitalState(phmm);
	O[1]=GenSymbol(phmm,q[1]);

	for(t=2;t<=T;t++)
	{
		q[t]=GenNextState(phmm,q[t-1]);
		O[t]=GenSymbol(phmm,q[t]);
	}
}

/***产生初始状态***/
int GenInitalState(HMM*phmm)
{
	double val,accum;
	int i,q_t;

	val=hmmgetrand();
	accum=0.0;
	q_t=phmm->N;
	for(i=1;i<=phmm->N;i++)
	{
		if(val<phmm->pi[i]+accum)
		{
			q_t=i;
			break;
		}
		else
		{
			accum+=phmm->pi[i];
		}
	}

	return q_t;
}

/***生产下一个状态***/
int GenNextState(HMM*phmm,int q_t)
{
	double val,accum;
	int j,q_next;

	val=hmmgetrand();
	accum=0.0;
	q_next=phmm->N;
	for(j=1;j<=phmm->N;j++)
	{
		if(val<phmm->A[q_t][j]+accum)
		{
			q_next=j;
			break;
		}
		else
			accum+=phmm->A[q_t][j];
	}

	return q_next;
}

/***生成观察值***/
int GenSymbol(HMM*phmm,int q_t)
{
	double val,accum;
	int j,o_t;

	val=hmmgetrand();
    accum=0.0;
	o_t=phmm->M;
	for(j=1;j<=phmm->M;j++)
	{
		if(val<phmm->B[q_t][j]+accum)
		{
			o_t=j;
			break;
		}
		else
			accum+=phmm->B[q_t][j];
	}

	return o_t;
}

/***读取序列***/
void ReadSequence(FILE*fp,int*pT,int**pO)
{
	int*O;
	int i;

	fscanf(fp,"T=%d\n",pT);
	O=ivector(1,*pT);
	for(i=1;i<=*pT;i++)
		fscanf(fp,"%d",&O[i]);
	*pO=O;
}

/***保存序列***/
void PrintSequence(FILE*fp,int T,int *O)
{
	int i;

	fprintf(fp,"T=%d\n",T);
	for(i=1;i<=T;i++)
		fprintf(fp,"%d",O[i]);
	printf("\n");
}