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

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 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 #ifdef HAVE_CONFIG_H
  24 #include <config.h>
  25 #endif
  26
  27 #include <gr_types.h>
  28 #include <qa_gri_fir_filter_with_buffer_scc.h>
  29 #include <gri_fir_filter_with_buffer_scc.h>
  30 #include <string.h>
  31 #include <iostream>
  32 #include <cmath>
  33 #include <cppunit/TestAssert.h>
  34 #include <random.h>
  35 #include <malloc16.h>
  36 #include <string.h>
  37
  38 typedef short           i_type;
  39 typedef gr_complex      o_type;
  40 typedef gr_complex      tap_type;
  41 typedef gr_complex      acc_type;
  42
  43 using std::vector;
  44
  45 #define ERR_DELTA       (1e-5)
  46
  47 #define NELEM(x) (sizeof (x) / sizeof (x[0]))
  48
  49 static float
  50 uniform ()
  51 {
  52   return 2.0 * ((float) random () / RANDOM_MAX - 0.5);  // uniformly (-1, 1)
  53 }
  54
  55 static void
  56 random_shorts (short *buf, unsigned n)
  57 {
  58   for (unsigned i = 0; i < n; i++)
  59     buf[i] = (short) rint (uniform () * 16384);
  60 }
  61
  62 static void
  63 random_complex (gr_complex *buf, unsigned n)
  64 {
  65   for (unsigned i = 0; i < n; i++){
  66     float re = rint (uniform () * 32767);
  67     float im = rint (uniform () * 32767);
  68     buf[i] = gr_complex (re, im);
  69   }
  70 }
  71
  72 static o_type
  73 ref_dotprod (const i_type input[], const tap_type taps[], int ntaps)
  74 {
  75   acc_type      sum = 0;
  76   for (int i = 0; i < ntaps; i++) {
  77     sum += (float)input[i] * taps[i];
  78   }
  79
  80   return sum;
  81 }
  82
  83 void
  84 qa_gri_fir_filter_with_buffer_scc::t1 ()
  85 {
  86   test_decimate(1);
  87 }
  88
  89 void
  90 qa_gri_fir_filter_with_buffer_scc::t2 ()
  91 {
  92   test_decimate(2);
  93 }
  94
  95 void
  96 qa_gri_fir_filter_with_buffer_scc::t3 ()
  97 {
  98   test_decimate(5);
  99 }
 100
 101 //
 102 // Test for ntaps in [0,9], and input lengths in [0,17].
 103 // This ensures that we are building the shifted taps correctly,
 104 // and exercises all corner cases on input alignment and length.
 105 //
 106 void
 107 qa_gri_fir_filter_with_buffer_scc::test_decimate (unsigned int decimate)
 108 {
 109   const int     MAX_TAPS        = 9;
 110   const int     OUTPUT_LEN      = 17;
 111   const int     INPUT_LEN       = MAX_TAPS + OUTPUT_LEN;
 112
 113   // Mem aligned buffer not really necessary, but why not?
 114   i_type       *input = (i_type *)malloc16Align(INPUT_LEN * sizeof(i_type));
 115   i_type       *dline = (i_type*)malloc16Align(INPUT_LEN * sizeof(i_type));
 116   o_type        expected_output[OUTPUT_LEN];
 117   o_type        actual_output[OUTPUT_LEN];
 118   tap_type      taps[MAX_TAPS];
 119
 120   srandom (0);  // we want reproducibility
 121   memset(dline, 0, INPUT_LEN*sizeof(i_type));
 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_shorts (input, INPUT_LEN);
 130       random_complex (taps, MAX_TAPS);
 131
 132       // compute expected output values
 133       memset(dline, 0, INPUT_LEN*sizeof(i_type));
 134       for (int o = 0; o < (int)(ol/decimate); o++){
 135         // use an actual delay line for this test
 136         for(int dd = 0; dd < (int)decimate; dd++) {
 137           for(int oo = INPUT_LEN-1; oo > 0; oo--)
 138             dline[oo] = dline[oo-1];
 139           dline[0] = input[decimate*o+dd];
 140         }
 141         expected_output[o] = ref_dotprod (dline, taps, n);
 142       }
 143
 144       // build filter
 145       vector<tap_type> f1_taps(&taps[0], &taps[n]);
 146       gri_fir_filter_with_buffer_scc *f1 = new gri_fir_filter_with_buffer_scc(f1_taps);
 147
 148       // zero the output, then do the filtering
 149       memset (actual_output, 0, sizeof (actual_output));
 150       f1->filterNdec (actual_output, input, ol/decimate, decimate);
 151
 152       // check results
 153       //
 154       // we use a sloppy error margin because on the x86 architecture,
 155       // our reference implementation is using 80 bit floating point
 156       // arithmetic, while the SSE version is using 32 bit float point
 157       // arithmetic.
 158
 159       for (int o = 0; o < (int)(ol/decimate); o++){
 160         CPPUNIT_ASSERT_COMPLEXES_EQUAL(expected_output[o], actual_output[o],
 161                                        abs (expected_output[o]) * ERR_DELTA);
 162       }
 163       delete f1;
 164     }
 165   }
 166   free16Align(input);
 167 }