//实现(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();
}
维特比译码算法
最新推荐文章于 2023-03-07 22:51:00 发布