git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/gen_interpolator_taps/gen_interpolator_taps.c

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 2, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #include <stdio.h>
  24 #include <unistd.h>
  25 #include <stdlib.h>
  26
  27 #define NSTEPS           10     // how many steps of mu are in the generated table
  28 #define MAX_NSTEPS      256
  29 #define NTAPS             8     // # of taps in the interpolator
  30 #define MAX_NTAPS       128
  31
  32 extern void initpt (double x[], int ntaps);
  33 extern double objective (double x[], int *ntaps);
  34 extern double global_mu;
  35 extern double global_B;
  36
  37 // fortran
  38 extern double prax2_ (double (fct)(double x[], int *ntaps),
  39                       double initv[], int *ntaps, double result[]);
  40
  41 static void
  42 usage (char *name)
  43 {
  44   fprintf (stderr, "usage: %s [-v] [-n <nsteps>] [-t <ntaps>] [-B <bw>]\n", name);
  45   exit (1);
  46 }
  47
  48 static void
  49 printline (double x[], int ntaps, int imu, int nsteps)
  50 {
  51   int   i;
  52
  53   printf ("  { ");
  54   for (i = 0; i < ntaps; i++){
  55     printf ("%12.5e", x[i]);
  56     if (i != ntaps - 1)
  57       printf (", ");
  58     else
  59       printf (" }, // %3d/%d\n", imu, nsteps);
  60   }
  61 }
  62
  63 int
  64 main (int argc, char **argv)
  65 {
  66   double        xx[MAX_NSTEPS+1][MAX_NTAPS];
  67   int           ntaps = NTAPS;
  68   int           nsteps = NSTEPS;
  69   int           i, j;
  70   double        result;
  71   double        step_size;
  72   int           c;
  73   int           verbose = 0;
  74
  75   global_B = 0.25;
  76
  77   while ((c = getopt (argc, argv, "n:t:B:v")) != EOF){
  78     switch (c){
  79     case 'n':
  80       nsteps = strtol (optarg, 0, 0);
  81       break;
  82
  83     case 't':
  84       ntaps = strtol (optarg, 0, 0);
  85       break;
  86
  87     case 'B':
  88       global_B = strtod (optarg, 0);
  89       break;
  90
  91     case 'v':
  92       verbose = 1;
  93       break;
  94
  95     default:
  96       usage (argv[0]);
  97       break;
  98     }
  99   }
 100
 101   if ((nsteps & 1) != 0){
 102     fprintf (stderr, "%s: nsteps must be even\n", argv[0]);
 103     exit (1);
 104   }
 105
 106   if (nsteps > MAX_NSTEPS){
 107     fprintf (stderr, "%s: nsteps must be < %d\n", argv[0], MAX_NSTEPS);
 108     exit (1);
 109   }
 110
 111   if ((ntaps & 1) != 0){
 112     fprintf (stderr, "%s: ntaps must be even\n", argv[0]);
 113     exit (1);
 114   }
 115
 116   if (nsteps > MAX_NTAPS){
 117     fprintf (stderr, "%s: ntaps must be < %d\n", argv[0], MAX_NTAPS);
 118     exit (1);
 119   }
 120
 121   if (global_B < 0 || global_B > 0.5){
 122     fprintf (stderr, "%s: bandwidth must be in the range (0, 0.5)\n", argv[0]);
 123     exit (1);
 124   }
 125
 126   step_size = 1.0/nsteps;
 127
 128   // the optimizer chokes on the two easy cases (0/N and N/N).  We do them by hand...
 129
 130   for (i = 0; i < ntaps; i++)
 131     xx[0][i] = 0.0;
 132   xx[0][ntaps/2] = 1.0;
 133
 134
 135   // compute optimal values for mu <= 0.5
 136
 137   for (j = 1; j <= nsteps/2; j++){
 138
 139     global_mu = j * step_size;          // this determines the MU for which we're computing the taps
 140
 141     // initialize X to a reasonable starting value
 142
 143     initpt (&xx[j][0], ntaps);
 144
 145     // find the value of X that minimizes the value of OBJECTIVE
 146
 147     result = prax2_ (objective, &xx[j][0], &ntaps, &xx[j][0]);
 148
 149     if (verbose){
 150       fprintf (stderr, "Mu: %10.8f\t", global_mu);
 151       fprintf (stderr, "Objective: %g\n", result);
 152     }
 153   }
 154
 155   // now compute remaining values via symmetry
 156
 157   for (j = 0; j < nsteps/2; j++){
 158     for (i = 0; i < ntaps; i++){
 159       xx[nsteps - j][i] = xx[j][ntaps-i-1];
 160     }
 161   }
 162
 163   // now print out the table
 164
 165   printf ("\
 166 /*\n\
 167  * This file was machine generated by gen_interpolator_taps.\n\
 168  * DO NOT EDIT BY HAND.\n\
 169  */\n\n");
 170
 171
 172   printf ("static const int     NTAPS     = %4d;\n", ntaps);
 173   printf ("static const int     NSTEPS    = %4d;\n", nsteps);
 174   printf ("static const double  BANDWIDTH = %g;\n\n", global_B);
 175
 176   printf ("static const float taps[NSTEPS+1][NTAPS] = {\n");
 177   printf ("  //    -4            -3            -2            -1             0             1             2             3                mu\n");
 178
 179
 180   for (i = 0; i <= nsteps; i++)
 181     printline (xx[i], ntaps, i, nsteps);
 182
 183   printf ("};\n\n");
 184
 185   return 0;
 186 }