git.gag.com Git - debian/gnuradio/blob - gnuradio-core/src/lib/filter/qa_gr_fir_fcc.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002,2010 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 3, 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  * FIXME.  This code is virtually identical to qa_gr_fir_?CC.cc
  25  *   Kludge up some kind of macro to handle the minor differences.
  26  */
  27
  28 #ifdef HAVE_CONFIG_H
  29 #include <config.h>
  30 #endif
  31
  32 #include <gr_types.h>
  33
  34 typedef float           i_type;
  35 typedef gr_complex      o_type;
  36 typedef gr_complex      tap_type;
  37 typedef gr_complex      acc_type;
  38
  39
  40 #include <qa_gr_fir_fcc.h>
  41 #include <gr_fir_fcc.h>
  42 #include <gr_fir_util.h>
  43 #include <string.h>
  44 #include <iostream>
  45 #include <cmath>
  46 #include <gr_types.h>
  47 #include <cppunit/TestAssert.h>
  48 #include <random.h>
  49 #include <string.h>
  50
  51 using std::vector;
  52
  53 #define ERR_DELTA       (1e-5)
  54
  55 #define NELEM(x) (sizeof (x) / sizeof (x[0]))
  56
  57
  58 //
  59 // typedef for something logically "pointer to constructor".
  60 // there may be a better way, please let me know...
  61 //
  62 typedef gr_fir_fcc* (*fir_maker_t)(const std::vector<tap_type> &taps);
  63
  64
  65 static float
  66 uniform ()
  67 {
  68   return 2.0 * ((float) random () / RANDOM_MAX - 0.5);  // uniformly (-1, 1)
  69 }
  70
  71 static void
  72 random_input (i_type *buf, unsigned n)
  73 {
  74   for (unsigned i = 0; i < n; i++)
  75     buf[i] = (i_type) rint (uniform () * 32767);
  76 }
  77
  78 static void
  79 random_complex (gr_complex *buf, unsigned n)
  80 {
  81   for (unsigned i = 0; i < n; i++){
  82     float re = rint (uniform () * 32767);
  83     float im = rint (uniform () * 32767);
  84     buf[i] = gr_complex (re, im);
  85   }
  86 }
  87
  88 static o_type
  89 ref_dotprod (const i_type input[], const tap_type taps[], int ntaps)
  90 {
  91   acc_type      sum = 0;
  92   for (int i = 0; i < ntaps; i++)
  93     sum += input[i] * taps[ntaps - i - 1];
  94
  95   return sum;
  96 }
  97
  98 //
  99 // Test for ntaps in [0,9], and input lengths in [0,17].
 100 // This ensures that we are building the shifted taps correctly,
 101 // and exercises all corner cases on input alignment and length.
 102 //
 103
 104 static void
 105 test_random_io (fir_maker_t maker)
 106 {
 107   const int     MAX_TAPS        = 9;
 108   const int     OUTPUT_LEN      = 17;
 109   const int     INPUT_LEN       = MAX_TAPS + OUTPUT_LEN;
 110
 111   i_type        input_raw[INPUT_LEN + 4*16/sizeof(i_type)];
 112   CPPUNIT_ASSERT(((intptr_t) input_raw & 0x3) == 0);
 113   memset(input_raw, 0, sizeof(input_raw));
 114   i_type       *input = &input_raw[2*16/sizeof(i_type)];
 115
 116   o_type        expected_output[OUTPUT_LEN];
 117   o_type        actual_output[OUTPUT_LEN];
 118   tap_type      taps[MAX_TAPS];
 119
 120
 121   srandom (0);  // we want reproducibility
 122
 123   for (int n = 0; n <= MAX_TAPS; n++){
 124     for (int ol = 0; ol <= OUTPUT_LEN; ol++){
 125
 126       // cerr << "@@@ n:ol " << n << ":" << ol << endl;
 127
 128       // build random test case
 129       random_input (input, INPUT_LEN);
 130       random_complex (taps, MAX_TAPS);
 131
 132       // compute expected output values
 133       for (int o = 0; o < ol; o++){
 134         expected_output[o] = ref_dotprod (&input[o], taps, n);
 135       }
 136
 137       // build filter
 138
 139       vector<tap_type> f1_taps (&taps[0], &taps[n]);
 140       gr_fir_fcc *f1 = maker (f1_taps);
 141
 142       // zero the output, then do the filtering
 143       memset (actual_output, 0, sizeof (actual_output));
 144       f1->filterN (actual_output, input, ol);
 145
 146       // check results
 147       //
 148       // we use a sloppy error margin because on the x86 architecture,
 149       // our reference implementation is using 80 bit floating point
 150       // arithmetic, while the SSE version is using 32 bit float point
 151       // arithmetic.
 152
 153       for (int o = 0; o < ol; o++){
 154         CPPUNIT_ASSERT_COMPLEXES_EQUAL(expected_output[o],
 155                                        actual_output[o],
 156                                        abs (expected_output[o]) * ERR_DELTA);
 157       }
 158
 159       delete f1;
 160     }
 161   }
 162 }
 163
 164 static void
 165 for_each (void (*f)(fir_maker_t))
 166 {
 167   std::vector<gr_fir_fcc_info>          info;
 168   gr_fir_util::get_gr_fir_fcc_info (&info);     // get all known fcc implementations
 169
 170   for (std::vector<gr_fir_fcc_info>::iterator p = info.begin ();
 171        p != info.end ();
 172        ++p){
 173
 174     std::cerr << " [" << p->name << "]";
 175     f (p->create);
 176   }
 177
 178   std::cerr << std::endl;
 179 }
 180
 181 void
 182 qa_gr_fir_fcc::t1 ()
 183 {
 184   for_each (test_random_io);
 185 }