berg_algs.c

有关信道估计和信道均衡的仿真程序

   if( is_cmplx ){
     mexCallMATLAB(num_out,plhs2,num_in,prhs,"rand"); 
     di = mxGetPr(plhs2[0]);
   }
  /* normalize dither to [-alpha,alpha) */
   for (k=0; k<N; k++){ dr[k] -= 0.5; dr[k] *= 2*alpha; }
   if( is_cmplx ){
     for (k=0; k<N; k++){ di[k] -= 0.5; di[k] *= 2*alpha; }
   }
  
  /* normalize stepsize according to alpha */
  mu = alpha*mu_without_alpha;

  /* allocate local memory */
  fr_cur = (double *)mxCalloc(N_f,sizeof(double));
  fi_cur = (double *)mxCalloc(N_f,sizeof(double));

  /* initialize equalizer */
  for (k=0; k<N_f; k++){
    fr_cur[k] = fr_init[k];
  }
  if( is_cmplx )
    for (k=0; k<N_f; k++){
      fi_cur[k] = fi_init[k];
    }

  /* adaptation routine */
  reg_ind = N_f-1;
  if( is_cmplx ){		/* if complex-valued... */
    for (i=0; i<N; ){
      reg_ind += 1+is_FSE;		/* update regressor index */

      /* calc equalizer output */
      yr_cur = 0;			
      yi_cur = 0;			
      for (k=0; k<N_f; k++){ 		
        yr_cur += rr[reg_ind-k]*fr_cur[k] - ri[reg_ind-k]*fi_cur[k];
        yi_cur += rr[reg_ind-k]*fi_cur[k] + ri[reg_ind-k]*fr_cur[k];
      }
      /* calc square-root CM error */
      srcme_cur = kappa-(yr_cur*yr_cur+yi_cur*yi_cur);	

      /* dithered decisions */
      if( yr_cur*srcme_cur + dr[i] > 0 )
        mu_yr_srcme = mu;		/* retain only sign of real err */
      else
        mu_yr_srcme = -mu;
      if( yi_cur*srcme_cur + di[i] > 0 )
        mu_yi_srcme = mu;		/* retain only sign of imag err */
      else
        mu_yi_srcme = -mu;

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
 	fr_cur[k] += rr[reg_ind-k]*mu_yr_srcme + ri[reg_ind-k]*mu_yi_srcme;
 	fi_cur[k] += rr[reg_ind-k]*mu_yi_srcme - ri[reg_ind-k]*mu_yr_srcme;
      }
      i++;				/* convenient to increment here */
      if( !(i % D_f) ) 
        for (k=0; k<N_f; k++){		/* store equalizer coefficients */
          fr[f_ind] = fr_cur[k];
          fi[f_ind++] = fi_cur[k];
        }
      if( !(i % D_e) ){
        yr[y_ind] = yr_cur;		/* store output */
        yi[y_ind] = yi_cur;
	srcme[y_ind++] = srcme_cur;	/* store srcme */
      }
    } 
  }
  else{				/* if real-valued... */
    for (i=0; i<N; ){
      reg_ind += 1+is_FSE;		/* update regressor index */

      yr_cur = 0;			
      for (k=0; k<N_f; k++){ 		/* calc equalizer output */
        yr_cur += rr[reg_ind-k]*fr_cur[k];
      }
      srcme_cur = kappa-yr_cur*yr_cur;	/* calc square-root CM error */

      /* dithered decision */
      if( yr_cur*srcme_cur + dr[i] > 0 )
        mu_yr_srcme = mu;		/* retain only sign of error */
      else
        mu_yr_srcme = -mu;

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
 	fr_cur[k] += rr[reg_ind-k] * mu_yr_srcme;
      }
      i++;				/* convenient to increment here */
      if( !(i % D_f) ) 
        for (k=0; k<N_f; k++){		/* store equalizer coefficients */
          fr[f_ind++] = fr_cur[k];
        }
      if( !(i % D_e) ){
        yr[y_ind] = yr_cur;		/* store ouput */
	srcme[y_ind++] = srcme_cur;	/* store srcme */
      }
    } 
  }

  /* store final outputs */
  for (k=0; k<N_f; k++){
    fr[N_f*(len_f-1)+k] = fr_cur[k];
  }
  yr[len_y-1] = yr_cur;
  srcme[len_y-1] = srcme_cur;
  if( is_cmplx ){
    for (k=0; k<N_f; k++){
      fi[N_f*(len_f-1)+k] = fi_cur[k];
    }
    yi[len_y-1] = yi_cur;
  }
}/* end of dsecma() */




