encode.c

用于LDPC编码译码的仿真实现。包括随机生成校验矩阵、由校验矩阵产生生成矩阵、编码、加随机噪声、译码等内容。原作者是老外

/* ENCODE.C - Encode message blocks. */

/* Copyright (c) 2000, 2001 by Radford M. Neal 
 *
 * Permission is granted for anyone to copy, use, or modify this program 
 * for purposes of research or education, provided this copyright notice 
 * is retained, and note is made of any changes that have been made. 
 *
 * This program is distributed without any warranty, express or implied.
 * As this program was written for research purposes only, it has not been
 * tested to the degree that would be advisable in any important application.
 * All use of this program is entirely at the user's own risk.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "rand.h"
#include "alloc.h"
#include "mod2sparse.h"
#include "rcode.h"
#include "enc.h"
#include "encode.h"



/* MAIN PROGRAM. */

void encode
(
  char *pchk_file,
  char *gen_file,
  int *data_source,
  int *data_encoded
)
{
  int *chks;
  int i;

  /* Look at arguments. */
  /* Read parity check file */

  read_pchk(pchk_file);       //读校验矩阵，设置全局变量，H、M、N

  if (N<=M)
  { fprintf(stderr,
 "Can't encode if number of bits (%d) not greater than number of checks (%d)\n",
      N,M);
    exit(1);
  }

  /* Read generator matrix file. */

  read_gen(gen_file,0,0);    //读生成矩阵，设置全局变量cols、rows、L、U

  chks = chk_alloc (M, sizeof *chks);

  /* Encode successive blocks. */

  /* Compute encoded block. */
  sparse_encode (data_source, data_encoded);   //编码data-source

  /* Check that encoded block is a code word. H与编码比特向量相乘，结果存入chks*/
  mod2sparse_mulvec (H, data_encoded, chks);

  for (i = 0; i<M; i++) 
  { if (chks[i]==1)
     { fprintf(stderr,"Output block is not a code word!  (Fails check %d)\n",i);
       abort(); 
     }
  }

  mod2sparse_free(H); //lxh
  mod2sparse_free(L); //lxh
  mod2sparse_free(U); //lxh
 // free(cols);//lxh
 // free(rows);//lxh
  free(chks);
}



void huffman_tree(double *p, int size, HFTREE T[])
{
	int i,j;
	int lc,rc;
	double minwt1,minwt2;

	for(i=0;i<2*size-1;i++) //置初值
	{
		T[i].parent=-1;
		T[i].lchild=-1;
		T[i].rchild=-1;
		if (i<size)
			T[i].wt=p[i];
		else
			T[i].wt=0.0;
	}

	for (i=0;i<size-1;i++)
	{
		lc=0;
		rc=0;
		minwt1=1;
		minwt2=1;
		for(j=0;j<size+i;j++)
			if((T[j].wt<minwt2)&&(T[j].parent==-1))
			{minwt1=minwt2;rc=lc;minwt2=T[j].wt; lc=j;}
			else if((T[j].wt<minwt1)&&(T[j].parent==-1))
			{ minwt1=T[j].wt;rc=j;}

		T[lc].parent=size+i;
		T[rc].parent=size+i;
		T[size+i].wt=minwt1+minwt2;
		T[size+i].lchild=lc;
		T[size+i].rchild=rc;
	}
}
/*========================================================*/
//codelength为各码字的码长

void huffman(double *p,int size,int *huffman_code,int codelength[],HFTREE T[])
{
	int i,j,k,temp,totollength=0;

    huffman_tree(p,size,T);

	//Huffman编码
	for(i=0;i<size;i++)
	{
		j=i;

		while(j!=2*size-2)
			{
			k=T[j].parent;
			
			if(T[k].lchild==j)
				*(huffman_code+totollength)=1;
			else
				*(huffman_code+totollength)=0;

		    j=T[j].parent;
			codelength[i]++;
			totollength++;
			}

			for(j=0;j<codelength[i]/2;j++)
			{
				temp=*(huffman_code+totollength-codelength[i]+j);
				*(huffman_code+totollength-codelength[i]+j)=*(huffman_code+totollength-j-1);
				*(huffman_code+totollength-j-1)=temp;
			}
	}

}

//信源编码器
/*
* 函数介绍：对长度为symbollength的（0,1,...,size-1）的信源符号，信源编码输出(0,1)bit
* 输入参数：symbol,信源输出符号；symbollength，信源序列的长度
            code,对各信源的0，1编码,codelength，各信源的编码长度,size,符号数
* 输出参数：data，函数所产生对信源符号信源编码后的输出（0,1）序列
* 返回值：  无
*/
 
