切比雪夫逼近法计算FIR系数
程序来自网络，可以计算普通的FIR滤波器，也可以计算微分器和希尔伯特变换器。
感觉该代码的主流程和《数字信号处理 C语言程序集》中的相似，但书中的代码我抄写下来编译后计算结果明显是错的（很可能是我抄错，可我已经比较了两遍，放弃）。
PS：十分反感使用l（lima）作变量，书中打印的l和1根本分不清。程序中将l全部替换为lima
/* TAB P VER 068 $Id: remez.c,v 1.3 1998/06/15 09:04:37 egil Exp $ * * Remez exchange algorithm * * the core of this program is converted from * an original Fortran source. see below for details. * * adaption for the PC and addition of FFT by: * egil kvaleberg * husebybakken 14a * 0379 oslo, norway * Email: * egil@kvaleberg.no * Web: * http://www.kvaleberg.com/ * * Bugs: * converted from Fortran, and no attempt has been made of hiding the fact * the user interface is pretty dismal, the error handling horrible, but * at least you should get your moneys worth */ #include <stdio.h> #include <string.h> #include <math.h> #include <unistd.h> /* *----------------------------------------------------------------------- * MAIN PROGRAM: FIR LINEAR PHASE FILTER DESIGN PROGRAM * * AUTHORS: JAMES H. MCCLELLAN * * DEPARTMENT OF ELECTRICAL ENGINEERING AND COMPUTER SCIENCE * MASSACHUSETTS INSTITUTE OF TECHNOLOGY * CAMBRIDGE, MASS. 02139 * * THOMAS W. PARKS * DEPARTMENT OF ELECTRICAL ENGINEERING * RICE UNIVERSITY * HOUSTON, TEXAS 77001 * * LAWRENCE R. RABINER * BELL LABORATORIES * MURRAY HILL, NEW JERSEY 07974 * * INPUT: * hz-- Sampling frequency (1.0 gives unity frequency operation) * nfilt-- FILTER LENGTH * jtype-- TYPE OF FILTER * 1 = MULTIPLE PASSBAND/STOPBAND FILTER * 2 = DIFFERENTIATOR * 3 = HILBERT TRANSFORM FILTER * nbands-- NUMBER OF BANDS * lgrid-- GRID DENSITY, WILL BE SET TO 16 UNLESS * SPECIFIED OTHERWISE BY A POSITIVE CONSTANT. * * edge(2*nbands)-- BANDEDGE ARRAY, LOWER AND UPPER EDGES FOR EACH BAND * WITH A MAXIMUM OF 10 BANDS. * * fx(nbands)-- DESIRED FUNCTION ARRAY (OR DESIRED SLOPE IF A * DIFFERENTIATOR) FOR EACH BAND. * * wtx(nbands)-- WEIGHT FUNCTION ARRAY IN EACH BAND. FOR A * DIFFERENTIATOR, THE WEIGHT FUNCTION IS INVERSELY * PROPORTIONAL TO F. * * SAMPLE INPUT DATA SETUP: * 1 * 32 * 1 * 3 * 0 * 0.0 0.1 0.2 0.35 0.425 0.5 * 0.0 1.0 0.0 * 10.0 1.0 10.0 * * THIS DATA SPECIFIES A LENGTH 32 BANDPASS FILTER WITH * STOPBANDS 0 TO 0.1 AND 0.425 TO 0.5, AND PASSBAND FROM * 0.2 TO 0.35 WITH WEIGHTING OF 10 IN THE STOPBANDS AND 1 * IN THE PASSBAND. THE GRID DENSITY DEFAULTS TO 16. * THIS IS THE FILTER IN FIGURE 10. * * THE FOLLOWING INPUT DATA SPECIFIES A LENGTH 32 FULLBAND * DIFFERENTIATOR WITH SLOPE 1 AND WEIGHTING OF 1/F. * THE GRID DENSITY WILL BE SET TO 20. * 1 * 32 * 2 * 1 * 20 * 0 0.5 * 1.0 * 1.0 * *----------------------------------------------------------------------- */ /* * program maximum values * tune as required */ #ifndef MAXSIZE #define MAXSIZE 4096 /* NOTE: max 128 on a 16-bit machine */ #endif #define MAXBANDS 10 #define MINFFT 128 #define ITRMAX 25 /* maximum number of iterations */ #define PI 3.1415926 void remez(double *dev, double des[], double grid[], double edge[], double wt[], int ngrid, int nbands, int iext[], double alpha[], int nfcns); void error(void); void toomuch(void); void ouch(int niter); double eff(double freq, double *fx, int lband, int jtype); double wate(double freq, double *fx, double *wtx, int lband, int jtype); double gee(int k, int n, double *grid, double *x, double *y, double *ad); double d(int k, int n, int m, double *x); double dofft(double h[], int nfcns, int nfilt, int neg, int nodd, double hz); void plotparam(double edge[],double fx[],int nbands, double min_db, double hz); #define DIMSIZE (MAXSIZE/2 + 2) #define WRKSIZE (16 * DIMSIZE) #define MINFILT 2 /* was 3, but 2 seems to work (although pretty useless) */ char file_coef[256]; char file_fft[256]; char file_param[256]; int main() { int jtype = 1; int nbands = 0; int nfilt = 0; int lgrid = 16; int nz; int neg, nodd, nm1; int j, jb, k, kup, lima, lband; double delf, change, fup, temp; static double edge[MAXBANDS+MAXBANDS+1],fx[MAXBANDS+1], wtx[MAXBANDS+1],deviat[MAXBANDS+1]; /* * THE FOLLOWING PARAMETERS ARE PASSED TO REMEZ */ int ngrid; double dev; static double h[DIMSIZE+1], des[WRKSIZE+1], grid[WRKSIZE+1], wt[WRKSIZE+1]; static int iext[DIMSIZE+1]; int nfcns; static double alpha[DIMSIZE+1]; char buf[256]; double hz = 1.0; /* * THE PROGRAM IS SET UP FOR A MAXIMUM LENGTH OF 128, BUT * THIS UPPER LIMIT CAN BE CHANGED BY REDIMENSIONING THE * ARRAYS iext, ad, alpha, x, y, h TO BE MAXSIZE/2 + 2. * THE ARRAYS des, grid, AND wt MUST DIMENSIONED 16(MAXSIZE/2 + 2). */ /* * PROGRAM INPUT SECTION */ printf("sampling frequency (1 for unity frequency)\n"); printf("Freq: "); fgets(buf,sizeof(buf),stdin); sscanf(buf,"%lf",&hz); printf("\n"); if (hz == 0.0) hz = 1.0; printf("number of sampling points (%d..%d)\n",MINFILT,MAXSIZE); printf("NFILT: "); fgets(buf,sizeof(buf),stdin); sscanf(buf,"%d",&nfilt ); printf("\n"); if (nfilt < MINFILT) error(); if (nfilt > MAXSIZE) toomuch(); printf("filter type: 1=bandpass (etc.), 2=differentiator, 3=Hilbert\n"); printf("JTYPE: "); fgets(buf,sizeof(buf),stdin); sscanf(buf,"%d",&jtype ); printf("\n"); if (jtype <= 0 || jtype > 3) error(); if (jtype==1) printf("filter bands, e.g. 2 for simple low/highpass, 3 for simple bandpass\n"); else printf("filter bands, normally %d\n",nbands=1); printf("NBANDS: "); fgets(buf,sizeof(buf),stdin); sscanf(buf,"%d", &nbands); printf("\n"); if (nbands <= 0) error(); if (nbands > MAXBANDS) error(); /* * GRID DENSITY IS ASSUMED TO BE 16 UNLESS SPECIFIED */ lgrid = 16; printf("grid interpolation (default %d)\n",lgrid); printf("LGRID: "); fgets(buf,sizeof(buf),stdin); sscanf(buf,"%d", &lgrid ); printf("\n"); if (lgrid <= 0) error(); jb = 2 * nbands; printf("%sfrequency band edges",hz==1.0 ? "relative ":"absolute "); if (jtype==1) printf(" (low and high) for each band"); printf(", %s0.0 to %.1f%s\n",hz==1.0?"scale ":"",0.5*hz,hz==1.0?"":"Hz"); printf("ENTER EDGES: "); for(j = 1; j <= jb; j++) { printf("band %d, %s: ",(j+1)/2,(j&1)?"low":"high"); scanf("%lf",&edge[j]); edge[j] /= hz; } printf("\n"); if (jtype==1) printf("desired gain for each band, scale 0.0 to 1.0 (or more)\n"); else if (jtype==2) printf("desired slope for each band, normally 1\n"); else printf("desired gain for each band, normally 1\n"); printf("ENTER FX: "); for(j = 1; j <= nbands; j++) { printf("band %d gain: ",j); scanf("%lf",&fx[j]); } printf("\n"); if (jtype==1) printf("desired ripple weight for each band (e.g. 1 for passband, 10 for stopband)\n"); else printf("desired ripple weight for each band (normally 1)\n"); printf("ENTER WTX: "); for(j = 1; j <= nbands; j++) { printf("band %d ripple weight: ",j); scanf("%lf",&wtx[j]); } printf("\n"); { char *p; printf("\nresult file for filter coefficients: "); fgets(file_coef,sizeof(file_coef),stdin); if (file_coef[0]=='\n') /* hack, to flush previous scanf() */ fgets(file_coef,sizeof(file_coef),stdin); p = strchr(file_coef,'\n'); if (p != NULL) *p = '\0'; printf("\nplot result file for FFT result: "); fgets(file_fft,sizeof(file_fft),stdin); p = strchr(file_fft,'\n'); if (p != NULL) *p = '\0'; printf("\nplot result file for parameters: "); fgets(file_param,sizeof(file_param),stdin); p = strchr(file_param,'\n'); if (p != NULL) *p = '\0'; } printf("\n"); neg = 1; if (jtype == 1) neg = 0; nodd = nfilt % 2; nfcns = nfilt / 2; if (nodd == 1 && neg == 0) nfcns = nfcns + 1; /* * 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) { if (edge[1] < delf) grid[1] = delf; } j = 1; lima = 1; lband = 1; /* * CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT * FUNCTION ON THE GRID */ for (;;) { fup = edge[lima + 1]; do { temp = grid[j]; des[j] = eff(temp,fx,lband,jtype); wt[j] = wate(temp,fx,wtx,lband,jtype); if (++j > WRKSIZE) toomuch(); /* too many points, or too dense grid */ grid[j] = temp + delf; } while (grid[j] <= fup); grid[j-1] = fup; des[j-1] = eff(fup,fx,lband,jtype); wt[j-1] = wate(fup,fx,wtx,lband,jtype); lband++; lima += 2; if (lband > nbands) break; grid[j] = edge[lima]; } ngrid = j - 1; if (neg == nodd) { if (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) { for(j = 1; j <= ngrid; j++) { change = cos(PI*grid[j]); des[j] = des[j] / change; wt[j] = wt[j] * change; } } } else { if (nodd != 1) { for(j = 1; j <= ngrid; j++) { change = sin(PI*grid[j]); des[j] = des[j] / change; wt[j] = wt[j] * change; } } else { for(j = 1; j <= ngrid; j++) { change = sin((2*PI) * grid[j]); des[j] = des[j] / change; wt[j] = wt[j] * change; } } } /*XX*/ temp = (double)(ngrid-1) / (double)nfcns; for(j = 1; j <= nfcns; j++) { iext[j] = (int)((j-1)*temp) + 1; /* round? !! */ } iext[nfcns+1] = ngrid; nm1 = nfcns - 1; nz = nfcns + 1; /* * CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM */ remez(&dev, des, grid, edge, wt, ngrid, nbands, iext, alpha, nfcns); /* * CALCULATE THE IMPULSE RESPONSE. */ if (neg <= 0) { if (nodd != 0) { for(j = 1; j <= nm1; j++) h[j] = 0.5 * alpha[nz-j]; h[nfcns] = alpha[1]; } else { 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 { if (nodd != 0) { 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; } else { 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]; } } /*XX*/ printf("\n**********************************************************\n"); printf(" FINITE IMPULSE RESPONSE (FIR)\n"); printf(" LINEAR PHASE DIGITAL FILTER DESIGN\n"); printf(" REMEZ EXCHANGE ALGORITHM\n"); if (jtype == 1) printf(" BANDPASS FILTER\n"); if (jtype == 2) printf(" DIFFERENTIATOR\n"); if (jtype == 3) printf(" HILBERT TRANSFORMER\n"); printf(" FILTER LENGTH = %3d\n",nfilt); printf(" ***** IMPULSE RESPONSE *****\n"); for(j = 1; j <= nfcns; j++) { k = nfilt + 1 - j; if (neg == 0) printf(" h(%2d) = %13lf = h(%2d)\n",j,h[j],k); else printf(" h(%2d) = %13lf = - h(%2d)\n",j,h[j],k); } if (neg == 1 && nodd == 1) printf(" h(%2d) = %13lf\n",nz,0.0); /* * print result to file */ if (file_coef[0]) { FILE *f= fopen(file_coef,"w"); if (f != NULL) { for (j=1; j <= nfcns; ++j) fprintf(f,"%13lf\n",h[j]); if (neg == 1 && nodd == 1) fprintf(f,"%13lf\n",0.0); /* corr. by Robert.Rossmair@t-online.de */ for (j=nfcns; j >= 1; --j) fprintf(f,"%13lf\n",neg == 0 ? h[j] : -h[j]); fclose(f); } } for (k=1; k <= nbands; k += 4) { kup = k + 3; if (kup > nbands) kup = nbands; printf("\n "); for(j = k; j <= kup; j++)printf(" BAND%3d ",j); printf("\nLOWER BAND EDGE "); for(j = k; j <= kup; j++)printf("%13lf%s ", edge[2*j-1]*hz,hz==1.0 ? "":"Hz"); printf("\nUPPER BAND EDGE "); for(j = k; j <= kup; j++)printf("%13lf%s ", edge[2*j]*hz,hz==1.0 ? "":"Hz"); if (jtype != 2) printf("\nDESIRED VALUE "); else printf("\nDESIRED SLOPE "); for(j = k; j <= kup; j++)printf("%13lf ", fx[j]); printf("\nWEIGHTING "); for(j = k; j <= kup; j++)printf("%13lf ", wtx[j]); for(j = k; j <= kup; j++)deviat[j] = dev / wtx[j]; printf("\nDEVIATION "); for(j = k; j <= kup; j++)printf("%13lf ", deviat[j]); if (jtype == 1) { printf("\nDEVIATION IN dB "); for(j = k; j <= kup; j++)printf("%13lf ", 20.0 * log10(deviat[j])); } } for(j = 1; j <= nz; j++) grid[j] = grid[iext[j]]; printf("\n\nEXTREMAL FREQUENCIES (MAXIMA OF THE ERROR CURVE)\n"); for(j = 1; j <= nz; j++) printf("%13lf%s%c",grid[j]*hz,hz==1.0 ? "":"Hz",(j%5) ? ' ':'\n'); printf("\n************************************************************\n"); temp = dofft(h,nfcns,nfilt,neg,nodd,hz); plotparam(edge,fx,nbands,temp,hz); return 0; } /* *----------------------------------------------------------------------- * FUNCTION: eff * FUNCTION TO CALCULATE THE DESIRED MAGNITUDE RESPONSE * AS A FUNCTION OF FREQUENCY. * AN ARBITRARY FUNCTION OF FREQUENCY CAN BE * APPROXIMATED IF THE USER REPLACES THIS FUNCTION * WITH THE APPROPRIATE CODE TO EVALUATE THE IDEAL * MAGNITUDE. NOTE THAT THE PARAMETER FREQ IS THE * VALUE OF NORMALIZED FREQUENCY NEEDED FOR EVALUATION. *----------------------------------------------------------------------- */ double eff(double freq, double *fx, int lband, int jtype) { if (jtype != 2) return fx[lband]; else return fx[lband] * freq; } /* *----------------------------------------------------------------------- * FUNCTION: wate * FUNCTION TO CALCULATE THE WEIGHT FUNCTION AS A FUNCTION * OF FREQUENCY. SIMILAR TO THE FUNCTION eff, THIS FUNCTION CAN * BE REPLACED BY A USER-WRITTEN ROUTINE TO CALCULATE ANY * DESIRED WEIGHTING FUNCTION. *----------------------------------------------------------------------- */ double wate(double freq, double *fx, double *wtx, int lband, int jtype) { if (jtype != 2) return wtx[lband]; if (fx[lband] >= 0.0001) return wtx[lband] / freq; return wtx[lband]; } /* *----------------------------------------------------------------------- * SUBROUTINE: error * THIS ROUTINE WRITES AN ERROR MESSAGE IF AN * ERROR HAS BEEN DETECTED IN THE INPUT DATA. *----------------------------------------------------------------------- */ void error(void) { printf("Error in input data.\n"); _exit(1); } /* *----------------------------------------------------------------------- * SUBROUTINE: toomuch * THIS ROUTINE WRITES AN ERROR MESSAGE IF TOO * MANY SAMPLE POINTS OR TOO DENSE GRID IS ENTERED. *----------------------------------------------------------------------- */ void toomuch(void) { printf("Too many sample points or too dense grid.\n"); _exit(1); } /* *----------------------------------------------------------------------- * SUBROUTINE: remez * THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM * FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS * FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE * ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE * DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE * GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE * EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV * ERROR BY DETERMINING THE BEST LOCATION OF THE EXTREMAL * FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES * THE COEFFICIENTS OF THE BEST APPROXIMATION. *----------------------------------------------------------------------- */ void remez(double *dev, double des[], double grid[], double edge[], double wt[], int ngrid, int nbands, int iext[], double alpha[], int nfcns) /* dev, iext, alpha are output types */ /* des, grid, edge, wt, ngrid, nbands, nfcns are input types */ { int k, k1, kkk, kn, knz, klow, kup, nz, nzz, nm1; int cn; int j, jchnge, jet, jm1, jp1; int lima, luck, nu, nut, nut1, niter; double ynz, comp, devl, gtemp, fsh, y1, err, dtemp, delf, dnum, dden; double aa, bb, ft, xe, xt; static double a[DIMSIZE+1], p[DIMSIZE+1], q[DIMSIZE+1]; static double ad[DIMSIZE+1], x[DIMSIZE+1], y[DIMSIZE+1]; devl = -1.0; nz = nfcns+1; nzz = nfcns+2; niter = 0; do { L100: iext[nzz] = ngrid + 1; ++niter; if (niter > ITRMAX) break; printf("ITERATION %2d: ",niter); for(j = 1; j <= nz; j++) x[j] = cos(grid[iext[j]]*(2*PI)); jet = (nfcns-1) / 15 + 1; for(j = 1; j <= nz; j++) ad[j] = d(j,nz,jet,x); dnum = 0.0; dden = 0.0; k = 1; for(j = 1; j <= nz; j++) { lima = iext[j]; dnum += ad[j] * des[lima]; dden += (double)k * ad[j] / wt[lima]; k = -k; } *dev = dnum / dden; printf("DEVIATION = %lg\n",*dev); nu = 1; if ( (*dev) > 0.0 ) nu = -1; (*dev) = -(double)nu * (*dev); k = nu; for(j = 1; j <= nz; j++) { lima = iext[j]; y[j] = des[lima] + (double)k * (*dev) / wt[lima]; k = -k; } if ( (*dev) <= devl ) { /* finished */ ouch(niter); break; } devl = (*dev); jchnge = 0; k1 = iext[1]; knz = iext[nz]; klow = 0; nut = -nu; j = 1; /* * SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION */ L200: if (j == nzz) ynz = comp; if (j >= nzz) goto L300; kup = iext[j+1]; lima = iext[j]+1; nut = -nut; if (j == 2) y1 = comp; comp = (*dev); if (lima >= kup) goto L220; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L220; comp = (double)nut * err; L210: if (++lima >= kup) goto L215; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L215; comp = (double)nut * err; goto L210; L215: iext[j++] = lima - 1; klow = lima - 1; ++jchnge; goto L200; L220: --lima; L225: if (--lima <= klow) goto L250; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) > 0.0) goto L230; if (jchnge <= 0) goto L225; goto L260; L230: comp = (double)nut * err; L235: if (--lima <= klow) goto L240; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L240; comp = (double)nut * err; goto L235; L240: klow = iext[j]; iext[j] = lima+1; ++j; ++jchnge; goto L200; L250: lima = iext[j]+1; if (jchnge > 0) goto L215; L255: if (++lima >= kup) goto L260; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L255; comp = (double)nut * err; goto L210; L260: klow = iext[j++]; goto L200; L300: if (j > nzz) goto L320; if (k1 > iext[1] ) k1 = iext[1]; if (knz < iext[nz]) knz = iext[nz]; nut1 = nut; nut = -nu; lima = 0; kup = k1; comp = ynz*(1.00001); luck = 1; L310: if (++lima >= kup) goto L315; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L310; comp = (double) nut * err; j = nzz; goto L210; L315: luck = 6; goto L325; L320: if (luck > 9) goto L350; if (comp > y1) y1 = comp; k1 = iext[nzz]; L325: lima = ngrid+1; klow = knz; nut = -nut1; comp = y1*(1.00001); L330: if (--lima <= klow) goto L340; err = (gee(lima,nz,grid,x,y,ad)-des[lima]) * wt[lima]; if (((double)nut*err-comp) <= 0.0) goto L330; j = nzz; comp = (double) nut * err; luck = luck + 10; goto L235; L340: if (luck == 6) goto L370; for(j = 1; j <= nfcns; j++) { iext[nzz-j] = iext[nz-j]; } iext[1] = k1; goto L100; L350: kn = iext[nzz]; for(j = 1; j <= nfcns; j++) iext[j] = iext[j+1]; iext[nz] = kn; goto L100; L370: ; } while (jchnge > 0); /* * CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION * USING THE INVERSE DISCRETE FOURIER TRANSFORM */ nm1 = nfcns - 1; fsh = 1.0e-06; gtemp = grid[1]; x[nzz] = -2.0; cn = 2*nfcns - 1; delf = 1.0/cn; lima = 1; kkk = 0; /*XX*/ #if 0 if (grid[1] < 0.01 && grid[ngrid] > 0.49) kkk = 1; #else if (edge[1] == 0.0 && edge[2*nbands] == 0.5) kkk = 1; #endif if (nfcns <= 3) kkk = 1; if (kkk != 1) { dtemp = cos((2*PI)*grid[1]); dnum = cos((2*PI)*grid[ngrid]); aa = 2.0/(dtemp-dnum); bb = -(dtemp+dnum)/(dtemp-dnum); } for(j = 1; j <= nfcns; j++) { ft = (j - 1) * delf; xt = cos((2*PI)*ft); if (kkk != 1) { xt = (xt-bb)/aa; #if 0 /*XX* ckeck up !! */ xt1 = sqrt(1.0-xt*xt); ft = atan2(xt1,xt)/(2*PI); #else ft = acos(xt)/(2*PI); #endif } L410: xe = x[lima]; if (xt > xe) goto L420; if ((xe-xt) < fsh) goto L415; ++lima; goto L410; L415: a[j] = y[lima]; goto L425; L420: if ((xt-xe) < fsh) goto L415; grid[1] = ft; a[j] = gee(1,nz,grid,x,y,ad); L425: if (lima > 1) lima = lima-1; } grid[1] = gtemp; dden = (2*PI) / cn; for(j = 1; j <= nfcns; j++) { dtemp = 0.0; dnum = (j-1) * dden; if (nm1 >= 1) { for(k = 1; k <= nm1; k++) { dtemp += a[k+1] * cos(dnum*k); } } alpha[j] = 2.0 * dtemp + a[1]; } for(j = 2; j <= nfcns; j++) alpha[j] *= 2.0 / cn; alpha[1] /= cn; if (kkk != 1) { p[1] = 2.0*alpha[nfcns]*bb+alpha[nm1]; p[2] = 2.0*aa*alpha[nfcns]; q[1] = alpha[nfcns-2]-alpha[nfcns]; for(j = 2; j <= nm1; j++) { if (j >= nm1) { aa *= 0.5; bb *= 0.5; } p[j+1] = 0.0; for(k = 1; k <= j; k++) { a[k] = p[k]; p[k] = 2.0 * bb * a[k]; } p[2] += a[1] * 2.0 *aa; jm1 = j - 1; for(k = 1; k <= jm1; k++) p[k] += q[k] + aa * a[k+1]; jp1 = j + 1; for(k = 3; k <= jp1; k++) p[k] += aa * a[k-1]; if (j != nm1) { for(k = 1; k <= j; k++)q[k] = -a[k]; q[1] += alpha[nfcns - 1 - j]; } } for(j = 1; j <= nfcns; j++) alpha[j] = p[j]; } if (nfcns <= 3) alpha[nfcns+1] = alpha[nfcns+2] = 0.0; } /* *----------------------------------------------------------------------- * FUNCTION: d * FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION * COEFFICIENTS FOR USE IN THE FUNCTION gee. *----------------------------------------------------------------------- */ double d(int k, int n, int m, double *x) { int j, lima; double q, retval; retval = 1.0; q = x[k]; for(lima = 1; lima <= m; lima++) { for (j = lima; j <= n; j += m) { if (j != k) retval *= 2.0 * (q - x[j]); } } return 1.0 / retval; } /* *----------------------------------------------------------------------- * FUNCTION: gee * FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE * LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM *----------------------------------------------------------------------- */ double gee(int k, int n, double *grid, double *x, double *y, double *ad) { int j; double p,c,d,xf; d = 0.0; p = 0.0; xf = cos((2*PI) * grid[k]); for(j = 1; j <= n; j++) { c = ad[j] / (xf - x[j]); d += c; p += c * y[j]; } return p/d; } /* *----------------------------------------------------------------------- * SUBROUTINE: ouch * WRITES AN ERROR MESSAGE WHEN THE ALGORITHM FAILS TO * CONVERGE. THERE SEEM TO BE TWO CONDITIONS UNDER WHICH * THE ALGORITHM FAILS TO CONVERGE: (1) THE INITIAL * GUESS FOR THE EXTREMAL FREQUENCIES IS SO POOR THAT * THE EXCHANGE ITERATION CANNOT GET STARTED, OR * (2) NEAR THE TERMINATION OF A CORRECT DESIGN, * THE DEVIATION DECREASES DUE TO ROUNDING ERRORS * AND THE PROGRAM STOPS. IN THIS LATTER CASE THE * FILTER DESIGN IS PROBABLY ACCEPTABLE, BUT SHOULD * BE CHECKED BY COMPUTING A FREQUENCY RESPONSE. *----------------------------------------------------------------------- */ void ouch(int niter) { printf("Failure to converge after %d iterations.\n", niter); printf("Design maybe correct anyways, but verify with FFT.\n"); } /*************************************************************************/ /* * FFT section */ struct complex { double r; double i; }; static struct complex a[MAXSIZE+1]; void fft(int m); /* * prepare for FFT * return minimum dB for plotting */ double dofft(double h[], int nfcns, int nfilt, int neg, int nodd, double hz) { int j; int n,m; int k; double d; double db; double min_db; double fq; double deg; FILE *f; for (j = 1; j <= nfcns; ++j) { k = nfilt + 1 - j; a[j].r = h[j]; a[j].i = 0.0; if (neg == 0) a[k].r = h[j]; else a[k].r = -h[j]; a[k].i = 0.0; } if (neg == 1 && nodd == 1) { a[nfcns+1].r = 0.0; a[nfcns+1].i = 0.0; } /* select fourier size */ m = 1; for (;;) { n = 2 << (m-1); if (n >= MAXSIZE) break; if (n >= MINFFT && n >= nfilt) break; ++m; } /* pad with zeros */ for (j = nfilt+1; j <= n; ++j) { a[j].r = 0.0; a[j].i = 0.0; } #if 0 printf("\n\ndeBUG: h[n]\n"); for (j = 1; j <= n; ++j) { printf("%13lf%c",a[j].r,(j%5) ? ' ':'\n'); } #endif fft(m); /* * print result to screen */ printf("\nRESULTING FREQUENCY RESPONSE\n"); printf("%12s%s %12s %8s %8s \n", "FREQ.", hz==1.0 ? "":" ", "MAGNITUDE","","PHASE"); min_db = 0.0; for (j = 1; j <= (n/2 + 1); ++j) { /* frequency */ fq = (1.0*(j-1))/n; /* magnitude */ d = sqrt(a[j].r*a[j].r + a[j].i*a[j].i); if (d == 0.0) { db = -999.0; deg = 0.0; } else { db = 20.0 * log10(d); if (db < min_db) min_db = db; /* phase, subtracting fixed delay through filter */ deg = atan2(a[j].i,a[j].r) * 360.0 / (2*PI); deg -= (360.0 * (nfilt-1) / 2.0) * fq; while (deg <= -180.0) deg += 360.0; } printf("%12.4lf%s %12.8lf %8.3lf dB %8.3lf deg\n", fq*hz, hz==1.0 ? "":"Hz", d, db, deg); } printf("\n************************************************************\n"); /* make min_db suitable for plot scaling */ min_db = (double)(10 * ((((int)min_db) / 10) - 1)); /* * print result to file * suitable for gnuplot */ if (file_fft[0]) { f = fopen(file_fft,"w"); if (f != NULL) { for (j = 1; j <= (n/2 + 1); ++j) { /* frequency */ fq = (1.0*(j-1))/n; /* magnitude */ d = sqrt(a[j].r*a[j].r + a[j].i*a[j].i); if (d == 0.0) db = min_db; else db = 20.0 * log10(d); fprintf(f,"%12.4lf, %8.3lf\n", fq*hz, db); } fclose(f); } } return min_db; } /* * output input parameters for bandpass filter */ void plotparam(double edge[],double fx[],int nbands,double min_db,double hz) { int j; double d; double db; double fq; FILE *f; /* * print result to file * suitable for gnuplot */ if (file_param[0]) { f = fopen(file_param,"w"); if (f != NULL) { for (j = 1; j <= 2*nbands; ++j) { /* frequency */ fq = edge[j]; /* magnitude */ d = fx[(j+1)/2]; if (d == 0.0) db = min_db; else db = 20.0 * log10(d); fprintf(f,"%12.4lf, %8.3lf\n", fq*hz, db); } fclose(f); } } } /* * fast fourier, 2^m * decimation in time, radix 2, in place * n = 2**m, e.g. 4 -> 16 */ void fft(int m) { struct complex u,w,t; int n; int nv2; int i,j,k; int lima,le,le1; int ip; n = 2 << (m-1); nv2 = n/2; j = 1; /* swap input values */ for (i = 1; i < n; ++i) { if (i < j) { t.r = a[j].r; t.i = a[j].i; a[j].r = a[i].r; a[j].i = a[i].i; a[i].r = t.r; a[i].i = t.i; } k = nv2; while (k < j) { j -= k; k >>= 1; } j += k; } /* loop all m-stages */ for (lima = 1; lima <= m; ++lima) { le = 2 << (lima-1); le1 = le/2; u.r = 1.0; u.i = 0.0; w.r = cos(PI/le1); w.i = sin(PI/le1); /* inverse fft: -sin */ /* loop all butterflies within stage */ for (j = 1; j <= le1; ++j) { /* inner loop */ for (i = j; i <= n; i += le) { ip = i + le1; t.r = a[ip].r * u.r - a[ip].i * u.i; t.i = a[ip].r * u.i + a[ip].i * u.r; a[ip].r = a[i].r - t.r; a[ip].i = a[i].i - t.i; a[i].r = a[i].r + t.r; a[i].i = a[i].i + t.i; } t.r = u.r * w.r - u.i * w.i; t.i = u.r * w.i + u.i * w.r; u.r = t.r; u.i = t.i; } } /* inverse fft: multiply result by 1/n? */ }
以下是使用上述代码计算一个128阶hilbert变换器的过程
sampling frequency (1 for unity frequency) Freq: 1 number of sampling points (2..4096) NFILT: 128 filter type: 1=bandpass (etc.), 2=differentiator, 3=Hilbert JTYPE: 3 filter bands, normally 1 NBANDS: 1 grid interpolation (default 16) LGRID: 16 relative frequency band edges, scale 0.0 to 0.5 ENTER EDGES: band 1, low: 0.01 band 1, high: 0.5 desired gain for each band, normally 1 ENTER FX: band 1 gain: 1 desired ripple weight for each band (normally 1) ENTER WTX: band 1 ripple weight: 1 result file for filter coefficients: coe_hilbert plot result file for FFT result: fft_hilbert plot result file for parameters: param_hilbert ITERATION 1: DEVIATION = 2.1678e-05 ITERATION 2: DEVIATION = 0.000763858 ITERATION 3: DEVIATION = 0.00482411 ITERATION 4: DEVIATION = 0.00672199 ITERATION 5: DEVIATION = 0.00681893 ITERATION 6: DEVIATION = 0.00681904 ********************************************************** FINITE IMPULSE RESPONSE (FIR) LINEAR PHASE DIGITAL FILTER DESIGN REMEZ EXCHANGE ALGORITHM HILBERT TRANSFORMER FILTER LENGTH = 128 ***** IMPULSE RESPONSE ***** h( 1) = 0.003780 = - h(128) h( 2) = 0.000768 = - h(127) h( 3) = 0.000846 = - h(126) h( 4) = 0.000929 = - h(125) h( 5) = 0.001016 = - h(124) h( 6) = 0.001109 = - h(123) h( 7) = 0.001206 = - h(122) h( 8) = 0.001310 = - h(121) h( 9) = 0.001419 = - h(120) h(10) = 0.001536 = - h(119) h(11) = 0.001659 = - h(118) h(12) = 0.001787 = - h(117) h(13) = 0.001920 = - h(116) h(14) = 0.002062 = - h(115) h(15) = 0.002219 = - h(114) h(16) = 0.002373 = - h(113) h(17) = 0.002541 = - h(112) h(18) = 0.002717 = - h(111) h(19) = 0.002901 = - h(110) h(20) = 0.003096 = - h(109) h(21) = 0.003300 = - h(108) h(22) = 0.003515 = - h(107) h(23) = 0.003740 = - h(106) h(24) = 0.003978 = - h(105) h(25) = 0.004228 = - h(104) h(26) = 0.004492 = - h(103) h(27) = 0.004768 = - h(102) h(28) = 0.005060 = - h(101) h(29) = 0.005369 = - h(100) h(30) = 0.005693 = - h(99) h(31) = 0.006036 = - h(98) h(32) = 0.006399 = - h(97) h(33) = 0.006783 = - h(96) h(34) = 0.007191 = - h(95) h(35) = 0.007623 = - h(94) h(36) = 0.008083 = - h(93) h(37) = 0.008572 = - h(92) h(38) = 0.009095 = - h(91) h(39) = 0.009654 = - h(90) h(40) = 0.010254 = - h(89) h(41) = 0.010899 = - h(88) h(42) = 0.011596 = - h(87) h(43) = 0.012352 = - h(86) h(44) = 0.013173 = - h(85) h(45) = 0.014070 = - h(84) h(46) = 0.015056 = - h(83) h(47) = 0.016144 = - h(82) h(48) = 0.017354 = - h(81) h(49) = 0.018706 = - h(80) h(50) = 0.020233 = - h(79) h(51) = 0.021971 = - h(78) h(52) = 0.023971 = - h(77) h(53) = 0.026299 = - h(76) h(54) = 0.029052 = - h(75) h(55) = 0.032359 = - h(74) h(56) = 0.036420 = - h(73) h(57) = 0.041532 = - h(72) h(58) = 0.048181 = - h(71) h(59) = 0.057205 = - h(70) h(60) = 0.070187 = - h(69) h(61) = 0.090518 = - h(68) h(62) = 0.127019 = - h(67) h(63) = 0.212023 = - h(66) h(64) = 0.636559 = - h(65) BAND 1 LOWER BAND EDGE 0.010000 UPPER BAND EDGE 0.500000 DESIRED VALUE 1.000000 WEIGHTING 1.000000 DEVIATION 0.006819 EXTREMAL FREQUENCIES (MAXIMA OF THE ERROR CURVE) 0.010000 0.011953 0.016836 0.023184 0.030508 0.037832 0.045156 0.052480 0.060293 0.068105 0.075918 0.083730 0.091543 0.099355 0.107168 0.114980 0.122793 0.130605 0.138418 0.146230 0.154043 0.161855 0.169668 0.177480 0.185293 0.193105 0.200918 0.209219 0.217031 0.224844 0.232656 0.240469 0.248281 0.256094 0.263906 0.271719 0.279531 0.287832 0.295645 0.303457 0.311270 0.319082 0.326895 0.334707 0.342520 0.350332 0.358633 0.366445 0.374258 0.382070 0.389883 0.397695 0.405508 0.413320 0.421133 0.429434 0.437246 0.445059 0.452871 0.460684 0.468496 0.476309 0.484121 0.491934 0.500000 ************************************************************ RESULTING FREQUENCY RESPONSE FREQ. MAGNITUDE PHASE 0.0000 0.00000000 -999.000 dB 0.000 deg 0.0078 0.92539060 -0.673 dB -90.000 deg 0.0156 0.99481965 -0.045 dB -90.000 deg 0.0234 1.00679832 0.059 dB -90.000 deg 0.0312 0.99376102 -0.054 dB -90.000 deg 0.0391 1.00559361 0.048 dB -90.000 deg 0.0469 0.99493877 -0.044 dB -90.000 deg 0.0547 1.00472872 0.041 dB -90.000 deg 0.0625 0.99554824 -0.039 dB -90.000 deg 0.0703 1.00426273 0.037 dB -90.000 deg 0.0781 0.99586118 -0.036 dB -90.000 deg 0.0859 1.00406315 0.035 dB -90.000 deg 0.0938 0.99597764 -0.035 dB -90.000 deg 0.1016 1.00400849 0.035 dB -90.000 deg 0.1094 0.99598559 -0.035 dB -90.000 deg 0.1172 1.00403654 0.035 dB -90.000 deg 0.1250 0.99592936 -0.035 dB -90.000 deg 0.1328 1.00411503 0.036 dB -90.000 deg 0.1406 0.99583264 -0.036 dB -90.000 deg 0.1484 1.00422683 0.037 dB -90.000 deg 0.1562 0.99570818 -0.037 dB -90.000 deg 0.1641 1.00436228 0.038 dB -90.000 deg 0.1719 0.99556292 -0.039 dB -90.000 deg 0.1797 1.00451671 0.039 dB -90.000 deg 0.1875 0.99539952 -0.040 dB -90.000 deg 0.1953 1.00468968 0.041 dB -90.000 deg 0.2031 0.99521379 -0.042 dB -90.000 deg 0.2109 1.00488182 0.042 dB -90.000 deg 0.2188 0.99503678 -0.043 dB -90.000 deg 0.2266 1.00504343 0.044 dB -90.000 deg 0.2344 0.99487692 -0.045 dB -90.000 deg 0.2422 1.00520353 0.045 dB -90.000 deg 0.2500 0.99471555 -0.046 dB -90.000 deg 0.2578 1.00536680 0.046 dB -90.000 deg 0.2656 0.99454929 -0.047 dB -90.000 deg 0.2734 1.00553733 0.048 dB -90.000 deg 0.2812 0.99437160 -0.049 dB -90.000 deg 0.2891 1.00571784 0.050 dB -90.000 deg 0.2969 0.99421325 -0.050 dB -90.000 deg 0.3047 1.00585343 0.051 dB -90.000 deg 0.3125 0.99408105 -0.052 dB -90.000 deg 0.3203 1.00598440 0.052 dB -90.000 deg 0.3281 0.99395013 -0.053 dB -90.000 deg 0.3359 1.00611622 0.053 dB -90.000 deg 0.3438 0.99381642 -0.054 dB -90.000 deg 0.3516 1.00625342 0.054 dB -90.000 deg 0.3594 0.99367200 -0.055 dB -90.000 deg 0.3672 1.00637348 0.055 dB -90.000 deg 0.3750 0.99358347 -0.056 dB -90.000 deg 0.3828 1.00645887 0.056 dB -90.000 deg 0.3906 0.99349934 -0.057 dB -90.000 deg 0.3984 1.00654260 0.057 dB -90.000 deg 0.4062 0.99341536 -0.057 dB -90.000 deg 0.4141 1.00662753 0.057 dB -90.000 deg 0.4219 0.99332845 -0.058 dB -90.000 deg 0.4297 1.00671737 0.058 dB -90.000 deg 0.4375 0.99326665 -0.059 dB -90.000 deg 0.4453 1.00674866 0.058 dB -90.000 deg 0.4531 0.99323654 -0.059 dB -90.000 deg 0.4609 1.00677829 0.059 dB -90.000 deg 0.4688 0.99320708 -0.059 dB -90.000 deg 0.4766 1.00680785 0.059 dB -90.000 deg 0.4844 0.99317707 -0.059 dB -90.000 deg 0.4922 1.00683912 0.059 dB -90.000 deg 0.5000 0.99317854 -0.059 dB -90.000 deg ************************************************************
浙公网安备 33010602011771号