remez.c

dsp AD公司ADSP21的代码,里面有FFT FIR IIR EQULIZER G722_21F 等可以在项目中直接应用的代码.此代码的来源是ADI公司自己出版的书籍,此书在美国购得

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

/**************************************************************************

REMEZ.C - PROGRAM TO MAKE FIR FILTER COEFFICIENTS USING
          REMEZ EXCHANGE FIR FILTER DESIGN PROGRAM

 Translated into C.
 Based on FORTRAN program by Parks McClellan as
 listed in Prentice Hall book by Rabiner & Gold
 "Theory and Application of Digital Signal Processing"
 and BASIC translation by N. Loy in Prentice Hall book
 "An Engineer's Guide to FIR Digital Filters"

INPUTS: NUMBER OF COEFFICIENTS, SAMPLING RATE, CUTOFF FREQUENCY.

OUTPUTS: FILTER COEFFICIENTS, DSP DATA FILE WITH MAGNITUDE

*************************************************************************/

int get_int(char *title_string,int low_limit,int up_limit);
double get_float(char *title_string,double low_limit,double up_limit);

#define PI  (3.141592653589793)
#define PI2 (6.283185307179586)
#define NFMAX (128)  /* maximum filter length */
#define EPS (1.e-30)  /* small number */
#define ITRMAX (25)  /* maximum number of iterations */
#define NBMAX  (10)  /* maximum number of bands */
#define MAG_LENGTH 512    /* size of magnitude vector to write */

int nfcns,ngrid;
int iext[(NFMAX/2+2)+1];
double alpha[(NFMAX/2+2)+1],des[16*(NFMAX/2+2)+1];
double grid[16*(NFMAX/2+2)+1];
double wt[16*(NFMAX/2+2)+1];
double dev;
double ad[(NFMAX/2+2)+1],x[(NFMAX/2+2)+1],y[(NFMAX/2+2)+1];

