home *** CD-ROM | disk | FTP | other *** search
/ Enigma Amiga Life 109 / EnigmaAmiga109CD.iso / software / musica / resample / src / rtest.c < prev   
C/C++ Source or Header  |  1999-10-31  |  3KB  |  98 lines
  1. /*************************************************************************\
  2.    "rtest.c": test program for the Remez Exchange algorithm ("remez.c"),
  3.                       creates a low pass FIR filter
  4.   With option '-p', a gnuplot file of the frequency response is generated.
  5. \*************************************************************************/
  6.  
  7. #include <math.h>
  8. #include <stdio.h>
  9. #include "remez.h"
  10.  
  11. float fpass, fstop;     /* describe a low pass filter */
  12. float stopwt;           /* extra weight for the stop band */
  13. extern double grid[];
  14. extern int iext[];      /* extremal frquencies */
  15.  
  16. int main( int argc, char* argv[] )
  17.     {
  18.     int n, j, plot_it;
  19.     float fr, pi, sum, err;
  20.     char filename[] = "filter.dat";
  21.     FILE *outfile;
  22.  
  23.     plot_it = (argc>1 && *argv[1] == '-' && *(argv[1]+1) == 'p');
  24.  
  25.     /* input and setup section */
  26.     printf("filter length (max. %d): ",NFMAX); scanf("%d", &n);
  27.     printf("pass and stop frequency (0 < fp < fs < 0.5): ");
  28.     scanf("%f %f", &fpass, &fstop);
  29.     printf("extra weight for errors in the stop band (>= 1): ");
  30.     scanf("%f", &stopwt);
  31.  
  32.     remez_talk = 1;
  33.     if (makefilter(n) > REMEZ_UNSURE)
  34.         return 10;
  35.  
  36.     /* output section */
  37.     printf("FIR filter designed by remez exchange algorithm.\n");
  38.     printf("filter length = %d\n", n);
  39.     for (j = 0; j<n; j++)
  40.         printf("%11.8f # h(%2d)\n", h[j], j+1);
  41.     if (plot_it)
  42.         outfile = fopen(filename, "w");
  43.     else
  44.         outfile = NULL;
  45.     if (outfile != NULL)
  46.         {
  47.         printf("writing frequency response to file '%s'\n", filename);
  48.         pi = 4 * atan(1);
  49.         for (fr = 0; fr < 0.5; fr += 0.002)
  50.             {
  51.             sum = 0;
  52.             for (j = 0; j<n; j++)
  53.                 sum += h[j] * cos(pi*fr * (n-1-2*j));
  54.             fprintf(outfile, "%f %e \n", fr, fabs(sum));
  55.             }
  56.         fclose(outfile);
  57.         }
  58.     printf("extremal frequencies:\n");
  59.     for (j = 0; 2*j < n; j++)
  60.         {
  61.         fr = grid[iext[j]];
  62.         sum = extreme(fr);
  63.         err = sum - desire(fr);
  64.         printf("%.3f: %8.5f, error: %8.5f (%.0f dB)\n",
  65.              fr, sum, err, 20*log10(fabs(err)));
  66.         }
  67.     }
  68.  
  69. float desire(float freq)
  70. /* Calculates the desired magnitude response as a function of frequency,
  71.  * in this implementation that of a low pass.
  72.  * Note that the parameter is the normalized frequency, where a value of
  73.  * 1.0 denotes the sampling rate (i. e. only values 0..0.5 will occur).
  74.  */
  75.     {
  76.     if (freq > fstop)
  77.         return 0;
  78.     else
  79.         return 1;
  80.     }
  81.  
  82. float weight(float freq)
  83. /* Calculates the weight function, i. e. how strict are we to follow the
  84.  * magnitude response specified by eff()? It is a good idea to put more
  85.  * weight on stop bands than on pass bands, and there had really better
  86.  * be transition bands (with a weight of zero) between them.
  87.  */
  88.     {
  89.     if (freq > fstop)
  90.         return stopwt;
  91.     else if (freq > fpass)
  92.         return 0;
  93.     else
  94.         return 1;
  95.     }
  96.  
  97.  
  98.