turbode.cpp

wcdma协议中的turbo码链路仿真

/*Turbo码译码程序，2006-08-29。使用Log-MAP译码算法，在译码中不使用尾比特*/

#include<math.h>
//#include<time.h>
//#include<stdlib.h>
#include<iostream.h>
#include <fstream.h> 
 
ofstream myf("d:\\WCDMA.txt", ios::trunc);

//variable
int tmp;
const int block_size=40;  //码块长度
int niter=10;  //迭代的次数
double * rx_data;  //接收到的数据
double SNR_db;
double SNR;
double rate=(double)1/3;
double *rx_s1;
double *rx_s2;
double *rx_p1;
double *rx_p2;  //接收到的信息比特1,信息比特2，校验比特1，校验比特2,信息比特2和信息比特1相同
double *ext1,*ext2;  //译码器1输出的外信息，译码器2输出的外信息
//int *de_data;//译码输出的数据
int de_data[block_size];
double **alpha,**beta;  //前项路径度量，反向路径度量
double **gama;  //分支度量
double Lc; //信道置信度 Lc=4*(Es/N0)
int const Infty_1=32767; //表示无穷大,用于限制alpha，beta
int const Infty_2=8191;  //表示无穷大，用于限制gamma

int *intlv_array;  //交织数组
int *de_intlv_array;  //解交织数组

//function
void decode_init();  //初始化
void turbo_decode();  //译码调度函数
void log_map_ext(double *rx_s,double *rx_p,double *ext_in,double *ext_out);  //Log-MAP译码算法,计算外信息
void log_map_llr(double *rx_s,double *rx_p,double *ext_in,double *llr);  //Log-MAP译码算法，计算对数似然比
int max(int a,int b); //比较a和b的最大值并返回
double max(double a,double b); //比较a和b的最大值并返回
int gcd(int a, int b);  //判断a和b是否互质，返回1:互质；0：不互质
double max_star(double a, double b);  //计算max star 运算 max*(a,b)=max(a,b)+log(1+exp(-abs(a-b))
void creat_interleaver();  //产生交织数组和解交织数组
void interleaver(int *data);  //交织
void deinterleaver(int *data);  //解交织
void interleaver(double *data);  //交织
void deinterleaver(double *data);  //解交织

void decision(double *indata);  //判决
void freemem();  //释放内存


void print(int * data,int len)
{
	for(int i=0;i<len;i++)
		cout<<data[i]<<" ";
	cout<<endl;
}

void print(double * data,int len)
{
	for(int i=0;i<len;i++)
		cout<<data[i]<<"  ";
	cout<<endl;
}

void zero(int *data,int len)
{
	for(int i=0;i<len;i++)
		data[i]=0;
}
void zero(double *data,int len)
{
	for(int i=0;i<len;i++)
		data[i]=0;
}

void rsc(int * input,int * output,int len)
{
	int reg[3];
	zero(reg,3);
	
	int feedback[3] = {0,0,0};
	for(int i=0;i<len;i++)
	{
		feedback[2] = reg[2] ^ reg[1];
		feedback[0] = feedback[2] ^ input[i];
		feedback[1] = feedback[0] ^ reg[0];
		output[i] =  feedback[1] ^ reg[2];
		reg[2] = reg[1];
		reg[1] = reg[0];
		reg[0] = feedback[0];
	}
}

void encode(int * input)
{
	rx_data = new double[3*block_size];
//	print(input,block_size);
//	int * y =new int[block_size];
	int y[block_size];
	zero(y,block_size);
	rsc(input,y,block_size);
//	print(y,block_size);
//	int * inter = new int[block_size];
	int inter[block_size];
	zero(inter,block_size);
	//init();
	for(int i=0;i<block_size;i++)
		inter[i] = input[intlv_array[i]];
//	interleaver(input);
	
//	int *y1 = new int[block_size];
	int y1[block_size];
	zero(y1,block_size);
	rsc(inter,y1,block_size);
	print(inter,block_size);
		for(i=0;i<block_size;i++)
		{
			rx_data[3*i] = input[i];
			rx_data[3*i+1] = y[i];
			rx_data[3*i+2] = y1[i];
		}
	//	print(rx_data,3*block_size);
		for(i=0;i<3*block_size;i++)
		{
			if(rx_data[i]==0)
				rx_data[i]=-1;
		}
}


