berg_cmegd.c

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

}/* end of cmegd_cPAM() */




/*********************** 
 * cmplx-QAM CME-GD *
 ***********************
 *
 * the cmplx-channel (rotationally invariant noise) cmplx-source CME-GD 
 * 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;       % for constellation plot...
 *     h = C*f_cur;
 *     hTh = norm(h)^2;
 *     fTf = norm(f_cur)^2;
 *     srcme_cur = sqrt( (kappa-2)*sum(abs(h).^4) + 2*hTh^2 ...
 *	 + 2*sigma2_n^2*fTf^2 + 4*sigma2_n*hTh*fTf...
 *       - 2*kappa*(hTh + sigma2_n*fTf) + kappa^2 );
 *     dJdf = (kappa-2)*C'*(h.*(abs(h).^2))... 
 *       + (2*hTh + 2*sigma2_n*fTf - kappa)*(C'*h + sigma2_n*f_cur);
 *     f_cur = f_cur - mu*dJdf;
 *     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 cmegd_cQAM(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, int is_cmplx,
	 int len_f, int len_y, int N_h,
	 double sigma2_n, double *Cr, double *Ci)
{
  double *fr_cur, *fi_cur, *hr_cur, *hi_cur, *h2_cur, *f2_cur, *gradr, *gradi;
  double yr_cur, yi_cur, srcme_cur;
  double sqnrm_h, sqnrm_f, sum_h4, tmp;
  int    i, k, n, 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));
  f2_cur = (double *)mxCalloc(N_f,sizeof(double));
  hr_cur = (double *)mxCalloc(N_h,sizeof(double));
  hi_cur = (double *)mxCalloc(N_h,sizeof(double));
  h2_cur = (double *)mxCalloc(N_h,sizeof(double));
  gradr = (double *)mxCalloc(N_f,sizeof(double));
  gradi = (double *)mxCalloc(N_f,sizeof(double));

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

  /* adaptation routine */
  reg_ind = N_f-1;
  for (i=0; i<N; ){
    reg_ind += 1+is_FSE;		/* update regressor index */
    i++;				/* convenient to increment here */

    sqnrm_f = 0; 
    for (k=0; k<N_f; k++){ 		/* calc |f|^2 and ||f||^2 */
      f2_cur[k] = fr_cur[k]*fr_cur[k] + fi_cur[k]*fi_cur[k]; 
      sqnrm_f += f2_cur[k];
    }
    sqnrm_h = 0; 
    for (n=0; n<N_h; n++){ 		/* calc h, hth, |h|^2, & ||h||^2 */
      hr_cur[n] = 0;
      hi_cur[n] = 0;
      for (k=0; k<N_f; k++){ 
        hr_cur[n] += Cr[k*N_h+n]*fr_cur[k] - Ci[k*N_h+n]*fi_cur[k]; 
        hi_cur[n] += Cr[k*N_h+n]*fi_cur[k] + Ci[k*N_h+n]*fr_cur[k]; 
      } 
      h2_cur[n] = hr_cur[n]*hr_cur[n] + hi_cur[n]*hi_cur[n];
      sqnrm_h += h2_cur[n];
    }

    /* calculate and store output (y) and square-root CM error (srcme) */
    if( !(i % D_e) ){
      yr[y_ind] = 0;			
      yi[y_ind] = 0;			
      for (k=0; k<N_f; k++){ 		
        yr[y_ind] += rr[reg_ind-k]*fr_cur[k] - ri[reg_ind-k]*fi_cur[k];
        yi[y_ind] += rr[reg_ind-k]*fi_cur[k] + ri[reg_ind-k]*fr_cur[k];
      }
      sum_h4 = 0; for( n=0; n<N_h; n++ ){ sum_h4 += h2_cur[n]*h2_cur[n]; }
      srcme[y_ind++] = sqrt(kappa*kappa - 2*kappa*(sqnrm_h + sigma2_n*sqnrm_f)
  		+ 2*sigma2_n*sigma2_n*sqnrm_f*sqnrm_f 
        	+ 2*sqnrm_h*(sqnrm_h + 2*sigma2_n*sqnrm_f)
        	+ (kappa-2)*sum_h4);
    }

    /* calc CME gradient */
    tmp = 2*sqnrm_h + 2*sigma2_n*sqnrm_f - kappa;
    for (k=0; k<N_f; k++){
      gradr[k] = sigma2_n*fr_cur[k];
      gradi[k] = sigma2_n*fi_cur[k];
      for (n=0; n<N_h; n++){
        gradr[k] += Cr[k*N_h+n]*hr_cur[n] + Ci[k*N_h+n]*hi_cur[n];
        gradi[k] += Cr[k*N_h+n]*hi_cur[n] - Ci[k*N_h+n]*hr_cur[n];
      }
      gradr[k] *= tmp;
      gradi[k] *= tmp;
      for (n=0; n<N_h; n++){
        gradr[k] += (kappa-2)*h2_cur[n]*
		(hr_cur[n]*Cr[k*N_h+n] + hi_cur[n]*Ci[k*N_h+n]);
        gradi[k] += (kappa-2)*h2_cur[n]*
		(hi_cur[n]*Cr[k*N_h+n] - hr_cur[n]*Ci[k*N_h+n]);
      }
    }

    /* update equalizer coefficients */
    for (k=0; k<N_f; k++){		
      fr_cur[k] -= mu*gradr[k];
      fi_cur[k] -= mu*gradi[k];
    }

    /* store equalizer coefficients */
    if( !(i % D_f) ) 
      for (k=0; k<N_f; k++){		
        fr[f_ind] = fr_cur[k];
        fi[f_ind++] = fi_cur[k];
      }
  } 

  /* store final equalizer */
  for (k=0; k<N_f; k++){
    fr[N_f*(len_f-1)+k] = fr_cur[k];
    fi[N_f*(len_f-1)+k] = fi_cur[k];
  }

  /* store final y and srcme */
  yr[len_y-1] = 0;
  yi[len_y-1] = 0;
  for (k=0; k<N_f; k++){
    yr[len_y-1] += rr[reg_ind-k]*(fr_cur[k]+mu*gradr[k])
			- ri[reg_ind-k]*(fi_cur[k]+mu*gradi[k]);
    yi[len_y-1] += rr[reg_ind-k]*(fi_cur[k]+mu*gradi[k])
			+ ri[reg_ind-k]*(fr_cur[k]+mu*gradr[k]);
  }
  sum_h4 = 0; for( n=0; n<N_h; n++ ){ sum_h4 += h2_cur[n]*h2_cur[n]; }
  srcme[len_y-1] = sqrt(kappa*kappa - 2*kappa*(sqnrm_h + sigma2_n*sqnrm_f)
  		+ 2*sigma2_n*sigma2_n*sqnrm_f*sqnrm_f 
        	+ 2*sqnrm_h*(sqnrm_h + 2*sigma2_n*sqnrm_f)
        	+ (kappa-2)*sum_h4);

}/* end of cmegd_cQAM() */




