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