复件 turbode.cpp

wcdma协议中的turbo码链路仿真

int max(int a,int b)
{
	return a>=b?a:b;
}

double max(double a,double b)
{
	return a>=b?a:b;
}

void log_map_ext(double *rx_s,double *rx_p,double *ext_in,double *ext_out)
{	
	int i,j;
	//calculate gama, gama = 0.5*ext*Xs + 0.5*Lc*Ys*Xs + 0.5*Lc*Yp*Xp
	//Xs:发送的信息位，Ys：接收到的信息位rx_s，Xp:发送的信息位，Yp：接收到的信息位rx_p
	for(i=0;i<block_size;i++)
	{
		gama[i][0]=-0.5*ext_in[i]-0.5*Lc*rx_s[i]-0.5*Lc*rx_p[i]; //gama(00)Xs=-1，Xp=-1
		gama[i][1]=-0.5*ext_in[i]-0.5*Lc*rx_s[i]+0.5*Lc*rx_p[i]; //gama(01)Xs=-1，Xp=1
		gama[i][2]=0.5*ext_in[i]+0.5*Lc*rx_s[i]-0.5*Lc*rx_p[i];  //gama(10)Xs=1，Xp=-1
		gama[i][3]=0.5*ext_in[i]+0.5*Lc*rx_s[i]+0.5*Lc*rx_p[i];  //gama(11)Xs=1，Xp=1
		for(j=0;j<4;j++)  //防止算法溢出
		{
			if(gama[i][j]>Infty_2)
			{
				gama[i][j]=Infty_2;
			}
			else if(gama[i][j]<-Infty_2)
			{
				gama[i][j]=-Infty_2;
			}
		}
	}

	//计算 alpha
	double temp1,temp2;
	double *temp_max=new double[block_size];
	temp_max[0]=0;
	for(i=1;i<block_size;i++)
	{
		//alpha(0)
		temp1=alpha[i-1][0]+gama[i-1][0];
		temp2=alpha[i-1][1]+gama[i-1][3];
		alpha[i][0]=max_star(temp1,temp2); 
		//alpha(4)
		temp1=alpha[i-1][0]+gama[i-1][3];
		temp2=alpha[i-1][1]+gama[i-1][0];
		alpha[i][4]=max_star(temp1,temp2);
		//alpha(1)
		temp1=alpha[i-1][2]+gama[i-1][2];
		temp2=alpha[i-1][3]+gama[i-1][1];
		alpha[i][1]=max_star(temp1,temp2);
		//alpha(5)
		temp1=alpha[i-1][2]+gama[i-1][1];
		temp2=alpha[i-1][3]+gama[i-1][2];
		alpha[i][5]=max_star(temp1,temp2);
		//alpha(2)
		temp1=alpha[i-1][4]+gama[i-1][1];
		temp2=alpha[i-1][5]+gama[i-1][2];
		alpha[i][2]=max_star(temp1,temp2);
		//alpha(6)
		temp1=alpha[i-1][4]+gama[i-1][2];
		temp2=alpha[i-1][5]+gama[i-1][1];
		alpha[i][6]=max_star(temp1,temp2);
		//alpha(3)
		temp1=alpha[i-1][6]+gama[i-1][3];
		temp2=alpha[i-1][7]+gama[i-1][0];
		alpha[i][3]=max_star(temp1,temp2);
		//alpha(7)
		temp1=alpha[i-1][6]+gama[i-1][0];
		temp2=alpha[i-1][7]+gama[i-1][3];
		alpha[i][7]=max_star(temp1,temp2);

		//归一化
		temp_max[i]=alpha[i][0];
		for(j=1;j<8;j++)
		{
			if(alpha[i][j]>temp_max[i])
			{
				temp_max[i]=alpha[i][j];
			}
		}
		for(j=0;j<8;j++)
		{
			alpha[i][j]=alpha[i][j]-temp_max[i];
			if (alpha[i][j]>Infty_1)
			{
				alpha[i][j]=Infty_1;
			}
			else if(alpha[i][j]<-Infty_1)
			{
				alpha[i][j]=-Infty_1;
			}
		}

	}

	//计算 beta
	for(i=block_size-1;i>0;i--)
	{
		//beta(0)
		temp1=beta[i][0]+gama[i][0];
		temp2=beta[i][4]+gama[i][3];
		beta[i-1][0]=max_star(temp1,temp2);
		//beta(1)
		temp1=beta[i][0]+gama[i][3];
		temp2=beta[i][4]+gama[i][0];
		beta[i-1][1]=max_star(temp1,temp2);
		//beta(2)
		temp1=beta[i][1]+gama[i][2];
		temp2=beta[i][5]+gama[i][1];
		beta[i-1][2]=max_star(temp1,temp2);
		//beta(3)
		temp1=beta[i][1]+gama[i][1];
		temp2=beta[i][5]+gama[i][2];
		beta[i-1][3]=max_star(temp1,temp2);
		//beta(4)
		temp1=beta[i][2]+gama[i][1];
		temp2=beta[i][6]+gama[i][2];
		beta[i-1][4]=max_star(temp1,temp2);
		//beta(5)
		temp1=beta[i][2]+gama[i][2];
		temp2=beta[i][6]+gama[i][1];
		beta[i-1][5]=max_star(temp1,temp2);
		//beta(6)
		temp1=beta[i][3]+gama[i][3];
		temp2=beta[i][7]+gama[i][0];
		beta[i-1][6]=max_star(temp1,temp2);
		//beta(7)
		temp1=beta[i][3]+gama[i][0];
		temp2=beta[i][7]+gama[i][3];
		beta[i-1][7]=max_star(temp1,temp2);
	
		//归一化
		for(j=0;j<8;j++)
		{
			beta[i-1][j]=beta[i-1][j]-temp_max[i-1];
			if (beta[i-1][j]>Infty_1)
			{
				beta[i-1][j]=Infty_1;
			}
			else if(beta[i-1][j]<-Infty_1)
			{
				beta[i-1][j]=-Infty_1;
			}
		}
		
	}
	delete []temp_max;
	//计算输出的外信息
	int pm_b_1[8]={4,0,1,5,6,2,3,7};
	int pm_b_0[8]={0,4,5,1,2,6,7,3};
	int bm_1[8]={3,3,2,2,2,2,3,3};
	int bm_0[8]={0,0,1,1,1,1,0,0};
	double temp_star1,temp_star0;

	for(i=0;i<block_size;i++)
	{
		temp1=0;
		temp2=0;
		temp_star1=alpha[i][0]+beta[i][4]+gama[i][3];
		temp_star0=alpha[i][0]+beta[i][0]+gama[i][0];
		for(j=1;j<8;j++)
		{
			//计算等于1的
			temp1=temp_star1;
			temp2=alpha[i][j]+beta[i][pm_b_1[j]]+gama[i][bm_1[j]];
			temp_star1=max_star(temp1,temp2);

			//计算等于0的
			temp1=temp_star0;
			temp2=alpha[i][j]+beta[i][pm_b_0[j]]+gama[i][bm_0[j]];
		    temp_star0=max_star(temp1,temp2);
		}
		//计算 LLR
		ext_out[i]=temp_star1-temp_star0;
		//计算外信息
		ext_out[i]=ext_out[i]-Lc*rx_s[i]-ext_in[i];
	}
}

