studyviterbi.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;
      }
    }
  }
}