12.cpp

viterbi算法的c实现

///viterbi decoding///

#include "iostream.h"
#include "stddef.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"

#define SIZE 254

struct strtrellis
{int nextstate;               //下一状态
 int out;                     //校验位 
 int curstate;                //当前状态
 int info;                    //信息位
 };

struct strtrellis trellis[16][2];
int code[SIZE];
int receivecode1[SIZE];
int receivecode2[SIZE];
int coded[SIZE];
int finalcode[SIZE];
int outcode[16][SIZE];
double matric[SIZE][16];
double transcode1[SIZE];
double transcode2[SIZE];

int punc[94]={0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0, //23
	       1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0, // 22
		   0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0, // 20
		   1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,  //21
		   0,0,0,0,0,0,0,0};//8

double noise()   
{ 
	int i;
	double x;
    x=0;
	for(i=0;i<12;i++)
		x+=(double)rand()/RAND_MAX;
	x=x-6.0;
	return x;
}

void createtrellis()
{
	int state;
    int q,i,remain,s;
    int vector[4];
    for(i=0;i<2;i++)  //输入为0，1两种情况
	{
 	   for(state=0;state<16;state++) //共16个状态 m为4？？
	   {
	       q=state;
		   for(s=0;s<4;s++)  //将状态存入寄存器
		   {
			  remain=q%2;
			  q=(int)q/2;
			  vector[s]=remain;
		   }
 		   trellis[state][i].curstate = state;  //当前状态
 		   trellis[state][i].info = i;   //输入即信息位
 		   q=vector[0];
 		   vector[0]=vector[1];
 		   vector[1]=vector[2];
 		   vector[2]=vector[3];
 		   vector[3]=(q + i + vector[0])%2;
 		   s = vector[3]*8 + vector[2]*4 + vector[1]*2 + vector[0];  //下一状态
 	       trellis[state][i].nextstate = s;  
 		   trellis[state][i].out = (q + vector[3]+vector[2]+vector[0])%2;  //输出校验位
	   }
	}
}

void encoder()
{
	int i,p;
    int vector[4];
    for(i=0;i<4;i++)  //初始化状态寄存器，均置0
	   vector[i]=0;
    for(i=0;i<SIZE-4;i++)
	{
 	   code[i]=rand()%2; //随机输入码字
	   p=vector[0];
	   vector[0]=vector[1];
	   vector[1]=vector[2];
	   vector[2]=vector[3];
       vector[3]=(code[i]+vector[0]+p)%2;
	   coded[i]=(vector[3]+vector[2]+vector[0]+p)%2; //输出校验位
	}
 //termination
 for(i=SIZE-4;i<SIZE;i++)  //输出最后四段校验位
 {
   p=vector[0];
   vector[0]=vector[1];
   vector[1]=vector[2];
   vector[2]=vector[3];

   vector[3]=0;  //无输入时
   code[i]=(p+vector[0])%2;
   coded[i]=(vector[3]+vector[2]+vector[0]+p)%2;
 }
}

void channel(double squre)
{int i;
 double noise();
 for(i=0;i<SIZE;i++)
 {
	transcode1[i]=1-2*code[i]+sqrt(squre)*noise(); //信息位
	transcode2[i]=1-2*coded[i]+sqrt(squre)*noise(); //校验位
   }
   
  /* for(i=0;i<94;i++)
 {
  if(punc[i]==0)
		  {  
		   transcode2[160+i]=sqrt(squre)*noise();
		  }
       	
  }*/
 

}

void harddecoder()
{
	int i;
    for(i=0;i<SIZE;i++)
	{
		if(transcode1[i]>0)
		    receivecode1[i]=0;
	    else
		    receivecode1[i]=1;
	    if(transcode2[i]>0)
		    receivecode2[i]=0;
	    else
		  receivecode2[i]=1;
    }

/*for(i=0;i<94;i++)
 {
  if(punc[i]==0)
		  {  
		   receivecode2[160+i]=rand()%2;
		  }
       	
  }*/

}

