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