维特比译码算法

//实现（4，1，6）卷积码的维特比译码源程序，采用了最大似然算法
//介绍了软判决维特比译码算法过程的三个步骤：初始化，度量更新和回溯译码
#include<stdio.h>
#define DATA_LENGTH  4000000
#define TRACEBACK_LENGTH  48   //译码回溯深度，一般为m的5-10倍，m为寄存器个数
#define TB_DATA_OUT  16        //回溯输出译码长度，即每次回溯输出16bits译码信息

unsigned long Trans_table[2*TRACEBACK_LENGTH]={0};//一个单元32位，每两个存储单元存储一个时刻的64状态幸存信息
unsigned long Tran;  //存储每个时刻的分支选择路径信息值，32位
int bef[64]={0};    //前64状态
int aft[64]={0};    //后64状态

void BFLY_A(int dp, int dm, int *p1, int *p2, int *p3, int select);  //用于计算路径度量
void Tran_exchange();      //用于调整Tran中第2，3位的次序
void Trace_back(int *pout); //回溯
int Max_state(int *pbef);   //用于找到具有最大值所对应的状态点
int Min(int *pbef);         //求出64个状态的最小值

void Viterbi()
{
    //采用软判决，16电平量化，量化范围：-8到7，理想情况下1为-8，0为7
   
    int t0,t1,t2,t3;  //分支度量因子
    int n,i,j;
    int st,mini;
   
    int input[4]={0};  //keep input data
    int output[TB_DATA_OUT]={0};  //keep output data
    int *pb,*pa1,*pa2;  //pb指向前状态，pa1,pa2指向后状态
    int *p4;            //存储路径信息
   
    FILE *fp;
    FILE *fq;
   
    if((fp=fopen("quantizer_BPSK0dB.bin","r"))==NULL)
    {
        printf("cannot not open file!/n");   //file that provides input data
        exit(1);
    }   
   
    if((fq=fopen("viterbi_data_out.bin","wb"))==NULL)  //file that keep output data
    {
        printf("cannot not open file viterbi_data_out!/n");
        exit(1);
    }
   
    bef[0]=1000;  //从0状态开始，任意赋一个初值
   
    for(n=0;n<DATA_LENGTH/TB_DATA_OUT;n++)  //回溯次数，由处理数据长度定
    {
        if(n==0)
            p4=Trans_table;
        else if(n>=3)
            p4=Trans_table+2*(TRACEBACK_LENGTH-TB_DATA_OUT);
           
        for(j=0;j<TB_DATA_OUT;j++)  //每次循环对应一个时刻
        {
            pb=bef;
            pa1=aft;
            pa2=aft+2;
           
            fread(input,sizeof(int),4,fp);  //read input data
           
            //4个不同的分支度量因子，对应前4个蝶形
            //他们对应的因子为：0000，0110，0010，0100
            t0=input[0] + input[1] + input[2] + input[3];
            t1=input[0] - input[1] - input[2] + input[3];
            t0=input[0] + input[1] - input[2] + input[3];
            t0=input[0] - input[1] + input[2] + input[3];
//            printf("t0=%d,t1=%d,t2=%d,t3=%d/n",t0,t1,t2,t3);

      //每时刻64个状态，32个蝶形，两个蝶形为一组，共循环16次
      Tran=0;
      for(i=0;i<16;i++)
      {
          //16个蝶形的执行顺序是state:0-1-1-0--2-3-3-2--4-5-5-4--6-7-7-6
          if(i==0||i==3)
              st=0;
          else if(i==1||i==2)
              st=1;
          else if(i==4||i==7)
              st=2;
          else if(i==5||i==6)
              st=3;
          else if(i==8||i==11)
              st=4;
          else if(i==9||i==10)
              st=5;
          else if(i==12||i==15)
              st=6;
          else
              st=7;
             
          switch(st)
          {
              case 0:
                  BFLY_A(t0, t1, pb, pa1, pa2, 0);
                  pa1++;
                  pa2++;
                  BFLY_A(t0, t1, pb, pa1, pa2, 1);
                  //以下调整Tran中第2，3位的次序
                  Tran_exchange();
                  break;
              case 1:
                  BFLY_A(-t2, -t3, pb, pa1, pa2, 0);
                  pa1++;
                  pa2++;
                  BFLY_A(-t2, -t3, pb, pa1, pa2, 1);
                  //以下调整Tran中第2，3位的次序
                  Tran_exchange();
                  break;
              case 2:
                  BFLY_A(-t2, -t3, pb, pa1, pa2, 1);
                  pa1++;
                  pa2++;
                  BFLY_A(-t2, -t3, pb, pa1, pa2, 0);
                  //以下调整Tran中第2，3位的次序
                  Tran_exchange();
                  break;
              case 3:
                  BFLY_A(t0, t1, pb, pa1, pa2, 1);
                  pa1++;
                  pa2++;
                  BFLY_A(t0, t1, pb, pa1, pa2, 0);
                  //以下调整Tran中第2，3位的次序
                  Tran_exchange();
                  break;
              case 4:
                  if(i==8)
                  {
                      *(p4++)=Tran;  //前32位存储一次
                      Tran=0;
                  }
                  BFLY_A(t1, t0, pb, pa1, pa2, 1);
                  pa1++;
                  pa2++;
                  BFLY_A(t1, t0, pb, pa1, pa2, 0);
                  Tran_exchange();
                  break;
              case 5:
                  BFLY_A(-t3, -t2, pb, pa1, pa2, 1);
                  pa1++;
                  pa2++;
                  BFLY_A(-t3, -t2, pb, pa1, pa2, 0);
                  Tran_exchange();
                  break;
              case 6:
                  BFLY_A(-t3, -t2, pb, pa1, pa2, 0);
                  pa1++;
                  pa2++;
                  BFLY_A(-t3, -t2, pb, pa1, pa2, 1);
                  Tran_exchange();
                  break;
              case 7:
                  BFLY_A(t1, t0, pb, pa1, pa2, 0);
                  pa1++;
                  pa2++;
                  BFLY_A(t1, t0, pb, pa1, pa2, 1);
                  Tran_exchange();
                  break;
              default:
                  printf("state error!/n");
                  break;
          }
          pa1=pa1+3;  //进入下一组蝶形运算
          pa2=pa2+3;
          pb=pb+2;            
      }//end of i
      *(p4++)=Tran;  //存储分支路径，后32位存储一次
     
      for(i=0;i<64;i++)
           bef[i]=aft[i];
           
        }//end of j
       
        //48时刻后，每16时刻回溯一次，输出16bits译码信息
        if(n>=2)
        {
            Trace_back(output);
           
            //将输出写到文件中
            fwrite(output,sizeof(int),TB_DATA_OUT,fq);
           
            //为了防止溢出，每个状态减去最小值
            mini=Min(bef);
            for(i=0;i<64;i++)
            {
                bef[i]=bef[i]-mini;
            }
           
            //Trans_table信息移位
            for(i=0;i<2*(TRACEBACK_LENGTH-TB_DATA_OUT);i++)
            {
               Trans_table[i]=Trans_table[i+32];
            }
        }
       
        //清零输出
        for(j=0;j<TB_DATA_OUT;j++)
            {output[j]=0;}
    }//end of n
    fclose(fq);
    fclose(fp);
}

