rx_ldpc_dec.cpp

模拟cmmb系统中ldpc译码的部分

/************************************/
/*函数说明：
/*功能:LDPC编码
/*输出参数：
/*rx_ldpcdecoded_bit：经过LDPC解码后的bit流
/*输入参数：
/*rx_ldpc_llr：bit解交织后的初始化似然比；
/*rx_bitdeinterved_bit_len：经过bit解交织后的初始化似然比的长度；
/*ldpc_rate:LDPC码率；
/*备注：采用bp算法
/************************************/

/*检验译码是否正确，c为0则正确*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rx_ldpc_dec.h"
#include "alloc.h"
#include "mod2sparse.h"

int check
( mod2sparse *H,	/* Parity check matrix */
  char *dblk,		/* Guess for codeword */
  char *pchk		/* Place to store parity checks */
)
{
	int m, i, c;
	m = mod2sparse_rows(H);  
	mod2sparse_mulvec (H, dblk, pchk);		/*    pchk=H x dblk   */
	c = 0;
	for (i = 0; i<m; i++) 
	{ 
		c += pchk[i];
	}
	return c;
}



char *initprp
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  char *dblk		/* Place to store decoding */
)
{ 
	mod2entry *e;
	int n;
	int j;

	n = mod2sparse_cols(H);

	for (j = 0; j<n; j++)
	{ 
		for (e = mod2sparse_first_in_col(H,j);   /* Returns the first entry in column j of H */
			!mod2sparse_at_end(e);   			 /* Returns 1 if e is a special entry obtained 
                                						by moving past the end, returns 0 otherwise */
			e = mod2sparse_next_in_col(e))
		{
			e->pr = lratio[j];                                    /*probabilities*/
			e->lr = 0;							/*Likelihood ratios*/
		}
    dblk[j] = lratio[j]<0;
	}
	return dblk;
}

  
char *iterprp
( mod2sparse *H,
  double *lratio,	
  char *dblk
)
{
  double pr, temp;
  double min[2];  
  mod2entry *e;
  int N, M;
  int i, j,neg_num;
  double alpha=0.8;

  M = mod2sparse_rows(H);
  N = mod2sparse_cols(H);

  for (i = 0; i<M; i++)
  {
		neg_num=0;
		e = mod2sparse_first_in_row(H,i);
		if(e->pr<0){
    		neg_num++;	
    	}
    min[0]=fabs(e->pr);
    e = mod2sparse_next_in_row(e);
		if(e->pr<0){
    		neg_num++;	
    	}
    min[1]=fabs(e->pr);
    if(min[0]<min[1])   //最小值存在min[1],次小值存在min[0]
    	{
    		temp=min[0];
    		min[0]=min[1];
    		min[1]=temp;
    	}
    for (e = mod2sparse_next_in_row(e);
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_row(e))
    { 
    	if(e->pr<0){
    		neg_num++;	
    	}	
    	if(fabs(e->pr)<min[0]){
    			min[0]=fabs(e->pr);
    			if(min[0]<min[1])   //最小值存在min[1],次小值存在min[0]
    				{
    					temp=min[0];
    					min[0]=min[1];
    					min[1]=temp;
    					}
    		
    		}
    	
    }
    for (e = mod2sparse_last_in_row(H,i);
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))
    { 
    	if(neg_num%2==0){
    		if(e->pr>=0){
    			if(e->pr!=min[1])
    				e->lr=min[1]*alpha;
    			else
    				e->lr=min[0]*alpha;
    			}
    		else{
    			if(fabs(e->pr)!=min[1])
    				e->lr=-min[1]*alpha;
    			else
    				e->lr=-min[0]*alpha;
    			}
    	}
    	else{
    		if(e->pr>=0){
    			if(e->pr!=min[1])
    				e->lr=-min[1]*alpha;
    			else
    				e->lr=-min[0]*alpha;
    			}
    		else{
    			if(fabs(e->pr)!=min[1])
    				e->lr=min[1]*alpha;
    			else
    				e->lr=min[0]*alpha;
    			}
    	}
    				
    }
  }

  /* Recompute probability ratios.  Also find the next guess based on the
     individually most likely values. */

  for (j = 0; j<N; j++)
  { pr = lratio[j];
    for (e = mod2sparse_first_in_col(H,j);
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_col(e))
    { 
      pr += e->lr;
    }
    dblk[j] = pr<0;    
    for (e = mod2sparse_last_in_col(H,j);
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_col(e))
    { 
		e->pr = pr-e->lr;
    }
  }
  return dblk;  
}   


