berg_algs.c

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

/******************************************************************************
 *									      *
 * berg_algs.c - adaptation routines for the BERGulator (called by ui_traj.m) *
 *									      *
 * written by: 	Phil "Bert" Schniter 					      *
 *		Blind Equalization Research Group 			      *
 * 		Cornell University					      *
 *		9/18/97							      *
 *									      *
 * Copyright 1997-1998 Phil Schniter					      *
 *									      *
 ******************************************************************************/

#include "mex.h"
#include <math.h>

/* to be consistent with the user interface algorithm ordering... */
#define	 CMA 	1
#define	 NCMA 	2
#define	 SECMA  3
#define	 SRCMA  4
#define	 SSCMA  5
#define	 DSECMA 6
#define	 WSGA   7
#define	 CMAGD  8


/******************** 
 * standard CMA 2-2 *
 ********************
 *
 * the CMA 2-2 routine replaces the following matlab code:
 *
 * --------------------------------------------------
 * f_cur = f_init;
 * for i = D_u:D_u: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'*(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 cma(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, int D_u)
{
  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)*D_u;	/* 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 */
 	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 += D_u;				/* 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)*D_u;	/* 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 */
 	fr_cur[k] += rr[reg_ind-k] * mu_yr_srcme;
      }
      i += D_u;				/* 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 cma() */




/********************** 
 * normalized CMA 2-2 *
 **********************
 *
 * the NCMA routine replaces the following matlab code:
 *
 * --------------------------------------------------
 * f_cur = f_init;
 * alpha = 0.001;
 * 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/(alpha+norm(reg)^2)*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 ncma(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 alpha)
{
  double *fr_cur, *fi_cur;
  double yr_cur, yi_cur, srcme_cur, mu_yr_srcme, mu_yi_srcme, sqnrm_reg;
  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 and squared regressor norm */
  reg_ind = N_f-1;
  sqnrm_reg = 0;
  for (k=0; k<N_f; k++){
    fr_cur[k] = fr_init[k];
    sqnrm_reg += rr[reg_ind-k]*rr[reg_ind-k];
  }
  if( is_cmplx )
    for (k=0; k<N_f; k++){
      fi_cur[k] = fi_init[k];
      sqnrm_reg += ri[reg_ind-k]*ri[reg_ind-k];
    }

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

      /* update squared regressor norm */
      sqnrm_reg += rr[reg_ind]*rr[reg_ind] + ri[reg_ind]*ri[reg_ind]
	           - rr[reg_ind-N_f]*rr[reg_ind-N_f] 
                   - ri[reg_ind-N_f]*ri[reg_ind-N_f];
      if( is_FSE ){
        sqnrm_reg += rr[reg_ind-1]*rr[reg_ind-1] + ri[reg_ind-1]*ri[reg_ind-1]
  	             - rr[reg_ind-N_f-1]*rr[reg_ind-N_f-1]
  	             - ri[reg_ind-N_f-1]*ri[reg_ind-N_f-1];
      }
      /* 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/(alpha + sqnrm_reg);
      mu_yi_srcme = mu*yi_cur*srcme_cur/(alpha + sqnrm_reg);

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

      /* update squared regressor norm */
      sqnrm_reg += rr[reg_ind]*rr[reg_ind] 
	           - rr[reg_ind-N_f]*rr[reg_ind-N_f];
      if( is_FSE ){
        sqnrm_reg += rr[reg_ind-1]*rr[reg_ind-1] 
  	             - rr[reg_ind-N_f-1]*rr[reg_ind-N_f-1];
      }
      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/(alpha + sqnrm_reg);

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




/************************ 
 * signed-error CMA 2-2 *
 ************************
 *
 * the SE-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*reg'*...
 *      (sign(real(y_cur)*srcme_cur)+j*sign(imag(y_cur)*srcme_cur));
 *   %f_cur = f_cur + mu*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 secma(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 */
 	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 */