void viterbidecoder()   //判决方法
{
	int i,j;
    int state,m0,m1,s0,s1,s[16];
    double m[16][5];
    int statenum;
    int curstate,curmatric;
    int mcode[16][4];

//initialize
 for(i=0;i<16;i++)
  	for(j=0;j<4;j++)
   	mcode[i][j]=0;
 for(i=1;i<16;i++)
	matric[0][i]=10000;    //？？？？？
 matric[0][0]=0;
 for(i=1;i<SIZE;i++)
	for(j=0;j<16;j++)
   	matric[i][j]=0;
 for (i=0;i<16;i++)
	for(j=0;j<SIZE;j++)
   	outcode[i][j]=0;

 for(i=0;i<SIZE-4;i++)  
 {
  	for(state=0;state<16;state++)    /* s0-s15每一个状态都有可能有输入0和输入1两条路径到达，比较度量值m0,m1大小，将大的舍去 */
	{
   	    m0=0; m1=0;           /*m0,m1纪录到达当前状态输入为0和为1时的d*/
   	    for(j=0;j<16;j++)     //当前状态的循环

		{                   /*  寻找下一状态为state当输入为0、1的当前状态  */
			if(trellis[j][0].nextstate==state)   //如果输入的序列是0的当前状态
				s0=j;
			if(trellis[j][1].nextstate==state)  //如果输入是1的当前状态
				s1=j;
		}
      
	    for(int k=0;k<94;k++)   
		{
			if(punc[k+160]!=0)
	  
			{
				if(trellis[s0][0].out!=receivecode2[i])    //判断当前状态时s0，输入是0，输出是0
        	       m0++;
                if(trellis[s0][0].info!=receivecode1[i])   //判断当前状态s0，输入是0，信息序列是否为0
                	m0++;
                if(trellis[s1][1].out!=receivecode2[i])    //判断当前状态s1时
                	m1++;
                if(trellis[s1][1].info!=receivecode1[i])
                	m1++;
			}
	        else
			{
				if(trellis[s0][0].info!=receivecode1[i])   //判断当前状态s0，输入是0，信息序列是否为0
                	m0++;
				if(trellis[s1][1].info!=receivecode1[i])
                   	m1++;
			}
		}  
	  
	    m0=m0+matric[i][s0];
        m1=m1+matric[i][s1];   //存储错误的比特数
        if(m0<m1)
		{
			matric[i+1][state]=m0;
            outcode[state][i]=0;    //判断输出在某一个状态输出时0还是1，并存储
		}
       else
	   {
      	   matric[i+1][state]=m1;
           outcode[state][i]=1;
	   } 
    }
 }

 statenum=1; 
 for(state=0;state<16;state++)
	for(j=0;j<5;j++)
		m[state][j]=0;
 for(i=SIZE-1;i>SIZE-5;i--)   /*归零前的四个输入*/
 {
	if(statenum==1)    /*statenum用来标记需要检查多少个前一状态，比如i = size -1 使归零，则 i = size -2有两个状态*/
	{
   	 statenum=statenum*2;
       for(state=0;state<16;state++)      /*最后一位输入1 0时归零*/
	   {
       	 if(trellis[state][0].nextstate==0) 
			 s[0]=state;
         if(trellis[state][1].nextstate==0) 
			 s[1]=state;
       }
       if(trellis[s[0]][0].info!=receivecode1[i])
       	m[s[0]][1]++;
       if(trellis[s[0]][0].out!=receivecode2[i])
       	m[s[0]][1]++;
       if(trellis[s[1]][1].info!=receivecode1[i])
       	m[s[1]][1]++;
       if(trellis[s[1]][1].out!=receivecode2[i])
       	m[s[1]][1]++;
       mcode[s[0]][SIZE-i-1]=0;    //如果已经到归零
       mcode[s[1]][SIZE-i-1]=1;
   }
   else
   {
   	for(j=0;j<statenum;j++)        /*第一次循环时statenum = 2对应 i = SIZE -2*/
	{
         for(state=0;state<16;state++)       //还没有到归零
		 {
       		if(trellis[state][0].nextstate==s[j]) /*寻找SIZE - 3的state使得下一状态为s[j]*/ 
				s0=state;                         /*第一次循环时的s[j]对应SIZE - 1的state*/
         	if(trellis[state][1].nextstate==s[j]) 
				s1=state;
         }

         m[s0][SIZE-i]=m[s[j]][SIZE-i-1];  /*第一次循环将size - 1的d作为size - 2 的初值d*/  
         if(trellis[s0][0].info!=receivecode1[i])   //m 存储的是最后几个归零状态的信息
         	m[s0][SIZE-i]++;
       	 if(trellis[s0][0].out!=receivecode2[i])
         	m[s0][SIZE-i]++;

         m[s1][SIZE-i]=m[s[j]][SIZE-i-1];
       	 if(trellis[s1][1].info!=receivecode1[i])
         	m[s1][SIZE-i]++;
       	 if(trellis[s1][1].out!=receivecode2[i])
         	m[s1][SIZE-i]++;

         s[j]=s0;
         s[j+statenum]=s1;
         mcode[s[j]][SIZE-i-1]=0;          //归零的几位的信息
         mcode[s[j+statenum]][SIZE-i-1]=1;
     }
      statenum=statenum*2;
   }
 }

 curstate=0;
 for(state=0;state<16;state++)
	matric[SIZE-1][state]=matric[SIZE-4][state]+m[state][4];
 curmatric=matric[SIZE-1][0];
 for(state=1;state<16;state++)
	if(curmatric>matric[SIZE-1][state])
	{
		curmatric=matric[SIZE-1][state];
        curstate=state;
    }
 s0=curstate;     //当前状态
 for(i=SIZE-4;i<SIZE;i++)
 {
	 finalcode[i]=mcode[s0][SIZE-i-1];    //最后四位的信息位
     s0=trellis[s0][mcode[s0][SIZE-i-1]].nextstate;    //对应的最后四位输入序列
 } 
 s0=curstate;
 for(i=SIZE-5;i>=0;i--)
 {
	finalcode[i]=outcode[s0][i];      //当前状态为s0的前输出校验序列
   for(state=0;state<16;state++)
   {
   	if(trellis[state][outcode[s0][i]].nextstate==s0)
	{
      	s0=state;
         break;
    }
   }
 }
}