void BFLY_A(int dp, int dm, int *p1, int *p2, int *p3, int select)
{
    int a,b,c,d;
   
    if(select==0)
    {
        a=*p1+dp;
        b=*(p1+1)+dm;
        c=*(p1+32)-dp;
        d=*(p1+33)-dm;
    }
    else if(select==1)
    {
        a=*p1-dp;
        b=*(p1+1)-dm;
        c=*(p1+32)+dp;
        d=*(p1+33)+dm;
    }
    else
    {
        printf("Parameter input errror!/n");
        exit(1);
    }
   
    Tran=Tran<<1;
   
    *p2=a;
    *p3=b;
    if(c>a)  //c路径标记为1，a路径标记为0
    {
        Tran=Tran|0x01;
        *p2=c;
    }
   
    Tran=Tran<<1;
   
    if(d>b)   //d路径标记为1，b路径标记为0
    {
        Tran=Tran|0x01;
        *p3=d;
    }
}

void Trace_back(int *pout)
{
    int *p4,*pt;
    int state;  //标识状态
    int i,num;
    int tb_data[TRACEBACK_LENGTH]={0}; //回溯数据，其最高16位为回溯译码输出
   
    p4=Trans_table+2*TRACEBACK_LENGTH; //移到数组最后
    pt=tb_data+TRACEBACK_LENGTH;    //移到数组最后，从后往前填充输出数据
    state=Max_state(bef);    //找到此时刻包含最大值的状态
   
    for(i=0;i<TRACEBACK_LENGTH;i++)
    {
        p4=p4-2;  //两个单元为一个整体，前一个单元处理前32状态
        pt--;    //输出数据
       
        //下面分两部分操作，前一部分处理前32状态，后一部分处理后32状态
        if(state>=0&&state<32)
        {
            if(state&0x01)  //只有state为奇数时，输入才可能为1
                 *pt=1;
            else
                 *pt=0;
                 
            num=31-state;
            state=state>>1;   //回到前一时刻
            if((*p4>>num)&0x01)
                state=state+32;
        }
        else if(state>=32&&state<64)
        {
            p4++;
            if(state&0x01)   //只有state为奇数时，输入才可能为1
                 *pt=1;
          else
                 *pt=0;
                 
            num=63-state;
            state=state>>1;   //回到前一时刻
            if((*p4>>num)&0x01)
                state=state+32;   
            p4--;
        }
    }//end of i
   
    //输出最高位16bits译码信息
    pt=tb_data;
    for(i=0;i<TB_DATA_OUT;i++)
    {
        *(pout++)=*(pt++);
    }
}

void Tran_exchange()
{
    int x,y;
    x=((Tran&0x04)?1:0);  //取出Tran的第3位
    y=((Tran&0x02)?1:0);  //取出Tran的第2位
    if(x!=y)  //如果不等则交换
    {
        if(x>y)
            Tran=Tran-2;  //Tran的第2，3位互换
        else
            Tran=Tran+2;  //Tran的第2，3位互换
    }
}

int Max_state(int *pbef)  //包含最大值的状态
{
    int i;
    int a,max,max_state;
   
    max=*pbef;
    max_state=0;
    for(i=1;i<64;i++)
    {
        a=*(++pbef);
        if(a>max)
        {
            max=a;
            max_state=i;
        }
    }
    return max_state;
}

int Min(int *pbef)
{
    int i;
    int a;
    int mini;
   
    mini=*pbef;
    for(i=1;i<64;i++)
    {
        a=*(++pbef);
        if(a<mini)
            mini=a;
    }
    return mini;
}

void main()
{
    Vitervi();
}