背景
此为我的本科毕设卷积码的viterbi译码，是我开始学编程之后的第一个作品，由于年代久远，仅当记录
源码
编码
//encode.h
#include <stdlib.h>
#include <time.h>
#include <math.h>
#define PI 3.14
#define length 10
#define M 3
#define N 3
void gensignal(int* origin,int mode)//产生信号
{//生成随机二进制数作为编码前的信息序列
if(mode==0) for(int i=0;i<length;i++) origin[ i]=0;//生成一位信息位
if(mode==1) for(int i=0;i<length;i++) origin[ i]=1;//生成一位信息位
if(mode==2) for(int i=0;i<length;i++) origin[ i]=i&1;//生成一位信息位
if(mode==3)
{
for(int i=0;i<length;i++) origin[ i]=rand()&1;//生成一位信息位
}
}
void encode(int encodeori[length][N],int* origin,int* generate)//编码器
{//生成卷积码,码率为1/N,一个码元有N比特
int shifter;//模拟左移移位寄存器当前状态[b(M-1),b(M-2),b(M-3),...,b(1)],最大表示16位
int i,j,k,temp=shifter;
for(k=0;k<length;k++)
{//由此信息位生成N bit码元
temp=(temp<<1)^origin[k];//模拟移位寄存器左移
for(i=0;i<N;i++)
{//生成此卷积码码元的1bit
int sum=0,temp1=temp&generate[ i];
encodeori[k][ i]=0;
for(j=0;j<M;j++) sum^=(temp1>>j)&1;
//实现卷积码编码电路输入输出关系,乘矩阵的列后自模2加
encodeori[k][ i]=(encodeori[k][ i]<<1)^sum;//组合成卷积码码元
}
shifter=temp;
}
if((encodeori[0][0]!=encodeori[0][1])||(encodeori[0][1]!=encodeori[0][2])||(encodeori[0][0]!=encodeori[0][2])) cout<<"\\\\\\";
}
void PSK(int encodeori[length][N],int PSKmaker[length][N])//PSK调制(1-->-1;0-->1)
{
for(int i=0;i<length;i++)
for(int j=0;j<N;j++)
{
if(encodeori[ i][j]) PSKmaker[ i][j]=-1;
else PSKmaker[ i][j]=1;
}
}
double gaosimaker(double DX)
{
double rand1,rand2,result;
rand1=(double)(rand())/RAND_MAX;//均匀分布
rand2=(double)(rand())/RAND_MAX;//均匀分布
result=sqrt((-2)*log(rand1))*cos(2*PI*rand2);
//变换抽样,均匀分布变为正态分布
return DX*result;
}
void addgaosi(double send[length][N],int PSKmaker[length][N],double DX)
{
for(int i=0;i<length;i++)
for(int j=0;j<N;j++)
send[ i][j]=PSKmaker[ i][j]+gaosimaker(DX);
}
int result(int x)
{
if(x=1) return 1;
else if(x>1)
{
return result(x-1)+(int)pow(10,x-1);
}
else return 65535;
}
解码
//decode.h
#include <math.h>
#include <iostream.h>
#include <stdlib>
#define length 10
#define M 3
#define N 3
#define statenum (int)pow(2,M-1)
struct fencenode{
int currentstate;//左移寄存器当前状态(...000,...001,...010,...)
int incode[2];//输入数据后先后到达此点的两个点分别生成的卷积码
int outcode[2];//outcode[0]:输入0后生成的卷积码;outcode[1]:输入1后生成的卷积码
fencenode* in[2];//输入数据后到达此点的两个点
fencenode* out[2];//out[0]:输入0后从此点到达的点;out[1]:输入1后从此点到达的点
};
struct survive{
int hamingdist;//累积汉明距离
int sufferstate[length+1];//历经状态
int outcode[length];//对应的原信号
};
fencenode* initstatediag(int* generate);
int dist(int dist1,int* dist2);
int distcomp(int dist1,int dist2,int* dist0,int ham1,int ham2);
void copysur(survive* mysurvivenew,survive* survive,int i);
void copysur1(survive* surv1,int m,survive* surv2,int n,int i);
survive* findmin(survive* surv,int i);
void demodulate(int demodule[length][N],double receive[length][N])//解调2PSK信号
{//硬判决,判决门限0(小于0-->1,大于0-->0)
for(int i=0;i<length;i++)
{
for(int j=0;j<N;j++)
{
if(receive[ i][j]<0) demodule[ i][j]=1;
else demodule[ i][j]=0;
}
}
}
survive* viterbidecoder(int demodule[length][N],int* recover,int* generate)
{
int i,j;
survive* mysurvive=new survive[statenum];//最多有状态数(statenum)个幸存路径
survive* mysurvivenew=new survive[statenum];//最多有状态数(statenum)个幸存路径
fencenode* myfence=initstatediag(generate);
for(i=0;i<statenum;i++)//对幸存路径的初始化
{
mysurvive[ i].hamingdist=0;//设置汉明初值
mysurvive[ i].sufferstate[0]=0;//初始状态为0
mysurvivenew[ i].hamingdist=0;
mysurvivenew[ i].sufferstate[0]=0;
for(j=0;j<length;j++)
{
mysurvive[ i].sufferstate[j+1]=0;
mysurvivenew[ i].outcode[j]=0;
mysurvive[ i].sufferstate[j+1]=0;
mysurvivenew[ i].outcode[j]=0;
}
}
for(i=1;i<=length;i++)
{//生成第i个时间段内的幸存路线片段
if(pow(2,i)<=statenum)//对应篱笆图左部下落区
{
for(j=0;j<statenum;j++)
{//生成第j状态对应的幸存路线
int temp1=(j>>(M-i))&3,temp2=(j>>(M-1-i))&1;
mysurvive[j].sufferstate[ i]=myfence[temp1].out[temp2]->currentstate;
mysurvive[j].hamingdist+=dist(myfence[temp1].outcode[temp2],demodule[i-1]);
mysurvive[j].outcode[i-1]=temp2;
}
}
else
{//此篱笆图区域内每个状态必由2条路径到达,取汉明距离较小的为此时间段内此状态的幸存路径
copysur(mysurvivenew,mysurvive,i);
//准备由前一次结果重新排序幸存队列,将mysurvive赋给mysurvivenew
for(j=0;j<statenum;j++)
{
fencenode* node1=myfence[j].in[0];fencenode* node2=myfence[j].in[1];//找到转移至状态j的两个状态的地址
int code1=myfence[j].incode[0],code2=myfence[j].incode[1];
int k=distcomp(code1,code2,demodule[i-1],mysurvivenew[node1->currentstate].hamingdist,mysurvivenew[node2->currentstate].hamingdist);
if(0==k)
{//node1的距离小
copysur1(mysurvive,j,mysurvivenew,node1->currentstate,i);//赋予新路径
mysurvive[j].sufferstate[ i]=j;
mysurvive[j].hamingdist+=dist(code1,demodule[i-1]);
mysurvive[j].outcode[i-1]=j&1;
}
else if(1==k)
{//node2的距离小
copysur1(mysurvive,j,mysurvivenew,node2->currentstate,i);//赋予新路径
mysurvive[j].sufferstate[ i]=j;
mysurvive[j].hamingdist+=dist(code2,demodule[i-1]);
mysurvive[j].outcode[i-1]=j&1;
}
else
{
if(rand()<RAND_MAX/2)
{
copysur1(mysurvive,j,mysurvivenew,node1->currentstate,i);//赋予新路径
mysurvive[j].sufferstate[ i]=j;
mysurvive[j].hamingdist+=dist(code1,demodule[i-1]);
mysurvive[j].outcode[i-1]=j&1;
}
else
{
copysur1(mysurvive,j,mysurvivenew,node2->currentstate,i);//赋予新路径
mysurvive[j].sufferstate[ i]=j;
mysurvive[j].hamingdist+=dist(code2,demodule[i-1]);
mysurvive[j].outcode[i-1]=j&1;
}
}
}
}
}
survive* thebest=findmin(mysurvive,statenum);//找出statenum条路径中汉明距离最小的
for(i=0;i<statenum;i++)
{
for(j=0;j<length;j++) recover[j]=thebest->outcode[j];
}
delete []mysurvive;
delete []mysurvivenew;
delete []myfence;
return mysurvive;
}
fencenode* initstatediag(int* generate)
{
fencenode* trellisdiagram=new fencenode[statenum];//网格图雏形
for(int i=0;i<statenum;i++) trellisdiagram[ i].currentstate=i;
trellisdiagram[0].in[0]=&trellisdiagram[0];trellisdiagram[0].incode[0]=0;
trellisdiagram[0].in[1]=&trellisdiagram[2];trellisdiagram[0].incode[1]=3;
trellisdiagram[0].out[0]=&trellisdiagram[0];trellisdiagram[0].outcode[0]=0;
trellisdiagram[0].out[1]=&trellisdiagram[1];trellisdiagram[0].outcode[1]=7;
trellisdiagram[1].in[0]=&trellisdiagram[0];trellisdiagram[1].incode[0]=7;
trellisdiagram[1].in[1]=&trellisdiagram[2];trellisdiagram[1].incode[1]=4;
trellisdiagram[1].out[0]=&trellisdiagram[2];trellisdiagram[1].outcode[0]=1;
trellisdiagram[1].out[1]=&trellisdiagram[3];trellisdiagram[1].outcode[1]=6;
trellisdiagram[2].in[0]=&trellisdiagram[1];trellisdiagram[2].incode[0]=1;
trellisdiagram[2].in[1]=&trellisdiagram[3];trellisdiagram[2].incode[1]=2;
trellisdiagram[2].out[0]=&trellisdiagram[0];trellisdiagram[2].outcode[0]=3;
trellisdiagram[2].out[1]=&trellisdiagram[1];trellisdiagram[2].outcode[1]=4;
trellisdiagram[3].in[0]=&trellisdiagram[1];trellisdiagram[3].incode[0]=6;
trellisdiagram[3].in[1]=&trellisdiagram[3];trellisdiagram[3].incode[1]=5;
trellisdiagram[3].out[0]=&trellisdiagram[2];trellisdiagram[3].outcode[0]=2;
trellisdiagram[3].out[1]=&trellisdiagram[3];trellisdiagram[3].outcode[1]=5;
return trellisdiagram;
}
int dist(int dist1,int* dist2)//dist1:二进制数,dist2:十进制数组,宽均为N
{
int sum=0;
switch(4*dist2[0]+2*dist2[1]+dist2[2])
{
case 0:
switch(dist1)
{
case 0:sum=0;break;
case 1: case 2: case 4:sum=1;break;
case 3: case 5: case 6:sum=2;break;
case 7:sum=3;break;
}
break;
case 1:
switch(dist1)
{
case 1:sum=0;break;
case 0: case 3: case 5:sum=1;break;
case 2: case 4: case 7:sum=2;break;
case 6:sum=3;break;
}
break;
case 2:
switch(dist1)
{
case 2:sum=0;break;
case 0: case 3: case 6:sum=1;break;
case 1: case 4: case 7:sum=2;break;
case 5:sum=3;break;
}
break;
case 3:
switch(dist1)
{
case 3:sum=0;break;
case 1: case 2: case 7:sum=1;break;
case 0: case 5: case 6:sum=2;break;
case 4:sum=3;break;
}
break;
case 4:
switch(dist1)
{
case 4:sum=0;break;
case 0: case 5: case 6:sum=1;break;
case 1: case 2: case 7:sum=2;break;
case 3:sum=3;break;
}
break;
case 5:
switch(dist1)
{
case 5:sum=0;break;
case 1: case 4: case 7:sum=1;break;
case 0: case 3: case 6:sum=2;break;
case 2:sum=3;break;
}
break;
case 6:
switch(dist1)
{
case 6:sum=0;break;
case 2: case 4: case 7:sum=1;break;
case 0: case 3: case 5:sum=2;break;
case 1:sum=3;break;
}
break;
case 7:
switch(dist1)
{
case 7:sum=0;break;
case 3: case 5: case 6:sum=1;break;
case 1: case 2: case 4:sum=2;break;
case 0:sum=3;break;
}
break;
default:break;
}
/*
int temp,sum=0,i;//按位异或
for(i=0;i<N;i++)
{
temp=((dist1>>i)&1)^dist2[N-1-i];
sum=sum+temp;
}//求汉明距离*/
return sum;
}
int distcomp(int dist1,int dist2,int* dist0,int ham1,int ham2)
{
if((dist(dist1,dist0)+ham1)<(dist(dist2,dist0)+ham2)) return 0;
else if((dist(dist1,dist0)+ham1)>(dist(dist2,dist0)+ham2)) return 1;
else return 2;
}
void copysur(survive* mysurvivenew,survive* mysurvive,int i)
{//将路径数组mysurvive中每个路径的前i个数据赋给mysurvivene中的每个路径
for(int j=0;j<statenum;j++)
{
for(int l=1;l<i;l++)//找到不同路径的源头
{
if(mysurvivenew[j].sufferstate[l]!=mysurvive[j].sufferstate[l]) break;
}
mysurvivenew[j].hamingdist=mysurvive[j].hamingdist;
for(int k=l;k<i;k++)
{
mysurvivenew[j].sufferstate[k]=mysurvive[j].sufferstate[k];
mysurvivenew[j].outcode[k-1]=mysurvive[j].outcode[k-1];
}
}
}
void copysur1(survive* surv1,int m,survive* surv2,int n,int i)
{//将路径surv2的前i个数据赋给surv1,变成新的surv1
for(int l=1;l<i;l++)//找到不同路径的源头
{
if(surv1[m].sufferstate[l]!=surv2[n].sufferstate[l]) break;
}
surv1[m].hamingdist=surv2[n].hamingdist;
for(int k=l;k<i;k++)
{
surv1[m].sufferstate[k]=surv2[n].sufferstate[k];
surv1[m].outcode[k-1]=surv2[n].outcode[k-1];
}
}
survive* findmin(survive* surv,int i)//在i个路径中找汉明距离最小的
{
int j,k=0,best=0,smaller=surv[0].hamingdist;
for(j=0;j<i;j++)
{
if(surv[j].hamingdist<smaller) {smaller=surv[j].hamingdist;best=j;}
/*
else if(surv[j].hamingdist=smaller)
{
if(rand()&1) best=j;
}*/
}
return &surv[best];
}
int pe(int* origin,int* recover)
{
int i,sum=0;
for(i=0;i<length;i++)
if(origin[ i]!=recover[ i]) return 1;
return 0;
}
主模块
//main.cpp
//卷积码(N,1,M)的viterbi译码
#include <iostream.h>
#include "encode.h"
#include "decode.h"
#include <stdlib.h>
#include <math.h>
#include <iomanip.h>
#define length 10
#define M 3
#define N 3
#define statenum (int)pow(2,M-1)
int generate[N]={1,5,7};//001,101,111
int origin[length];
//encodeori[length][N]由信息位生成卷积码,码率为1/N,一个码元有N比特
int encodeori[length][N];
//PSKmaker[length][N]将卷积码调制成PSK信号
int PSKmaker[length][N];
//send[length][N]加入高斯白噪声后的信号
double send[length][N];
//从信道里接收的信号
int demodule[length][N];
//解码后的信号;由statenum个幸存路径筛选得出
int recover[length];
int viterbikey(double DX,int mode)
{
gensignal(origin,mode);
encode(encodeori,origin,generate);//产生信号
PSK(encodeori,PSKmaker);//PSK调制
addgaosi(send,PSKmaker,DX);//加高斯白噪声
demodulate(demodule,send);//解调
survive* sur=viterbidecoder(demodule,recover,generate);
int PE=pe(origin,recover);
return PE;
}
double viterbi()
{
int mode=3,sum=0,i,j,k;
double DX=1.0;
for(i=0;i<100;i++)
{
if(viterbikey(DX,mode)==1)
{
sum=sum+1;
cout<<"origin: ";
for(j=0;j<length;j++) cout<<origin[j];
cout<<endl<<"encodeori: ";
for(j=0;j<length;j++)
{
for(k=0;k<N;k++) cout<<encodeori[j][k];
cout<<" ";
}
cout<<endl<<"PSKmaker: ";
for(j=0;j<length;j++)
{
for(k=0;k<N;k++) cout<<;PSKmaker[j][k];
cout<<" ";
}
cout<<endl<<"send ";
for(j=0;j<length;j++)
{
for(k=0;k<N;k++) cout<<setprecision(2)<<send[j][k]<<" ";
cout<<" ";
}
cout<<endl<<"demodule: ";
for(j=0;j<length;j++)
{
for(k=0;k<N;k++) cout<<demodule[j][k];
cout<<" ";
}
cout<<endl<<"recover: ";
for(j=0;j<length;j++) cout<<recover[j];
cout<<endl<<endl;
}
}
return 1.0*sum/1000;
}
void main()
{
unsigned long seed;
time((long*)&seed);
srand(seed);
cout<<"误码率为： "<<viterbi()<<endl;
}
matlab版
length=1000;M=3;N=3;statenum=2^(M-1);error=zeros(1,101);
for xyz=1:1:101
DX=(xyz-1)/10;
origin=randi([0,1],1,length);generate=[0,1,1;0,0,1;1,1,1];
encodeori=zeros(length,N);PSKmaker=zeros(length,N);
send=zeros(length,N);demodule=zeros(length,N);recover=zeros(1,length);
shift=zeros(length,M);
shift(1,M)=origin(1);
for i=2:1:length
shift(i,:)=[shift(i-1,2:end),origin(i)];
end
encodeori=mod(shift*generate,2);
PSKmaker=ones(size(encodeori))-encodeori*2;
send=PSKmaker+DX*randn(size(PSKmaker));
demodule=send<0;
shift=zeros(statenum,1,2*M-3);
for i=1:1:statenum
j=i-1;k=0;
while j~=0
shift(i,1,end-k)=mod(j,2);
j=bitshift(j,-1);
k=k+1;
end
end
arraytonum=ones(M-1,1);
for i=M-1:-1:2
arraytonum(i-1)=2*arraytonum(i);
end
shift1=squeeze(shift(:,:,end-M+2:end));
myfence.incode(1).value=mod([zeros(statenum,1),shift1]*generate,2);
myfence.incode(2).value=mod([ones(statenum,1),shift1]*generate,2);
mysurvive.hamingdist=zeros(statenum,length);
mysurvive.sufferstate=zeros(statenum,length+1,M-1);
mysurvive.outcode=zeros(statenum,length);
for i=2:1:M
mysurvive.sufferstate(:,i,:)=shift(:,1,(i-1):(i+M-3));
end
mysurvive.outcode(:,1:M-1)=shift1;
mysurvive.sufferstate1=cat(3,zeros(statenum,M-1,1),mysurvive.sufferstate(:,2:M,:));
mysurvive.hamingdist(:,1)=sum(mod(squeeze(mysurvive.sufferstate1(:,1,:))*generate,2)~=(ones(statenum,1)*demodule(1,:)),2)';
for i=2:1:M-1
mysurvive.hamingdist(:,i)=sum(mod(squeeze(mysurvive.sufferstate1(:,i,:))*generate,2)~=(ones(statenum,1)*demodule(i,:)),2)+mysurvive.hamingdist(:,i-1);
end
for i=M:1:length
mysurvivenew=mysurvive;
for j=1:1:statenum
mysurvive.outcode(j,i)=mod(j-1,2);
mysurvive.sufferstate(j,i+1,:)=shift1(j,:);
state1=[0,shift1(j,1:end-1)]*arraytonum+1;state2=[1,shift1(j,1:end-1)]*arraytonum+1;
dist1=sum(myfence.incode(1).value(j,:)~=demodule(i,:))+mysurvivenew.hamingdist(state1,i-1);
dist2=sum(myfence.incode(2).value(j,:)~=demodule(i,:))+mysurvivenew.hamingdist(state2,i-1);
if (dist1<dist2)||((dist1==dist2)&&(randi([0,1])==0))
mysurvive.outcode(j,1:i-1)=mysurvivenew.outcode(state1,1:i-1);
mysurvive.sufferstate(j,2:i,:)=mysurvivenew.sufferstate(state1,2:i,:);
mysurvive.hamingdist(j,1:i-1)=mysurvivenew.hamingdist(state1,1:i-1);
mysurvive.hamingdist(j,i)=dist1;
else
mysurvive.outcode(j,1:i-1)=mysurvivenew.outcode(state2,1:i-1);
mysurvive.sufferstate(j,2:i,:)=mysurvivenew.sufferstate(state2,2:i,:);
mysurvive.hamingdist(j,1:i-1)=mysurvivenew.hamingdist(state2,1:i-1);
mysurvive.hamingdist(j,i)=dist2;
end
end
end
thebest=min(find(mysurvive.hamingdist(:,end)==min(mysurvive.hamingdist(:,end))));
recover=mysurvive.outcode(thebest,:);
error(xyz)=sum(recover~=origin)/length;
end
plot(linspace(0,10,101),error)