Adding ability for FIR filter with internal buffer to decimate.
[debian/gnuradio] / gnuradio-core / src / lib / filter / qa_gri_fir_filter_with_buffer_fff.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_fff.h>
29 #include <gri_fir_filter_with_buffer_fff.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 float   o_type;
40 typedef float   tap_type;
41 typedef float   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   return sum;
70 }
71
72 void
73 qa_gri_fir_filter_with_buffer_fff::t1 ()  
74 {
75   test_decimate(1);
76 }
77
78 void
79 qa_gri_fir_filter_with_buffer_fff::t2 ()
80 {
81   test_decimate(2);
82 }
83
84 void
85 qa_gri_fir_filter_with_buffer_fff::t3 ()
86 {
87   test_decimate(5);
88 }
89
90 //
91 // Test for ntaps in [0,9], and input lengths in [0,17].
92 // This ensures that we are building the shifted taps correctly,
93 // and exercises all corner cases on input alignment and length.
94 //
95 void
96 qa_gri_fir_filter_with_buffer_fff::test_decimate(unsigned int decimate)
97 {
98   const int     MAX_TAPS        = 9;
99   const int     OUTPUT_LEN      = 17;
100   const int     INPUT_LEN       = MAX_TAPS + OUTPUT_LEN;
101
102   // Mem aligned buffer not really necessary, but why not?
103   i_type       *input = (i_type *)malloc16Align(INPUT_LEN * sizeof(i_type));
104   i_type       *dline = (i_type*)malloc16Align(INPUT_LEN * sizeof(i_type));
105   o_type        expected_output[OUTPUT_LEN];
106   o_type        actual_output[OUTPUT_LEN];
107   tap_type      taps[MAX_TAPS];
108
109   srandom (0);  // we want reproducibility
110   memset(dline, 0, INPUT_LEN*sizeof(i_type));
111
112   for (int n = 0; n <= MAX_TAPS; n++){
113     for (int ol = 0; ol <= OUTPUT_LEN; ol++){
114
115       // cerr << "@@@ n:ol " << n << ":" << ol << endl;
116
117       // build random test case
118       random_floats (input, INPUT_LEN);
119       random_floats (taps, MAX_TAPS);
120
121       // compute expected output values
122       memset(dline, 0, INPUT_LEN*sizeof(i_type));
123       for (int o = 0; o < (int)(ol/decimate); o++){
124         // use an actual delay line for this test
125         for(int dd = 0; dd < (int)decimate; dd++) {
126           for(int oo = INPUT_LEN-1; oo > 0; oo--)
127             dline[oo] = dline[oo-1];
128           dline[0] = input[decimate*o+dd];
129         }
130         expected_output[o] = ref_dotprod (dline, taps, n);
131       }
132
133       // build filter
134       vector<tap_type> f1_taps(&taps[0], &taps[n]);
135       gri_fir_filter_with_buffer_fff *f1 = new gri_fir_filter_with_buffer_fff(f1_taps);
136
137       // zero the output, then do the filtering
138       memset (actual_output, 0, sizeof (actual_output));
139       f1->filterNdec (actual_output, input, ol/decimate, decimate);
140
141       // check results
142       //
143       // we use a sloppy error margin because on the x86 architecture,
144       // our reference implementation is using 80 bit floating point
145       // arithmetic, while the SSE version is using 32 bit float point
146       // arithmetic.
147       
148       for (int o = 0; o < (int)(ol/decimate); o++){
149         CPPUNIT_ASSERT_DOUBLES_EQUAL(expected_output[o], actual_output[o],
150                                      fabsf (expected_output[o]) * ERR_DELTA);
151       }
152       delete f1;
153     }
154   }
155   free16Align(input);
156 }