void turbo_decode()  //控制译码的函数
{

	interleaver(rx_s2);  //对输入第二级译码器的数据进行交织运算
//	cout<<"rx_s2:";
//	print(rx_s2,block_size);
	int i;  //用于控制迭代译码次数的循环变量
	
	for(i=0;i<niter-1;i++)  //首先迭代N-1次，在译码器之间传递外信息
	{
		log_map_ext(rx_s1,rx_p1,ext2,ext1);  //第一级译码器译码，输出外信息
		interleaver(ext1);  //对第一级译码器输出的外信息进行交织，然后输入到第二级译码器
		log_map_ext(rx_s2,rx_p2,ext1,ext2);  //第二级译码器译码，输出外信息
		deinterleaver(ext2);  //对第二级译码器输出的外信息解交织，然后输入到第一级译码器
	}
	//最后一次译码
	log_map_ext(rx_s1,rx_p1,ext2,ext1);  //第一级译码器译码，输出外信息
	interleaver(ext1);   //对第一级译码器输出的外信息进行交织，然后输入到第二级译码器
	log_map_llr(rx_s2,rx_p2,ext1,ext2); //第二级译码器译码，输出对数似然比
	deinterleaver(ext2);  //将第二级译码器输出的对数似然比解交织
/*	for(i=0;i<block_size;i++)
	{
		cout<<ext2[i]<<" ";
	}
	cout<<endl;*/

	
	decision(ext2);  //对第二级译码器输出的对数似然比进行硬判决，作为最后的译码输出

	
	for (i=0;i<block_size;i++ )
	{
		de_data[i]=(int)ext2[i];
	}
		
	freemem(); 	
}

void decode_init()
{
	int i=0,j=0;
    rx_s1=new double[block_size];
	zero(rx_s1,block_size);
	rx_s2=new double[block_size];
	zero(rx_s2,block_size);
	rx_p1=new double[block_size];
	zero(rx_p1,block_size);
	rx_p2=new double[block_size];
	zero(rx_p2,block_size);
/*			for(i=0;i<3*block_size;i++)
			{
					cout<<rx_data[i]<<"  ";
			}
			cout<<endl;*/
	for(i=0;i<3*block_size;i=i+3)
	{
		rx_s1[j]=rx_data[i];
		rx_s2[j]=rx_data[i];
		rx_p1[j]=rx_data[i+1];
		rx_p2[j]=rx_data[i+2];
		j++;
	}
//cout<<"rx_s2:";	print(rx_s2,block_size);
/*	cout<<"rx_s2:";
	print(rx_s1,block_size);
	cout<<endl;
	cout<<"rx_p1:";
	print(rx_p1,block_size);
	cout<<endl;
	cout<<"rx_p2:";
	print(rx_p2,block_size);
	cout<<endl;*/

	//初始化其它
	ext1=new double[block_size]; 
	ext2=new double[block_size];
	zero(ext1,block_size);
	zero(ext2,block_size);

	alpha=new double*[block_size];
	beta=new double*[block_size];
	for(i=0;i<block_size;i++)  //创建alpha,beta矩阵
	{
		alpha[i]=new double[8];
		beta[i]=new double[8];
	}
	for(i=0;i<block_size;i++)
	{
		for(j=0;j<8;j++)
		{
			alpha[i][j]=0;
			beta[i][j]=0;
		}
	}

	alpha[0][0]=0;  //初始化alpha	
	for(i=1;i<8;i++)
	{
		alpha[0][i]=-Infty_1;
	}
	
	
	for(i=0;i<8;i++)	//初始化beta
	{
		beta[block_size-1][i]=-Infty_1;
	}
	
	gama=new double*[block_size];
	for(i=0;i<block_size;i++)
	{
		gama[i]=new double[4];//00 01 10 11
	}
	for(i=0;i<block_size;i++)
	{
		for(j=0;j<4;j++)
			gama[i][j]=0;
	}
	//Lc=4*pow(10,SNR/10);  //Lc=4*(EsN0)
//	Lc=4*0.5*pow(10,SNR_db/10);//SNR=(long double)pow(10.0,SNR_db/10.0);
	Lc=4*rate*SNR;
//	cout<<"snr"<<SNR<<endl;
//	cout<<"lc"<<Lc<<endl;

//	de_data=new int[block_size];
//	int de_data[block_size];
	zero(de_data,block_size);	
}

