dec.c

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

  

}

void bubble(double a[],int size,int d[] )
{

	double temp;
	int i,j,k;


	for(i=1;i<size;i++)
	{
		for(j=0;j<size-i;j++)
			if(a[j]<a[j+1])
			{
				temp=a[j];
				a[j]=a[j+1];
				a[j+1]=temp;

				k=d[j];
				d[j]=d[j+1];
				d[j+1]=k;

			}
	}
}

void ibubble(int a[],int size,int d[] )
{

	int temp;
	int i,j,k;


	for(i=1;i<size;i++)
	{
		for(j=0;j<size-i;j++)
			if(a[j]<a[j+1])
			{
				temp=a[j];
				a[j]=a[j+1];
				a[j+1]=temp;

				k=d[j];
				d[j]=d[j+1];
				d[j+1]=k;

			}
	}
}


/* DO ONE ITERATION OF PROBABILITY PROPAGATION. */

void iterprp_min_sum
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  int *dblk,		/* Place to store decoding */
  double *LLR
)
{
  double pr, dl;
  mod2entry *e;
  int N, M;
  int i, j,temp;

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

   /* Recompute likelihood ratios. */

  for (i = 0; i<M; i++)
  { e = mod2sparse_first_in_row(H,i);
	dl = e->pr;
	e->lr=0x7fffffffff;
    for (e = mod2sparse_next_in_row(e);                                   /* Find the second entry in a row */
         !mod2sparse_at_end(e);             /* See if we've reached the end     */
         e = mod2sparse_next_in_row(e))    /*Move from one entry to the right direction*/
    { e->lr = dl;
	  temp = (dl*e->pr)>0 ? 1 :-1;
	  dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	  dl*=temp;
	  
    }

    dl = 0x7fffffffff;
    for (e = mod2sparse_last_in_row(H,i); /* Find the last entry in a row */
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))   /*Move from one entry to the left direction*/
    { 
	  temp = (dl*e->lr)>0 ? 1 :-1;
	  e->lr= (fabs(dl)<fabs(e->lr)) ? fabs(dl):fabs(e->lr);
	  e->lr *=temp;

	  temp = (dl*e->pr)>0 ? 1 :-1;
	  dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	  dl *=temp;
    }
  }

 /* 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); /* Find the first entry in a column */
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_col(e))   /*Move from one entry to the down direction*/
    { e->pr = pr;
      pr += e->lr;
    }

    LLR[j]=pr;//lxh
	dblk[j] =(pr>=0);
    pr = 0;
    for (e = mod2sparse_last_in_col(H,j); /* Find the last entry in a column */
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_col(e))  /* Find the up entry in a column */
    { e->pr += pr;
/*      if (isnan(e->pr)) 
      { e->pr = 1;
      }
*/      pr += e->lr;
    }
  }

}


void iterprp_min_sum_constant
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  int *dblk,		/* Place to store decoding */
  double *LLR,
  double cons
)
{
  double pr, dl,modify;
  mod2entry *e;
  int N, M;
  int i, j,temp;

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

   /* Recompute likelihood ratios. */

  for (i = 0; i<M; i++)
  { e = mod2sparse_first_in_row(H,i);
	dl = e->pr;
	e->lr=0x7fffffffff;
    for (e = mod2sparse_next_in_row(e);                                   /* Find the second entry in a row */
         !mod2sparse_at_end(e);             /* See if we've reached the end     */
         e = mod2sparse_next_in_row(e))    /*Move from one entry to the right direction*/
    { e->lr = dl;
	 temp = (dl*e->pr)>0 ? 1 :-1;
	 if ((fabs(dl-e->pr))>2*(fabs(dl+e->pr)) && (fabs(dl+e->pr)<2))
		 modify=cons;
	 else if ((fabs(dl+e->pr))>2*(fabs(dl-e->pr)) && (fabs(dl-e->pr)<2))
		 modify=-cons;
	 else 
		 modify=0;
	 /*验算是否和BP算法一致*/
	// modify = log(1+exp(-(fabs(dl+e->pr))))-log(1+exp(-(fabs(dl-e->pr))));

	  dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	  dl*=temp;
	  dl +=modify;
	  
    }

    dl = 0x7fffffffff;
    for (e = mod2sparse_last_in_row(H,i); /* Find the last entry in a row */
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))   /*Move from one entry to the left direction*/
    { 
	 temp = (dl*e->lr)>0 ? 1 :-1;
	 if ((fabs(dl-e->lr))>2*(fabs(dl+e->lr)) && (fabs(dl+e->lr)<2))
		 modify=cons;
	 else if ((fabs(dl+e->lr))>2*(fabs(dl-e->lr)) && (fabs(dl-e->pr)<2))
		 modify=-cons;
	 else 
		 modify=0;

		//	 modify = log(1+exp(-(fabs(dl+e->lr))))-log(1+exp(-(fabs(dl-e->lr))));

	  e->lr= (fabs(dl)<fabs(e->lr)) ? fabs(dl):fabs(e->lr);
	  e->lr *=temp;
  	  e->lr +=modify;


	  temp = (dl*e->pr)>0 ? 1 :-1;
	   if ((fabs(dl-e->pr))>2*(fabs(dl+e->pr)) && (fabs(dl+e->pr)<2))
		 modify=cons;
	 else if ((fabs(dl+e->pr))>2*(fabs(dl-e->pr)) && (fabs(dl-e->pr)<2))
		 modify=-cons;
	 else 
		 modify=0;
	// modify = log(1+exp(-(fabs(dl+e->pr))))-log(1+exp(-(fabs(dl-e->pr))));


	  dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	  dl *=temp;
  	  dl +=modify;

    }
  }

 /* 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); /* Find the first entry in a column */
         !mod2sparse_at_end(e);
         e = mod2sparse_next_in_col(e))   /*Move from one entry to the down direction*/
    { e->pr = pr;
      pr += e->lr;
    }

    LLR[j]=pr;//lxh
	dblk[j] =(pr>=0);
    pr = 0;
    for (e = mod2sparse_last_in_col(H,j); /* Find the last entry in a column */
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_col(e))  /* Find the up entry in a column */
    { e->pr += pr;
/*      if (isnan(e->pr)) 
      { e->pr = 1;
      }
*/      pr += e->lr;
    }
  }

}

