viterbi.c

vitebi算法的程序

#include<stdio.h>
#include<math.h>
#include"mex.h"
#define Initial 1000
#define Pow2_Mem_Max 1024
#define Sur_Max 60
#define Pow2_Rate_Max 64
const double Gen_2[2]={369,491},Gen_3[3]={367,435,457};

void my_vit(double *In_Bit,double *Out_Sym,double *Gen,int Rate,int Mem_Len,int Sur_Len,int Fram_Len)
{
 int i,j,k,m,sl,Xor_Self,Tmpr,flag;
 static int Stat_Mac[Pow2_Mem_Max];/*[1<<Mem_Len]*/
 int Tmp_Bit[Sur_Max][Pow2_Mem_Max];/*[Sur_Len*(1<<Mem_Len)]*/
 double Min_Tmp;
 double Tmp[Pow2_Rate_Max];/*[(1<<Rate)]*/
 double Sum_Tmp1[Pow2_Mem_Max];/*[1<<Mem_Len]*/
 static double Sum_Tmp[Pow2_Mem_Max];/*[(1<<Mem_Len)]*/
/*Initial part*/
 for(i=0;i<(1<<Mem_Len);i++)
 {Stat_Mac[i]=0;
  for(j=0;j<Rate;j++)
   {k=Xor_Self=0;
    Tmpr = (i<<1)&(int)Gen[j];
    while(Tmpr>>k)
     Xor_Self ^=  (Tmpr>>k++)&1;  
     Stat_Mac[i] += Xor_Self << j;
   }     
 }

 /*Processing*/
 for(sl=0;sl<=Fram_Len/Sur_Len;sl++)
 {
   if((1+sl)*Sur_Len<=Fram_Len) m=Sur_Len;
   else {m=Fram_Len%Sur_Len;if(m==0) break;}
   for(i=0;i<m;i++)
    {for(j=0;j<(1<<Rate);j++)
      {Tmp[j]=0;
       for(k= 0;k<Rate;k++)
       Tmp[j] += ((((j>>k)&1)<<1)-1)*In_Bit[Rate*(i+sl*Sur_Len)+k];
      }
     for(k=0;k<(1<<Mem_Len);k++)
         Sum_Tmp1[k]=Sum_Tmp[k];
     for(k=0;k<((1<<Mem_Len)>>1);k++)
      for(j=0;j<2;j++)
      {  Tmp_Bit[i][(k<<1)+j] =(Sum_Tmp1[k]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k]])<
       (Sum_Tmp1[k+((1<<Mem_Len)>>1)]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k+((1<<Mem_Len)>>1)]])?0:1;
         Sum_Tmp[(k<<1)+j]  = (Tmp_Bit[i][(k<<1)+j]==0)?(Sum_Tmp1[k]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k]]):
       (Sum_Tmp1[k+((1<<Mem_Len)>>1)]+Tmp[((1<<Rate)-1)*j-((j<<1)-1)*Stat_Mac[k+((1<<Mem_Len)>>1)]]);
      }
    }
   Min_Tmp=Sum_Tmp[0];
   for(k=1;k<(1<<Mem_Len);k++)
    Min_Tmp=(Min_Tmp<=Sum_Tmp[k])?Min_Tmp:Sum_Tmp[k];
   for(k=0;k<(1<<Mem_Len);k++)
   if(Min_Tmp==Sum_Tmp[k])
    {flag=k;
     for(j=m-1;j>=0;j--)
      {flag+=Tmp_Bit[j][flag]*(1<<Mem_Len);
       if(sl*Sur_Len+j<Fram_Len-Mem_Len) Out_Sym[sl*Sur_Len+j] =(double)(flag&1);
       flag >>= 1;
      }
    }
  }
}    
void mexFunction(int nlhs,mxArray *plhs[], int nrhs,const mxArray *prhs[])
{ double *Input,*Output,*Gen,Tmp,Tmp1[6];
  double Mem_Len,Sur_Len;
  int k,m,n,Gm,Gn,Chg;
  Gen=Tmp1;  
  Input=mxGetPr(prhs[0]); /*(Length,Rate)*/
  m=mxGetM(prhs[0]);
  n=mxGetN(prhs[0]);
  
  if(m>n) {Chg=n;n=m;m=Chg;}  
  Mem_Len=mxGetScalar(prhs[1]);
  Sur_Len=mxGetScalar(prhs[2]);
  
  Gm=mxGetM(prhs[3]);
  Gn=mxGetN(prhs[3]);

  if((Gm*Gn==1))
    {Tmp=mxGetScalar(prhs[3]);
     if(Tmp==2.0) for(k=0;k<2;k++) Tmp1[k]=Gen_2[k];
     if(Tmp==3.0) for(k=0;k<3;k++) Tmp1[k]=Gen_3[k];
     Gn=(int)Tmp;
    }
  else Gen=mxGetPr(prhs[3]);     
  
  if(Gm>Gn) Gn=Gm; /*Gn is Rate*/
        
  if(nrhs<4) 
    mexErrMsgTxt("The number of input variable Must be at least four!\nYou shoud input in the following formation:\n    convol(In_Symbol,Memory_Length,Gen_Polynomial)\nGood Luck! More detail,see convol.c");
  
  plhs[0]= mxCreateDoubleMatrix(n/Gn-(int)Mem_Len,m,mxREAL);
 
  Output=mxGetPr(plhs[0]);
for(k=0;k<m;k++)
  my_vit(Input+k*n,Output+k*(n/Gn-(int)Mem_Len),Gen,Gn,(int)Mem_Len,(int)Sur_Len,n/Gn);
 }