main()
{
    void rerror(),remez();
    double d(),gee();

    int i,j,k,l;
    int keyboard, example1, example2, example3, example4;
    int jtype,kup,lband,lgrid,nfilt,nbands;
    int jb,neg,nodd,nm1,nz;
    double delf,fup,temp,change,tempflt;
    char *cptr,*trailer;
    float *htotal,*hin;
    double *mag;

    double deviat[NBMAX+1],wtx[NBMAX+1];
    double edge[2*NBMAX+1],fx[NBMAX+1];

    static double h[(NFMAX/2+2)+1];

    FILE *chan;

keyboard = 0;
example1 = 0;
example2 = 0;
example3 = 0;
example4 = 0;

printf("  ...   REMEZ EXCHANGE FIR FILTER DESIGN PROGRAM ... \n\n");

   printf("1:  EXAMPLE1 --- LOWPASS FILTER\n");
   printf("2:  EXAMPLE2 --- BANDPASS FILTER\n");
   printf("3:  EXAMPLE3 --- DIFFERENTIATOR\n");
   printf("4:  EXAMPLE4 --- HILBERT TRANSFORMER\n");
   printf("5:  KEYBOARD --- GET INPUT PARAMETERS FROM KEYBOARD\n");
   i = get_int("selection",1,5);

   switch(i){
      case 1:  example1 = 1;
               break;
      case 2:  example2 = 1;
               break;
      case 3:  example3 = 1;
               break;
      case 4:  example4 = 1;
               break;
      case 5:  keyboard = 1;
               break;
   }

if(keyboard){     /* GET DATA FROM KEYBOARD */

/* GET INPUT DATA FROM USER */
nfilt = get_int("number of coefficients",3,NFMAX);
printf("\nFilter types are:  1=Bandpass, 2=Differentiator, 3=Hilbert\n");
jtype = get_int("filter type",1,3);
nbands = get_int("number of bands",1,NBMAX);

lgrid = 16;

jb = 2 * nbands;
printf("Now inputting edge (corner) frequencies for %d band edges\n",jb);
edge[0]=0.0;
for(i=1;i<=jb;i++){
  char temps[80];
  sprintf(temps,"edge frequency for edge (corner) # %d",i);
  edge[i] = get_float(temps,edge[i-1],0.5);
}

for(i=1;i<=nbands;i++){
  char temps[80];
  sprintf(temps,"gain of band # %d",i);
  fx[i] = get_float(temps,0.0,1.e3);
  sprintf(temps,"weight of band # %d",i);
  wtx[i] = get_float(temps,0.1,100.0);
}
}  /* end if(keyboard) */

else if(example1){      /* LOWPASS FILTER */
nfilt = 24;
jtype = 1;
nbands = 2;
lgrid = 16;
edge[1] = 0;
edge[2] = 0.08;
edge[3] = 0.16;
edge[4] = 0.5;
fx[1] = 1;
fx[2] = 0;
wtx[1] = 1;
wtx[2] = 1;
}

else if(example2){         /* BANDPASS FILTER */
nfilt = 24;
jtype = 1;
nbands = 3;
lgrid = 16;
edge[1] = 0;
edge[2] = 0.1;
edge[3] = 0.2;
edge[4] = 0.35;
edge[5] = 0.425;
edge[6] = 0.5;
fx[1] = 0;
fx[2] = 1;
fx[3] = 0;
wtx[1] = 2;
wtx[2] = 1;
wtx[3] = 10;
}

else if(example3){          /* DIFFERENTIATOR */
nfilt = 32;
jtype = 2;
nbands = 1;
lgrid = 16;
edge[1] = 0;
edge[2] = 0.5;
fx[1] = 1;
wtx[1] = 1;
}

else if(example4){         /* HILBERT TRANSFORMER */
nfilt = 20;
jtype = 3;
nbands = 1;
lgrid = 16;
edge[1] = 0.05;
edge[2] = 0.5;
fx[1] = 1;
wtx[1] = 1;
}
/* END INPUT DATA SECTION */

   if(nfilt > NFMAX){
      printf("Error: # coeff > %d\n...now exiting to system...\n",NFMAX);
      exit(1);
   }

   if(nfilt < 3) rerror("Error: # coeff < 3");
   if(nbands <= 0) nbands = 1;
   if(lgrid <= 0) lgrid = 16;

   printf("#coeff = %d\nType = %d\n#bands =  %d\nGrid = %d\n",
      nfilt,jtype,nbands,lgrid);
   for(i=1;i<=2*nbands;i++) printf("E[%d] = %.2lf\n",i,edge[i]);
   for(i=1;i<=nbands;i++) printf("Gain, wt[%d] = %.2lf  %.2lf\n",i,fx[i],wtx[i]);
   putchar('\n');

   neg = 1;
   if(jtype == 1) neg = 0;
   nodd  = nfilt % 2;
   nfcns = nfilt / 2;
   if(nodd == 1 && neg == 0) nfcns++;

/* SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
/* IS (FILTER LENGTH + 1)*GRID DENSITY/2 */
   grid[1] = edge[1];
   delf = lgrid * nfcns;
   delf = 0.5/delf;
   if(neg != 0 && edge[1] < delf) grid[1] = delf;

/* CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
/* FUNCTION ON THE GRID */
j = 1;
l = 1;
lband = 1;
while(1){   /* loop until break */
   fup = edge[l+1];
   do{
      temp = grid[j];
      des[j] = fx[lband];
      if(jtype == 2) des[j] *= temp;
      wt[j] = wtx[lband];
      if(jtype == 2 && fx[lband] >= 0.0001) wt[j] /= temp;
      j++;
      grid[j] = temp + delf;
   }while(grid[j] <= fup);

   grid[j-1] = fup;
   des[j-1] = fx[lband];
   if(jtype == 2) des[j-1] *= fup;
   wt[j-1] = wtx[lband];
   if(jtype == 2 && fx[lband] >= 0.0001) wt[j-1] /= fup;
   lband++;
   l += 2;
   if(lband > nbands) break;
   grid[j] = edge[l];
}  /* END while(1) */

ngrid = j-1;
if(neg == nodd && grid[ngrid] > (0.5 - delf)) ngrid--;

/* SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
/* TO THE ORIGINAL PROBLEM */
if(neg <= 0){
   if(nodd == 1);    /* DO NOTHING */
   else{
      for(j=1;j <= ngrid;j++){
         change = cos(PI * grid[j]);
         wt[j] *= change;
         if(change == 0.0) change = EPS;
         des[j] /= change;
      }
   }
}
else{
   if(nodd == 1){
      for(j=1;j <= ngrid;j++){
         change = sin(PI2 * grid[j]);
         wt[j] *= change;
         if(change == 0.0) change = EPS;
         des[j] /= change;
      }
   }
   else{
      for(j=1;j <= ngrid;j++){
         change = sin(PI * grid[j]);
         wt[j] *= change;
         if(change == 0.0) change = EPS;
         des[j] /= change;
      }
   }
}

/* INITIAL GUESS FOR THE EXTREMAL FREQUENCIES -- EQUALLY
/* SPACED ALONG THE GRID */
temp = (ngrid - 1)/nfcns;
for(j=1;j<=nfcns;j++) iext[j] = (j-1) * temp + 1;
iext[nfcns+1] = ngrid;
nm1 = nfcns - 1;
nz = nfcns + 1;

/* CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM */
remez(edge,nbands);

if(neg <= 0){
   if(nodd == 0){
      h[1] = 0.25 * alpha[nfcns];
      for(j=2;j<=nm1;j++)
         h[j] = 0.25 * (alpha[nz - j] + alpha[nfcns + 2 - j]);
      h[nfcns] = 0.5 * alpha[1] + 0.25 * alpha[2];
   }
   else{
      for(j=1;j<=nm1;j++) h[j] = 0.5 * alpha[nz - j];
      h[nfcns] = alpha[1];
   }
}
else{
   if(nodd == 0){
      h[1] = 0.25 * alpha[nfcns];
      for(j=2;j<=nm1;j++)
         h[j] = 0.25 * (alpha[nz - j] - alpha[nfcns + 2 - j]);
      h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[2];
   }
   else{
      h[1] = 0.25 * alpha[nfcns];
      h[2] = 0.25 * alpha[nm1];
      for(j=3;j<=nm1;j++)
         h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns + 3 - j]);
      h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[3];
      h[nz] = 0.0;
   }
}

