max-log-map.cpp

max-log-map算法的C语言的实现

//Turbo码：max-log-map 译码算法 
#include "CTC.h"	//	包含头文件
//using namespace std;
//ofstream ofileExtri ( "Extri.txt", ios::app );
//ofstream ofileLLR ( "LLR.txt", ios::app );
/********************************************************************
功能:	为每一个待译码bit计算其LLR (Max-log-Map) 
Input:	r_s[] 接收的系统位bit失量
		r_p1[]，r_p2[] 分别为接收的校验位bit失量.
		prior[]:先验概率　(计算Gama时，每次迭代后先验概率都改变)
		sgmaSquare:	不同SNR时值不同.
Output: LLR[].	
*********************************************************************/
void computeLLR( double *r_s, double *r_p1, double *prior, double sgmaSquare, double *LLR )
//void computeLLR( double *r_s, double *r_p1, double *prior, double A, double *LLR )
{
	int i,j;
	double tempx,tempy;
	double sum_temp1, sum_temp0;


	//alfa[]和beta[]以交织长度为单位，再以s0,s1..,s7为次序依次存储
	//其中gama先存input=1后存input=0;	结果：由当前状态和input决定
	double *gama = new double [2 * n_states * interLength];//2*8*interLength；
	double *alfa = new double [n_states * (interLength+1)];//8*(interLength+1)；
	double *beta = new double [n_states * (interLength+1)];//8*(interLength+1)；

	double *tempmax = new double [ interLength+1 ];

	for ( i = 0; i < interLength+1; i++ )
		tempmax[i] = 0;
		//tempmax[i]=minusInfinity;

	//Initialization of alfa and beta，只对t=0和t=interLength时初试化
	//零状态时
	alfa[0] = 0;//1
	beta[interLength] = log(1.0/8.0);//1;//0
	//非零状态时
	for ( i = 1; i < n_states; i++ )
	{
		alfa[i * (interLength+1)] = minusInfinity;//0
		beta[i * (interLength+1) + interLength] = log(1.0/8.0);//0.125;// log(1/8);//minusInfinity;
		//ofileLLR << alfa[i * (interLength+1)] << beta[i * (interLength+1) + interLength]<<endl;
	}// End-Initialization of alfa and beta 

	//exit(1);
	
	//当初始状态分别为s0,s1...s7, and input = 1; 数组中元素为output,参考Trellis图
	int Trellis_1[n_states][2] = {{1,1}, {1,1}, {1,-1}, {1,-1}, {1,-1}, {1,-1}, {1,1}, {1,1}};
	//input = 0; 数组中元素为output
	int Trellis_0[n_states][2] = {{-1,-1}, {-1,-1}, {-1,1}, {-1,1}, 
	{-1,1}, {-1,1}, {-1,-1}, {-1,-1}};
	
	// compute gama,for t = 1,2,...,interLength; s = s0, s2,...,s7;
	// t 为时刻
	for (i = 0; i < interLength; i++)			
	{
		for (j = 0; j < n_states; j++)		//分别为s0,s1...s7
		{
			// 起始状态，and input=1,
			gama[(j*interLength*2) + 2*i+0] =( r_s[i]*Trellis_1[j][0] 
				+ r_p1[i]*Trellis_1[j][1] )/sgmaSquare  + 0.5*prior[i]; 
			//gama[(j*interLength*2) + 2*i+0] = exp( 0.5*prior[i] + 
				//( r_s[i]*Trellis_1[j][0] + r_p1[i]*Trellis_1[j][1] )/sgmaSquare );

			// 起始状态，and input=0
			//gama[(j*interLength*2) + 2*i+1] = exp( -0.5*prior[i] + 
				//( r_s[i]*Trellis_0[j][0] + r_p1[i]*Trellis_0[j][1] )/sgmaSquare );
			gama[(j*interLength*2) + 2*i+1] =( r_s[i]*Trellis_0[j][0] 
				+ r_p1[i]*Trellis_0[j][1] )/sgmaSquare -0.5*prior[i];
			//ofileLLR << gama[(j*interLength*2) + 2*i+0] << endl;
			
			//if(i==1000) exit(1);
			//if(gama[(j*interLength*2) + 2*i+0] < 1e-300)
			//{
			//	ofileLLR << "i = " <<i <<"j = " <<j<<",error";
			//	exit(1);
			//}
		
		}
	}//End-compute gama
    		
	//compute Alfa, for t = 1,2,...,interLength; s = s0, s2,...,s7
	//?-计算alfa似乎只需：alfa[0]到alfa[interLength-1]	
	int Forward_trellis[n_states][2] = { {0,1}, {2,3}, {4,5}, {6,7}, {0,1}, {2,3}, {4,5}, {6,7} };
	//第1,3个为状态，第2,4个元素：1表示input=0, 0表示input = 1 
	int mapA[n_states][4] = {{0,1,1,0}, {2,0,3,1}, {4,1,5,0}, {6,0,7,1}, {0,0,1,1}, {2,1,3,0}, 
	{4,0,5,1}, {6,1,7,0}};

	for (i = 1; i < interLength + 1; i++) //共interLength个
	{
		for (j = 0; j < n_states; j++)		//分别为s0,s1...s7
		{
			tempx = alfa[ Forward_trellis[j][0]*(interLength+1) + (i-1) ] 
				+ gama[ mapA[j][0]*2*interLength + 2*(i-1) + mapA[j][1] ];		
			tempy = alfa[ Forward_trellis[j][1]*(interLength+1) + (i-1) ] 
				+ gama[ mapA[j][2]*2*interLength + 2*(i-1) + mapA[j][3] ];

			//tempx = alfa[ Forward_trellis[j][0]*(interLength+1) + (i-1) ] * 
			//	gama[ mapA[j][0]*2*interLength + 2*(i-1) + mapA[j][1] ];		
			//tempy = alfa[ Forward_trellis[j][1]*(interLength+1) + (i-1) ] * 
			//	gama[ mapA[j][2]*2*interLength + 2*(i-1) + mapA[j][3] ];

			//alfa[j*(interLength+1) + i] = tempx + tempy;

			//if ( exp(tempx + tempy) < 1e-100 )
			//	alfa[j*(interLength+1) + i] = minusInfinity;
			//else
				
			alfa[j*(interLength+1) + i] = tempx > tempy ? tempx : tempy;
						
			//ofileLLR << alfa[j*(interLength+1) + i] << endl;
			//if(i==1000) 
			
		}
		/*-------------------------------------------添加的
		double sumalfa=0;
		for (j = 0; j < n_states; j++)
		{
			sumalfa += alfa[j*(interLength+1) + i];
		}
		for (j = 0; j < n_states; j++)	
		{	
			alfa[j*(interLength+1) + i] /= sumalfa;


		}*/
		
			//{
			//	ofileLLR << "i = " <<i <<"j = " <<j<<",error";
			//	exit(1);
		//--------------------------------------------

	}//End- compute Alfa
	//exit(1);

	////End- compute Beta, for t = 1,2,...,interLength; s = s0, s2,...,s7;
	int Bacward_trellis[n_states][2] = {{0,4}, {0,4}, {1,5}, {1,5}, {2,6}, {2,6}, {3,7}, {3,7}};
	int mapB[n_states][2] = {{1,0}, {0,1}, {0,1}, {1,0}, {1,0}, {0,1}, {0,1}, {1,0}};//1对应input = 0

	for ( i = interLength - 1; i >= 0; i-- ) //共interLength个
	{
		for (j = 0; j < n_states; j++)		//分别为s0,s1...s7
		{
			tempx = beta[ Bacward_trellis[j][0]*(interLength+1) + (i+1) ] 
				+ gama[ j*2*interLength + 2*i + mapB[j][0] ];				
			tempy = beta[ Bacward_trellis[j][1]*(interLength+1) + (i+1) ] 
				+ gama[ j*2*interLength + 2*i + mapB[j][1] ];

			//tempx = beta[ Bacward_trellis[j][0]*(interLength+1) + (i+1) ] * 
				//gama[ j*2*interLength + 2*i + mapB[j][0] ];				
			//tempy = beta[ Bacward_trellis[j][1]*(interLength+1) + (i+1) ] * 
				//gama[ j*2*interLength + 2*i + mapB[j][1] ];

			//beta[j*(interLength+1) + i] = tempx + tempy; 

		
			/*ofileLLR << beta[ Bacward_trellis[j][1]*(interLength+1) + (i+1) ]<<
				","<<gama[ j*2*interLength + 2*i + mapB[j][0] ] << endl;
			if(i==3160-5) exit(1); */
			
			//if ( exp(tempx + tempy) < 1e-100 )
				//beta[j*(interLength+1) + i] = minusInfinity;
			//else
				
			beta[j*(interLength+1) + i] = tempx > tempy ? tempx : tempy; 
			//ofileLLR <<beta[j*(interLength+1) + i]<<endl;

		}
		
		/*-------------------------------------------添加的
		double sumbeta=0;
		for (j = 0; j < n_states; j++)		
			 sumbeta += beta[j*(interLength+1) + i];
		

		//-------------------------------------------
		for (j = 0; j < n_states; j++)
		{
			beta[j*(interLength+1) + i] /= sumbeta;
			//ofileLLR << beta[j*(interLength+1) + i] << endl;
			//if(i==1000) exit(1);
		}*/

	}//End-compute Beta
	//exit(1);

	//compute LLR
	double temp1[n_states];	
	double temp0[n_states];
	int map1[n_states] = {4,0,1,5,6,2,3,7};
	int map0[n_states] = {0,4,5,1,2,6,7,3};
	for ( i = 0; i < interLength; i++ ) //共interLength个
	{
		sum_temp1 = 0;
		sum_temp0 = 0;
		for ( j = 0; j < n_states; j++ )	//分别为s0,s1...s7
		{
			temp1[j] = alfa[j*(interLength+1) + i] + gama[j*2*interLength + (2*i+0)] 
			+ beta[map1[j]*(interLength+1) + i+1]; //input=1
			temp0[j] = alfa[j*(interLength+1) + i] + gama[j*2*interLength + (2*i+1)] 
			+ beta[map0[j]*(interLength+1) + i+1];
			//temp1[j] = alfa[j*(interLength+1) + i] * gama[j*2*interLength + (2*i+0)] 
			//	* beta[map1[j]*(interLength+1) + i+1]; //input=1
			//temp0[j] = alfa[j*(interLength+1) + i] * gama[j*2*interLength + (2*i+1)] 
			//	* beta[map0[j]*(interLength+1) + i+1]; //input=0
		
			//ofileLLR << alfa[j*(interLength+1) + i]<<","<< gama[j*2*interLength + (2*i+0)]
				//<<","<< beta[map1[j]*(interLength+1) + i+1]<< endl;
		}
		/*for ( j = 0; j < n_states; j++ )
		{
			sum_temp1 += temp1[j];
			sum_temp0 += temp0[j];
		}*/
		//if( sum_temp1<1e-300)
			//LLR[i] = minusInfinity;
		//else if (sum_temp0<1e-300) 
		//	LLR[i] = -minusInfinity;
		//else 
			
		//LLR[i] = log(sum_temp1/sum_temp0);
		
		LLR[i] = getMax( temp1, n_states ) - getMax( temp0, n_states );
		//ofileLLR << LLR[i] << endl; 
		//exit(1);
		//
	}
	//ofileLLR << endl;
	//exit(1);
	delete []tempmax;

	delete []gama;
	delete []alfa;
	delete []beta;
}

