dec_ldpc.cpp

CMMB中的LDPC译码程序

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
//// 本程序用于CMMB标准中码率R=1/2的LDPC码译码的误码率评估
//// R=1/2的LDPC码，码长n＝9216，k=4608，
//// 校验矩阵H为规则矩阵，每行重量nrou = 6，每列重量ngamma = 3
//// 校验矩阵H每行中元素1的位置在协议的附录中给出
//// =====================================================
////////////////////////////////////////////////////////////////////////////////////////
double dec_LDPC_MSA(float* fdata_in,int* f_orig_data,int nH_rowdef[4608][6],int nH_coldef[9216][3],int nrowH,int ncolH,int nrou,int ngamma,int niter_max)
{
/////////////////////////////////////////////////////////
//输入参数
//输入待译码数据y,
//噪声能量dsigma2,
//H矩阵每行不为1的列H_rowdef,
//H矩阵每列不为1的行H_coldef,
//H矩阵的列数nrowH,
//H矩阵的行数ncolH,
//最大校验顶点度, 每行中1的个数nrou,
//最大码元顶点度, 每列中1的个数ngamma,
//最大循环次数niter_max
//////////////////////////////////////////////////////////
	
int i, j, k, l;
double LLR_pi[9216] = {0};
double LLR_qi[9216] = {0};
int result_z[9216] = {0};
double * LLR_rji = (double *)malloc(4608*9216*sizeof(double));
double * LLR_qji = (double *)malloc(4608*9216*sizeof(double));
int s;	
double sij;
double Aij;
double biterror[30] = {0};
int pos_nozero;
int j_col;
int i_row;
int col_nor_k = 0;
int row_nor_k = 0;
double ber;


for (i = 0;i <= 4608*9216-1;i++)
{
	*(LLR_rji+i) = 0;
	*(LLR_qji+i) = 0;
}

for (i = 0;i<= 9216-1;i++)
{
	for (j = 0;j <= 2;j++)
	{
		nH_coldef[i][j] = nH_coldef[i][j] - 1;
	}
}


for (i = 0; i <= 9216-1; i++)
{
	LLR_pi[i] = fdata_in[i];
	
	if (LLR_pi[i] > 0)
	{
       result_z[i] = 1;
	}
	else
	{
		result_z[i] = 0;
	}
}


for (i = 1; i <= niter_max; i++)
{
	for (j=0;j<=4608-1;j++)
	{
		s = 0;
		for (k=0;k<= nrou-1;k++)
		{
			pos_nozero = nH_rowdef[j][k];
			s = s + result_z[pos_nozero];
			s = s%2;
		}
		if (s == 1)
		{
			break;
		}
	}

	if (s == 1)
	{
	    for (j=0;j<=9216-1;j++)
		{
           for (k=0;k <= ngamma-1;k++)
		   {
			   j_col = nH_coldef[j][k];
			  *(LLR_qji + j_col*9216 + j) = LLR_pi[j];//LLR_qji[j_col][j]=LLR_pi[j];
			   for (l=0;l<= ngamma-1;l++)
			   {
				   if (l != k)
				   {
                      col_nor_k = nH_coldef[j][l];
                     *(LLR_qji + j_col*9216 + j) = *(LLR_qji + j_col*9216 + j) + *(LLR_rji + col_nor_k*9216 + j); // LLR_qji[j_col][j] = LLR_qji[j_col][j] + LLR_rji[nH_coldef[j][l]][j];
				   }
			   }
		   }
		}

		for (j=0;j<= 4608-1;j++)
		{
			for (k=0;k<= nrou-1;k++)
			{
				i_row  =nH_rowdef[j][k];
                Aij = pow(10,100);
                sij = 1.0;
				for (l=0;l <= nrou-1;l++)
				{
					if (l != k)
					{
                        row_nor_k = nH_rowdef[j][l];
                        if (*(LLR_qji+j*9216+row_nor_k) < 0)
						{ 
                            sij = -sij;
							if ((-*(LLR_qji+j*9216+row_nor_k)) < Aij)
							{
								Aij = -*(LLR_qji+j*9216+row_nor_k);//if (abs(LLR_qji(j,i_temp))<Aij);Aij = abs(LLR_qji(j,i_temp));
							}
                        }
						else
						{
							if (*(LLR_qji+j*9216+row_nor_k) < Aij)
							{
								Aij = *(LLR_qji+j*9216+row_nor_k);//if (abs(LLR_qji(j,i_temp))<Aij);Aij = abs(LLR_qji(j,i_temp));
							}
						}
					}
				}
				*(LLR_rji+j*9216+i_row) = pow(-1,nrou)*sij*Aij;  //LLR_rji[j][i_row] = ((-1)^nrou)*sij*fnphi_quan(Aij); 
			}
		}

         
        for (j=0;j<=9216-1;j++)
		{
			LLR_qi[j]=0;
		}


		for (j=0;j<=9216-1;j++)
		{
			for (k=0;k<=4608-1;k++)
			{
				LLR_qi[j]=LLR_qi[j] + *(LLR_rji+k*9216+j);//LLR_qi[j]=LLR_qi[j]+LLR_rji[k][j];
			}
			LLR_qi[j]= LLR_qi[j]+ LLR_pi[j];
			
		}

		for (j=0;j<=9216-1;j++)
		{	
            if (LLR_qi[j] > 0)
			{
                result_z[j]=1;
			}    
            else
			{
                result_z[j]=0;
			}

			if (result_z[j] != f_orig_data[j])
			{
				biterror[i-1] = biterror[i-1] + 1;
			} 
		}

		
	}
	else
	{
		for (j=0;j<=9216-1;j++)
		{	
            if (LLR_qi[j] > 0)
			{
                result_z[j]=1;
			}    
            else
			{
                result_z[j]=0;
			}
			
			if (result_z[j] != f_orig_data[j])
			{
				biterror[i-1] = biterror[i-1] + 1;
			} 
		}
		break;
	}
	printf("%lf \n",biterror[i-1]/9216);
}

ber = biterror[niter_max-1]/9216;
return ber;
free(LLR_rji);
free(LLR_qji);
}
/*----------------------------------------------------------------------------
FUNCTION main()
RETURNS  NONE

  PARAMETERS NONE.
		
		  PURPOSE
		  CMMB Simulation Platform main file.
		  
			WRITTEN BY
			CS
			
			  DATE:
			  2008-2-16  
----------------------------------------------------------------------------*/

