1dc869ef768502d759674240669fbe7c2c999047
[debian/gnuradio] / gnuradio-core / src / lib / filter / qa_gri_fir_filter_with_buffer_fsf.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_fsf.h>
29 #include <gri_fir_filter_with_buffer_fsf.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 float   i_type;
39 typedef short   o_type;
40 typedef float   tap_type;
41 typedef int     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_floats (float *buf, unsigned n)
57 {
58   for (unsigned i = 0; i < n; i++)
59     buf[i] = (float) rint (uniform () * 32767);
60 }
61
62 static o_type
63 ref_dotprod (const i_type input[], const tap_type taps[], int ntaps)
64 {
65   acc_type      sum = 0;
66   for (int i = 0; i < ntaps; i++) {
67     sum += input[i] * taps[i];
68   }
69       
70   return sum;
71 }
72
73 //
74 // Test for ntaps in [0,9], and input lengths in [0,17].
75 // This ensures that we are building the shifted taps correctly,
76 // and exercises all corner cases on input alignment and length.
77 //
78
79 void
80 qa_gri_fir_filter_with_buffer_fsf::t1 ()  
81 {
82   const int     MAX_TAPS        = 9;
83   const int     OUTPUT_LEN      = 17;
84   const int     INPUT_LEN       = MAX_TAPS + OUTPUT_LEN;
85
86   // Mem aligned buffer not really necessary, but why not?
87   i_type       *input = (i_type *)malloc16Align(INPUT_LEN * sizeof(i_type));
88   i_type       *dline = (i_type*)malloc16Align(INPUT_LEN * sizeof(i_type));
89   o_type        expected_output[OUTPUT_LEN];
90   o_type        actual_output[OUTPUT_LEN];
91   tap_type      taps[MAX_TAPS];
92
93   srandom (0);  // we want reproducibility
94   memset(dline, 0, INPUT_LEN*sizeof(i_type));
95
96   for (int n = 0; n <= MAX_TAPS; n++){
97     for (int ol = 0; ol <= OUTPUT_LEN; ol++){
98
99       // cerr << "@@@ n:ol " << n << ":" << ol << endl;
100
101       // build random test case
102       random_floats (input, INPUT_LEN);
103       random_floats (taps, MAX_TAPS);
104
105       // compute expected output values
106       memset(dline, 0, INPUT_LEN*sizeof(i_type));
107       for (int o = 0; o < ol; o++){
108         // use an actual delay line for this test
109         for(int oo = INPUT_LEN-1; oo > 0; oo--)
110           dline[oo] = dline[oo-1];
111         dline[0] = input[o];
112         expected_output[o] = ref_dotprod (dline, taps, n);
113       }
114
115       // build filter
116       vector<tap_type> f1_taps(&taps[0], &taps[n]);
117       gri_fir_filter_with_buffer_fsf *f1 = new gri_fir_filter_with_buffer_fsf(f1_taps);
118
119       // zero the output, then do the filtering
120       memset (actual_output, 0, sizeof (actual_output));
121       f1->filterN (actual_output, input, ol);
122
123       // check results
124       //
125       // we use a sloppy error margin because on the x86 architecture,
126       // our reference implementation is using 80 bit floating point
127       // arithmetic, while the SSE version is using 32 bit float point
128       // arithmetic.
129       
130       for (int o = 0; o < ol; o++){
131         CPPUNIT_ASSERT_DOUBLES_EQUAL(expected_output[o], actual_output[o],
132                                      abs (expected_output[o]) * ERR_DELTA);
133       }
134       delete f1;
135     }
136   }
137   free16Align(input);
138 }