iracode.txt

非规则重复累积码(IRA)在AWGN信道的编译码C程序

码率为0.5的IRA码在AWGN信道下的仿真程序
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#define maxfi 50 //信息位长度
#define N 7928  //交织位长度
#define a 8  //IRA码校验节点参数
#define codelth N/a+maxk
#define r N/a  //校验位个数
#define ENN 5  //信噪比个数
#define framenum 1000  //帧个数
#define MAX_iter 100  //译码循环总次数
int fin[maxfi],y[N/a+maxk],u2[maxk];
float rate,sigma;
float fi[maxk];
int alpha[7928];  //交织器
unsigned ui[maxk];
void encod(int ui[]);
void decoder();
void permute();
float caltan(float aa);
int sign(float num);
float normal(float mean,float variance);
main()
{
   float EBN0db[ENN]={1,4,1.5,1.6,1.7,1.8};
   float berr=0,en;
   int i,temp,frame,EN,err;
   float li[maxfi],s,Pb[ENN],receiver[codelth];


/////////////////////////////////////////////////PAGE2
   FILE*fperm;
//初始化
   for(i=0;i<maxfi;i++)
        {
         li[i]=0;
         fi[i]=0;
        }
   li[3]=0.252744;
   li[11]=0.081476;
   li[12]=0.327162;
   li[46]=0.184589;
   li[48]=0.154029;
   s=0;
   for(i=2;i<maxfi;i++)
      s+=li[i]/i;
   for(i=2;i<maxfi;i++)
      fi[i]=li[i]/(i*s);
   rate=(float)maxk/(float)(maxk+N/a);  //码率
   //读交织器数据
   permute();
   //information node parameter
   for(i=0;i<maxfi;i++)
      fin[i]=0;
   fin[3]=668;
   fin[11]=60;
   fin[12]=216;
   fin[46]=32;
   fin[48]=25;
//产生随机数据
   srand((unsigned)time(NULL));
   for(EN=0;EN<ENN;EN++))
   {
      berr=0;
      en=pow(10,EbN0db[EN]/10);
      sigma=1/sqrt(2*rate*en);
//generate input data;
      for(frame=0;frame<framenum;frame++)


/////////////////////////////////////////////////PAGE3
         {
            for(i=0;i<maxk;i++)
               {
                  temp=(int)(0.5+(float)rand()/0x7fff+1));
                  ui[i]=temp;
               }
//编码
            encode(ui);
//AWGN噪声信道
            for(i=0;i<codelth;i++)
                  receiver[i]=y[i]+sigma*normal(0,1);
//译码
            decoder(receiver);
            err=0;
            for(i=0;i<maxk;i++)
               {
                  if(ui[i]!=u2[i])
                      {
                          err=err+1;
                      }
               }
            berr=berr+err;
         }
//end of frame
      Pb[EN]=berr/(maxk*framenum);
      printf("%f ",Pb[EN]);
      }
//end of EN
      if((fperm=fopen("pb.txt","w"))!=NULL)
         {
           for(EN=0;EN<ENN;EN++)
           s=fprintf(fperm,"%f\n",Pb[EN]);
           fclose(fperm);
   }

void encod(unsigned ui[])
   {
      int i,j,jj,s,k;

/////////////////////////////////////////////////PAGE4

      unsigned v[N],c[N],x[N/a+1];
//repeat
      s=0;
      k=0;
      for(i=3;i<maxfi;i++)
      {
      for(j=0;j<fin[i];j++)
         {
         for(jj=0;jj<i;jj++)
            {
            v[s]=ui[k];
            s++;
            }
            k++;
         }
      }
//permute
      for(i=0;i<N;i++)
                c[i]=v[alpha[i]];
      //output data
      x[0]=0;
      for(j=1;j<r+1;j++)
         {
         s=0;
         for(i=0;i<a;i++)
         s+=c[(j-1)*a+i];
         x[j]=(x[j-1]+s)%2;
         }
      for(i=0;i<maxk;i++)
         y[i]=2*ui[i]-1;
      for(i=0;i<r;i++)
         y[i+maxk]=2*x[i+1]-1;
   }