int gcd(int a, int b)
{
    int v, u, t, q;
    u = abs(a);
    v = abs(b);
    while (v>0) 
   {
      q = u / v;
      t = u - v*q;
      u = v;
      v = t;
    }
    return(u);
}

void creat_interleaver()
{
	intlv_array=new int[block_size];  //交织和解交织数组初始化
	de_intlv_array=new int[block_size];

    const int MAX_INTERLEAVER_SIZE = 5114;
    const int MIN_INTERLEAVER_SIZE = 40;
    int K;  
    int R;  
    int C;  
    int p=0;  
    int v=0;  
    int *s; 
    int *q; 
    int *r; 
    int *T; 
    int **U; 
    int i, j, qj, temp, row, col, index, count;
    int interleaver_size = block_size;

    if(interleaver_size > MAX_INTERLEAVER_SIZE)
    {
	cout<<"td-scdma_turbo_interleaver_sequence: The interleaver size is too large"<<endl;
    }	
    if(interleaver_size < MIN_INTERLEAVER_SIZE)
    {
		cout<<"td-scdma_turbo_interleaver_sequence: The interleaver size is too small"<<endl;
	}
    		 
      K = interleaver_size;
      static int primes[55] = {2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257};
      static int roots[55] = {0,0,0,3,2,2,3,2,5,2,3,2,6,3,5,2,2,2,2,7,5,3,2,3,5,2,5,2,6,3,3,2,3,2,2,6,5,2,5,2,2,2,19,5,2,3,2,3,2,6,3,7,7,6,3};
      int primes_length  = 55;
      if ((K>=40) && (K<=159)) 
      {
	      R = 5;
      } 
      else if ( ((K>=160)&&(K<=200)) || ((K>=481)&&(K<=530)) ) 
      		{
			R = 10;
      		} 
      	    else {
			R = 20;
     		 }
      if ((K>=481)&&(K<=530)) 
      {
		p = 53;
		v = 2;
		C = p;
      } 
      else 
      {
		for (i=0; i<primes_length; i++) 
		{
		  if ( ( double(primes[i]+1) - double(K)/double(R) ) >= 0.0 ) 
		  {
		    p = primes[i];
		    v = roots[i];
		    break;
		  }
		}
	if ( ( double(p) - double(K)/double(R) ) >= 0.0 ) 
	{
	  if ( ( double(p) - 1.0 - double(K)/double(R) ) >= 0.0 ) 
	  {
	    C = p-1;
	  } 
	  else 
	  {
	    C = p;
	  }
	} 
	else 
	{
	  C = p+1;
	}
	  }
      
	  s = new int[p-1];
      s[0] = 1;
      for (i=1; i<=(p-2); i++) 
      {
		  s[i] = (int) fmod((v * s[i-1]), p);
      }
      q = new int[R];
      q[0] = 1;
      for (j=1; j<=(R-1); j++) 
      {
		  for (i=0; i<primes_length; i++) 

		  {
			  qj = primes[i];	  
			  if ( (qj>6) && (qj>q[j-1]) ) 	 
			  {	   
				  if (gcd((int) qj,(int) (p-1)) == 1) 	  
				  {
					  q[j] = qj;
					  break;
				  }
			  }
		  }     
	  }
      static int Pat1[20] = {19,9,14,4,0,2,5,7,12,18,10,8,13,17,3,1,16,6,15,11};
      static int Pat2[20] = {19,9,14,4,0,2,5,7,12,18,16,13,17,15,3,1,6,11,8,10};
      static int Pat3[10] = {9,8,7,6,5,4,3,2,1,0};
      static int Pat4[5] = {4,3,2,1,0};
      if (K>=3211)
	  {
		  T = Pat1;      
	  } 
	  else if (K>=3161)
	  {	
		  T = Pat2;     
	  } 
	  else if (K>=2481)
	  {
		  T = Pat1;
      } 
	  else if (K>=2281) 
	  {
		  T = Pat2;
      } 
	  else if (K>=531) 
	  {	
		  T = Pat1;
      } 
	  else if (K>=481)
	  {	
		  T = Pat3;
      }
	  else if (K>=201)
	  {
		  T = Pat1;
      } 
	  else if (K>=160) 
	  {
		  T = Pat3;
      }
	  else
	  {
		  T = Pat4;
      }

      r=new int[R];
      for (j=0; j<=(R-1); j++)
	  {
		  r[ T[j] ] = q[j];
      }
      U = new int*[R];
      for (i=0;i<R;i++) 
      {	
		  U[i] = new int[C];
      }

      if (C==p) 
	  {    	
		  for (j=0; j<=(R-1); j++) 
		  {	  
			  for (i=0; i<=(p-2); i++) 
			  {	   
				  U[j][i] = s[(int) fmod((int) (i*r[j]), (int) (p-1)) ];
			  }	  
			  U[j][p-1] = 0;	
		  }     
	  } 
	  else if (C==(p+1))
	  {	
		  for (j=0; j<=(R-1); j++)
		  {
			  for (i=0; i<=(p-2); i++) 
			  {	    
				  U[j][i] = s[(int) fmod((int) (i*r[j]),(int) (p-1)) ]; 	  
			  }	 
			  U[j][p-1] = 0;	 
			  U[j][p] = p;	
		  }	
		  if (K == (C*R))
		  {	  
			  temp = U[R-1][p];	 
			  U[R-1][p] = U[R-1][0];	  
			  U[R-1][0] = temp;
		  }     
	  } 
	  else if (C==(p-1)) 
	  {	
		  for (j=0; j<=(R-1); j++)
		  {
			  for (i=0; i<=(p-2); i++)
			  {
				  U[j][i] = s[(int) fmod((int) (i*r[j]),(int) (p-1)) ] - 1;	 
			  }	
		  }
      }
      
	  count = 0;    
	  for (i=0; i<C; i++)     
	  {	
		  for (j=0; j<R; j++) 	
		  {		
			  row = T[j];		
			  col = U[row][i];		
			  index = row * C + col;		
			  if (index < K) 		
			  {			
				  intlv_array[count] = index;			
				  count++;	    
			  }	
		  }
	  }
      for (i=0;i<K;i++)
      {
		  de_intlv_array[intlv_array[i]] = i;
	  }          
	  for (i=0;i<R;i++)       
	  {	
		  delete U[i];
      }
	  delete []s;
	  delete []q;
	  delete []r;
//	  delete []T;
	  delete []U;
}

void interleaver(int *data)
{
	int i;
	int *temp=new int[block_size];
	for(i=0;i<block_size;i++)
	{
		temp[i]=data[i];
	}
	for(i=0;i<block_size;i++)
	{
		data[i]=temp[intlv_array[i]];
	}
	delete []temp;
}

void interleaver(double *data)
{

	int i;
	double *temp=new double[block_size];
	for(i=0;i<block_size;i++)
	{
		temp[i]=data[i];
	}
	for(i=0;i<block_size;i++)
	{
		data[i]=temp[intlv_array[i]];
	}
	delete []temp;
}

void deinterleaver(int *data)
{
	int i;
	int *temp=new int[block_size];
	for(i=0;i<block_size;i++)
	{
		temp[i]=data[i];
	}
	for(i=0;i<block_size;i++)
	{
		data[i]=temp[de_intlv_array[i]];
	}
	delete []temp;
}

void deinterleaver(double *data)
{
	int i;
	double *temp=new double[block_size];
	for(i=0;i<block_size;i++)
	{
		temp[i]=data[i];
	}
	for(i=0;i<block_size;i++)
	{
		data[i]=temp[de_intlv_array[i]];
	}
	delete []temp;
}