int source_code(int *symbol,int symbollength,int *code,int codelength[],int size,int *data)
{
	int i,j,bit_length=1;
	int *codeoffset; //信源符号在编码符号code的起始位置

   	if ((codeoffset=(int *)malloc(sizeof(int)*size))==NULL)
	{
	  printf("\n fail to allocate memory of codeoffset \n");
	  exit(1);  
	}

	*codeoffset=0;
	for(i=1;i<size;i++)
		*(codeoffset+i)=*(codeoffset+i-1)+codelength[i-1];

    for(i=1;i<symbollength+1;i++)
	{
		for(j=0;j<codelength[symbol[i]];j++)
		{
			*(data+bit_length)=*(code+codeoffset[symbol[i]]+j);
			bit_length++;
		}
	}
	free(codeoffset);
	return(bit_length-1);

}


/*
生成 N=16200
R=1/2的 LDPC 码
*/
void encode_DVB(int *source_data, int source_length, int *encoded_data,int encoded_length)
{
	int i,j,info; // 循环变量
	int H_rows,H_cols,row,col,pos;
	int **DVB_H;
	int q=90,c;
	int *Parity;
	FILE *PPP;
	int Addr_Check[25][8]={{20, 712, 2386,6354, 4061, 1062, 5045,-1},{5158,-1},
	{21, 2543, 5748, 4822, 2348, 3089, 6328,-1},{5876,-1},
	{22, 926, 5701, 269, 3693, 2438, 3190,-1},{3507,-1},
	{23, 2802, 4520, 3577, 5324, 1091, 4667,-1},{4449,-1},
	{24, 5140, 2003, 1263, 4742, 6497, 1185,-1},{6202,-1},
	{0, 4046, 6934,-1},{1, 2855, 66,-1},{2, 6694, 212,-1},{3, 3439, 1158,-1},
	{4, 3850, 4422,-1},{5, 5924, 290,-1},{6, 1467, 4049,-1},{7, 7820, 2242,-1},
	{8, 4606, 3080,-1},{9, 4633, 7877,-1},{10, 3884, 6868,-1},{11, 8935, 4996,-1},
	{12, 3028, 764,-1},{13, 5988, 1057,-1},{14, 7411,-1}};
	
	if ((PPP=fopen( "test_pos.txt", "a+" )) == NULL)
    {
	  printf("\nError! Can not open file PPP.txt\n");
	  exit(1);
    }

	H_cols = 16200;
	H_rows = 8100;


	if ((DVB_H = (int **) malloc(H_rows*sizeof(int *))) == NULL)
	{
		printf("Not enough memory to allocate buffer.\n");
		exit(1);
	}
	for (i=0;i<H_rows;i++)
	{
		if ((DVB_H[i]= (int *) malloc((H_cols - H_rows)*sizeof(int))) == NULL)
		{
			printf("Not enough memory to allocate buffer.\n");
			exit(1);
		}
	}

	for(i=0;i<H_rows;i++)
		for(j=0;j<(H_cols-H_rows);j++)
			DVB_H[i][j]=0;
	

	
	if ((Parity = (int *)malloc( H_rows*sizeof(int)))==NULL)
	{
	  printf("\n fail to allocate memory of Parity \n");
	  exit(1);  
	}



	/*set initial value for check bits*/
	for (i=0;i<H_rows;i++)
		Parity[i]=0;

	for (info=0; info<H_cols-H_rows;info++)
	{
		row = info/360;
		col = 0;
		while (Addr_Check[row][col]!=-1)
		{
			pos = (int)fmod(info,360)*q; // the position of check to accummlator
			pos+=Addr_Check[row][col]; 
			pos = (int)fmod(pos, H_cols-H_rows);
			fprintf(PPP,"%d  ",pos);
			Parity[pos] ^= source_data[info];
			DVB_H[pos][info] ^= 1;
			col++;
		}
		fprintf(PPP,"\n");		
	}

	for (i=1;i<H_rows;i++)
	{
		Parity[i] ^= Parity[i-1];
		for(j=0;j<(H_cols-H_rows);j++)
			DVB_H[i][j] ^= DVB_H[i-1][j];

	}

	for (i=0;i<H_cols-H_rows;i++)
		encoded_data[i]=source_data[i];
	for (i=0;i<H_rows;i++)
		encoded_data[H_cols-H_rows+i]=Parity[i];

	/*check whether the encoded_data is a really code*/

	c=0;
	for (i=0;i<H_rows;i++)
	{
		for(j=0;j<(H_cols-H_rows);j++)
			c+=DVB_H[i][j]*source_data[j];
		c+=Parity[i];
		c=(int)fmod(c,2);
		if(c!=0)
			printf("%d",-1);
	}
  
	free(Parity);
	for(i=0; i<H_rows;i++)
	{
		free(DVB_H[i]);
	}
	free(DVB_H);



	fclose(PPP);

}