3 * Copyright 2010 Free Software Foundation, Inc.
5 * This file is part of GNU Radio
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)
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.
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.
27 #include <gri_fft_filter_ccc_sse.h>
32 #include <xmmintrin.h>
35 gri_fft_filter_ccc_sse::gri_fft_filter_ccc_sse (int decimation,
36 const std::vector<gr_complex> &taps)
37 : d_fftsize(-1), d_decimation(decimation), d_fwdfft(0), d_invfft(0)
39 d_xformed_taps = (gr_complex*)fftwf_malloc(1*sizeof(gr_complex));
43 gri_fft_filter_ccc_sse::~gri_fft_filter_ccc_sse ()
45 fftwf_free(d_xformed_taps);
52 print_vector_complex(const std::string label, const std::vector<gr_complex> &x)
55 for (unsigned i = 0; i < x.size(); i++)
56 std::cout << x[i] << " ";
63 * determines d_ntaps, d_nsamples, d_fftsize, d_xformed_taps
66 gri_fft_filter_ccc_sse::set_taps (const std::vector<gr_complex> &taps)
69 compute_sizes(taps.size());
71 d_tail.resize(tailsize());
72 for (i = 0; i < tailsize(); i++)
75 gr_complex *in = d_fwdfft->get_inbuf();
76 gr_complex *out = d_fwdfft->get_outbuf();
78 float scale = 1.0 / d_fftsize;
80 // Compute forward xform of taps.
81 // Copy taps into first ntaps slots, then pad with zeros
82 for (i = 0; i < d_ntaps; i++)
83 in[i] = taps[i] * scale;
85 for (; i < d_fftsize; i++)
88 d_fwdfft->execute(); // do the xform
90 // now copy output to d_xformed_taps
91 for (i = 0; i < d_fftsize; i++)
92 d_xformed_taps[i] = out[i];
97 // determine and set d_ntaps, d_nsamples, d_fftsize
100 gri_fft_filter_ccc_sse::compute_sizes(int ntaps)
102 int old_fftsize = d_fftsize;
104 d_fftsize = (int) (2 * pow(2.0, ceil(log(ntaps) / log(2))));
105 d_nsamples = d_fftsize - d_ntaps + 1;
108 fprintf(stderr, "gri_fft_filter_ccc_sse: ntaps = %d, fftsize = %d, nsamples = %d\n",
109 d_ntaps, d_fftsize, d_nsamples);
111 assert(d_fftsize == d_ntaps + d_nsamples -1 );
113 if (d_fftsize != old_fftsize){ // compute new plans
116 d_fwdfft = new gri_fft_complex(d_fftsize, true);
117 d_invfft = new gri_fft_complex(d_fftsize, false);
119 fftwf_free(d_xformed_taps);
120 d_xformed_taps = (gr_complex*)fftwf_malloc((d_fftsize)*sizeof(gr_complex));
125 gri_fft_filter_ccc_sse::filter (int nitems, const gr_complex *input, gr_complex *output)
129 int ninput_items = nitems * d_decimation;
131 for (int i = 0; i < ninput_items; i += d_nsamples){
133 memcpy(d_fwdfft->get_inbuf(), &input[i], d_nsamples * sizeof(gr_complex));
135 for (j = d_nsamples; j < d_fftsize; j++)
136 d_fwdfft->get_inbuf()[j] = 0;
138 d_fwdfft->execute(); // compute fwd xform
140 float *a = (float*)(d_fwdfft->get_outbuf());
141 float *b = (float*)(&d_xformed_taps[0]);
142 float *c = (float*)(d_invfft->get_inbuf());
144 __m128 x0, x1, x2, t0, t1, m;
145 m = _mm_set_ps(-1, 1, -1, 1);
146 for (j = 0; j < 2*d_fftsize; j+=4) { // filter in the freq domain
147 x0 = _mm_load_ps(&a[j]);
148 t0 = _mm_load_ps(&b[j]);
150 t1 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(3, 3, 1, 1));
151 t0 = _mm_shuffle_ps(t0, t0, _MM_SHUFFLE(2, 2, 0, 0));
152 t1 = _mm_mul_ps(t1, m);
154 x1 = _mm_mul_ps(x0, t0);
155 x2 = _mm_mul_ps(x0, t1);
157 x2 = _mm_shuffle_ps(x2, x2, _MM_SHUFFLE(2, 3, 0, 1));
158 x2 = _mm_add_ps(x1, x2);
160 _mm_store_ps(&c[j], x2);
163 d_invfft->execute(); // compute inv xform
165 // add in the overlapping tail
167 for (j = 0; j < tailsize(); j++)
168 d_invfft->get_outbuf()[j] += d_tail[j];
170 // copy nsamples to output
172 while (j < d_nsamples) {
173 *output++ = d_invfft->get_outbuf()[j];
176 dec_ctr = (j - d_nsamples);
179 memcpy(&d_tail[0], d_invfft->get_outbuf() + d_nsamples,
180 tailsize() * sizeof(gr_complex));
183 assert(dec_ctr == 0);