/***************************************** 
 * Weekrackody's signed Godard algorithm *
 *****************************************
 *
 * the WSGA routine replaces the following matlab code:
 *
 * --------------------------------------------------
 * f_cur = f_init;
 * for i = 1:N,
 *   reg = r(1, N_f+(1+is_FSE)*i:-1:(1+is_FSE)*i+1);
 *   y_cur = reg*f_cur;
 *   srcme_cur = kappa-abs(y_cur)^2;
 *   f_cur = f_cur + mu*reg'*sign(kappa-abs(real(y_cur))-abs(imag(y_cur)))*...
 *              (sign(real(y_cur))+j*sign(imag(y_cur)));
 *   f(:,floor((i-1)/D_f)+1) = f_cur;
 *   y(floor((i-1)/D_e)+1) = y_cur;
 *   srcme(floor((i-1)/D_e)+1) = srcme_cur;
 * end
 * --------------------------------------------------
 */
void wsga(double *fr, double *fi, double *yr, double *yi, double *srcme, 
         double *rr, double *ri, double *fr_init, double *fi_init,
	 double kappa, double mu,
         int N, int N_f, int D_e, int D_f, int is_FSE, bool is_cmplx,
	 int len_f, int len_y)
{
  double *fr_cur, *fi_cur;
  double yr_cur, yi_cur, srcme_cur, mu_yr_srcme, mu_yi_srcme, sgn_yr, sgn_yi;
  int    i, k, f_ind=0, y_ind=0, reg_ind; 

  /* allocate local memory */
  fr_cur = (double *)mxCalloc(N_f,sizeof(double));
  fi_cur = (double *)mxCalloc(N_f,sizeof(double));

  /* initialize equalizer */
  for (k=0; k<N_f; k++){
    fr_cur[k] = fr_init[k];
  }
  if( is_cmplx )
    for (k=0; k<N_f; k++){
      fi_cur[k] = fi_init[k];
    }

  /* adaptation routine */
  reg_ind = N_f-1;
  if( is_cmplx ){		/* if complex-valued... */
    for (i=0; i<N; ){
      reg_ind += 1+is_FSE;		/* update regressor index */

      /* calc equalizer output */
      yr_cur = 0;			
      yi_cur = 0;			
      for (k=0; k<N_f; k++){ 		
        yr_cur += rr[reg_ind-k]*fr_cur[k] - ri[reg_ind-k]*fi_cur[k];
        yi_cur += rr[reg_ind-k]*fi_cur[k] + ri[reg_ind-k]*fr_cur[k];
      }
      /* calc square-root CM error */
      srcme_cur = kappa-(yr_cur*yr_cur+yi_cur*yi_cur);	
  
      /* calculate WSGA update */ 
      if( yr_cur > 0 ) sgn_yr = 1.0; else sgn_yr = -1.0;
      if( yi_cur > 0 ) sgn_yi = 1.0; else sgn_yi = -1.0;
      if( kappa - sgn_yr*yr_cur - sgn_yi*yi_cur > 0 ){
        mu_yr_srcme = mu*sgn_yr;	
        mu_yi_srcme = mu*sgn_yi;	
      }
      else{
        mu_yr_srcme = -mu*sgn_yr;
        mu_yi_srcme = -mu*sgn_yi;	
      }
      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
 	fr_cur[k] += rr[reg_ind-k]*mu_yr_srcme + ri[reg_ind-k]*mu_yi_srcme;
 	fi_cur[k] += rr[reg_ind-k]*mu_yi_srcme - ri[reg_ind-k]*mu_yr_srcme;
      }
      i++;				/* convenient to increment here */
      if( !(i % D_f) ) 
        for (k=0; k<N_f; k++){		/* store equalizer coefficients */
          fr[f_ind] = fr_cur[k];
          fi[f_ind++] = fi_cur[k];
        }
      if( !(i % D_e) ){
        yr[y_ind] = yr_cur;		/* store output */
        yi[y_ind] = yi_cur;
	srcme[y_ind++] = srcme_cur;	/* store srcme */
      }
    } 
  }
  else{				/* if real-valued (same as SE-CMA)... */
    for (i=0; i<N; ){
      reg_ind += 1+is_FSE;		/* update regressor index */

      yr_cur = 0;			
      for (k=0; k<N_f; k++){ 		/* calc equalizer output */
        yr_cur += rr[reg_ind-k]*fr_cur[k];
      }
      srcme_cur = kappa-yr_cur*yr_cur;	/* calc square-root CM error */
      if( yr_cur*srcme_cur > 0 )
        mu_yr_srcme = mu;		/* retain only sign of error */
      else
        mu_yr_srcme = -mu;

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
 	fr_cur[k] += rr[reg_ind-k] * mu_yr_srcme;
      }
      i++;				/* convenient to increment here */
      if( !(i % D_f) ) 
        for (k=0; k<N_f; k++){		/* store equalizer coefficients */
          fr[f_ind++] = fr_cur[k];
        }
      if( !(i % D_e) ){
        yr[y_ind] = yr_cur;		/* store ouput */
	srcme[y_ind++] = srcme_cur;	/* store srcme */
      }
    } 
  }

  /* store final outputs */
  for (k=0; k<N_f; k++){
    fr[N_f*(len_f-1)+k] = fr_cur[k];
  }
  yr[len_y-1] = yr_cur;
  srcme[len_y-1] = srcme_cur;
  if( is_cmplx ){
    for (k=0; k<N_f; k++){
      fi[N_f*(len_f-1)+k] = fi_cur[k];
    }
    yi[len_y-1] = yi_cur;
  }
}/* end of wsga() */




