make-ldpc.c

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

/* MAKE-LDPC.C - Make a Low Density Parity Check code's parity check matrix. */
//生成LDPC的校验矩阵
/* 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 "make-ldpc.h"
#include "mod2sparse.h"
#include "rand.h"
#include "rcode.h"
#include "alloc.h"
#include "intio.h"
#include "open.h"

/* CREATE A SPARSE PARITY-CHECK MATRIX.  Of size M by N, stored in H. */

void make_ldpc_l(char *pchk_file, double *lamita, int Dlmax,double *lao,int Drmax, int seed)
{
  FILE *file;
	mod2entry *e, *f, *g, *h;
  int added, elim4;
  int i, j, k, t;
  int col_max;						//total check bits (1s) in all cols
  int *rand_col;
  int *col_num;
  int temp;
  double col_avrg = 0.0;
  double temps = 0;

  struct rand_sp_parmt *col_sp_num;
  int row_max;		//total check bits (1s) in all rows
  int *rand_row;
  int *row_num;
  double row_avrg = 0.0;
  struct rand_sp_parmt *row_sp_num;	//
  int row_temp, col_temp;
  int total_chk;				//total check bits (1s) in sp_matrix
  int no4cycle = 0;
  int cb = 3;
    /* Check for some problems. */
    
  if (cb>M)
  { 
	  fprintf(stderr,
      "Number of checks per bit (%d) is greater than total checks (%d)\n", cb, M);
	  exit(1);
  }

//  no4cycle = 1;

  if (cb==M && N>1 && no4cycle)
  { 
	  fprintf(stderr, "Can't eliminate cycles of length four with this many checks per bit\n");
	  exit(1);
  }


	rand_col = chk_alloc(N, sizeof *rand_col);
	col_num =  chk_alloc (Dlmax+1, sizeof *col_num);
	col_sp_num =chk_alloc(N,sizeof *col_sp_num);

	for (i=2; i<Dlmax+1; i++)
	{
		col_avrg += lamita[i]/i;
	}
	col_avrg = 1/col_avrg; 
		

	//not repeat random int number from 0 to N-1
	rand_int_num( rand_col, N, seed);		


	col_num[0] = 0;
	col_max = 0;
	for (i=1; i<Dlmax+1; i++)
	{
		col_num[i] = (int)(floor((N*col_avrg*lamita[i]/i)+0.5)); //度数为i节点的列数
		col_max += col_num[i]; //总列数
	}

    if ((temp = N-col_max)!=0)
	{
		i = 0;
		for (j=0; j<Dlmax+1; j++)
		{
			if (col_num[j] != 0) 
			{
				col_num[j]++;
				if ((++i) == temp) break;
			}
		}
	}


	i = 0;
	col_max = 0;
	for (j=0; j<Dlmax+1; j++)
	{
		temp = col_num[j];
		if (temp>0)
		{
			do
			{
				col_sp_num[i].chk_num = j;				//store chekc number
				col_sp_num[i].pos_no = rand_col[i];		//store col number	
				i++;	
 				temp --;
			} while (temp);	
			col_max += col_num[j]*j;
		}
	}


	free(rand_col);
	free(col_num);


	rand_row = chk_alloc(M, sizeof *rand_row);
	row_num = chk_alloc (Drmax+1, sizeof *row_num);
	row_sp_num = chk_alloc(M,sizeof *row_sp_num);
	

	for (i=2; i<Drmax+1; i++)
	{
		row_avrg += lao[i]/i;
	}
	row_avrg = 1.0/row_avrg;
		

	//not repeat random int number from 0 to M-1
	rand_int_num( rand_row, M, seed);		

	row_num[0] = 0;
	row_max = 0;
	for (i=1; i<Drmax+1; i++)
	{
		row_num[i] = (int)(floor((M*row_avrg*lao[i]/i)+0.5));
		row_max +=row_num[i];
	}


   if ((temp = M-row_max)!=0)
	{
		i = 0;
		for (j=0; j<Drmax+1; j++)
		{
			if (row_num[j] != 0) 
			{
				row_num[j]++;
				if ((++i) == temp) break;
			}
		}
	}


   	i = 0;
	row_max = 0;
	for (j=0; j<Drmax+1; j++)
	{
		temp = row_num[j];
		if (temp>0)
		{
			do
			{
				row_sp_num[i].chk_num = j;				//store chekc number
				row_sp_num[i].pos_no = rand_row[i];		//store row number	
				i++;	
				temp --;
			} while (temp);	
			row_max += row_num[j]*j;
		}
	}

	free(rand_row);
	free(row_num);


  total_chk = (col_max<=row_max) ? col_max : row_max;	//get the smaller number of check bits 
//  printf("col_max=%d row_max=%d total_chk=%d\n", col_max, row_max, total_chk);

  chk_pos = chk_alloc(total_chk, sizeof *chk_pos);

  temp = 0;
  row_temp = 0;
  do
  {
	  i = rand_int(M, seed);
	  if ((row_sp_num[i].chk_num))
	  {
		  chk_pos[row_temp].row = row_sp_num[i].pos_no;
		  row_temp ++;
		  row_sp_num[i].chk_num --;
	  }
  } while (row_temp != total_chk);

  temp = 0;
  col_temp = 0;
  do
  {
	  i = rand_int(N, seed);
	  if ((col_sp_num[i].chk_num))
	  {
		  chk_pos[col_temp].col = col_sp_num[i].pos_no;
		  col_temp ++;
		  col_sp_num[i].chk_num --;
	  }
  } while (col_temp != total_chk);

  free(col_sp_num);
  free(row_sp_num);

  H = mod2sparse_allocate(M,N);
  mod2sparse_clear(H);

  /* Create the initial version of the parity check matrix. */

// make sp martrix acording the position matrix
  for (i=0; i<total_chk; i++)
  {
	  if (mod2sparse_find(H, chk_pos[i].row, chk_pos[i].col))
	  {
		  if (move_ps( chk_pos[i].row, chk_pos[i].col, seed) == 0)	//move check position.
		  {
			  printf(" can not moving postion row=%d, col=%d\n", chk_pos[i].row, chk_pos[i].col);
			  exit(0);
		  }
	  }
	  else
		  mod2sparse_insert(H, chk_pos[i].row, chk_pos[i].col);
  }
  free(chk_pos);


 /* Add extra bits to avoid rows with less than two checks. */
//spn:
  added = 0;

  for (i = 0; i<M; i++)
  { e = mod2sparse_first_in_row(H,i);
    if (mod2sparse_at_end(e))
    { j = rand_int(N, seed);
      e = mod2sparse_insert(H,i,j);
      added += 1;
    }
    e = mod2sparse_first_in_row(H,i);
    if (mod2sparse_at_end(mod2sparse_next_in_row(e)) && N>1)
    { do 
      { j = rand_int(N, seed); 
      } while (j==mod2sparse_col(e));
      mod2sparse_insert(H,i,j);
      added += 1;
    }
  }


  if (added>0)
  { fprintf(stderr,
           "Added %d extra bit-checks to make row counts at least two\n",
           added);
  }

  /* Add extra bits to try to avoid problems with even column counts. */

  if (cb%2==0 && cb<M && N>1 && added<2)
  { int a;
    for (a = 0; added+a<2; a++)
    { do
      { i = rand_int(M, seed);
        j = rand_int(N, seed);
      } while (mod2sparse_find(H,i,j));
      mod2sparse_insert(H,i,j);
    }
    fprintf(stderr,
 "Added %d extra bit-checks to try to avoid problems from even column counts\n",
      a);
  }

  /* Eliminate cycles of length four, if asked, and if possible. */
  if (no4cycle)
  { 
    elim4 = 0;

    for (t = 0; t<10; t++) 
    { k = 0;
      for (j = 0; j<N; j++)
      { for (e = mod2sparse_first_in_col(H,j);
             !mod2sparse_at_end(e);
             e = mod2sparse_next_in_col(e))
        { for (f = mod2sparse_first_in_row(H,mod2sparse_row(e));
               !mod2sparse_at_end(f);
               f = mod2sparse_next_in_row(f))
          { if (f==e) continue;
            for (g = mod2sparse_first_in_col(H,mod2sparse_col(f));
                 !mod2sparse_at_end(g);
                 g = mod2sparse_next_in_col(g))
            { if (g==f) continue;
              for (h = mod2sparse_first_in_row(H,mod2sparse_row(g));
                   !mod2sparse_at_end(h);
                   h = mod2sparse_next_in_row(h))
              { if (mod2sparse_col(h)==j)
                { do
                  { i = rand_int(M,seed);
                  } while (mod2sparse_find(H,i,j));
                  mod2sparse_delete(H,e);
                  mod2sparse_insert(H,i,j);
                  elim4 += 1;
                  k += 1;
                  goto nextj;
                }
              }
            }
          }
        }
      nextj: ;
      }
      if (k==0) break;
    }

    if (elim4>0)
    { fprintf(stderr,
        "Eliminated %d cycles of length four by moving checks within column\n",
         elim4);
    }

    if (t==20) 
    { fprintf(stderr,
        "Couldn't eliminate all cycles of length four in 10 passes\n");
    }
  }
   file = open_file_std(pchk_file,"wb");
  if (file==NULL) 
  { fprintf(stderr,"Can't create parity check file: %s\n",file);
    exit(1);
  }

  intio_write(file,('P'<<8)+0x80);
  
  if (ferror(file) || !mod2sparse_write(file,H) || fclose(file)!=0)
  { fprintf(stderr,"Error writing to parity check file %s\n",file);
    exit(1);
  }

  fclose(file);
  mod2sparse_free(H);

 
}


