Imported Upstream version 3.2.2
[debian/gnuradio] / gnuradio-core / src / lib / filter / gr_fir_ccf_simd.cc
1 /* -*- c++ -*- */
2 /*
3  * Copyright 2002,2007 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 #include <gr_fir_ccf_simd.h>
27
28 #include <assert.h>
29 #include <malloc16.h>
30 #include <iostream>
31 #include <stdexcept>
32
33 using std::cerr;
34 using std::endl;
35
36 gr_fir_ccf_simd::gr_fir_ccf_simd ()
37   : gr_fir_ccf_generic ()
38 {
39   // cerr << "@@@ gr_fir_ccf_simd\n";
40
41   d_fcomplex_dotprod = 0;
42   
43   d_aligned_taps[0] = 0;
44   d_aligned_taps[1] = 0;
45   d_aligned_taps[2] = 0;
46   d_aligned_taps[3] = 0;
47 }
48
49 gr_fir_ccf_simd::gr_fir_ccf_simd (const std::vector<float> &new_taps)
50   : gr_fir_ccf_generic (new_taps)
51 {
52   // cerr << "@@@ gr_fir_ccf_simd\n";
53
54   d_fcomplex_dotprod = 0;
55   
56   d_aligned_taps[0] = 0;
57   d_aligned_taps[1] = 0;
58   d_aligned_taps[2] = 0;
59   d_aligned_taps[3] = 0;
60   set_taps (new_taps);
61 }
62
63 gr_fir_ccf_simd::~gr_fir_ccf_simd ()
64 {
65   free16Align (d_aligned_taps[0]);
66   free16Align (d_aligned_taps[1]);
67   free16Align (d_aligned_taps[2]);
68   free16Align (d_aligned_taps[3]);
69 }
70
71 void
72 gr_fir_ccf_simd::set_taps (const std::vector<float> &inew_taps)
73 {
74   gr_fir_ccf::set_taps (inew_taps);     // call superclass
75   const std::vector<float> new_taps = gr_reverse(inew_taps);
76
77   unsigned len = new_taps.size ();
78
79   // Make 4 copies of the coefficients, one for each data alignment
80   // Note use of special 16-byte-aligned version of calloc()
81   
82   for (unsigned i = 0; i < 4; i++){
83     free16Align (d_aligned_taps[i]);    // free old value
84
85     // this works because the bit representation of a IEEE floating point
86     // +zero is all zeros.  If you're using a different representation,
87     // you'll need to explictly set the result to the appropriate 0.0 value.
88     
89     d_aligned_taps[i] = (float *) calloc16Align (1 + (len + i - 1) / 4,
90                                                4 * sizeof (float));
91     if (d_aligned_taps[i] == 0){
92       // throw something...
93       cerr << "@@@ gr_fir_ccf_simd d_aligned_taps[" << i << "] == 0\n";
94     }
95
96     for (unsigned j = 0; j < len; j++)
97       d_aligned_taps[i][j+i] = new_taps[j];
98   }
99 }
100  
101 gr_complex 
102 gr_fir_ccf_simd::filter (const gr_complex input[])
103 {
104   if (ntaps () == 0)
105     return 0.0;
106
107   if (((intptr_t) input & 0x7) != 0)
108     throw std::invalid_argument("gr_complex must be 8-byte aligned");
109
110   // Round input data address down to 16 byte boundary
111   // NB: depending on the alignment of input[], memory
112   // before input[] will be accessed. The contents don't matter since 
113   // they'll be multiplied by zero coefficients. I can't conceive of any
114   // situation where this could cause a segfault since memory protection
115   // in the x86 machines is done on much larger boundaries.
116   
117   const gr_complex *ar = (gr_complex *)((unsigned long) input & ~15);
118
119   // Choose one of 4 sets of pre-shifted coefficients. al is both the
120   // index into d_aligned_taps[] and the number of 0 words padded onto
121   // that coefficients array for alignment purposes.
122
123   unsigned al = input - ar;
124
125   // call assembler routine to do the work, passing number of 2x4-float blocks.
126
127   // assert (((unsigned long) ar & 15) == 0);
128   // assert (((unsigned long) d_aligned_taps[al] & 15) == 0);
129
130   // cerr << "ar: " << ar << " d_aligned_taps[ar]: " << d_aligned_taps[al]
131   //  << " (ntaps() + al - 1)/2 + 1: " << (ntaps() + al -1) / 2 + 1 << endl;
132
133   float result[2];
134
135   // the trick here is to invert input and taps, and reuse FCC speedup
136   d_fcomplex_dotprod (d_aligned_taps[al], (float*)ar, (ntaps() + al - 1) / 2 + 1, result);
137
138   // cerr << "result = " << result[0] << " " << result[1] << endl;
139
140   return gr_complex(result[0], result[1]);
141 }