berg_algs.c

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

        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 secma() */




/**************************** 
 * signed-regressor CMA 2-2 *
 ****************************
 *
 * the SR-CMA 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*(sign(real(reg))+j*sign(imag(reg)))'*(y_cur*srcme_cur);
 *   %f_cur = f_cur + mu*(reg./abs(reg))'*(y_cur*srcme_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 srcma(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;
  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);	
      mu_yr_srcme = mu*yr_cur*srcme_cur;	
      mu_yi_srcme = mu*yi_cur*srcme_cur;	

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
        if (rr[reg_ind-k] > 0){
 	  fr_cur[k] += mu_yr_srcme;
 	  fi_cur[k] += mu_yi_srcme;
     	}
 	else{
 	  fr_cur[k] -= mu_yr_srcme;
 	  fi_cur[k] -= mu_yi_srcme;
	}
        if (ri[reg_ind-k] > 0){
 	  fr_cur[k] += mu_yi_srcme;
 	  fi_cur[k] -= mu_yr_srcme;
	}
 	else{
 	  fr_cur[k] -= mu_yi_srcme;
 	  fi_cur[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 */
      mu_yr_srcme = mu*yr_cur*srcme_cur;	

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
        if (rr[reg_ind-k] > 0)
 	  fr_cur[k] += mu_yr_srcme;
 	else
 	  fr_cur[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 srcma() */




/***************************************** 
 * signed-regressor signed-error CMA 2-2 *
 *****************************************
 *
 * the SS-CMA 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*(sign(real(reg))+j*sign(imag(reg)))'*...
 *     (sign(real(y_cur)*srcme_cur)+j*sign(imag(y_cur)*srcme_cur));
 *   %f_cur = f_cur+mu*(reg./abs(reg))'*(y_cur*srcme_cur)/abs(y_cur*srcme_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 sscma(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;
  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);	
      if( yr_cur*srcme_cur > 0 )
        mu_yr_srcme = mu;               /* retain only sign of real error */
      else
        mu_yr_srcme = -mu;
      if( yi_cur*srcme_cur > 0 )
        mu_yi_srcme = mu;               /* retain only sign of imag error */
      else
        mu_yi_srcme = -mu;

      for (k=0; k<N_f; k++){		/* update equalizer coefficients */
        if (rr[reg_ind-k] > 0){
 	  fr_cur[k] += mu_yr_srcme;
 	  fi_cur[k] += mu_yi_srcme;
     	}
 	else{
 	  fr_cur[k] -= mu_yr_srcme;
 	  fi_cur[k] -= mu_yi_srcme;
	}
        if (ri[reg_ind-k] > 0){
 	  fr_cur[k] += mu_yi_srcme;
 	  fi_cur[k] -= mu_yr_srcme;
	}
 	else{
 	  fr_cur[k] -= mu_yi_srcme;
 	  fi_cur[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 */
      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 */
        if (rr[reg_ind-k] > 0)
 	  fr_cur[k] += mu_yr_srcme;
 	else
 	  fr_cur[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 sscma() */




/********************************* 
 * dithered signed-error CMA 2-2 *
 *********************************
 *
 * the DSE-CMA routine replaces the following matlab code:
 *
 * --------------------------------------------------
 * dith_r = alpha*(2*rand(1,N)-1);
 * dith_i = alpha*(2*rand(1,N)-1);
 * mu = mu*alpha;
 * 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(real(y_cur)*srcme_cur+dith_r(i));
 *   if ~isreal(r),
 *     f_cur = f_cur + j*mu*reg'*sign(imag(y_cur)*srcme_cur+dith_i(i));
 *   end;
 *   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 dsecma(double *fr, double *fi, double *yr, double *yi, double *srcme, 
         double *rr, double *ri, double *fr_init, double *fi_init,
	 double kappa, double mu_without_alpha,
         int N, int N_f, int D_e, int D_f, int is_FSE, bool is_cmplx,
	 int len_f, int len_y, double alpha)
{
  double *fr_cur, *fi_cur;
  double yr_cur, yi_cur, srcme_cur, mu; 
  double mu_yr_srcme, mu_yi_srcme;
  int    i, k, f_ind=0, y_ind=0, reg_ind; 
  double *dr, *di;				/* dither records */
  int 	  num_in=2, num_out=1;			/* # i/o args for "rand" */
  mxArray *prhs[2], *plhs1[1], *plhs2[1];	/* arrays of matlab arrays */
  double  rand_in[2];				/* inputs to "rand": {N,1} */


  /* first create dither signal(s)... */
  /* set inputs to "rand" function */
   prhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);	/* make matlab array */
   prhs[1] = mxCreateDoubleMatrix(1,1,mxREAL);	/* make matlab array */
   rand_in[0] = N;				/* set numerical value */
   rand_in[1] = 1;				/* set numerical value */
   memcpy(mxGetPr(prhs[0]), rand_in, sizeof(double));  	/* copy values */
   memcpy(mxGetPr(prhs[1]), rand_in+1, sizeof(double));	/* copy values */
  /* call matlab's "rand" */
   mexCallMATLAB(num_out,plhs1,num_in,prhs,"rand");
   dr = mxGetPr(plhs1[0]);