adapt.bak

基于Matlab

     outFr[i]=F.v[i].re;
     outFi[i]=F.v[i].im;
   }
   for (i=0; i<Ng; i++) {
     outGr[i]=G.v[i].re;
     outGi[i]=G.v[i].im;
   }
}




/* 
 * adaptFSE - update fractionally spaced equalizer using gradient
 * search algorithm (either LMS, CMA, or DD).
 */
void
adaptFSE(Nfa,Nfb,Ng,Fa0r,Fa0i,Fb0r,Fb0i,G0r,G0i,
         L,ar,ai,rar,rai,rbr,rbi,yr,yi,xr,xi,ahatr,ahati,
         muf,mug,offset,delay,pdelay,
         constellation,M,sig,algorithm,regfilt,structure,
         outFar,outFai,outFbr,outFbi,outGr,outGi,outyr,outyi,
         outxr,outxi,outzr,outzi,outahatr,outahati,
         outefr,outefi,outegr,outegi,outphi)
int Nfa,Nfb,Ng;
double *Fa0r,*Fa0i;
double *Fb0r,*Fb0i;
double *G0r,*G0i;
int L;
double *ar,*ai;
double *rar,*rai;
double *rbr,*rbi;
double *yr,*yi;
double *xr,*xi;
double *ahatr,*ahati;
double muf,mug;
int offset;
int delay,pdelay;
int constellation,M;
double sig;
int algorithm;
int regfilt;
int structure;
double *outFar,*outFai;
double *outFbr,*outFbi;
double *outGr,*outGi;
double *outyr,*outyi;
double *outxr,*outxi;
double *outzr,*outzi;
double *outahatr,*outahati;
double *outefr,*outefi;
double *outegr,*outegi;
double *outphi;
{ 
   complexVec Fa;             /* even forward filter */ 
   complexVec Fb;             /* odd forward filter */ 
   complexVec G;              /* feedback filter */
   complex *a;                /* complex symbols */
   complex *ra;               /* complex even received data */
   complex *rb;               /* complex odd received data */
   complex *y;                /* complex equalizer output */
   complex *x;                /* filtered equalizer output */
   complex *z;                /* complex equalizer phase corrected output */
   complex *ahat;             /* complex symbol estimates */
   complex q;                 /* decisions for phase recovery */
   complex ef;                /* forward filter update error */
   complex eg;                /* feedback filter update error */
   complex corr;              /* phase corrrection */ 
   complex ep0,ep1;           /* current and past phase error */ 
   double epmod;              /* modulus of phase error */ 
   double phi,c;              /* phase estimate and correction */ 
   static double a1=1.6;      /* loop filter coefficients */
   static double a0=-0.9;     /* loop filter coefficient */
   double g;                  /* CMA constant */
   double m4,m2;              /* Constellation 4th and 2nd moments */
   int i;                     /* Tap index */
   int n;                     /* sample counter */

   /* Get constellation variance and define CMA constant */
   constellationStats(constellation,M,sig,&m4,&m2);
   g=m4;

   /* Initialize initial guess */
   a=(complex *)mxCalloc(L,sizeof(complex));
   ra=(complex *)mxCalloc(L,sizeof(complex));
   rb=(complex *)mxCalloc(L,sizeof(complex));
   y=(complex *)mxCalloc(L,sizeof(complex));
   x=(complex *)mxCalloc(L,sizeof(complex));
   z=(complex *)mxCalloc(L,sizeof(complex));
   ahat=(complex *)mxCalloc(L,sizeof(complex));
   ep0.re=1.0; ep0.im=0.0;
   ep1.re=1.0; ep1.im=0.0;
   c=0;
   phi=0;
   for (n=0; n<L; n++) {
      a[n].re=ar[n];
      a[n].im=ai[n];
      ra[n].re=rar[n];
      ra[n].im=rai[n];
      rb[n].re=rbr[n];
      rb[n].im=rbi[n];
      y[n].re=yr[n];
      y[n].im=yi[n];
      x[n].re=xr[n];
      x[n].im=xi[n];
      ahat[n].re=ahatr[n];
      ahat[n].im=ahati[n];
   }

   /* Initialize filters */
   Fa.len=Nfa;
   for (i=0; i<Nfa; i++) {
      Fa.v[i].re=Fa0r[i];
      Fa.v[i].im=Fa0i[i];
   }

   Fb.len=Nfb;
   for (i=0; i<Nfb; i++) {
      Fb.v[i].re=Fb0r[i];
      Fb.v[i].im=Fb0i[i];
   }

   G.len=Ng;
   for (i=0; i<Ng; i++) {
      G.v[i].re=G0r[i];
      G.v[i].im=G0i[i];
   }

   /* Run adaptation until output variance is below threshold */ 
   for (n=offset; n<L; n++) { 

      y[n]=zero;
      /* Compute equalizer output */
      for (i=0; i<Nfa; i++) {
         y[n].re=y[n].re+(Fa.v[i].re*ra[n-i].re-Fa.v[i].im*ra[n-i].im);
         y[n].im=y[n].im+(Fa.v[i].re*ra[n-i].im+Fa.v[i].im*ra[n-i].re);
      }
      for (i=0; i<Nfb; i++) {
         y[n].re=y[n].re+(Fb.v[i].re*rb[n-i].re-Fb.v[i].im*rb[n-i].im);
         y[n].im=y[n].im+(Fb.v[i].re*rb[n-i].im+Fb.v[i].im*rb[n-i].re);
      }
      if (structure==1)
         for (i=0; i<Ng; i++) {
            y[n].re=y[n].re-(G.v[i].re*y[n-pdelay-i].re-
                             G.v[i].im*y[n-pdelay-i].im);
            y[n].im=y[n].im-(G.v[i].re*y[n-pdelay-i].im+
                             G.v[i].im*y[n-pdelay-i].re);
         }
      else
         for (i=0; i<Ng; i++) {
            y[n].re=y[n].re-(G.v[i].re*ahat[n-pdelay-i].re-
                             G.v[i].im*ahat[n-pdelay-i].im);
            y[n].im=y[n].im-(G.v[i].re*ahat[n-pdelay-i].im+
                             G.v[i].im*ahat[n-pdelay-i].re);
         }

      /* Use regressor filtering if desired */
      x[n]=y[n-pdelay];
      if ((regfilt==1) & (structure==1)) 
         for (i=0; i<Ng; i++) {
            x[n].re=x[n].re-(G.v[i].re*x[n-pdelay-i].re-
                             G.v[i].im*x[n-pdelay-i].im);
            x[n].im=x[n].im-(G.v[i].re*x[n-pdelay-i].im+
                             G.v[i].im*x[n-pdelay-i].re);
         } 


      /* LMS */
      if (algorithm==1)  {
         ahat[n]=a[n-delay]; 
         ef.re=y[n].re-ahat[n].re; ef.im=y[n].im-ahat[n].im;
         eg=ef;

         /* Phase recovered output */
         z[n].re=y[n].re;
         z[n].im=y[n].im;
      }
      /* DD */
      else if (algorithm==2) {
         ahat[n]=quantize(y[n],constellation,M,sig,n); 
         ef.re=y[n].re-ahat[n].re; ef.im=y[n].im-ahat[n].im;
         eg=ef;

         /* Phase recovered output */
         z[n].re=y[n].re;
         z[n].im=y[n].im;
      }
      /* CMA */
      else if (algorithm==3) {
         ef.re=y[n].re*(y[n].re*y[n].re+y[n].im*y[n].im-g);
         ef.im=y[n].im*(y[n].re*y[n].re+y[n].im*y[n].im-g);
         eg=ef;

         /* Phase recovered output */
         corr.re=cos(phi); corr.im=-sin(phi);
         z[n].re=y[n].re*corr.re-y[n].im*corr.im;
         z[n].im=y[n].im*corr.re+y[n].re*corr.im;
         ahat[n]=quantize(z[n],constellation,M,sig,n); 

         q=quantize(z[n],constellation,M,sig,n); 

         ep1.re=z[n].re*q.re+z[n].im*q.im;
         ep1.im=-z[n].re*q.im+z[n].im*q.re;
         epmod=(double)sqrt(ep1.re*ep1.re+ep1.im*ep1.im); 
         ep1.re=ep1.re/epmod;
         ep1.im=ep1.im/epmod;

         c=c+a1*ep1.im+a0*ep0.im;
         ep0=ep1;
         phi=phi+c;
      }


      /* Update equalizers */
      for (i=0; i<Nfa; i++) {
         Fa.v[i].re=Fa.v[i].re-muf*(ef.re*ra[n-i].re+ef.im*ra[n-i].im);
         Fa.v[i].im=Fa.v[i].im-muf*(-ef.re*ra[n-i].im+ef.im*ra[n-i].re);
      }
      for (i=0; i<Nfb; i++) {
         Fb.v[i].re=Fb.v[i].re-muf*(ef.re*rb[n-i].re+ef.im*rb[n-i].im);
         Fb.v[i].im=Fb.v[i].im-muf*(-ef.re*rb[n-i].im+ef.im*rb[n-i].re);
      }
      if (structure==1)
         for (i=0; i<Ng; i++) {
            G.v[i].re=G.v[i].re+mug*(eg.re*x[n-i].re+
                                     eg.im*x[n-i].im);
            G.v[i].im=G.v[i].im+mug*(-eg.re*x[n-i].im+
                                      eg.im*x[n-i].re);
         }
      else
         for (i=0; i<Ng; i++) {
           G.v[i].re=G.v[i].re+mug*(eg.re*ahat[n-pdelay-i].re+
                                    eg.im*ahat[n-pdelay-i].im);
           G.v[i].im=G.v[i].im+mug*(-eg.re*ahat[n-pdelay-i].im+
                                     eg.im*ahat[n-pdelay-i].re);
         }

      /* Store output & error */
      outyr[n]=y[n].re;
      outyi[n]=y[n].im;
      outxr[n]=x[n].re;
      outxi[n]=x[n].im;
      outzr[n]=z[n].re;
      outzi[n]=z[n].im;
      outahatr[n]=ahat[n].re;
      outahati[n]=ahat[n].im;
      outefr[n]=ef.re;
      outefi[n]=ef.im; 
      outegr[n]=ef.re;
      outegi[n]=ef.im; 
      outphi[n]=phi;

      /* Check for instability */
      if (ef.re*ef.re+ef.im*ef.im>100000 | eg.re*eg.re+eg.im*eg.im>100000)
         mexErrMsgTxt("Adaptation blowing up! Stopping adaptation.");

   } 


   /* Store final filters */
   for (i=0; i<Nfa; i++) {
     outFar[i]=Fa.v[i].re;
     outFai[i]=Fa.v[i].im;
   }
   for (i=0; i<Nfb; i++) {
     outFbr[i]=Fb.v[i].re;
     outFbi[i]=Fb.v[i].im;
   }
   for (i=0; i<Ng; i++) {
     outGr[i]=G.v[i].re;
     outGi[i]=G.v[i].im;
   }
}