void log_map_llr(double *rx_s,double *rx_p,double *ext_in,double *llr)
{
	int i,j;
	//calculate gama, gama = 0.5*ext*Xs + 0.5*Lc*Ys*Xs + 0.5*Lc*Yp*Xp
	//Xs:发送的信息位，Ys：接收到的信息位rx_s，Xp:发送的信息位，Yp：接收到的信息位rx_p
	for(i=0;i<block_size;i++)
	{
		gama[i][0]=-0.5*ext_in[i]-0.5*Lc*rx_s[i]-0.5*Lc*rx_p[i]; //gama(00)Xs=-1，Xp=-1
		gama[i][1]=-0.5*ext_in[i]-0.5*Lc*rx_s[i]+0.5*Lc*rx_p[i]; //gama(01)Xs=-1，Xp=1
		gama[i][2]=0.5*ext_in[i]+0.5*Lc*rx_s[i]-0.5*Lc*rx_p[i];  //gama(10)Xs=1，Xp=-1
		gama[i][3]=0.5*ext_in[i]+0.5*Lc*rx_s[i]+0.5*Lc*rx_p[i];  //gama(11)Xs=1，Xp=1
		for(j=0;j<4;j++)  //防止算法溢出
		{
			if(gama[i][j]>Infty_2)
			{
				gama[i][j]=Infty_2;
			}
			else if(gama[i][j]<-Infty_2)
			{
				gama[i][j]=-Infty_2;
			}
		}
	}
	
	//计算 alpha
	double temp1,temp2;
	double *temp_max=new double[block_size];
	temp_max[0]=0;

	for(i=1;i<block_size;i++)
	{
		//alpha(0)
		temp1=alpha[i-1][0]+gama[i-1][0];
		temp2=alpha[i-1][1]+gama[i-1][3];
		alpha[i][0]=max_star(temp1,temp2);
		//alpha(4)
		temp1=alpha[i-1][0]+gama[i-1][3];
		temp2=alpha[i-1][1]+gama[i-1][0];
		alpha[i][4]=max_star(temp1,temp2);
		//alpha(1)
		temp1=alpha[i-1][2]+gama[i-1][2];
		temp2=alpha[i-1][3]+gama[i-1][1];
		alpha[i][1]=max_star(temp1,temp2);
		//alpha(5)
		temp1=alpha[i-1][2]+gama[i-1][1];
		temp2=alpha[i-1][3]+gama[i-1][2];
		alpha[i][5]=max_star(temp1,temp2);
		//alpha(2)
		temp1=alpha[i-1][4]+gama[i-1][1];
		temp2=alpha[i-1][5]+gama[i-1][2];
		alpha[i][2]=max_star(temp1,temp2);
		//alpha(6)
		temp1=alpha[i-1][4]+gama[i-1][2];
		temp2=alpha[i-1][5]+gama[i-1][1];
		alpha[i][6]=max_star(temp1,temp2);
		//alpha(3)
		temp1=alpha[i-1][6]+gama[i-1][3];
		temp2=alpha[i-1][7]+gama[i-1][0];
		alpha[i][3]=max_star(temp1,temp2);
		//alpha(7)
		temp1=alpha[i-1][6]+gama[i-1][0];
		temp2=alpha[i-1][7]+gama[i-1][3];
		alpha[i][7]=max_star(temp1,temp2);
		
		//归一化
		temp_max[i]=alpha[i][0];
		for(j=1;j<8;j++)
		{
			if(alpha[i][j]>temp_max[i])
			{
				temp_max[i]=alpha[i][j];
			}
		}
		for(j=0;j<8;j++)
		{
			alpha[i][j]=alpha[i][j]-temp_max[i];
			if (alpha[i][j]>Infty_1)
			{
				alpha[i][j]=Infty_1;
			}
			else if(alpha[i][j]<-Infty_1)
			{
				alpha[i][j]=-Infty_1;
			}
		}
	}

	//计算 beta
	for(i=block_size-1;i>0;i--)
	{
		//beta(0)
		temp1=beta[i][0]+gama[i][0];
		temp2=beta[i][4]+gama[i][3];
		beta[i-1][0]=max_star(temp1,temp2);
		//beta(1)
		temp1=beta[i][0]+gama[i][3];
		temp2=beta[i][4]+gama[i][0];
		beta[i-1][1]=max_star(temp1,temp2);
		//beta(2)
		temp1=beta[i][1]+gama[i][2];
		temp2=beta[i][5]+gama[i][1];
		beta[i-1][2]=max_star(temp1,temp2);
		//beta(3)
		temp1=beta[i][1]+gama[i][1];
		temp2=beta[i][5]+gama[i][2];
		beta[i-1][3]=max_star(temp1,temp2);
		//beta(4)
		temp1=beta[i][2]+gama[i][1];
		temp2=beta[i][6]+gama[i][2];
		beta[i-1][4]=max_star(temp1,temp2);
		//beta(5)
		temp1=beta[i][2]+gama[i][2];
		temp2=beta[i][6]+gama[i][1];
		beta[i-1][5]=max_star(temp1,temp2);
		//beta(6)
		temp1=beta[i][3]+gama[i][3];
		temp2=beta[i][7]+gama[i][0];
		beta[i-1][6]=max_star(temp1,temp2);
		//beta(7)
		temp1=beta[i][3]+gama[i][0];
		temp2=beta[i][7]+gama[i][3];
		beta[i-1][7]=max_star(temp1,temp2);

		//归一化
		for(j=0;j<8;j++)
		{
			beta[i-1][j]=beta[i-1][j]-temp_max[i-1];
			if (beta[i-1][j]>Infty_1)
			{
				beta[i-1][j]=Infty_1;
			}
			else if(beta[i-1][j]<-Infty_1)
			{
				beta[i-1][j]=-Infty_1;
			}
		}
		
	}
		delete []temp_max;

	//计算输出的对数似然比LLR
	int pm_b_1[8]={4,0,1,5,6,2,3,7};
	int pm_b_0[8]={0,4,5,1,2,6,7,3};
	int bm_1[8]={3,3,2,2,2,2,3,3};
	int bm_0[8]={0,0,1,1,1,1,0,0};
	double temp_star1,temp_star0;

	for(i=0;i<block_size;i++)
	{
		temp1=0;
		temp2=0;
		temp_star1=alpha[i][0]+beta[i][4]+gama[i][3];
		temp_star0=alpha[i][0]+beta[i][0]+gama[i][0];
		for(j=1;j<8;j++)
		{
			//计算等于1的
			temp1=temp_star1;
			temp2=alpha[i][j]+beta[i][pm_b_1[j]]+gama[i][bm_1[j]];
			temp_star1=max_star(temp1,temp2);
			//计算等于0的
			temp1=temp_star0;
			temp2=alpha[i][j]+beta[i][pm_b_0[j]]+gama[i][bm_0[j]];
			temp_star0=max_star(temp1,temp2);
		}
		llr[i]=temp_star1-temp_star0;
	}
}

