#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

//定義一些可變的參數(shù)
#define Length 4096         //輸入偽隨機(jī)碼長度
#define sigma 0.8          //noise的標(biāo)準(zhǔn)方差 
#define FrameNumber 1       //幀數(shù)
#define IterationTimes  6   //迭代次數(shù)
int NumOFGenerater=4;  //寄存器的長度，不包括輸入位
int generater1 = 0x1f;       // g1  十六進(jìn)制，生成器包括輸入的數(shù)據(jù)位
int generater2 = 0x11;       // g2          生成器最左邊位對應(yīng)輸入數(shù)據(jù)位


//產(chǎn)生任意均值與方差的n個正態(tài)分布隨機(jī)數(shù)
//函數(shù)的參數(shù)：mu和ro(μ和σ值)要事先給定。
void normal_noise(double *out,double mu,double ro,double seed,int number)
{
 int i;
 int M;
 M=65536;
 for ( i=0; i<number; i++)
 {
  double rand=0.0;
  for (int j=1; j<=12; j++)    //產(chǎn)生多個[0,1]區(qū)間內(nèi)均勻分布偽隨機(jī)數(shù)
  {
   seed=(seed)*2053+13849;
   seed=(unsigned long)seed%M;   //以65536為模對種子取余
   rand=rand+(seed)/(double)M;
  }
  out[i]=mu+ro*(rand-6.0);
 }

}

//RSC encoder
void RSC_SISO(int *pout,int *pin,int Length1,int state,int G1,int G2)
{

 int i,temp1,temp2,k,times,a;
 int regist=0; 
 
 for(i=0;i<Length1;i++)
 {
  temp1=(*pin)<<state;      
  regist=regist | temp1;
  
  temp2=regist & G1;
  times=0;
  for(k=0;k<state+1;k++)
  {
   if(temp2&(1<<k))
     times++;
  }
  a=times%2;
  
  temp1=1<<state;    //清第一位
  temp1=~temp1;
  regist=regist & temp1;

  temp1=a<<state;    //設(shè)置第一位
  regist=regist | temp1;

  temp2=regist & G2;
  times=0;
  for(k=0;k<state+1;k++)
  {
   if(temp2&(1<<k))
     times++;
  }
  *pout=times%2;
  pout++;


  //register 向低位移位
  regist=regist>>1;

  pin++;
 }

}


//deinterleaving algorith
//pi is the input array,pout is the output array;
template <class T>
void deinterleaving(T *pout,const T *pin,const int length)
{
 int i,j;
 int ir,jr,n,m,temp;
 static int P[8]={17,37,19,29,41,23,13,7};
 int N;
 N=(int)sqrt(length);
 typedef T* INT;
 INT *M;
 M=new INT[N];
 for (i=0;i<N;i++)
 {
  M[i]=new T[N];
 }

 for (i=0;i<N;i++)
  for (j=0;j<N;j++)
  {
   n=(N/2-1)*(i+j);
   ir=n%N;
   temp=(i+j)%8;
   m=(P[temp]*(j+1))-1;
   jr=m%N;
   M[ir][jr]=*pin;
   pin++;
  }

 for (j=0;j<N;j++)
  for (i=0;i<N;i++)
  {
   *pout=M[i][j];
   pout++;
  }

 for (i=0;i<N;i++)
 {
  delete[] M[i];  
 } 
 delete M;
}


//interleaving algorith
//pi is the input array,return a interleaved pointer
template <class T>
void interleaving(T *pout,const T *pin,const int length)
{ 
 int i,j;
 int ir,jr,n,m,temp; 
 static int P[8]={17,37,19,29,41,23,13,7};
 int N;
 N=(int)sqrt(length);
 typedef T* INT;
 INT *M;
 M=new INT[N];
 for (i=0;i<N;i++)
 {
  M[i]=new T[N];
 }

 for (j=0;j<N;j++)
  for (i=0;i<N;i++)
  {
   M[i][j]=*pin;
   pin++;
  }
 
 for (i=0;i<N;i++)
  for (j=0;j<N;j++)
  {
   n=(N/2-1)*(i+j);
   ir=n%N;
   temp=(i+j)%8;
   m=(P[temp]*(j+1))-1;
   jr=m%N;
   *pout=M[ir][jr];
   pout++;
  }

 for (i=0;i<N;i++)
 {
  delete[] M[i];  
 } 
 delete M;

}