int move_ps( int i, int j, int seed)
{
	mod2entry *t;
	int k, m, temp = 0, temp1 = 0;
	
n_r:if ((k = rand_int(N, seed))==j) 
		goto n_r;
	if ( mod2sparse_find(H, i, k))	
		// find i_th row and k_th col is empty 
	{
		if ((temp++) > (N*N))
			return (0);
		goto n_r;
	}

	if ((temp1++) > (N*N))
		return (0);
	temp = 0;

	for (m=0; m<M; m++)
	{
		if ( ((t = mod2sparse_find(H, m, k))!=0) &&
		// find m_th row and k_th col is not empty 

		( mod2sparse_find(H, m, j)==0))
		// find m_th row and j_th col is empty 
		goto m_r;
	}
	goto n_r;
m_r:mod2sparse_insert(H,i,k);
	mod2sparse_delete(H,t);
    mod2sparse_insert(H,m,j);

//	printf("e row=%d col=%d\n", i, j);
	return (1);
}


//Eliminate cycles of length four for lin jiaru, and if possible.

int eliminate( mod2entry *e, int j, int seed)
{
	mod2entry *t;
	int i, k, m, temp = 0, temp1 = 0, temp2 = 0;
	
	i = mod2sparse_row(e);
n_r:if ((k = rand_int(N, seed))==j) 
		goto n_r;
	if ( mod2sparse_find(H, i, k))	
		// find i_th row and k_th col is empty 
	{
		if ((temp++) > (N*N))
			return (0);
		goto n_r;
	}
	if ((temp1++) > (N*N))
		return (0);
	temp = 0;
	temp2 = 0;
	
m_r:m = rand_int(M, seed);
	if ( (t = mod2sparse_find(H, m, k))==0)
		// find m_th row and k_th col is not empty 
		goto m_r;
	if ( mod2sparse_find(H, m, j)!=0)
		// find m_th row and j_th col is empty 
	{
		if ((temp2++) > (N*N))
			goto n_r;
		goto m_r;
	}
/*
	for (m=0; m<M; m++)
	{
		if ( ((t = mod2sparse_find(H, m, k))!=0) &&
		// find m_th row and k_th col is not empty 

		( mod2sparse_find(H, m, j)==0))
		// find m_th row and j_th col is empty 
		goto m_r;
	}
	goto n_r;
	*/
	mod2sparse_delete(H,e);
    mod2sparse_insert(H,i,k);
	mod2sparse_delete(H,t);
    mod2sparse_insert(H,m,j);
	return (1);

//	printf(" row=%d col=%d,", i, j);
}



