HMM forward algorithm - c source code

This is HMM forward algorithm source code.
I made this code today during about half the day.

I wouldn't introduce HMM and forward algorithm on this page. Because there are so many good material in the web.

The source code need Txt file included A, B, pi and observation values.
After load the txt file, the code return likelihood value.

Refer tp the real excute screen.

txt file include...
number of state, number of observation, number of observation sequence
state matrix by 1 line
observation matrix by 1 line
observation sequence

This picture describes the result of forward algorithm.

C source code






#include < stdio.h>


double HMM_forwardAlgorithm(double* &pi, double** &A, double ** &B, int* &obS, int sN, int oN, int obs_N);

void main()
{
 ///////////////////////////////////////////////
 //Input HMM ramda
 int stateN, observeN, ob_seq_N;
 double *pi;
 double **A;
 double **B;
 int *obS;
 ////////////////////////////////////////////////

 FILE * fp;
 fp = fopen("HMM_Ramda2.txt", "r");
 fscanf(fp,"%d %d %d",&stateN, &observeN, &ob_seq_N);
 /////////////////////////////////////////////////
 //alloc buffer
 pi = new double[stateN];
 A = new double*[stateN];
 B = new double*[stateN];
 obS = new int[ob_seq_N];

 for(int i=0; i< stateN; ++i)
 {
  A[i] = new double[stateN];
  B[i] = new double[observeN];
 }
 /////////////////////////////////////////////////
 //read Data 
 //pi
 for(int i=0; i< stateN; ++i)
  fscanf(fp, "%lf", &(pi[i]) );
 //A
 for(int i=0; i< stateN; ++i)
  for(int j=0; j< stateN; ++j)
   fscanf(fp, "%lf", &(A[i][j]) );
 //B
 for(int i=0; i< stateN; ++i)
  for(int j=0; j< observeN; ++j)
   fscanf(fp, "%lf", &(B[i][j]) );

 //observe sequence
 for(int i=0; i< ob_seq_N; ++i)
  fscanf(fp, "%d", &(obS[i]));

 ///////////////////////////////////////////////////
 HMM_forwardAlgorithm(pi, A, B, obS, stateN, observeN, ob_seq_N);
 

 ////////////////////////////////////////////////
 fclose(fp);

 for(int i=0; i< stateN; ++i)
 {
  delete[] A[i];
  delete[] B[i];
 }
 delete[] A;
 delete[] B;
 delete[] pi;
 delete[] obS; 

}



double HMM_forwardAlgorithm(double* &pi, double** &A, double ** &B, int* &obS, int sN, int oN, int obs_N)
{
 int i,j,t;
 //make forward matrix
 double** fMtx = new double *[sN];
 for(i=0; i< sN; ++i)
  fMtx[i] = new double[obs_N];

 //first step
 for(i=0; i< sN; ++i)
 {
  fMtx[i][0] = pi[i] * B[i][ obS[0]-1 ]; 
 }

 //routine
 double sum;
 for(t=1; t< obs_N; ++t)
 {  
  for(j=0; j< sN; ++j)
  {   
   sum=0;
   for(i=0; i< sN; ++i)
   {
    sum += (fMtx[i][t-1]*A[i][j]);    
   }
   fMtx[j][t] = sum * B[j][ obS[t]-1];
  }
 }


 //final
 sum=0; 
 for(j=0; j< sN; ++j)
 {
  sum += (fMtx[j][obs_N-1]);
 } 

 ////////////////////////////////////////////////////////
 //report
 printf("state N=%d / Observe N=%d\n\n", sN, oN);

 printf("Initial Matrix\n");
 for(int i=0; i< sN; ++i)
  printf("%lf ", pi[i] );
 printf("\n\n");

 printf("A matrix\n");
 for(int i=0; i< sN; ++i)
 {
  for(int j=0; j< sN; ++j)
  {
   printf("%lf ", A[i][j]);
  }
  printf("\n");
 }
 printf("\n");
 printf("B matrix\n");
 for(int i=0; i< sN; ++i)
 {
  for(int j=0; j< oN; ++j)
  {
   printf("%lf ", B[i][j]);
  }
  printf("\n");
 }
 printf("\n");
 printf("Observation sequence\n");
 for(int i=0; i< obs_N; ++i)
  printf("%d ", obS[i]);

 printf("\n\n");


 printf("** Trellis matrix **\n");
 for(j=0; j< sN; ++j)
 {
  for(t=0; t< obs_N; ++t)
  {
   printf("%lf  ", fMtx[j][t]);

   if(t!=obs_N-1)
    printf("->  ");
  }
  printf("\n");
 }


 printf("\nLikelihood = %lf \n", sum);

 for(i=0; i< sN; ++i)
  delete[] fMtx[i];
 delete[] fMtx;

 return sum;
}

Comments

Popular posts from this blog

OpenCV Stitching example (Stitcher class, Panorama)

(OpenCV Study) Background subtractor MOG, MOG2, GMG example source code (BackgroundSubtractorMOG, BackgroundSubtractorMOG2, BackgroundSubtractorGMG)

Example source code of extract HOG feature from images, save descriptor values to xml file, using opencv (using HOGDescriptor )

Real-time N camera stitching Class.

8 point algorithm (Matlab source code) / The method to get the Fundamental Matrix and the Essential matrix

Optical Flow sample source code using OpenCV

Video Stabilization example source code, (using cvFindHomography, cvWarpPerspective functions in openCV)

(OpenCV Study) calcOpticalFlowFarneback example source code ( dense optical flow )

yuv422(YUYV) to RGB and RGB to yuv422(YUYV), (Using OpenCV and TBB)

(C, C++) TinyXML , xml read & write