//compute the variance of input data which be modulated by BPSK (+1,-1)
double computevariance(const double *pin,int length)
{
 int i;
 double mean,var,temp;
 temp=0;
 for (i=0;i<length;i++)
 {
  temp=temp+fabs(*(pin+i));
 }
 mean=temp/length;
 temp=0;
 for (i=0;i<length;i++)
 {
  temp=temp+pow(fabs(*(pin+i))-mean,2);
 }
 var=temp/length;
 return var;
}

//decode the data which use the MAP algorithm
//decoder_with_info(outdata,inx,iny,inz,Length)
void decoder_with_info(double *pout,const double *pinx,const double *piny,const double *pinz,const int Len,
        int state,int G1,int G2)
{
 int i,m,fm,k;
 double factor,temp0,temp1,out;
 double *temp;

 int NumOfState;
 NumOfState=1<<state;

 temp=new double[NumOfState];

 typedef double* DOUBLE;
 DOUBLE *a,*b;
 a=new DOUBLE[Len];
 b=new DOUBLE[Len+1];
 for(i=0;i<Len;i++)
 {
  a[i]=new double[NumOfState];
  b[i]=new double[NumOfState];
 }
 b[Len]=new double[NumOfState];

 if (a==NULL||b==NULL)
 {
  exit(1);
 }

 
 //initialize a[] and b[], a indicate forward path metric,b indicate backward path metric;
 for (i=1;i<NumOfState;i++)
 {
  a[0][i]=0;
 }
 a[0][0]=1;
 
 for (i=1;i<NumOfState;i++)
  b[Len][i]=1.0/Length;
 b[Len][0]=1;
 

 //compute the variance of input data
 double VAR=computevariance(pinx,Len);
 double Lc=2.0/VAR;

 //compute the gama,where gama indicate the current path metric
 double gama_with_info,gama_yk0,gama_yk1;
 int xs,xp; 
// int judge[16][16],c[16][16],d[16][16];

 typedef int* INT;
 INT *judge,*c,*d;
 judge=new INT[NumOfState];
 for(i=0;i<NumOfState;i++)
 {
  judge[i]=new int[NumOfState];  
 }
 c=new INT[NumOfState];
 for(i=0;i<NumOfState;i++)
 {
  c[i]=new int[NumOfState];  
 }
 d=new INT[NumOfState];
 for(i=0;i<NumOfState;i++)
 {
  d[i]=new int[NumOfState];  
 }

 int regist,indata,times,temp2,temp3;
 

 //matrix judge[i][j] judge the probability of the state i turn into the state j
 //the value of c[i][j] is the input data of encoder that tansform state i into state j
 //d[i][j]the output data when state i turn into state j

 for (int j=0;j<NumOfState;j++)
  for (i=0;i<NumOfState;i++)
  {
   judge[j][i]=0;
   c[j][i]=0;
   d[j][i]=0;
  }

 for (i=0;i<NumOfState;i++)
  for ( j=0;j<2;j++)
  {
   regist=i;
   temp2=j<<state;      
   regist=regist | temp2;
  
   temp3=regist & G1;
   times=0;
   for(k=0;k<state+1;k++)
   {
    if(temp3&(1<<k))
      times++;
   }
   indata=times%2;
   
   temp2=1<<state;    //清第一位
   temp2=~temp2;
   regist=regist & temp2;

   temp2=indata<<state;    //設(shè)置第一位
   regist=regist | temp2;
   
   temp3=regist & G2;     //編碼器輸出數(shù)據(jù)
   times=0;
   for(k=0;k<state+1;k++)
   {
    if(temp3&(1<<k))
      times++;
   }   

   //register 向低位移位,即下一個狀態(tài)
   regist=regist>>1;
   judge[i][regist]=1;
   c[i][regist]=j;
   d[i][regist]=times%2;

  }

 

 //compute the a[]
 for (k=1;k<Len;k++)
 {
  
  factor=0;
  for (m=0;m<NumOfState;m++)     //m--is the state of encoder C1 at time k
  {
   temp[m]=0;
      for (fm=0;fm<NumOfState;fm++)  //fm--is the state of encoder C1 at time k-1
   {
     if (judge[fm][m]==0)
      gama_with_info=0;
     else 
     {
      xs=2*c[fm][m]-1;
         xp=2*d[fm][m]-1;
      gama_with_info=exp(0.5*(xs*(*(pinz+k-1)+Lc*(*(pinx+k-1)))+Lc*(*(piny+k-1))*xp));
     }
       temp[m]=temp[m]+gama_with_info*a[k-1][fm];    
   }
   factor=factor+temp[m];
  }
  for (i=0;i<NumOfState;i++)
   a[k][i]=temp[i]/factor;   

 }

 //compute the value of b[]
 for (k=Len-1;k>0;k--)
 {
  factor=0;
  for (fm=0;fm<NumOfState;fm++)
  {
   temp[fm]=0;
   for (m=0;m<NumOfState;m++)
   {
    if (judge[fm][m]==0)
     gama_with_info=0;
    else
    {
     xs=2*c[fm][m]-1;
     xp=2*d[fm][m]-1;
     gama_with_info=exp(0.5*(xs*(*(pinz+k)+Lc*(*(pinx+k)))+Lc*(*(piny+k))*xp));
    }
    factor=factor+gama_with_info*a[k][fm];
    temp[fm]=temp[fm]+gama_with_info*b[k+1][m];
   }    
  }

  for (i=0;i<NumOfState;i++)
   b[k][i]=temp[i]/factor; 
 }

 //compute the value of output data
 for (k=1;k<=Len;k++)
 {
  
  temp0=0;
  temp1=0;
  for (m=0;m<NumOfState;m++)
   for (fm=0;fm<NumOfState;fm++)
   {
    if (judge[fm][m]==0)
    {
     gama_yk0=0;
     gama_yk1=0;
    }
    else 
    {
     xp=2*d[fm][m]-1;
     if(c[fm][m]==1) 
     {
      gama_yk0=0,gama_yk1=exp(0.5*Lc*(*(piny+k-1))*xp);
     }
     else 
     { 
      gama_yk1=0,gama_yk0=exp(0.5*Lc*(*(piny+k-1))*xp);
     }
    }
    temp1=temp1+gama_yk1*a[k-1][fm]*b[k][m];
    temp0=temp0+gama_yk0*a[k-1][fm]*b[k][m];
   }
  out=log(temp1/temp0);
  *(pout+k-1)=out;
 }
 for (i=0;i<Len;i++)
 {
  delete[] a[i];
  delete[] b[i];
 }
 delete[] b[Len];
 delete a;
 delete b;
}


