git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/gen_interpolator_taps/objective_fct.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 /*
  24  * generate MMSE FIR interpolation table values
  25  */
  26
  27 #include <math.h>
  28 #include <assert.h>
  29 #include "simpson.h"
  30
  31 #define MU      0.5                     /* the MU for which we're computing coeffs */
  32
  33 #define Ts      (1.0)                   /* sampling period              */
  34 #define B       (1.0/(4*Ts))            /* one-sided signal bandwidth   */
  35 //#define       B       (1.0/(8./3*Ts))         /* one-sided signal bandwidth   */
  36
  37 static unsigned global_n;
  38 static double  *global_h;
  39 double          global_mu = MU;
  40 double          global_B  = B;
  41
  42 /*
  43  * This function computes the difference squared between the ideal
  44  * interpolator frequency response at frequency OMEGA and the
  45  * approximation defined by the FIR coefficients in global_h[]
  46  *
  47  * See eqn (9-7), "Digital Communication Receivers", Meyr, Moeneclaey
  48  * and Fechtel, Wiley, 1998.
  49  */
  50
  51 static double
  52 integrand (double omega)
  53 {
  54   double real_ideal;
  55   double real_approx;
  56   double real_diff;
  57   double imag_ideal;
  58   double imag_approx;
  59   double imag_diff;
  60
  61   int    i, n;
  62   int    I1;
  63   double *h;
  64
  65   real_ideal = cos (omega * Ts * global_mu);
  66   imag_ideal = sin (omega * Ts * global_mu);
  67
  68   n = global_n;
  69   h = global_h;
  70   I1 = -(n / 2);
  71
  72   real_approx = 0;
  73   imag_approx = 0;
  74
  75   for (i = 0; i < n; i++){
  76     real_approx += h[i] * cos (-omega * Ts * (i + I1));
  77     imag_approx += h[i] * sin (-omega * Ts * (i + I1));
  78   }
  79
  80   real_diff = real_ideal - real_approx;
  81   imag_diff = imag_ideal - imag_approx;
  82
  83   return real_diff * real_diff + imag_diff * imag_diff;
  84 }
  85
  86 /*
  87  * Integrate the difference squared over all frequencies of interest.
  88  */
  89 double
  90 c_fcn (double *x, int n)
  91 {
  92   assert ((n & 1) == 0);        /* assert n is even */
  93   global_n = n;
  94   global_h = x;
  95   return qsimp (integrand, -2 * M_PI * global_B, 2 * M_PI * global_B);
  96 }
  97
  98 /* this is the interface expected by the calling fortran code */
  99
 100 double
 101 objective (double x[], int *ndim)
 102 {
 103   return c_fcn (x, *ndim);
 104 }
 105
 106 static double
 107 si (double x)
 108 {
 109   if (fabs (x) < 1e-9)
 110     return 1.0;
 111
 112   return sin(x) / x;
 113 }
 114
 115 /*
 116  * starting guess for optimization
 117  */
 118 void initpt (double x[], int ndim)
 119 {
 120   int i;
 121   for (i = 0; i < ndim; i++){
 122     x[i] = si (M_PI * ((double) (i - ndim/2) + global_mu));
 123   }
 124 }