viterbi416_new.c

实现了（4

//实现（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();
}