void main(void)
{
	int i,sum,j,a,b;
    double ber;

    void encoder();
    void createtrellis();
    void viterbidecoder();
    void channel(double squre);
    void harddecoder();
    double snr;
    double qure;
 
 FILE *fp;

 
//punc[94]=(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0, 1,0,1,0,0,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,0, 0,0,0,0,0,0,0,0);
 
 if((fp=fopen("result.dat","w+"))==NULL)
 {
 	printf("can't open file\n");
   	return;
 }
 else 
 {
 	for(snr=0.0;snr<10.0;snr+=0.5)
	{    
		b=160;
 		sum=0;
   	    ber=0;
   	    j=1;
 		qure=(1.0*448) / (2*250*pow(10,0.1*snr));
   	    fprintf(fp,"SNR=%1.2f\n",snr);
   	    cout<<"SNR="<<snr<<endl;
   	    while((sum<500)&&(ber<1e8))
		{
 			createtrellis();
  			encoder();
         /* for(i=0;i<94;i++)
		  { 
			  if(punc[i]==0)
		  {  
		    coded[160+i]=0.5;
		  }
       	
		  }*/
		  
		
			channel(qure);
 			harddecoder();
 			viterbidecoder();
			for(i=0;i<SIZE-4;i++)
 			    if(code[i]!=finalcode[i]) sum++;  //发送的码子的对应信息位不等于译出的码子的对应信息位
      	    ber=(double)sum/((SIZE-4)*j);
      	    j++;
		}
   	    fprintf(fp,"the BER of the decoder is: %e\n",ber);
   	    cout<<"the BER of the decoder is:"<<ber<<endl;
	}
    fclose(fp);
 }
}