/******************************************************************
功能:	得到一个序列中的最大值无素.
INPUT:data_seq - 
	      length - 序列长度.
RETURN VALUE: 返回最大值元素.
******************************************************************/
double getMax( double *dataSeq, int length )//本程序中为状态的个数,that is length=8=n_states;
{
	int i;
	double temp;
	temp = dataSeq[0];
	for ( i = 1; i < length; i++ )
	{
		if ( temp < dataSeq[i] )//
			temp = dataSeq[i];
	}

	return temp;
}
/******************************************************************
功能:	hard decisions
INPUT:data[] - 实际上为第二个译码器最后一次迭代后的对数似然比序列
OUPUT: result[]  0/1 sequence.
******************************************************************/
void decision( double *data, int *result )
{
	int i;
	for ( i = 0; i < interLength; i++ ) 
	{
		if( data[i] > 0) //if( data[i] >= 0)
			result[i] = 1;		
		else 
			result[i] = 0;
	}
}
/**************************************************************************************
功能:	计算外信息for each decoded bit.
Input:	LLR[]
		priorInfor[]先验信息
		r_s[]: 系统位信息
		sgmaSquare:	
Output: Extrinsic[] 由当前DECODER产生, 将作为先验信息用于下一个DECODER
****************************************************************************************/
//void getExtrinsic( double *LLR, double *priorInfor, double *r_s, double A, double *Extrinsic )
void getExtrinsic( double *LLR, double *priorInfor, double *r_s, double sgmaSquare, double *Extrinsic )
{
	int i = 0;
	for ( i = 0; i < interLength; i++ )
	{
		//add
		Extrinsic[i] = ( LLR[i] - r_s[i] * (2/sgmaSquare) - priorInfor[i] )*0.7;//
		//Extrinsic[i] =(LLR[i] - r_s[i] - priorInfor[i]) * 0.7;//
		//ofileExtri << LLR[i]<<","<< priorInfor[i] << endl; 
		/*if ( i==10 ) 
		{
			ofile1 << Extrinsic[i]<<","<<LLR[9]<<","<<priorInfor[9]<<endl;
		}*/
	}
	//exit(1);
	//ofileExtri << endl;
}