//main function
void main()
{
 clock_t start,finish;
 double duration;
 int j;
 int i,flag,temp;
 double err,err1,err2;
 double var,Lc;
 err=0;
 err1=0;
 err2=0;

 start=clock(); 
for (j=0;j<FrameNumber;j++)
{
 double *d2out,*aa,*aa2in,*d1in,*d2in,*d1out;
 double *noise,*d12,*d21,*out;
 d1out=new double[Length];
 d2out=new double[Length];
 aa=new double[Length];
 aa2in=new double[Length];
 d1in=new double[Length];
 d2in=new double[Length]; 
 d12=new double[Length];
 d21=new double[Length];
 out=new double[Length];

 int *a,*a2,*b1,*encoderout,*b2;
 double *decoderin;
 a=new int[Length];
 a2=new int[Length];
 b1=new int[Length];
 b2=new int[Length];
 
//產(chǎn)生偽隨機(jī)輸入序列
 srand( (unsigned)time( NULL ) );
 for (i=0;i<Length;i++)
 {
  flag=rand()%2;
  a[i]=flag;
 }

//編碼部分
// RSC code C1,a is the data ,b is the redundant code;
 RSC_SISO(b1,a,Length,NumOFGenerater,generater1,generater2);

//uniform interleaving of input data
 interleaving(a2,a,Length);
 //uniforminterleaving(a2,a,Length);

//RSC code C2
 RSC_SISO(b2,a2,Length,NumOFGenerater,generater1,generater2);
 delete [] a2;

//puncture
 encoderout=new int[Length];
 for (i=0;i<Length;i++)
 { 
  if (i%2==0) 
  {
   encoderout[i]=b1[i];
  }
  else 
  {
   encoderout[i]=b2[i];
  }
 }
 delete [] b1;
 delete [] b2;


//信道部分,加入AWGN
//產(chǎn)生噪聲,數(shù)學(xué)期望mean=0,標(biāo)準(zhǔn)方差 ro=1.5
 noise=new double[Length*2];
 double seed=(double)time(NULL);
 normal_noise(noise,0,sigma,seed,Length*2);

//加入高斯白噪聲,aa,bb
 decoderin=new double[Length];
 for (i=0;i<Length;i++)
 {
  decoderin[i]=(2*encoderout[i]-1)+noise[i];
  aa[i]=(2*a[i]-1)+noise[i+Length];
 }
 delete [] encoderout;
 delete [] noise;


//解碼部分
//解復(fù)用
 for (i=0;i<Length;i++)
 {
  if (i%2==0) 
  {
   d1in[i]=decoderin[i];
   d2in[i]=0;
  }
  else  
  {
   d1in[i]=0;
   d2in[i]=decoderin[i];
  }
 }
 delete [] decoderin;

 var=computevariance(aa,Length);
 cout<<"var="<<var<<endl;
 cout<<"Eb/No="<<10*log10(0.5/var)<<endl;
 cout<<endl;
 Lc=2.0/var;

//compute the BER without code
 for (i=0;i<Length;i++)
 {
  if (aa[i]>=0) 
  {
   temp=1;
  }
  else 
  { 
   temp=0;
  }
  //
  if (a[i]!=temp) 
  {
   err++;
  }
 }
 cout<<"未迭代譯碼前 誤碼數(shù)量為:"<<err<<"個"<<endl;
 cout<<"誤碼率為:"<<err/(Length*FrameNumber)<<endl;
 cout<<endl;
 cout<<endl;


//iteration decode 

 interleaving(aa2in,aa,Length);
 //uniforminterleaving(aa2in,aa,Length); 
 for (i=0;i<Length;i++)
  d21[i]=0;

    // iteration
 for (int j=0;j<IterationTimes;j++)
 {
     err2=0;
  decoder_with_info(d1out,aa,d1in,d21,Length,NumOFGenerater,generater1,generater2);//decoder1 decode
  interleaving(d12,d1out,Length);
  //uniforminterleaving(d12,d1out,Length); 
  decoder_with_info(d2out,aa2in,d2in,d12,Length,NumOFGenerater,generater1,generater2);//decoder2 decoder
  deinterleaving(d21,d2out,Length);
  //uniformdeinterleaving(d21,d2out,Length);
  for (i=0;i<Length;i++)
   d2out[i]=d2out[i]+Lc*aa2in[i]+d12[i];
  deinterleaving(out,d2out,Length);
 //uniformdeinterleaving(out,d2out,Length);

 //cout<<"This is the output of six iteration"<<endl; 
  for (i=0;i<Length;i++)
  {  
   if (out[i]>0) 
   {
    out[i]=1;
   }
   else 
   {
    out[i]=0;
   }
   if (a[i]!=out[i])  
   {
    err2++;
   }
  
  }
  cout<<"第"<<j+1<<"次迭代后:"<<endl;
  cout<<"誤碼數(shù)量："<<err2<<"個"<<endl;
  cout<<"誤碼率為："<<err2/(Length*FrameNumber)<<endl;
  cout<<endl;
 }
 
 delete[] a;
 delete[] d1out;
 delete[] d2out;
 delete[] aa;
 delete[] aa2in;
 delete[] d1in;
 delete[] d2in;
 delete[] d12;
 delete[] d21;
 delete[] out;

}
    //output the duration
 finish=clock();
 duration=(double)(finish-start)/CLOCKS_PER_SEC;
 cout<<"該過程花費時間為："<<endl;
 cout<<duration<<"秒"<<endl;

}