dec.c

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

	// 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;
    }
  }

}


/* DO ONE ITERATION OF PROBABILITY PROPAGATION. */
/* Efficient implementations of the sum-product algorithm for decoding LDPC codes.pdf */

void iterprp_min_sum_parallel
( 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;
  double modify;


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

   /* Recompute likelihood ratios. */

  for (i = 0; i<M; i++)
  { 
	dl=0x7fffffffff;

    for (e = mod2sparse_first_in_row(H,i);  // 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
    { 
		 temp = (dl*e->pr)>0 ? 1 :-1;
	  
	if ((fabs(dl+e->pr)<0.5) && (fabs(dl+e->pr)>=0))
		 modify=-0.5*fabs(dl+e->pr)+0.7;
	else if ((fabs(dl+e->pr)<1.6) && (fabs(dl+e->pr)>=0.5))
		 modify=-0.25*fabs(dl+e->pr)+0.575;
	else if ((fabs(dl+e->pr)<2.2) && (fabs(dl+e->pr)>=1.6))
		 modify=-0.125*fabs(dl+e->pr)+0.375;
	else if ((fabs(dl+e->pr)<3.2) && (fabs(dl+e->pr)>=2.2))
		 modify=-0.0625*fabs(dl+e->pr)+0.2375;
	else if ((fabs(dl+e->pr)<4.4) && (fabs(dl+e->pr)>=3.2))
		 modify=-0.03125*fabs(dl+e->pr)+0.1375;
	else 
		 modify=0;

	 if ((fabs(dl-e->pr)<0.5) && (fabs(dl-e->pr)>=0))
		 modify-=-0.5*fabs(dl-e->pr)+0.7;
	 else if ((fabs(dl-e->pr)<1.6) && (fabs(dl-e->pr)>=0.5))
		 modify-=-0.25*fabs(dl-e->pr)+0.575;
	 else if ((fabs(dl-e->pr)<2.2) && (fabs(dl-e->pr)>=1.6))
		 modify-=-0.125*fabs(dl-e->pr)+0.375;
	 else if ((fabs(dl-e->pr)<3.2) && (fabs(dl-e->pr)>=2.2))
		 modify-=-0.0625*fabs(dl-e->pr)+0.2375;
	 else if ((fabs(dl-e->pr)<4.4) && (fabs(dl-e->pr)>=3.2))
		 modify-=-0.03125*fabs(dl-e->pr)+0.1375;
	 else 
		 modify-=0;

	     dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	     dl*=temp;
	     dl+=modify;
			  
    }
/*	e = mod2sparse_first_in_row(H,i);
	dl = e->pr;
    for (e = mod2sparse_next_in_row(e);     
         !mod2sparse_at_end(e);             
         e = mod2sparse_next_in_row(e))    
    { dl = log(1+exp(dl+e->pr))-log(exp(dl)+exp(e->pr));
	}*/
   

    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*/
    { 
		if ((dl+e->pr)<-3)
			modify=0;
		else if (((dl+e->pr)<-0.68) && ((dl+e->pr)>=-3))
			modify=-0.25*(dl+e->pr)-0.75;
		else if (((dl+e->pr)<-0.27) && ((dl+e->pr)>=-0.68))
			modify=-2*(dl+e->pr)-1.94;
		else if (((dl+e->pr)<0.0) && ((dl+e->pr)>=-0.27))
			modify=-0.125*(dl+e->pr)-3.56;
		else if (((dl+e->pr)<0.15) && ((dl+e->pr)>=0))
			modify=16*(dl+e->pr)-4;
		else if (((dl+e->pr)<0.4) && ((dl+e->pr)>=0.15))
			modify=4*(dl+e->pr)-2.2;
		else if (((dl+e->pr)<1.3) && ((dl+e->pr)>=0.4))
			modify=2*(dl+e->pr)-1.4;
		else
			modify=(dl+e->pr)-0.1;

		if ((e->pr-dl)<-3)
			modify-=0;
		else if (((e->pr-dl)<-0.68) && ((e->pr-dl)>=-3))
			modify-=-0.25*(e->pr-dl)-0.75;
		else if (((e->pr-dl)<-0.27) && ((e->pr-dl)>=-0.68))
			modify-=-2*(e->pr-dl)-1.94;
		else if (((e->pr-dl)<0.0) && ((e->pr-dl)>=-0.27))
			modify-=-0.125*(e->pr-dl)-3.56;
		else if (((e->pr-dl)<0.15) && ((e->pr-dl)>=0))
			modify-=16*(e->pr-dl)-4;
		else if (((e->pr-dl)<0.4) && ((e->pr-dl)>=0.15))
			modify-=4*(e->pr-dl)-2.2;
		else if (((e->pr-dl)<1.3) && ((e->pr-dl)>=0.4))
			modify-=2*(e->pr-dl)-1.4;
		else
			modify-=(e->pr-dl)-0.1;
		
		e->lr=modify-dl;
	//	e->lr=log(fabs(exp((e->pr)+dl)-1))-log(fabs(exp((e->pr)-dl)-1))-dl;
		
    }
  }

 /* 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_lin
( mod2sparse *H,	/* Parity check matrix */
  double *lratio,	/* Likelihood ratios for bits */
  int *dblk,		/* Place to store decoding */
  double *LLR
)
{
  double pr, dl,dmax,dm;
  mod2entry *e;
  int N, M;
  int i, j, temp,sgn;
  double cons=0.8;


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

   /* Recompute likelihood ratios. */

 for (i = 0; i<M; i++)
  { 
	dl=0x7fffffffff;

    for (e = mod2sparse_first_in_row(H,i);  // 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
    { 
		 temp = (dl*e->pr)>0 ? 1 :-1;
	  
	 	 dl = (fabs(dl)<fabs(e->pr)) ? fabs(dl):fabs(e->pr);
	     dl*=temp;
		 dmax=(dmax>fabs(e->pr)) ? dmax:fabs(e->pr);
			  
    }
	
	dm=fabs(dl);
	for (e = mod2sparse_first_in_row(H,i);  // 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
    { 
	  if (dmax>fabs(e->pr) && (fabs(e->pr)>dm))
		  dmax=fabs(e->pr);
	  
    }


    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*/
    { 
		sgn=(e->pr)>0 ? 1 :-1;
		temp=(dl>0)?1:-1;
		if (fabs(dl)<fabs(e->pr))
		   e->lr=sgn*dl*0.77;
		else
			e->lr=sgn*temp*dmax*0.77;
		
		
    }
  }

 /********************/
 /*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);                                   
         !mod2sparse_at_end(e);            
         e = mod2sparse_next_in_row(e))    
    { 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); 
         !mod2sparse_at_end(e);
         e = mod2sparse_prev_in_row(e))   
    { 
	  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;
    }
  }

}

/* DO ONE ITERATION OF PROBABILITY PROPAGATION. */

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

  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 = fabs(e->pr);
	e->lr=0x7fffffffff;

	temp=mod2sparse_col(e);
	thea=dblk[temp];
    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;
	  dl = (dl<fabs(e->pr)) ? dl:fabs(e->pr);
      temp=mod2sparse_col(e);
	  thea+=dblk[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*/
    { 
	  e->lr= (dl<e->lr) ? dl:e->lr;
	 
	  temp=mod2sparse_col(e);
	  e->lr*=(-cons)*pow(-1,(int)(fmod(thea-dblk[temp],2)));

      dl = (dl<fabs(e->pr)) ? dl:fabs(e->pr);
    }
  }

 /* 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;
    }
  }

}