/*char *iterprp
( mod2sparse *H,	
  double *lratio,	
  char *dblk	
)
{
	double pr, dl, t;
	mod2entry *e;
	int n, m;
	int i, j;

	m = mod2sparse_rows(H);
	n = mod2sparse_cols(H);

	for (i = 0; i<m; i++)
	{ 
		dl = 1;
		for (e = mod2sparse_first_in_row(H,i);
			 !mod2sparse_at_end(e);
			 e = mod2sparse_next_in_row(e))
		{ 
			e->lr = dl;
			dl *= 2/(1+e->pr) - 1;
		}
		dl = 1;
		for (e = mod2sparse_last_in_row(H,i);
			 !mod2sparse_at_end(e);
			 e = mod2sparse_prev_in_row(e))
		{
			t = e->lr * dl;
			e->lr = (1-t)/(1+t);
			dl *= 2/(1+e->pr) - 1;
		}
	}
*/
  /* Recompute probability ratios.  Also find the next guess based on the
     individually most likely values. */

/*	for (j = 0; j<n; j++)
	{ 
		pr = lratio[j];
		for (e = mod2sparse_first_in_col(H,j);
			!mod2sparse_at_end(e);
			e = mod2sparse_next_in_col(e))
		{
			e->pr = pr;
			pr *= e->lr;
		}

	    dblk[j] = pr>=1;    
		pr = 1;
		for (e = mod2sparse_last_in_col(H,j);
			!mod2sparse_at_end(e);
			e = mod2sparse_prev_in_col(e))
		{
			e->pr *= pr;
			pr *= e->lr;
		}
	} 
	return dblk;  
}  */ 


/*sr3里存似然比，sr2里存码字, irs1里存行坐标，jc里存列坐标, dblk里返回译码结果*/
void dec( mod2sparse *H,double *lratio,char* dblk,int M)								
{ 
	char *pchk;
	int c,n;  
	int max_iter = 100; 
	int Num_OUT,Suc_OUT;   
	pchk = (char *)calloc(M,sizeof(char));/*检查译码是否退出时存储*/

	/*------------decode---------------*/    
	
	initprp(H,lratio,dblk);          /*初始化*/ 
	
	for (n = 0; ; n++)
	{	 
		c = check(H,dblk,pchk);		   
		if (n<max_iter && c==0)
		{		 		 	 			  		 	
			Suc_OUT=1;		   	 		
			Num_OUT=n+1;
			printf("%d,sucess\n",n);
			break; 				   	 		
		}
		   	 
		if (n==max_iter )
		{				   	 
	  		Suc_OUT=0;	  	  	
	   		Num_OUT=n;	   		 	
	   		break; 			   		 	
		}
		
		iterprp(H,lratio,dblk);/*译码*/	
	}
		  		
//	printf("success = %d,inum = %d\n",Suc_OUT,Num_OUT);/*输出译码信息，这个地方可能要改进*/

/* free the memory */	
	
	free(pchk);
	return;
}