void decoder(float receiver[])
   {
      float muc[N],mcu[N],myc[2][r],muc1[N],u1[N],baa[r],f[maxk];
      float s[maxk],b,ss,temp[2],miny[a+2],sal,div,al;

/////////////////////////////////////////////////PAGE5
      int i,j,k,pp,kkk,kj,err,mucs[N],mycs[2][r],sgn;
//初始化
//recerver data
b=(float)pow(sigma,2);
for(i=0;i<maxk;i++)
   {
   f[i]=2*receiver[i]/b;
   u2[i]=0;
   }
for(i=0;i<r;i++)
         baa[i]=2*receiver[maxk+i]/b;
for(i=0;i<r;i++)
   {
   myc[0][i]=0;
   mycs[0][i]=1;
   myc[1][i]=0;
   mycs[1][i]=1;
   mcy[0][i]=0;
   mcy[1][i]=0;
   }
for(i=0;i<N;i++)
   {
   mcu[i]=0;
   muc[i]=0;
   u1[i]=0;
   muc1[i]=0;
   mucs[i]=1;
   }
//****begin of decoding

for(k=0;k<MAX_iter;k++)
   {
//update my->c
   myc[0][r-1]=baa[r-1];
   mycs[0][r-1]=sign(baa[r-1]);
   for(i=0;i<r-1;i++)
      {
      myc[0][i]=baa[i]+mcy[0][i+1];

/////////////////////////////////////////////////PAGE6
      myc[1][i]=baa[i]+mcy[1][i];
      mycs[0][i]=sign(baa[i]+mcy[0][i+1]);
      mycs[1][i]=sign(baa[i]+mcy[1][i]);
      }
//update mu->c
   kk=0;
   pp=0;
   for(i=3;i<maxfi;i++)
   {
      for(j=0;j<fin[i];j++)
      {
      s[kk]=0;
      for(kj=0;kj<i;kj++)
         {
             s[kk]=u1[pp];
             pp++;
         }
      s[kk]=s[kk]+f[kk];
      if(s[kk]>0)
             u2[kk]=1;
      else
             u2[kk]=0;
      kk++;
      }
   }
   err=0;
   for(i=0;i<maxk;i++)
      {
      if(ui[i]!=u2[i])
           err=err+1;
      }
   if(err==0)
   break;
//get muc
   pp=0;
   kk=0;
   for(i=0;i<maxfi;i++)
      {
           for(j<fin[i];j++)
           {

/////////////////////////////////////////////////PAGE6
              for(kj=0;kj<i;kj++)
                 {
                    muc1[pp]=s[kk]-u1[pp];
                    pp++;
                 }
              kk++;
           }
      }
   for(i=0;i<N;i++)
   {
      muc[i]=muc1[alpha[i]];
      mucs[i]=sign(muc1[alpha[i]]);
   }
//update checknodes
//update mcu[0]
   ss=0;
   sgn=pow(-1,a+1);
   for(i=0;i<a;i++)
      {
      ss=ss+caltan(muc[i]);
      sgn=sgn*mucs[i];
      }
   ss=ss+caltan(myc[0][0]);
   sgn=sgn*mycs[0][0];
   for(i=0;i<a;i++)
      mcu[i]=sgn*mucs[i]*caltan(ss-caltan(muc[i]));
   for(i=1;i<r;i++)
   {
      ss=0;
      sgn=(int)pow(-1,a+2);
      for(j=0;j<a;j++)
        {
        ss=caltan(muc[i*a+j])+ss;
        sgn=sgn*mucs[i*a+j];
        }
      ss=ss+caltan(myc[1][i-1])+caltan(myc[0][i]);
      sgn=sgn*mycs[i][i-1]*mycs[0][i];
      for(j=0;j<a;j++)
        mcu[i*a+j]=sgn*mucs[i*a+j]*caltan(ss-caltan(muc[i*a+j]));

/////////////////////////////////////////////////PAGE7
      for(i=0;i<N;i++)
        u1[alpha[i]]=mcu[i];
//update mc->y
//cal mcy[1][0] 0-up 1-down
      ss=0;
      sgn=(int)pow(-1,a+1);
      for(i=0;i<a;i++)
        {
           ss=ss+caltan(muc[i]);
           sgn=sgn*mucs[i];
        }
        mcy[1][0]=sgn*caltan(ss);
        for(i=1;i<r;i++)
        {
           ss=0;sgn=(int)pow(-1,a+2);
           for(j=0;j<a;j++)
             {
                ss=ss+caltan(muc[i*a+j]);
                sgn=sgn*mucs[i*a+j];
             }
           mcy[0][i]=sgn*mycs[0][i]*caltan(ss+caltan(myc[0][i]));
           mcy[1][i]=sgn*mycs[1][i-1]*caltan(ss+caltan(myc[1][i-1]));
        }
  }
//end of k times decoding
}
//读交织器数据
void permute()
{
  FILE*fperm;
  int i,rr;
  if((fperm=fopen("s7928.txt","r"))!=NULL)
    {
       for(i=0;i<N;i++)
          rr=fscanf(fperm,"%d",&alpha[i]);
       fclose(fperm);

/////////////////////////////////////////////////PAGE8
    }
}
float normal(float mean,float variance)
{
  int i,NN;
  float x=0;
  NN=30;
  for(i=0;i<NN;i++)
     x+=(float)rand()/(RAND_MAX+1.0);
  return(mean+sqrt(variance)*(x-NN/2.0)/sqrt(NN/12.0));
}
float caltan(float aa)
{
  float tan1,aa1;
  aa1=fabs(aa);
  if(aa1==0)
    tan1=1000;
  else
    tan1=(float)log((exp(aa1)+1)/(exp(aa1)-1));
    if(fabs(tan1)>20)tan1=20;
  return tan1;
}
int sign(float num)
{
  int i;
  if(num>0)
    i=1;
  else
    if(num<0)
    i=-1;
  else
    i=0;
  return i;
}