/* CREATE A SPARSE PARITY-CHECK MATRIX.  Of size M by N, stored in H. */

void make_ldpc_zhj1(char *pchk_file, int v,int num_sub_matrix)
{
  FILE *file;
  int i, j, k;
  int row, col;
  int remainder;
  int **L;
  
  if ((L = (int **) malloc(num_sub_matrix*sizeof(int *))) == NULL)
  {
	printf("Not enough memory to allocate buffer.\n");
	exit(1);
  }
	for (i=0;i<num_sub_matrix;i++)
	{
		if ((L[i]= (int *) malloc(3*sizeof(int))) == NULL)
		{
			printf("Not enough memory to allocate buffer.\n");
			exit(1);
		}
	}

  if (num_sub_matrix>(int)(((v/6)-1)*0.5))
  {
	  fprintf(stderr,"can't construct so many matrices\n");
	  exit(1);
  }

  remainder=(int)(fmod(v,6));
  switch (remainder)
  {
  case 1:
		  for (i=0;i<num_sub_matrix;i++)
		  {
			  L[i][0]=2*i+1;
			  L[i][1]=5*(v/6)-4*i-1;
			  L[i][2]=(v/6)+2*i+1;
		  }
		  break;
  case 3:
		  for (i=0;i<num_sub_matrix;i++)
		  {
			  L[i][0]=2*i+1;
		      L[i][1]=5*(v/6)-4*i+1;
			  L[i][2]=(v/6)+2*i+1;
		  }
		  break;
  case 5:
		  for (i=0;i<num_sub_matrix;i++)
		  {
			  L[i][0]=2*i+1;
			  L[i][1]=5*(v/6)-4*i+2;
			  L[i][2]=(v/6)+2*i+2;
		  }
		  break;
  default: printf("can't construct for such 'v'\n");
  }

  for (i=0;i<num_sub_matrix;i++)
  {
	  L[i][2] = 0;
	  L[i][0] = (int) (fmod(L[i][2]+L[i][0],v));
	  L[i][1] = (int) (fmod(L[i][0]+L[i][1],v));
  }

  H = mod2sparse_allocate(M,N);