/********************** 
 * the MATLAB gateway *
 **********************/
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
 /**********************************************************************
  * input args: {N,N_f,D_f,D_e,is_FSE,kappa,mu,r,f_init,alg_num,misc}  *
  * output args: {f,y,srcme}                                           *
  **********************************************************************/

  bool	 is_cmplx;
  int 	 N, N_f, D_e, D_f, is_FSE, len_y, len_f, alg_num, D_u;
  double kappa, mu, alpha; 
  double *rr, *ri, *fr, *fi, *yr, *yi, *srcme, *fr_init, *fi_init; 

  /* check for the proper number of arguments */
  if(nrhs!=11)
    mexErrMsgTxt("Eleven inputs required.");
  if(nlhs!=3)
    mexErrMsgTxt("Three outputs required.");

  /* check for sizes on input arguments */
  if( mxGetM(prhs[0])!=1 || mxGetN(prhs[0])!=1 || !mxIsDouble(prhs[0]) )
    mexErrMsgTxt("First argument must be a scalar: N");
  if( mxGetM(prhs[1])!=1 || mxGetN(prhs[1])!=1 || !mxIsDouble(prhs[1]) )
    mexErrMsgTxt("Second argument must be a scalar: N_f");
  if( mxGetM(prhs[2])!=1 || mxGetN(prhs[2])!=1 || !mxIsDouble(prhs[2]) )
    mexErrMsgTxt("Third argument must be a scalar: D_f");
  if( mxGetM(prhs[3])!=1 || mxGetN(prhs[3])!=1 || !mxIsDouble(prhs[3]) )
    mexErrMsgTxt("Fourth argument must be a scalar: D_e");
  if( (int)mxGetScalar(prhs[4])!=0 && (int)mxGetScalar(prhs[4])!=1 )
    mexErrMsgTxt("Fifth argument must be in set {0,1}: is_FSE");
  if( mxGetM(prhs[5])!=1 || mxGetN(prhs[5])!=1 || !mxIsDouble(prhs[5]) )
    mexErrMsgTxt("Sixth argument must be a scalar: kappa");
  if( mxGetM(prhs[6])!=1 || mxGetN(prhs[6])!=1 || !mxIsDouble(prhs[6]) )
    mexErrMsgTxt("Seventh argument must be a scalar: mu");
  if( mxGetM(prhs[7])!=1 || !mxIsDouble(prhs[7]) )
    mexErrMsgTxt("Eighth argument must be a row vector: r");
  if( mxGetN(prhs[8])!=1 || !mxIsDouble(prhs[8]) )
    mexErrMsgTxt("Ninth argument must be a column vector: f_init");
  if( (int)mxGetScalar(prhs[9])<1 || (int)mxGetScalar(prhs[9])>8 )
    mexErrMsgTxt("Tenth argument must be an integer: 1 <= alg_num <= 8");

  /* retrieve scalars */ 
  N = (int)mxGetScalar(prhs[0]);
  N_f = (int)mxGetScalar(prhs[1]);
  D_f = (int)mxGetScalar(prhs[2]);
  D_e = (int)mxGetScalar(prhs[3]);
  is_FSE = (int)mxGetScalar(prhs[4]);
  kappa = (double)mxGetScalar(prhs[5]);
  mu = (double)mxGetScalar(prhs[6]);
  is_cmplx = mxIsComplex(prhs[7]);	/* true if r is complex */
  alg_num = (int)mxGetScalar(prhs[9]);

  /* set input pointer(s) */
  rr = mxGetPr(prhs[7]);
  fr_init = mxGetPr(prhs[8]);
  if( is_cmplx ){ 
    ri = mxGetPi(prhs[7]);
    if( mxIsComplex(prhs[8]) )		/* true if f_init is complex */
      fi_init = mxGetPi(prhs[8]);
    else
      fi_init = (double *)mxCalloc(N_f,sizeof(double));
  }

  /* determine decimated lengths */
  len_y = (N-1)/D_e + 1;
  len_f = (N-1)/D_f + 1;

  /* allocate memory for outputs */
  if( is_cmplx ){ 
    plhs[0] = mxCreateDoubleMatrix(N_f,len_f,mxCOMPLEX);	/* f */
    plhs[1] = mxCreateDoubleMatrix(1,len_y,mxCOMPLEX);   	/* y */
  }
  else{
    plhs[0] = mxCreateDoubleMatrix(N_f,len_f,mxREAL);    	/* f */
    plhs[1] = mxCreateDoubleMatrix(1,len_y,mxREAL);      	/* y */
  }
  plhs[2] = mxCreateDoubleMatrix(1,len_y,mxREAL); 		/* srcme */
  
  /* set output pointers */
  fr = mxGetPr(plhs[0]); 
  yr = mxGetPr(plhs[1]); 
  srcme = mxGetPr(plhs[2]); 
  if( is_cmplx ){ 
    fi = mxGetPi(plhs[0]); 
    yi = mxGetPi(plhs[1]);  
  }

  /* call algorithm subroutines */
  switch( alg_num ){
    case CMA: 
      D_u = (int)mxGetScalar(prhs[10]);
      cma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y, D_u); break;
    case NCMA: 
      alpha = (double)mxGetScalar(prhs[10]);
      ncma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y, alpha); break;
    case SECMA:
      secma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y); break;
    case SRCMA:
      srcma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y); break;
    case SSCMA:
      sscma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y); break;
    case DSECMA:
      alpha = (double)mxGetScalar(prhs[10]);
      dsecma(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y, alpha); break;
    case WSGA:
      wsga(fr, fi, yr, yi, srcme, rr, ri, fr_init, fi_init, kappa, 
     	  mu, N, N_f, D_e, D_f, is_FSE, is_cmplx, len_f, len_y); break;
    case CMAGD:
      mexErrMsgTxt("berg_algs.mexlx does not support CMA-GD."); break;
    default:
      mexErrMsgTxt("berg_algs.mexlx : You shouldn't be here!");
  }
}