/* 
 * addComplex - add complex numbers a and b 
 */ 
complex
addComplex(a,b)
complex a;
complex b;
{
   complex c;

   c.re=a.re+b.re;
   c.im=a.im+b.im;

   return c;
}


/* 
 * subComplex - subtract complex numbers a and b 
 */ 
complex
subComplex(a,b)
complex a;
complex b;
{
   complex c;

   c.re=a.re-b.re;
   c.im=a.im-b.im;

   return c;
}

 
/* 
 * multComplex - multiply complex numbers a and b 
 */ 
complex
multComplex(a,b)
complex a;
complex b;
{
   complex c;

   c.re=a.re*b.re-a.im*b.im;
   c.im=a.re*b.im+a.im*b.re;

   return c;
}


/* 
 * modulus - modulus |z| of z 
 */
double
modulus(z)
complex z;
{
   return sqrt(z.re*z.re + z.im*z.im);
}


/* 
 * conj - conjugate of a=(x,y) is (x,-y) 
 */
complex
conj(a)
complex a;
{
   complex c;

   c.re=a.re;
   c.im=-a.im;

   return c;
}





/*
 * quantize - make a symbol estimate.
 */
complex
quantize(z,constellation,M,sig,n)
complex z;
int constellation;
int M;
double sig;
int n;
{
   complex a;
   complex y;

   y.re=z.re*sig;
   y.im=z.im*sig;
   if (constellation==1)
      a=unpam(y,M);
   else if (constellation==2)
      a=unqam(y,M);
   else if (constellation==3) {
      if (n&1) {
         a=unpam(y,M);
      }
      else {
         y.re=y.im;
         y.im=0.0;
         a=unpam(y,M);
         a.im=a.re;
         a.re=0.0;
      }
   }
   a.re=a.re/sig;
   a.im=a.im/sig;


   return a;
}

     