void rx_ldpc_dec(double * rx_ldpc_llr, double ldpc_rate, char * rx_ldpcdecoded_bit) 
{	
	int i,j,k,M,K,iii;
	
	int *jcs1,*irs1;     //注意溢出，要动态分布

	char dblk[LDPC_CHECK_MATRIX_N];
	
	char temp[LDPC_CHECK_MATRIX_N];
	
	int order[LDPC_CHECK_MATRIX_N];
  
	mod2sparse *H;
	
	FILE *fp;			     
 	
	K = (int)(LDPC_CHECK_MATRIX_N*ldpc_rate);
	M = LDPC_CHECK_MATRIX_N-K;
	H = mod2sparse_allocate(M,LDPC_CHECK_MATRIX_N);    
	k=0;
	
	irs1 = (int *)calloc(LDPC_CHECK_MATRIX_nz ,sizeof(int));
	
	jcs1 = (int *)malloc(LDPC_CHECK_MATRIX_nz *sizeof(int *));
	
	/*对两种模式进行判断*/
	switch (M)
	{	
    case 4608:
    
    		/*读稀疏矩阵H*/    
    		
    fp=fopen("jcs1.dat","rb");	// H given by std
    	if(fp==NULL)
		{
			printf("jcs1 can not open the file!\n");
			return;
		}				
		fread(jcs1,sizeof(int),LDPC_CHECK_MATRIX_nz,fp);	
		fclose(fp);	
    		
				
		fp=fopen("irs1.dat","rb");    		
    	if(fp==NULL)
		{
			printf("irs1 can not open the file!\n");
			return;
		}				
		fread(irs1,sizeof(int),LDPC_CHECK_MATRIX_nz ,fp);				
		fclose(fp);
		
		fp = fopen("order12.dat","rb");  //读码字重排顺序
    			
		if(fp==NULL)
		{
			printf("can not open order12!\n");					
			return ;
		}									
		fread(order,sizeof(int),LDPC_CHECK_MATRIX_N,fp);					 	
		fclose(fp);

		for(i=0;i<LDPC_CHECK_MATRIX_nz ;i++)   
    		mod2sparse_insert(H,irs1[i],jcs1[i]);   /*矩阵的非零位置插入entry*/   
		
		/*-------对每一帧译码------*/   		

    dec(H,rx_ldpc_llr,dblk,M);  
		for(iii=0;iii<LDPC_CHECK_MATRIX_N;iii++)
			temp[iii] = dblk[order[iii]];
		for(j=0;j<K;j++)	
			rx_ldpcdecoded_bit[j]=temp[j+M]; /*译码结果存在rs_encoded->data中，用来输出*/

    break;
        
    case 2304: 		
    		
    	fp=fopen("jcs2.dat","rb");	// H given by std
    	if(fp==NULL)
		{
			printf("jcs2 can not open the file!\n");
			return;
		}				
		fread(jcs1,sizeof(int),LDPC_CHECK_MATRIX_nz,fp);	
		fclose(fp);	
    		
				
		fp=fopen("irs2.dat","rb");    		
    	if(fp==NULL)
		{
			printf("irs2 can not open the file!\n");
			return;
		}				
		fread(irs1,sizeof(int),LDPC_CHECK_MATRIX_nz ,fp);				
		fclose(fp);
		
		fp = fopen("order34.dat","rb");
    		
		if(fp==NULL)
		{
			printf("can not open order34!\n");					
			return ;
		}									
		fread(order,sizeof(int),LDPC_CHECK_MATRIX_N,fp);					 	
		fclose(fp);

		
		for(i=0;i<LDPC_CHECK_MATRIX_nz ;i++)   
    		mod2sparse_insert(H,irs1[i],jcs1[i]);   /*矩阵的非零位置插入entry*/ 


		/*对每一帧译码*/
    		
    		dec(H,rx_ldpc_llr,dblk,M);   			
    		for(iii=0;iii<LDPC_CHECK_MATRIX_N;iii++)
					temp[iii] = dblk[order[iii]];
    		for(j=0;j<K;j++)	
    			rx_ldpcdecoded_bit[j]=temp[j+M]; /*译码结果存在rs_encoded->data中，用来输出*/
		
        break;
        
    default:
        printf("LDPC_RATE PARAMETER ERROR IN LDPC_DEC");
        return;    
	}
	
	mod2sparse_free(H);
	free(irs1);
	free(jcs1);
}