void iterprp_min_sum_table
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  int *dblk,		/* Place to store decoding */
  double *LLR
)
{
  double pr, dl,modify;
  double table[8]={0.65,0.55,0.45,0.35,0.25,0.15,0.05,0.0};
  mod2entry *e;
  int N, M;
  int i, j,temp;

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

   /* Recompute likelihood ratios. */

  for (i = 0; i<M; i++)
  { e = mod2sparse_first_in_row(H,i);
	dl = e->pr;
	e->lr=0x7fffffffff;
    for (e = mod2sparse_next_in_row(e);                                   /* Find the second entry in a row */
         !mod2sparse_at_end(e);             /* See if we've reached the end     */
         e = mod2sparse_next_in_row(e))    /*Move from one entry to the right direction*/
    { e->lr = dl;
	 temp = (dl*e->pr)>0 ? 1 :-1;

	if ((fabs(dl+e->pr)<0.196) && (fabs(dl+e->pr)>=0))
		 modify=table[0];
	 else if ((fabs(dl+e->pr)<0.433) && (fabs(dl+e->pr)>=0.196))
		 modify=table[1];
	 else if ((fabs(dl+e->pr)<0.71) && (fabs(dl+e->pr)>=0.433))
		 modify=table[2];
	 else if ((fabs(dl+e->pr)<1.05) && (fabs(dl+e->pr)>=0.71))
		 modify=table[3];
	 else if ((fabs(dl+e->pr)<1.508) && (fabs(dl+e->pr)>=1.05))
		 modify=table[4];
	 else if ((fabs(dl+e->pr)<2.252) && (fabs(dl+e->pr)>=1.508))
		 modify=table[5];
	 else if ((fabs(dl+e->pr)<4.5) && (fabs(dl+e->pr)>=2.252))
		 modify=table[6];
	 else 
		 modify=table[7];

	 if ((fabs(dl-e->pr)<0.196) && (fabs(dl-e->pr)>=0))
		 modify-=table[0];
	 else if ((fabs(dl-e->pr)<0.433) && (fabs(dl-e->pr)>=0.196))
		 modify-=table[1];
	 else if ((fabs(dl-e->pr)<0.71) && (fabs(dl-e->pr)>=0.433))
		 modify-=table[2];
	 else if ((fabs(dl-e->pr)<1.05) && (fabs(dl-e->pr)>=0.71))
		 modify-=table[3];
	 else if ((fabs(dl-e->pr)<1.508) && (fabs(dl-e->pr)>=1.05))
		 modify-=table[4];
	 else if ((fabs(dl-e->pr)<2.252) && (fabs(dl-e->pr)>=1.508))
		 modify-=table[5];
	 else if ((fabs(dl-e->pr)<4.5) && (fabs(dl-e->pr)>=2.252))
		 modify-=table[6];
	 else 
		 modify-=table[7];


	 /*演算是否和BP算法一致*/
	// modify = log(1+exp(-(fabs(dl+e->pr))))-log(1+exp(-(fabs(dl-e->pr))));

	  dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	  dl*=temp;
	  dl +=modify;
	  
    }

    dl = 0x7fffffffff;
    for (e = mod2sparse_last_in_row(H,i); /* Find the last entry in a row */
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))   /*Move from one entry to the left direction*/
    { temp = (dl*e->lr)>0 ? 1 :-1;
	 if ((fabs(dl+e->lr)<0.196) && (fabs(dl+e->lr)>=0))
		 modify=table[0];
	 else if ((fabs(dl+e->lr)<0.433) && (fabs(dl+e->lr)>=0.196))
		 modify=table[1];
	 else if ((fabs(dl+e->lr)<0.71) && (fabs(dl+e->lr)>=0.433))
		 modify=table[2];
	 else if ((fabs(dl+e->lr)<1.05) && (fabs(dl+e->lr)>=0.71))
		 modify=table[3];
	 else if ((fabs(dl+e->lr)<1.508) && (fabs(dl+e->lr)>=1.05))
		 modify=table[4];
	 else if ((fabs(dl+e->lr)<2.252) && (fabs(dl+e->lr)>=1.508))
		 modify=table[5];
	 else if ((fabs(dl+e->lr)<4.5) && (fabs(dl+e->lr)>=2.252))
		 modify=table[6];
	 else 
		 modify=table[7];

	 if ((fabs(dl-e->lr)<0.196) && (fabs(dl-e->lr)>=0))
		 modify-=table[0];
	 else if ((fabs(dl-e->lr)<0.433) && (fabs(dl-e->lr)>=0.196))
		 modify-=table[1];
	 else if ((fabs(dl-e->lr)<0.71) && (fabs(dl-e->lr)>=0.433))
		 modify-=table[2];
	 else if ((fabs(dl-e->lr)<1.05) && (fabs(dl-e->lr)>=0.71))
		 modify-=table[3];
	 else if ((fabs(dl-e->lr)<1.508) && (fabs(dl-e->lr)>=1.05))
		 modify-=table[4];
	 else if ((fabs(dl-e->lr)<2.252) && (fabs(dl-e->lr)>=1.508))
		 modify-=table[5];
	 else if ((fabs(dl-e->lr)<4.5) && (fabs(dl-e->lr)>=2.252))
		 modify-=table[6];
	 else 
		 modify-=table[7];


		//	 modify = log(1+exp(-(fabs(dl+e->lr))))-log(1+exp(-(fabs(dl-e->lr))));

	  e->lr= (fabs(dl)<fabs(e->lr)) ? fabs(dl):fabs(e->lr);
	  e->lr *=temp;
  	  e->lr +=modify;


	  temp = (dl*e->pr)>0 ? 1 :-1;

	 if ((fabs(dl+e->pr)<0.196) && (fabs(dl+e->pr)>=0))
		 modify=table[0];
	 else if ((fabs(dl+e->pr)<0.433) && (fabs(dl+e->pr)>=0.196))
		 modify=table[1];
	 else if ((fabs(dl+e->pr)<0.71) && (fabs(dl+e->pr)>=0.433))
		 modify=table[2];
	 else if ((fabs(dl+e->pr)<1.05) && (fabs(dl+e->pr)>=0.71))
		 modify=table[3];
	 else if ((fabs(dl+e->pr)<1.508) && (fabs(dl+e->pr)>=1.05))
		 modify=table[4];
	 else if ((fabs(dl+e->pr)<2.252) && (fabs(dl+e->pr)>=1.508))
		 modify=table[5];
	 else if ((fabs(dl+e->pr)<4.5) && (fabs(dl+e->pr)>=2.252))
		 modify=table[6];
	 else 
		 modify=table[7];

	 if ((fabs(dl-e->pr)<0.196) && (fabs(dl-e->pr)>=0))
		 modify-=table[0];
	 else if ((fabs(dl-e->pr)<0.433) && (fabs(dl-e->pr)>=0.196))
		 modify-=table[1];
	 else if ((fabs(dl-e->pr)<0.71) && (fabs(dl-e->pr)>=0.433))
		 modify-=table[2];
	 else if ((fabs(dl-e->pr)<1.05) && (fabs(dl-e->pr)>=0.71))
		 modify-=table[3];
	 else if ((fabs(dl-e->pr)<1.508) && (fabs(dl-e->pr)>=1.05))
		 modify-=table[4];
	 else if ((fabs(dl-e->pr)<2.252) && (fabs(dl-e->pr)>=1.508))
		 modify-=table[5];
	 else if ((fabs(dl-e->pr)<4.5) && (fabs(dl-e->pr)>=2.252))
		 modify-=table[6];
	 else 
		 modify-=table[7];