/********************** 
 * 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, is_cmplx, *
  *              sigma2_n, C}  						   *
  * output args: {f, y, srcme}                                             *
  **************************************************************************/

  int 	 N, N_f, D_e, D_f, is_FSE, len_y, len_f, is_cmplx, N_h;
  double kappa, mu, alpha, sigma2_n; 
  double *rr, *ri, *fr, *fi, *yr, *yi, *srcme, *fr_init, *fi_init, *Cr, *Ci; 

  /* check for the proper number of arguments */
  if(nrhs!=12)
    mexErrMsgTxt("Twelve 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])<0 || (int)mxGetScalar(prhs[9])>2 )
    mexErrMsgTxt("Tenth argument must be an integer: 0 <= is_cmpx <= 2");
  if( mxGetM(prhs[10])!=1 || mxGetN(prhs[10])!=1 || !mxIsDouble(prhs[10]) )
    mexErrMsgTxt("Eleventh argument must be a scalar: sigma2_n");
  if( !mxIsDouble(prhs[11]) )
    mexErrMsgTxt("Twelvth argument must be a matrix: C");
  

  /* 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]);
  sigma2_n = (double)mxGetScalar(prhs[10]);
  N_h = mxGetM(prhs[11]);

  /* check sizes */
  if( N_f != mxGetN(prhs[11]) )
    mexErrMsgTxt("Column dimension on C is not N_f.");
  if( N_f != mxGetM(prhs[8]) )
    mexErrMsgTxt("Row dimension on f is not N_f.");

  /* complex-valued?: 0=real, 1=real source, cmplx chan, 2=cmplx */
  is_cmplx = (int)mxGetScalar(prhs[9]);

  /* set input pointer(s) */
  rr = mxGetPr(prhs[7]);
  fr_init = mxGetPr(prhs[8]);
  Cr = mxGetPr(prhs[11]);
  if( is_cmplx ){			/* when cmplx alg: */
    if( mxIsComplex(prhs[7]) ){			/* cmplx r */
      ri = mxGetPi(prhs[7]);
      if( mxIsComplex(prhs[11]) )		/* ...cmplx C */
        Ci = mxGetPi(prhs[11]);
      else					/* ...real C */
        Ci = (double *)mxCalloc(N_f*N_h,sizeof(double));
    }
    else{					/* real r */
      ri = (double *)mxCalloc(N,sizeof(double));
      if( mxIsComplex(prhs[11]) )		/* ...cmplx C */
        mexErrMsgTxt("Real data can't come from a cmplx channel!");
      else					/* ...real C */
        Ci = (double *)mxCalloc(N_f*N_h,sizeof(double));
    }
    if( mxIsComplex(prhs[8]) )			/* ...cmplx f_init */
      fi_init = mxGetPi(prhs[8]);
    else					/* ...real f_init */
      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( is_cmplx ){
    case 0: 
      cmegd_rPAM(fr, yr, srcme, rr, fr_init, kappa, mu, N, N_f, D_e, 
	  D_f, is_FSE, is_cmplx, len_f, len_y, N_h, sigma2_n, Cr); 
      break;
    case 1: 
      cmegd_cPAM(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, N_h, sigma2_n, 
	  Cr, Ci); 
      break;
    case 2: 
      cmegd_cQAM(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, N_h, sigma2_n, 
	  Cr, Ci); 
      break;
    default:
      mexErrMsgTxt("berg_cmegd.mexlx : You shouldn't be here!");
  }
}