www.pudn.com > 7894509matlab.rar > ldpc_generate.c, change:2000-05-09,size:8815b


/* Generate H for  LDPC*/  
/* Based on sparse.c by Matt Davey <mcdavey@mrao.cam.ac.uk> with his permission*/  
#include <math.h>  
#include "mex.h"  
  
/* Input Arguments: parameters*/  
#define	M_IN	prhs[0]          /* number of parity checks */  
#define	N_IN	prhs[1]          /* blocklength */   
#define	T_IN	prhs[2]          /* mean column weight */  
#define	Q_IN	prhs[3]          /* GF base  */  
#define	SEED_IN	prhs[4]          /* seed for random generator */  
  
  
/* Output Arguments: matrices*/  
#define	H_OUT	plhs[0]  
  
  
void mexFunction(  
                 int nlhs,       mxArray *plhs[],  
                 int nrhs, const mxArray *prhs[]  
		 )  
{  
  short	**M_list, **N_list, *M_target;  
  double *pp,*sr,*s,*ss;  
  int N,M,q,i,j,k,nzmax,*irs,*jcs,l,  
	  tr,tm,tc,done,redo,tmp,regime,tr_max,t_max, m_low;  
  float t;  
  long seed;  
  void adjust(int *, int *, int *, int *);  
  unsigned int K2,M2;  
  char c;  
  mxArray *arg_in[2], *arg_out[1]; /* to call rand generator of Matlab*/  
  
  /* Check for proper number of arguments */  
  if (nrhs != 5) {  
    mexErrMsgTxt("GENERATE requires five input arguments.");  
  } else if (nlhs > 1) {  
    mexErrMsgTxt("GENERATE requires one output argument.");  
  }  
  
  pp = mxGetPr(N_IN);    N = (int) (*pp);  
  pp = mxGetPr(M_IN);    M = (int) (*pp);  
  pp = mxGetPr(T_IN);    t = (float) (*pp);  
  pp = mxGetPr(Q_IN);    q = (int) (*pp);  
  pp = mxGetPr(SEED_IN); seed = (int) (*pp);  
    
  
  arg_in[0] = mxCreateString("state");  
  arg_in[1] = mxCreateDoubleMatrix(1, 1, mxREAL);  
  s = mxGetPr(arg_in[1]);  
  s[0] = seed; /* this will be used to call rand*/  
  
  /* initialize random generator */  
  mexCallMATLAB(0, NULL, 2, arg_in, "rand"); /* rand('state',seed) */  
  s[0] = 1; /* from now on use s to store "1"*/  
  
  
/* I have no idea about the details of the following - please  
ask the author - :).  igor*/  
  
  
  
  
  
/* Generate some sparse matrices for error-correcting codes.  Supply  
 * the following parameters:  
 * 	N: Blocklength  
 * 	M: Number of parity checks  
 * 	t: Mean column weight.  
 *  
 * If t<3, then we generate weight 2 columns systematically in the  
 * form of blocks of identity matrices to reduce probability of  
 * getting short cycle lengths.  Can't have more than M weight 2  
 * columns, though.  
 *  
 * Having generated any weight 2 columns, we fill the rest.  Calculate  
 * number of columns to fill with floor(t) and number with ceiling(t).  
 * Find mean row weight and calculate number of rows to fill with  
 * floor(r) and ceiling(r).  
 *  
 * We fill rows as follows:  
 * (Regime 0): Fill so that rows contain <= tr ones until M-tm rows  
 * 	    contain tr ones.  
 * (Regime 1): Fill so that rows contain <= (tr+1) ones until tm rows  
 * 	    contain (tr+1) ones.  
 * (Regime 2): Fill remaining rows containing < tr ones until done.  
 */  
  
   
  /* N=6112;M=4512;t=2.3;*/  
  
       
  t_max=(int)ceil(t);  
  M_list=(short **)mxMalloc(N*sizeof(short *));  
  M_target=(short *)mxMalloc(M*sizeof(short *));  
  N_list=(short **)mxMalloc(M*sizeof(short *));  
  for(i=0;i<N;i++){  
    M_list[i]=(short *)mxMalloc((t_max+1)*sizeof(short));  
  }  
  i=0;  
  /* Do we have any weight 2 columns?  
   *  
   * If so, do this first.  Remember that M mightn't have a large  
   * power of two as a divisor, so might need to find some M'<M to use  
   * as unit length.  
   */  
  K2=0;  
  if(t<3){  
    K2=ceil((double)N*(3-t));  
	if(K2>M)  
		mexErrMsgTxt("GENERATE: Can't have more than M weight 2 columns.");  
    j=2;  
    done=0;  
    for(i=0;!done;i++){  
      M2=floor((double)M/(double)j);  
      if((M2*(j-1))>=K2) done=1;  
      j*=2;  
    }  
    M2*=(j/4);  
  }  
  /*   
   * i contains number of identity blocks we'll need.... */  
  tr=((short)floor((double)(t*N)/(double)M));  
  /* Now we want `tr' to be final minimum row weight, `tr_max' to be  
   * final maximum row weight, `tm' to be number of rows which will  
   * have weight greater than `tr'.  'tc' will be a running count of  
   * how many rows we still have to fill up to weight `tr'.  Once we  
   * hit this many, we can start overfilling rows.  
   */  
  if (i>tr){  
    tr_max=i;  
    /* If identity blocks make overheavy rows, we need to calculate  
     * the minimum row weight.  
     */  
    done=0;  
    k=1;  
    j=floor((double)t*N)-2*K2; /* Number of ones left to distribute */  
    for(i=0;!done;i++){  
      /* (M-2*M2) rows will be empty after identity blocks */  
      j-=((M-2*M2)+(2*M2*(k-1))/k);  
      if(j<0) {  
	done=1;  
      }  
      else {  
	k*=2;  
      }  
    }  
    tr=i-1;  
    tm=M+j;  
  }  
  else {  
    /* This is easier! */  
    tr_max=tr+1;  
    tm=(int)floor((((double)t*N)/(double)M-tr)*M +0.5);  
  }  
  tc=M-tm;  
  for(i=0;i<M;i++){  
    N_list[i]=(short *)mxMalloc((tr_max+1)*sizeof(short));    
  }  
  for(i=0;i<M;i++){  
    N_list[i][0]=0;  
  }  
  
  regime=0;  
  /* Generate weight 2 columns.  First create two identity matrices on  
   * top of each other, then two 1/2 size ones in the lower rows of  
   * the matrix, and so on.  
   *  
   * j:   length of current identity block  
   * k:   base of this block  
   * i:   current column position  
   */  
  j=M2;  
  k=0;  
  i=0;  
  while(i<K2){  
    for(;(i-k)<j && i<K2;i++){  
      M_list[i][0]=2;  
      M_list[i][1]=i;  
      M_list[i][2]=i+j;  
      N_list[i][0]++;  
      if(N_list[i][0]==tr) adjust(&tm,&tr,&tc,®ime);  
      N_list[i][N_list[i][0]]=i;  
      N_list[i+j][0]++;  
      if(N_list[i+j][0]==tr) adjust(&tm,&tr,&tc,®ime);  
      N_list[i+j][N_list[i+j][0]]=i;  
    }  
    k=i;  
    j/=2;  
  }  
   
 /* Now fill the unsystematic columns, ensuring weight per row as even as poss. */  
  i=K2;  
  if(K2==0){  
    /* Fill low weight columns */  
    for(i=0;i<(int)(N*(t_max-t)+0.5);i++){  
      for(k=1;k<=(int)floor(t);k++){  
	done=0;  
	do {  
      mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand"); /* ss = rand(1) */  
	  ss = mxGetPr(arg_out[0]);  
	  j=(short)floor(M*ss[0]);  
      mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand");  
	  ss = mxGetPr(arg_out[0]);  
	  if((ss[0])<(1-(double)N_list[j][0]/(double)tr)) {  
	    done=1;  
	    for(l=1;l<k;l++) if(j==M_list[i][l]) done=0;  
	  }  
	} while(!done);  
	N_list[j][0]++;  
	N_list[j][N_list[j][0]]=i;  
	if(N_list[j][0]==tr) adjust(&tm,&tr,&tc,®ime);  
	M_list[i][k]=j;  
      }  
      M_list[i][0]=k-1;  
    }  
  }  
  redo=1;  
  for(;i<N;i++){  
    fprintf(stderr,"%d\r",i);  
    for(k=1;k<=t_max;k++){  
      done=0;  
      do {  
	/* find the lowest weight rows, and fill one of them */  
	if(redo){  
	  l=tr_max;  
	  for(j=0;j<M;j++) if (N_list[j][0]<l) l=N_list[j][0];  
	  m_low=0;  
	  for(j=0;j<M;j++) if (N_list[j][0]==l) {M_target[m_low]=j; m_low++;}  
	}  
	mexCallMATLAB(1, arg_out,1 , &arg_in[1], "rand");  
	ss = mxGetPr(arg_out[0]);  
	j=M_target[tmp=(short)floor(m_low*ss[0])];  
	/*	if(ss[0]<(1-(double)N_list[j][0]/(double)tr)) {*/  
	  done=1;  
	  for(l=1;l<k;l++) if(j==M_list[i][l]) done=0;  
	  if(done==1){  
	    if(m_low==1) redo=1;  
	    else {  
	      for(;tmp<(m_low-1);tmp++) M_target[tmp]=M_target[tmp+1];  
	      m_low--;  
	      redo=0;  
	    }  
	  }  
	  /*	}*/  
	    
      } while(!done);  
      N_list[j][0]++;  
      N_list[j][N_list[j][0]]=i;  
      if(N_list[j][0]==tr) adjust(&tm,&tr,&tc,®ime);  
      M_list[i][k]=j;  
    }  
    M_list[i][0]=k-1;  
  }  
  tr=((short)ceil((double)(3*N-K2)/(double)M));  
    
  
  for(i=0;i<M;i++){  
    mxFree(N_list[i]);  
  }  
  mxFree(N_list);  
  
  
/* done generating H matrix */  
  
  
  
  
    
    /* Allocate space for sparse matrix */  
    nzmax=0; for(j=0 ; j<N ; j++) nzmax +=M_list[j][0];  
	mexPrintf("%d \n",nzmax);  
    /* NOTE: The maximum number of non-zero elements cannot be less  
       than the number of columns in the matrix. */  
    if (N>nzmax){  
	nzmax=N;  
    }  
    plhs[0] = mxCreateSparse(M,N,nzmax,mxREAL);  
    sr  = mxGetPr(plhs[0]);  
    irs = mxGetIr(plhs[0]);  /* row */  
    jcs = mxGetJc(plhs[0]);  /* column */  
      
    /* Copy nonzeros */  
    k = 0;   
    for (j=0; (j<N ); j++) {  
	jcs[j] = k;  
	for (i=1; (i<=M_list[j][0] ); i++) {  
		sr[k] = 1;  
		irs[k] = M_list[j][i];  
		k++;	      
	}  
    }  
    jcs[N] = k;  
  
  
  
  for(i=0;i<N;i++){  
    mxFree(M_list[i]);  
  }  
  mxFree(M_list);  
  
  return;  
}  
  
  
  
void adjust(int *tm, int *tr, int *tc, int *regime){  
  switch(*regime){  
  case 0:  
    (*tc)--;  
    if((*tc)==0){  
      *regime=1;  
      (*tr)++;  
    }  
    break;  
  case 1:  
    (*tm)--;  
    if((*tm)==0){  
      *regime=2;  
      (*tr)--;  
    }  
    break;  
  }  
}