void decision(double *indata)
{
	int i;
	for(i=0;i<block_size;i++)
	{
		if(indata[i]>=0)
		{
			indata[i]=1;
		}
		else
		{
			indata[i]=0;
		}
	}
}

double max_star(double a,double b)
{
	double temp_max,temp_diff,lut;
	temp_max=max(a,b);
	temp_diff=fabs(a-b);
    double	temp_bit=2^5;

	lut=log(1+exp(-temp_diff/temp_bit))*temp_bit;//temp_bit有什么用？
	return temp_max+lut;
}
	
/*void random(int * data,int len)
{
	srand((int)time(NULL));
	for(int i=0;i<len;i++)
		data[i] = rand()%2;
}*/


double random()
{
	long z,k;
	static long s1=12345L;
	static long s2=12345473464L;
	
	k= s1 / 53668L;
	s1 = 40014L * (s1 - k*53668L) - k*12211L;
	if (s1<0)
		s1 = s1 + 2147483563L;
	k = s2 / 52774;
	s2 = 40692L * (s2 - k*52774L) - k*3791L;
	if (s2<0)
        s2 = s2 + 2147483399L;
	z=s1 - s2;
	if (z<1)
		z = z + 2147483562L;
	return (double) z / (double) 2147483563.0;
}

short generate_binary_source()
{
	if(random()<.5) return 0;
	else return 1;
}