void main()
{
	FILE *fp;
	int i,j;
	int nrowH = 4608;
	int ncolH = 9216;
    int nrou = 6;
	int ngamma = 3;
	int nR = 1/2;
	double nSNRb = 1.4; 
	int niter_max = 30;
	float fdata_in[9216] = {0};
	int f_orig_data[9216] = {0};
	double dsigma2 = 1/(2*0.5*(pow(10 ,(1.4/10))));
	int nH_rowdef[4608][6] = {0};
	int nH_coldef[9216][3] = {0};
	double derrorrate = {0};

	// 打开仿真/实测数据源
	//////////////////////////////////////////////////////////////////
	
	if ( (fp = fopen("dec_ldpc_data_in.txt","r")) == NULL)//读取待译码接收序列
	{
		printf("Open file \"dec_ldpc_data_in.txt\" error!");
			  exit(1);
	}
	else
	{
		for (i = 0; i <= ncolH - 1; i++)
		{
			fscanf(fp,"%f ",&fdata_in[i]);
			//printf("%f \n",fdata_in[i]);
		}
		fclose(fp);
	}
	//////////////////////////////////////////////////////////////////
	
	if ( (fp = fopen("ldpc_origin_data_in.txt","r")) == NULL)//读取原始数据序列
	{
		printf("Open file \"ldpc_origin_data_in.txt\" error!");
		exit(1);
	}
	else
	{
		for (i = 0; i <= ncolH - 1; i++)
		{
			fscanf(fp,"%d ",&f_orig_data[i]);
			//printf("%f \n",fdata_in[i]);
		}
		fclose(fp);
	}
	/////////////////////////////////////////////////////////////////////////////////		  
		  
	if ( (fp = fopen("H_12_t.txt","r")) == NULL)//读取H矩阵每行不为1的列
	{
		printf("Open file \"H_12_t.txt\" error!");
		exit(1);
	}
	else
	{
		for (i = 0; i <= nrowH - 1; i++)
		{
			for (j=0;j<=nrou - 1;j++)
			fscanf(fp,"%d",&nH_rowdef[i][j]);			
		}
		fclose(fp);//
	}
		  
 //////////////////////////////////////////////////////////////////////////		  
	if ( (fp = fopen("H_12_z.txt","r")) == NULL)//读取H矩阵每列不为1的行
	{
		printf("Open file \"H_12_z.txt\" error!");
		exit(1);
	}
	else
	{
		for (i = 0; i <= ncolH - 1; i++)
		{
			for (j=0;j<=ngamma - 1;j++)
			fscanf(fp,"%d",&nH_coldef[i][j]);			
		}
	    fclose(fp);
	}
////////////////////////////////////////////////////////////////////////////////////
	derrorrate = dec_LDPC_MSA(fdata_in,f_orig_data,nH_rowdef,nH_coldef,nrowH, ncolH,nrou,ngamma,niter_max);
	printf("f%",derrorrate);
}