/*
 * constellationStats - 4th and 2nd moments of constellation
 */
void
constellationStats(constellation,M,sig,m4,m2)
int constellation;
int M;
double sig;
double *m4;
double *m2;
{
   int i,j;
   complex a;
   int rootM;

   if (constellation==1) {
      /* Compute 4th moment */
      *m4=0;
      for (i=1; i<=M/2; i++) {
         *m4=*m4+2.0/(double)M*pow((double)(2*i-1)/sig,4.0);
      }

      /* Compute 2nd moment */
      *m2=((double)M*M-1.0)/3.0;
   }
   else if (constellation==2) {
      /* Compute 4th moment */
      *m4=0;
      rootM=sqrt((double) M);
      for (i=1; i<=rootM/2; i++)  {
         for (j=1; j<=rootM/2; j++)  {
            a.re=(2*i-1)/sig;
            a.im=(2*j-1)/sig;
            *m4=*m4+4.0/(double)M*pow(modulus(a),4.0);
         }
      }

      /* Compute 2nd moment */
      *m2=(2.0*((double)M-1.0))/3.0;
   }
}






/*
 * unpam - decode an M-PAM symbol to nearest neighbor.
 */
complex
unpam(z,M)
complex z;
int M;
{
   complex a=zero;

   if (modulus(z)>M-1.0) {
      a.re=sign(z.re)*((double)M-1.0);
   }
   else {
      a.re=sign(z.re)*(round((modulus(z)+1.0)/2.0)*2.0-1.0);
   }

   return a;
}


/*
 * unqam - make a M-QAM estimate.
 */
complex
unqam(z,M)
complex z;
int M;
{
   complex zre=zero;
   complex zim=zero;
   complex a,are,aim;
   int rootM;

   rootM=sqrt(M);

   zre.re=z.re;
   are=unpam(zre,rootM);
   zim.re=z.im;
   aim=unpam(zim,rootM);

   a.re=are.re;
   a.im=aim.re;

   return a;
}


/*
 * sign - signum function
 */
double
sign(x)
double x;
{
   double y;
 
   if (x>0.0)
      y=1.0;
   else
      y=-1.0;

   return y;
}


/*
 * round - round to nearest integer
 */
double
round(x)
double x;
{
   double y;

   if (x>((double)floor(x)+(double)ceil(x))/2.0)
      y=(double)ceil(x);
   else
      y=(double)floor(x);

   return y;
}