/* Noise generator
parameters : sigma2 = variance (global variable)
generate a gaussian distributed sequence with zero mean
*/
float AWGN(double sigma2)
{
	static int iset=0;
	static double gset;
	double fac,r,v1,v2;
	{	
		do
		{
			v1=2.0*random()-1.0;
			v2=2.0*random()-1.0;
			r=v1*v1+v2*v2;
		}
		while (r>=1.0 || r == 0.0);
		fac=sqrt(-2.0*sigma2*log(r)/r);
		gset = v1 * fac;
		iset = 1;
		return (float)( v2*fac);
	}
}




void main()
{				
	int *data = new int[block_size];
	int i;
	SNR_db=50.0;
	for(SNR_db=0.0;SNR_db<=5.0;SNR_db++)
	{  
		int err=0;
		int number=0; 
		SNR=(long double)pow(10.0,SNR_db/10.0);     
		double sigma=1/(2*rate*SNR);
		do
		{		
			zero(data,block_size);
			for(i=0;i<block_size;i++)
			{
				data[i]=generate_binary_source();
			}
//			random(data,block_size);
			creat_interleaver();
//			print(data,block_size);			
			encode(data);
			//	print(rx_data,3*block_size);
			
/*			for(i=0;i<3*block_size;i++)
			{				
				rx_data[i]=(double)(rx_data[i]+AWGN(sigma));				
			}
*/
/*			for(i=0;i<3*block_size;i++)
			{
					cout<<rx_data[i]<<"  ";
			}
			cout<<endl;*/
			
			decode_init();
			turbo_decode();	

			//cout<<"输出：";
//		    print(de_data,block_size);

			for(i=0;i<block_size;i++)
			{
				err+=data[i] ^ de_data[i];
			}
//			cout<<err<<"  "<<endl;
			number+=block_size;
		}while(number<=10000);				
		
		double err_ratio=(long double)err/number;
		cout<<"SNR_db ="<<SNR_db<<"-------"<<err_ratio<<endl;
		ofstream myf("d:\\WCDMA.txt", ios::ate);
		myf<<SNR_db<<"\t"<<err_ratio<<endl;
	}
	
	for(SNR_db=5.0;SNR_db<=6.0;SNR_db+=0.1)
	{  
		int err=0;
		int number=0; 
		SNR=(long double)pow(10.0,SNR_db/10.0);
		double sigma=1/(2*rate*SNR);
		do
		{		
			zero(data,block_size);
			for(i=0;i<block_size;i++)
			{
				data[i]=generate_binary_source();
			}

			creat_interleaver();
		
			encode(data);
			
			for(i=0;i<3*block_size;i++)
			{				
				rx_data[i]=(double)(rx_data[i]+AWGN(sigma));				
			}
			
			decode_init();
			turbo_decode();	

			for(i=0;i<block_size;i++)
			{
				err+=data[i] ^ de_data[i];
			}
			number+=block_size;
		}while(number<=10000);				
		
		double err_ratio=(long double)err/number;
		cout<<"SNR_db ="<<SNR_db<<"-------"<<err_ratio<<endl;
		ofstream myf("d:\\WCDMA.txt", ios::ate);
		myf<<SNR_db<<"\t"<<err_ratio<<endl;
	}
	delete []data;
	cin>>tmp;
}

void freemem()
{
	int i;
	delete []rx_s1;
	delete []rx_s2;
	delete []rx_p1;
	delete []rx_p2;
	delete []ext1;
	delete []ext2;
	for(i=0;i<block_size;i++)
	{
		delete alpha[i];
		delete beta[i];
		delete gama[i];
	}
	delete alpha;
	delete beta;
	delete gama;
	delete []intlv_array;
	delete []de_intlv_array;
}