home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_300 / 308_01 / rmezgoto.c < prev    next >
Text File  |  1990-06-17  |  19KB  |  730 lines
  1. /* FILE RMEZGOTO.C */
  2. /* 7 October 1989 */
  3. /* Remez Exchange FIR Filter Design Program
  4. /* Translated into C by Bob Briggs
  5. /* Based on FORTRAN program by Parks McClellan as
  6. /* listed in Prentice Hall book by Rabiner & Gold
  7. /* "Theory and Application of Digital Signal Processing"
  8. /* and BASIC translation by N. Loy in Prentice Hall book
  9. /* "An Engineer's Guide to FIR Digital Filters"
  10. /* Writes file "impulse.rmz"
  11. */
  12.  
  13. #include "stdio.h"
  14. #include "math.h"
  15. #define PI  (3.141592653589793)
  16. #define PI2 (6.283185307179586)
  17. #define NFMAX (128)  /* maximum filter length */
  18. #define EPS (1e-24)  /* small number */
  19. #define ITRMAX (25)  /* maximum number of iterations */
  20. #define NBMAX  (10)  /* maximum number of bands */
  21.  
  22. int nfcns,ngrid;
  23. int iext[(NFMAX/2+2)+1];
  24. double alpha[(NFMAX/2+2)+1],des[16*(NFMAX/2+2)+1];
  25. double grid[16*(NFMAX/2+2)+1];
  26. double wt[16*(NFMAX/2+2)+1];
  27. double dev;
  28. double ad[(NFMAX/2+2)+1],x[(NFMAX/2+2)+1],y[(NFMAX/2+2)+1];
  29.  
  30. main(){
  31. void rerror(),exit(),remez();
  32. int i,j,k,l;
  33. int keyboard, example1, example2, example3, example4;
  34. int jtype,kup,lband,lgrid,nfilt,nbands,result;
  35. int jb,neg,nodd,nm1,nz;
  36. double delf,fup,temp,change;
  37. double d(),gee();
  38. double deviat[NBMAX+1],wtx[NBMAX+1];
  39. double edge[2*NBMAX+1],fx[NBMAX+1],h[(NFMAX/2+2)+1];
  40. FILE *chan;
  41.  
  42. keyboard = 0;
  43. example1 = 0;
  44. example2 = 0;
  45. example3 = 0;
  46. example4 = 0;
  47.  
  48. jtype = 0;
  49.  
  50. printf("  ...   REMEZ EXCHANGE FIR FILTER DESIGN PROGRAM ... \n");
  51. putchar('\n');
  52.  
  53. result = 0;
  54. while(!result){
  55.    printf("Select 1 to 5:   ---  then press ENTER\n");
  56.    printf("1:  EXAMPLE1 --- LOWPASS FILTER\n");
  57.    printf("2:  EXAMPLE2 --- BANDPASS FILTER\n");
  58.    printf("3:  EXAMPLE3 --- DIFFERENTIATOR\n");
  59.    printf("4:  EXAMPLE4 --- HILBERT TRANSFORMER\n");
  60.    printf("5:  KEYBOARD --- GET INPUT PARAMETERS FROM KEYBOARD\n");
  61.    putchar('\n');
  62.    result = scanf("%d",&i);
  63.  
  64.    switch(i){
  65.       case 1:  example1 = 1;
  66.                break;
  67.       case 2:  example2 = 1;
  68.                break;
  69.       case 3:  example3 = 1;
  70.                break;
  71.       case 4:  example4 = 1;
  72.                break;
  73.       case 5:  keyboard = 1;
  74.                break;
  75.       default: result = 0;
  76.    }
  77. }
  78.  
  79. if(keyboard){     /* GET DATA FROM KEYBOARD */
  80. printf("Filter types are:  1=Bandpass, 2=Differentiator, 3=Hilbert\n");
  81.  
  82. /* GET INPUT DATA FROM USER */
  83. printf("Input: # coeff, Filter type, # bands, Grid density\n");
  84. printf("Use tab, space, or return to separate each input ...\n");
  85. result = scanf("%d%d%d%d",&nfilt,&jtype,&nbands,&lgrid);
  86.  
  87. putchar('\n');
  88. putchar('\n');
  89. printf("Band and edge (corner) definition:\n");
  90. putchar('\n');
  91. printf("Edge #:  1   2   3     4   5           6\n");
  92. printf("Band #:    1        2           3\n");
  93. printf("        |        -------\n");
  94. printf("        |       *       *\n");
  95. printf("    Mag |      *         *\n");
  96. printf("        |     *           *\n");
  97. printf("        |-----             -------------\n");
  98. printf("        |_______________________________\n");
  99. printf("                 Frequency\n");
  100.  
  101. jb = 2 * nbands;
  102. printf("Now inputting edge (corner) frequencies for %d band edges\n",jb);
  103. for(i=1;i<=jb;i++){
  104.   printf("Input edge frequency for edge (corner) # %d:\n",i);
  105.   result = scanf("%lf",&edge[i]);   /* lf = long float = double */
  106. }
  107.  
  108. for(i=1;i<=nbands;i++){
  109.   printf("Input gain, weight of band # %d:",i);
  110.   result = scanf("%lf%lf",&fx[i],&wtx[i]);
  111.   putchar('\n');
  112. }
  113. }  /* end if(keyboard) */
  114.  
  115. else if(example1){      /* LOWPASS FILTER */
  116. nfilt = 24;
  117. jtype = 1;
  118. nbands = 2;
  119. lgrid = 16;
  120. edge[1] = 0;
  121. edge[2] = 0.08;
  122. edge[3] = 0.16;
  123. edge[4] = 0.5;
  124. fx[1] = 1;
  125. fx[2] = 0;
  126. wtx[1] = 1;
  127. wtx[2] = 1;
  128. }
  129.  
  130. else if(example2){         /* BANDPASS FILTER */
  131. nfilt = 32;
  132. jtype = 1;
  133. nbands = 3;
  134. lgrid = 16;
  135. edge[1] = 0;
  136. edge[2] = 0.1;
  137. edge[3] = 0.2;
  138. edge[4] = 0.35;
  139. edge[5] = 0.425;
  140. edge[6] = 0.5;
  141. fx[1] = 0;
  142. fx[2] = 1;
  143. fx[3] = 0;
  144. wtx[1] = 10;
  145. wtx[2] = 1;
  146. wtx[3] = 10;
  147. }
  148.  
  149. else if(example3){          /* DIFFERENTIATOR */
  150. nfilt = 32;
  151. jtype = 2;
  152. nbands = 1;
  153. lgrid = 16;
  154. edge[1] = 0;
  155. edge[2] = 0.5;
  156. fx[1] = 1;
  157. wtx[1] = 1;
  158. }
  159.  
  160. else if(example4){         /* HILBERT TRANSFORMER */
  161. nfilt = 20;
  162. jtype = 3;
  163. nbands = 1;
  164. lgrid = 16;
  165. edge[1] = 0.05;
  166. edge[2] = 0.5;
  167. fx[1] = 1;
  168. wtx[1] = 1;
  169. }
  170. /* END INPUT DATA SECTION */
  171.  
  172.    if(nfilt > NFMAX){
  173.       printf("Error: # coeff > %d\n...now exiting to system...\n",NFMAX);
  174.       exit(1);
  175.    }
  176.  
  177.    if(nfilt < 3) rerror("Error: # coeff < 3");
  178.    if(nbands <= 0) nbands = 1;
  179.    if(lgrid <= 0) lgrid = 16;
  180.  
  181.    printf("#coeff = %d\nType = %d\n#bands =  %d\nGrid = %d\n",
  182.       nfilt,jtype,nbands,lgrid);
  183.    for(i=1;i<=2*nbands;i++) printf("E[%d] = %.2lf\n",i,edge[i]);
  184.    for(i=1;i<=nbands;i++) printf("Gain, wt[%d] = %.2lf  %.2lf\n",i,fx[i],wtx[i]);
  185.    putchar('\n');
  186.  
  187.    if(jtype == 0) rerror("Filter type = 0 not valid\n");
  188.  
  189.    neg = 1;
  190.    if(jtype == 1) neg = 0;
  191.    nodd  = nfilt % 2;
  192.    nfcns = nfilt / 2;
  193.    if(nodd == 1 && neg == 0) nfcns++;
  194.  
  195. /* SET UP THE DENSE GRID.  THE NUMBER OF POINTS IN THE GRID
  196. /* IS (FILTER LENGTH + 1)*GRID DENSITY/2 */
  197.    grid[1] = edge[1];
  198.    delf = lgrid * nfcns;
  199.    delf = 0.5/delf;
  200.    if(neg == 0) goto lable135;
  201.    if(edge[1] < delf) grid[1] = delf;
  202. lable135:   j = 1;
  203.    l = 1;
  204.    lband = 1;
  205. lable140:  fup = edge[l+1];
  206. lable145:  temp = grid[j];
  207.  
  208. /* CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT
  209. /* FUNCTION ON THE GRID */
  210.    des[j] = fx[lband];
  211.    if(jtype == 2) des[j] *= temp;
  212.    wt[j] = wtx[lband];
  213.    if(jtype == 2 && fx[lband] >= 0.0001) wt[j] /= temp;
  214.    j++;
  215.    grid[j] = temp + delf;
  216.    if(grid[j] > fup) goto lable150;
  217.    goto lable145;
  218. lable150:   grid[j-1] = fup;
  219.    des[j-1] = fx[lband];
  220.    if(jtype == 2) des[j-1] *= fup;
  221.    wt[j-1] = wtx[lband];
  222.    if(jtype == 2 && fx[lband] >= 0.0001) wt[j-1] /= fup;
  223.    lband++;
  224.    l += 2;
  225.    if(lband > nbands) goto lable160;
  226.    grid[j] = edge[l];
  227.    goto lable140;
  228. lable160: ngrid = j-1;
  229.    if(neg != nodd) goto lable165;
  230.    if(grid[ngrid] > (0.5 - delf)) ngrid--;
  231. lable165:  ;
  232.  
  233. /* SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
  234. /* TO THE ORIGINAL PROBLEM */
  235.    if(neg > 0) goto lable180;
  236.    if(nodd == 1) goto lable200;    /* DO NOTHING */
  237.    for(j=1;j <= ngrid;j++){
  238.       change = cos(PI * grid[j]);
  239.       if(change == 0) change = EPS;
  240.       des[j] = des[j]/change;
  241.       wt[j] *= change;
  242.    }
  243.    goto lable200;
  244. lable180:   if(nodd == 1) goto lable190;
  245.       for(j=1;j <= ngrid;j++){
  246.          change = sin(PI * grid[j]);
  247.          if(change == 0) change = EPS;
  248.          des[j] = des[j]/change;
  249.          wt[j] *= change;
  250.       }
  251.    goto lable200;
  252. lable190:   for(j=1;j <= ngrid;j++){
  253.                change = sin(PI2 * grid[j]);
  254.                if(change == 0) change = EPS;
  255.                des[j] = des[j]/change;
  256.                wt[j] *= change;
  257.             }
  258.  
  259. /* INITIAL GUESS FOR THE EXTREMAL FREQUENCIES -- EQUALLY
  260. /* SPACED ALONG THE GRID */
  261. lable200:   temp = (ngrid - 1)/nfcns;
  262. for(j=1;j<=nfcns;j++) iext[j] = (j-1) * temp + 1;
  263. iext[nfcns+1] = ngrid;
  264. nm1 = nfcns - 1;
  265. nz = nfcns + 1;
  266.  
  267. /* CALL THE REMEZ EXCHANGE ALGORITHM TO DO THE APPROXIMATION PROBLEM */
  268. remez(edge,nbands);
  269.  
  270. if(neg <= 0){
  271.    if(nodd == 0){
  272.       h[1] = 0.25 * alpha[nfcns];
  273.       for(j=2;j<=nm1;j++)
  274.          h[j] = 0.25 * (alpha[nz - j] + alpha[nfcns + 2 - j]);
  275.       h[nfcns] = 0.5 * alpha[1] + 0.25 * alpha[2];
  276.    }
  277.    else{
  278.       for(j=1;j<=nm1;j++) h[j] = 0.5 * alpha[nz - j];
  279.       h[nfcns] = alpha[1];
  280.    }
  281. }
  282. else{
  283.    if(nodd == 0){
  284.       h[1] = 0.25 * alpha[nfcns];
  285.       for(j=2;j<=nm1;j++)
  286.          h[j] = 0.25 * (alpha[nz - j] - alpha[nfcns + 2 - j]);
  287.       h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[2];
  288.    }
  289.    else{
  290.       h[1] = 0.25 * alpha[nfcns];
  291.       h[2] = 0.25 * alpha[nm1];
  292.       for(j=3;j<=nm1;j++)
  293.          h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns + 3 - j]);
  294.       h[nfcns] = 0.5 * alpha[1] - 0.25 * alpha[3];
  295.       h[nz] = 0.0;
  296.    }
  297. }
  298.  
  299. /* PROGRAM OUTPUT SECTION */
  300.  
  301. putchar('\n');
  302. for(j=1;j<=70;j++) putchar('*');
  303. putchar('\n');
  304. printf("                        FINITE IMPULSE RESPONSE (FIR)\n");
  305. printf("                        LINEAR PHASE DIGITAL FILTER DESIGN\n");
  306. printf("                        REMEZ EXCHANGE ALGORITHM\n");
  307. switch(jtype){
  308.    case 1: printf("                        BANDPASS FILTER\n");
  309.            break;
  310.    case 2: printf("                        DIFFERENTIATOR\n");
  311.            break;
  312.    case 3: printf("                        HILBERT TRANSFORMER\n");
  313.            break;
  314. }
  315. printf("              FILTER LENGTH = %d\n",nfilt);
  316. printf("              ***** IMPULSE RESPONSE *****\n");
  317. for(j=1;j<=nfcns;j++){
  318.    k = nfilt + 1 - j;
  319.    if(neg == 0)
  320.       printf("                   H(%3d) = %15.9e = H(%4d)\n",j,h[j],k);
  321.    if(neg == 1)
  322.       printf("                   H(%3d) = %15.9e = -H(%4d)\n",j,h[j],k);
  323. }
  324. if(neg == 1 && nodd == 1) printf("H(%3d) = 0.0\n",nz);
  325. putchar('\n');
  326. for(k=1;k<=nbands;k+=4){
  327.    kup = k + 3;
  328.    if(kup > nbands) kup = nbands;
  329.    for(j=1;j<=23;j++) putchar(' ');
  330.    for(j=k;j<=kup;j++) printf("BAND%2d         ",j);
  331.    putchar('\n');
  332.    printf(" LOWER BAND EDGE");
  333.    for(j=k;j<=kup;j++) printf("%15.8f",edge[2*j-1]);
  334.    putchar('\n');
  335.    printf(" UPPER BAND EDGE");
  336.    for(j=k;j<=kup;j++) printf("%15.8f",edge[2*j]);
  337.    putchar('\n');
  338.    if(jtype != 2){
  339.       printf(" DESIRED VALUE  ");
  340.       for(j=k;j<=kup;j++) printf("%15.8f",fx[j]);
  341.       putchar('\n');
  342.    }
  343.    if(jtype == 2){
  344.       printf(" DESIRED SLOPE  ");
  345.       for(j=k;j<=kup;j++) printf("%15.8f",fx[j]);
  346.       putchar('\n');
  347.    }
  348.    printf(" WEIGHTING      ");
  349.    for(j=k;j<=kup;j++) printf("%15.8f",wtx[j]);
  350.    putchar('\n');
  351.    for(j=k;j<=kup;j++) deviat[j] = dev/wtx[j];
  352.    printf(" DEVIATION      ");
  353.    for(j=k;j<=kup;j++) printf("%15.8f",deviat[j]);
  354.    putchar('\n');
  355.    if(jtype != 1) continue;
  356.    for(j=k;j<=kup;j++) deviat[j] = 20.0 * log10(deviat[j]);
  357.    printf(" DEVIATION IN DB");
  358.    for(j=k;j<=kup;j++) printf("%15.8f",deviat[j]);
  359.    putchar('\n');
  360. }
  361. putchar('\n');
  362. printf(" EXTREMAL FREQUENCIES\n");
  363. putchar(' ');
  364. for(j=1;j<=nz;j++){
  365.    printf("%12.7f",grid[iext[j]]);
  366.    if(!(j%5)){ putchar('\n');putchar(' ');}  /* CR every 5 columns */
  367. }
  368. putchar('\n');
  369. for(j=1;j<=70;j++) putchar('*');
  370. putchar('\n');
  371.  
  372. /* NOW WRITE DATA TO FILE "IMPULSE.RMZ" */
  373. chan = fopen("impulse.rmz","wb");  /* opens file for write binary */
  374. if(chan == 0){printf("Cannot open file %s\n","impulse.rmz");exit();}
  375. for(j=1;j<=nfcns;j++){
  376.    fprintf(chan,"%f ",h[j]);
  377.    if(!(j%4)){fputc('\r',chan);fputc('\n',chan);}  /* CRLF every 4 values */
  378. }
  379. if(neg == 1 && nodd == 1) fprintf(chan,"0.0 ");
  380. for(j=nfcns;j>=1;j--){
  381.    if(neg == 0) fprintf(chan,"%f ",h[j]);
  382.    if(neg == 1) fprintf(chan,"%f ",-h[j]);
  383.    if(!((nfcns-j+1)%4)){fputc('\r',chan);fputc('\n',chan);}
  384.      /* CRLF every 4 values */
  385. }
  386. fclose(chan);
  387.  
  388. }   /* END main() */
  389.  
  390. /* REMEZ EXCHANGE ALGORITHM */
  391. void remez(edge,nbands) double edge[];int nbands;{
  392. int j,k,l;
  393. int jchnge,jet,jm1,jp1;
  394. int k1,kkk,klow,kn,knz,kup;
  395. int loop1,luck;
  396. int niter,nm1,nu,nut,nut1,nz,nzz;
  397. int out1,out2,out3;
  398. double cn,fsh,gtemp;
  399. double delf,tmp;
  400. double devl,dtemp,dnum,dden,err;
  401. double comp,y1,ynz,aa,bb,ft,xt,xe;
  402. double d(),gee();
  403. double a[67],p[67],q[67];
  404.  
  405.    devl = -1.0;
  406.    nz = nfcns + 1;
  407.    nzz = nfcns + 2;
  408.    niter = 0;
  409.    comp = 0;
  410.  
  411.    printf("Iteration");
  412. lable100:
  413.    iext[nzz] = ngrid + 1;
  414.    niter++;
  415.    if(niter > ITRMAX) goto lable400;
  416.    printf(" %d",niter);
  417.    if(!(niter%20)) putchar('\n');   /* newline every 20 iterations */
  418.    for(j=1;j<=nz;j++){
  419.       dtemp = cos(PI2 * grid[iext[j]]);
  420.       x[j] = dtemp;
  421.    }
  422.    jet = (int)((nfcns - 1)/15 +1);
  423.    for(j=1;j<=nz;j++) ad[j] = d(j,nz,jet);
  424.    dnum = 0.0;
  425.    dden = 0.0;
  426.    k=1;
  427.    for(j=1;j<=nz;j++){
  428.       l=iext[j];
  429.       dnum += ad[j] * des[l];
  430.       if(wt[l] == 0.0) wt[l] = EPS;
  431.       dden += k*ad[j]/wt[l];
  432.       k = -k;
  433.    }
  434.    dev = dnum/dden;
  435.    nu = 1;
  436.    if(dev > 0.0) nu = -1;
  437.    dev *= -nu;
  438.    k = nu;
  439.    for(j=1;j<=nz;j++){
  440.       l = iext[j];
  441.       if(wt[l] == 0.0) wt[l] = EPS;
  442.       dtemp = k*dev/wt[l];
  443.       y[j] = des[l] + dtemp;
  444.       k = -k;
  445.    }
  446.    if(dev > devl) goto lable150;
  447.       printf("   ********* FAILURE TO CONVERGE **********\n");
  448.       printf("Probable cause is machine rounding error\n");
  449.       printf("The impulse response may be correct\n");
  450.       printf("Check with frequency response\n");
  451.    goto lable400;
  452. lable150:   devl = dev;
  453.    jchnge = 0;
  454.    k1 = iext[1];
  455.    knz = iext[nz];
  456.    klow = 0;
  457.    nut = -nu;
  458.    j = 1;
  459.  
  460. /* SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION */
  461.  
  462. lable200:
  463.    if(j == nzz) ynz = comp;
  464.    if(j >= nzz) goto lable300;
  465.    kup = iext[j+1];
  466.    l = iext[j] + 1;
  467.    nut = - nut;
  468.    if(j == 2) y1 = comp;
  469.    comp = dev;
  470.    if(l >= kup) goto lable220;
  471.    err = (gee(l,nz) - des[l]) * wt[l];
  472.    dtemp = nut * err - comp;
  473.    if(dtemp <= 0.0) goto lable220;
  474.    comp = nut * err;
  475. lable210:
  476.    l++;
  477.    if(l >= kup) goto lable215;
  478.    err = (gee(l,nz) - des[l]) * wt[l];
  479.    dtemp = nut * err - comp;
  480.    if(dtemp <= 0.0) goto lable215;
  481.    comp = nut * err;
  482.    goto lable210;
  483. lable215:
  484.    iext[j] = l - 1;
  485.    j++;
  486.    klow = l - 1;
  487.    jchnge++;
  488.    goto lable200;
  489. lable220:
  490.    l--;
  491. lable225:
  492.    l--;
  493.    if(l <= klow) goto lable250;
  494.    err = (gee(l,nz) - des[l]) * wt[l];
  495.    dtemp = nut * err - comp;
  496.    if(dtemp > 0.0) goto lable230;
  497.    if(jchnge <= 0) goto lable225;
  498.    goto lable260;
  499. lable230:
  500.    comp = nut * err;
  501. lable235:
  502.    l--;
  503.    if(l <= klow) goto lable240;
  504.    err = (gee(l,nz) - des[l]) * wt[l];
  505.    dtemp = nut * err - comp;
  506.    if(dtemp <= 0.0) goto lable240;
  507.    comp = nut * err;
  508.    goto lable235;
  509. lable240:   klow = iext[j];
  510.    iext[j] = l+1;
  511.    j++;
  512.    jchnge++;
  513.    goto lable200;
  514. lable250:
  515.    l = iext[j] + 1;
  516.    if(jchnge > 0) goto lable215;
  517. lable255:
  518.    l++;
  519.    if(l >= kup) goto lable260;
  520.    err = (gee(l,nz) - des[l]) * wt[l];
  521.    dtemp = nut * err - comp;
  522.    if(dtemp <= 0.0) goto lable255;
  523.    comp = nut * err;
  524.    goto lable210;
  525. lable260:
  526.    klow = iext[j];
  527.    j++;
  528.    goto lable200;
  529. lable300:
  530.    if(j > nzz) goto lable320;
  531.    if(k1 > iext[1]) k1 = iext[1];
  532.    if(knz < iext[nz]) knz = iext[nz];
  533.    nut1 = nut;
  534.    nut = -nu;
  535.    l = 0;
  536.    kup = k1;
  537.    comp = ynz * (1.00001);
  538.    luck = 1;
  539. lable310:
  540.    l++;
  541.    if(l >= kup) goto lable315;
  542.    err = (gee(l,nz) - des[l]) * wt[l];
  543.    dtemp = nut * err - comp;
  544.    if(dtemp <= 0.0) goto lable310;
  545.    comp = nut * err;
  546.    j = nzz;
  547.    goto lable210;
  548. lable315:
  549.    luck = 6;
  550.    goto lable325;
  551. lable320:
  552.    if(luck > 9) goto lable350;
  553.    if(comp > y1) y1 = comp;
  554.    k1 = iext[nzz];
  555. lable325:
  556.    l = ngrid + 1;
  557.    klow = knz;
  558.    nut = -nut1;
  559.    comp = y1 * (1.00001);
  560. lable330:
  561.    l--;
  562.    if(l <= klow) goto lable340;
  563.    err = (gee(l,nz) - des[l]) * wt[l];
  564.    dtemp = nut * err - comp;
  565.    if(dtemp <= 0.0) goto lable330;
  566.    j = nzz;
  567.    comp = nut * err;
  568.    luck += 10;
  569.    goto lable235;
  570. lable340:
  571.    if(luck == 6) goto lable370;
  572.    for(j=1;j<=nfcns;j++) iext[nzz - j] = iext[nz - j];
  573.    iext[1] = k1;
  574.    goto lable100;
  575. lable350:
  576.    kn = iext[nzz];
  577.    for(j=1;j<=nfcns;j++) iext[j] = iext[j+1];
  578.    iext[nz] = kn;
  579.    goto lable100;
  580. lable370:
  581.    if(jchnge > 0) goto lable100;
  582.  
  583. /* CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION
  584. /* USING THE INVERSE DISCRETE FOURIER TRANSFORM */
  585.  
  586. lable400:
  587.    nm1 = nfcns - 1;
  588.    fsh = 1.0e-6;
  589.    gtemp = grid[1];
  590.    x[nzz] = -2.0;
  591.    cn = 2 * nfcns - 1;
  592.    delf = 1.0 / cn;
  593.    l = 1;
  594.    kkk = 0;
  595.    if(edge[1] < 0.01 && edge[2*nbands] > 0.49) kkk = 1;
  596.    if(nfcns <= 3) kkk = 1;
  597. if(kkk != 1) {
  598.    dtemp = cos(PI2 * grid[1]);
  599.    dnum = cos(PI2 * grid[ngrid]);
  600.    tmp = dtemp - dnum;
  601.    if(tmp == 0.0) tmp = EPS;
  602.    aa = 2.0/tmp;
  603.    bb = -(dtemp + dnum)/tmp;
  604. }
  605.  
  606. for(j=1;j<=nfcns;j++){
  607.    ft = (j - 1) * delf;
  608.    xt = cos(PI2 * ft);
  609.    if(kkk != 1){
  610.       xt = (xt - bb)/aa;
  611.       ft = acos(xt)/PI2;
  612.    }
  613.    out1 = 0;
  614.    out2 = 0;
  615.    while(1){         /* LOOP 9 */
  616.       xe = x[l];
  617.       if(xt > xe){
  618.          out1 = 1;
  619.          break;
  620.       }
  621.       if((xe - xt) < fsh){
  622.          out2 = 1;
  623.          break;
  624.       }
  625.       l++;
  626.    }     /* END OF while LOOP 9 */
  627.    if(out1){
  628.       if((xt - xe) < fsh) a[j] = y[l];
  629.       else{
  630.          grid[1] = ft;
  631.          a[j] = gee(1,nz);
  632.       }
  633.    }
  634.    if(out2) a[j] = y[l];
  635.    if(l > 1) l--;
  636. }  /* END for() LOOP */
  637.  
  638. grid[1] = gtemp;
  639. dden = PI2/cn;
  640.  
  641. for(j=1;j<=nfcns;j++){
  642.    dtemp = 0.0;
  643.    dnum = (j-1) * dden;
  644.    if(nm1 >= 1) for(k=1;k<=nm1;k++)
  645.       dtemp += a[k+1] * cos(dnum * k);
  646.    alpha[j] = 2 * dtemp + a[1];
  647. }  /* END for() LOOP */
  648.  
  649. for(j=2;j<=nfcns;j++) alpha[j] *= 2/cn;
  650. alpha[1] /= cn;
  651. if(kkk != 1){
  652.    p[1] = 2.0 * alpha[nfcns] * bb + alpha[nm1];
  653.    p[2] = 2.0 * alpha[nfcns] * aa;
  654.    q[1] = alpha[nfcns - 2] - alpha[nfcns];
  655.    for(j=2;j<=nm1;j++){
  656.       if(j >= nm1){
  657.          aa *= 0.5;
  658.          bb *= 0.5;
  659.       }
  660.       p[j+1] = 0.0;
  661.       for(k=1;k<=j;k++){
  662.          a[k] = p[k];
  663.          p[k] = 2.0 * bb * a[k];
  664.       }
  665.       p[2] += 2.0 * aa * a[1];
  666.       jm1 = j-1;
  667.       for(k=1;k<=jm1;k++) p[k] += q[k] + aa * a[k+1];
  668.       jp1 = j + 1;
  669.       for(k=3;k<=jp1;k++) p[k] += aa * a[k-1];
  670.       if(j != nm1){
  671.          for(k=1;k<=j;k++) q[k] = - a[k];
  672.          q[1] += alpha[nfcns - 1 - j];
  673.       }
  674.    }     /* END OF for() LOOP */
  675.  
  676.    for(j=1;j<=nfcns;j++) alpha[j] = p[j];
  677. }     /* END OF if(kkk != 1) */
  678.  
  679. if(nfcns <= 3){
  680.    alpha[nfcns + 1] = 0.0;
  681.    alpha[nfcns + 2] = 0.0;
  682. }
  683. }     /* END REMEZ EXCHANGE ALGORITHM */
  684.  
  685. /* FUNCTION TO CALCULATE THE LAGRANGE INTERPOLATION
  686. /* COEFFICIENTS FOR USE IN THE FUNCTION GEE */
  687. double d(k,n,m) int k,n,m;{
  688.    int j,el;
  689.    double de,q;
  690.  
  691.    de = 1.0;
  692.    q = x[k];
  693.  
  694.    for(el=1;el<=m;el++)
  695.       for(j=el;j<=n;j+=m)
  696.          if(j != k) de *= 2.0 * (q - x[j]);
  697.  
  698.    if(de == 0.0) de = EPS;
  699.    return(1.0/de);
  700. }
  701.  
  702. /* FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE
  703. /* LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM */
  704. double gee(k,n) int k,n;{
  705.    int j;
  706.    double c,de,p,xf;
  707.  
  708.    de = 0.0;
  709.    p = 0.0;
  710.    xf = cos(PI2*grid[k]);
  711.  
  712.    for(j=1;j<=n;j++){
  713.       c = xf - x[j];
  714.       if(c == 0.0) c = EPS;
  715.       c = ad[j]/c;
  716.       de += c;
  717.       p += c * y[j];
  718.    }
  719.    if(de == 0.0) de = EPS;
  720.    return(p/de);
  721. }
  722.  
  723. void rerror(error_text) char error_text[];{
  724.    void exit();
  725.  
  726.    printf("%s\n",error_text);
  727.    printf("...now exiting to system...\n");
  728.    exit(1);
  729. }
  730.