/* PROGRAM OUTPUT SECTION */

putchar('\n');
for(j=1;j<=70;j++) putchar('*');
putchar('\n');
printf("                        FINITE IMPULSE RESPONSE (FIR)\n");
printf("                        LINEAR PHASE DIGITAL FILTER DESIGN\n");
printf("                        REMEZ EXCHANGE ALGORITHM\n");
switch(jtype){
   case 1: printf("                        BANDPASS FILTER\n");
           break;
   case 2: printf("                        DIFFERENTIATOR\n");
           break;
   case 3: printf("                        HILBERT TRANSFORMER\n");
           break;
}
printf("              FILTER LENGTH = %d\n",nfilt);
printf("              ***** IMPULSE RESPONSE *****\n");
for(j=1;j<=nfcns;j++){
   k = nfilt + 1 - j;
   if(neg == 0)
      printf("                   H(%3d) = %17.9e = H(%4d)\n",j,h[j],k);
   if(neg == 1)
      printf("                   H(%3d) = %17.9e = -H(%4d)\n",j,h[j],k);
}
if(neg == 1 && nodd == 1)
      printf("                   H(%3d) =  0.0\n",nz);   /* nz = nfcns + 1 */
putchar('\n');
for(k=1;k<=nbands;k+=4){
   kup = k + 3;
   if(kup > nbands) kup = nbands;
   for(j=1;j<=23;j++) putchar(' ');
   for(j=k;j<=kup;j++) printf("BAND%2d         ",j);
   putchar('\n');
   printf(" LOWER BAND EDGE");
   for(j=k;j<=kup;j++) printf("%15.8f",edge[2*j-1]);
   putchar('\n');
   printf(" UPPER BAND EDGE");
   for(j=k;j<=kup;j++) printf("%15.8f",edge[2*j]);
   putchar('\n');
   if(jtype != 2){
      printf(" DESIRED VALUE  ");
      for(j=k;j<=kup;j++) printf("%15.8f",fx[j]);
      putchar('\n');
   }
   if(jtype == 2){
      printf(" DESIRED SLOPE  ");
      for(j=k;j<=kup;j++) printf("%15.8f",fx[j]);
      putchar('\n');
   }
   printf(" WEIGHTING      ");
   for(j=k;j<=kup;j++) printf("%15.8f",wtx[j]);
   putchar('\n');
   for(j=k;j<=kup;j++) deviat[j] = dev/wtx[j];
   printf(" DEVIATION      ");
   for(j=k;j<=kup;j++) printf("%15.8f",deviat[j]);
   putchar('\n');
   if(jtype != 1) continue;
   for(j=k;j<=kup;j++) {
    if(fx[j] > 1e-5) deviat[j] = 20.0 * log10(deviat[j]);
    else  deviat[j] = -40.0 * log10(1.0-deviat[j]);
   }
   printf(" DEVIATION IN DB");
   for(j=k;j<=kup;j++) printf("%15.8f",deviat[j]);
   putchar('\n');
}
putchar('\n');
printf(" EXTREMAL FREQUENCIES\n");
putchar(' ');
for(j=1;j<=nz;j++){
   printf("%12.7f",grid[iext[j]]);
   if(!(j%5)){ putchar('\n');putchar(' ');}  /* CR every 5 columns */
}
putchar('\n');

  htotal = (float *) calloc(2*MAG_LENGTH, sizeof(float));
  if(!htotal){
    printf("\n  Unable to allocate output memory.\n");
    exit(1);
  }
  hin = (float *) calloc(2*MAG_LENGTH, sizeof(float));
  if(!hin){
    printf("\n  Unable to allocate output memory.\n");
    exit(1);
  }

k=0;
for(j=1;j<=nfcns;j++){
   htotal[k++] = h[j];
}
if(neg == 1 && nodd == 1) htotal[k++] = 0.0;
if(nodd == 1) j=nfcns-1; else j=nfcns;

for(;j>=1;j--){
   if(neg == 0) htotal[k++] = h[j];
   if(neg == 1) htotal[k++] = -h[j];
}
for(i=0; i < k ; i++) hin[i]=htotal[i];

    printf("\nFIR coefficients written to text file COEF.DAT\n");
    chan=fopen("coef.dat","wt");
    for(i = 0 ; i < k ; i++) {
        fprintf(chan,"\n%15.10f",htotal[i]);
    }

}   /* END main() */

/* REMEZ EXCHANGE ALGORITHM */
void remez(edge,nbands) double edge[];int nbands;{
int j,k,l;
int jchnge,jet,jm1,jp1;
int k1,kkk,klow,kn,knz,kup;
int loop1,luck;
int niter,nm1,nu,nut,nut1,nz,nzz;
int out1,out2,out3;
double cn,fsh,gtemp;
double delf,tmp;
double devl,dtemp,dnum,dden,err;
double comp,y1,ynz,aa,bb,ft,xt,xe;
double d(),gee();
static double a[67],p[67],q[67];

   devl = -1.0;
   nz = nfcns + 1;