SDRAngel  4.11.5
Developer docs for <a href="https://github.com/f4exb/sdrangel">SDRangel<\a>, an Open Source Qt5 / OpenGL 3.0+ SDR and signal analyzer frontend to various hardware.
kissfft.h
Go to the documentation of this file.
1 #ifndef INCLUDE_KISSFFT_H
2 #define INCLUDE_KISSFFT_H
3 
4 #include <complex>
5 #include <vector>
6 
7 /*
8 Copyright (c) 2003-2010 Mark Borgerding
9 
10 All rights reserved.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions
14 are met:
15 
16  * Redistributions of source code must retain the above copyright
17  notice, this list of conditions and the following disclaimer.
18  * Redistributions in binary form must reproduce the above copyright
19  notice, this list of conditions and the following disclaimer in the
20  documentation and/or other materials provided with the distribution.
21  * Neither the author nor the names of any contributors may be used to
22  endorse or promote products derived from this software without
23  specific prior written permission.
24 
25 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
29 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
30 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
32 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
33 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
34 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
35 THE POSSIBILITY OF SUCH DAMAGE.
36 */
37 
38 namespace kissfft_utils {
39 
40  template<typename T_scalar, typename T_complex>
41  struct traits {
42  typedef T_scalar scalar_type;
43  typedef T_complex cpx_type;
44  void fill_twiddles(std::complex<T_scalar>* dst, int nfft, bool inverse)
45  {
46  T_scalar phinc = (inverse ? 2 : -2) * acos((T_scalar)-1) / nfft;
47  for(int i = 0; i < nfft; ++i)
48  dst[i] = exp(std::complex<T_scalar>(0, i * phinc));
49  }
50 
51  void prepare(std::vector<std::complex<T_scalar> >& dst, int nfft, bool inverse, std::vector<int>& stageRadix, std::vector<int>& stageRemainder)
52  {
53  _twiddles.resize(nfft);
54  fill_twiddles(&_twiddles[0], nfft, inverse);
55  dst = _twiddles;
56 
57  //factorize
58  //start factoring out 4's, then 2's, then 3,5,7,9,...
59  int n = nfft;
60  int p = 4;
61  do {
62  while(n % p) {
63  switch(p) {
64  case 4:
65  p = 2;
66  break;
67  case 2:
68  p = 3;
69  break;
70  default:
71  p += 2;
72  break;
73  }
74  if(p * p > n)
75  p = n;// no more factors
76  }
77  n /= p;
78  stageRadix.push_back(p);
79  stageRemainder.push_back(n);
80  } while(n > 1);
81  }
82  std::vector<cpx_type> _twiddles;
83 
84  const cpx_type twiddle(int i)
85  {
86  return _twiddles[i];
87  }
88  };
89 
90 } // namespace
91 
92 template<typename T_Scalar, typename T_Complex, typename T_traits = kissfft_utils::traits<T_Scalar, T_Complex> >
93 class kissfft {
94 public:
95  typedef T_traits traits_type;
97  typedef typename traits_type::cpx_type cpx_type;
98 
100  {
101  }
102 
103  kissfft(int nfft, bool inverse, const traits_type & traits = traits_type()) :
104  _nfft(nfft), _inverse(inverse), _traits(traits)
105  {
106  _traits.prepare(_twiddles, _nfft, _inverse, _stageRadix, _stageRemainder);
107  }
108 
109  void configure(int nfft, bool inverse, const traits_type & traits = traits_type())
110  {
111  _twiddles.clear();
112  _stageRadix.clear();
113  _stageRemainder.clear();
114 
115  _nfft = nfft;
116  _inverse = inverse;
117  _traits = traits;
118  _traits.prepare(_twiddles, _nfft, _inverse, _stageRadix, _stageRemainder);
119  }
120 
121  void transform(const cpx_type* src, cpx_type* dst)
122  {
123  kf_work(0, dst, src, 1, 1);
124  }
125 
126 private:
127  void kf_work(int stage, cpx_type* Fout, const cpx_type* f, size_t fstride, size_t in_stride)
128  {
129  int p = _stageRadix[stage];
130  int m = _stageRemainder[stage];
131  cpx_type * Fout_beg = Fout;
132  cpx_type * Fout_end = Fout + p * m;
133 
134  if(m == 1) {
135  do {
136  *Fout = *f;
137  f += fstride * in_stride;
138  } while(++Fout != Fout_end);
139  } else {
140  do {
141  // recursive call:
142  // DFT of size m*p performed by doing
143  // p instances of smaller DFTs of size m,
144  // each one takes a decimated version of the input
145  kf_work(stage + 1, Fout, f, fstride * p, in_stride);
146  f += fstride * in_stride;
147  } while((Fout += m) != Fout_end);
148  }
149 
150  Fout = Fout_beg;
151 
152  // recombine the p smaller DFTs
153  switch(p) {
154  case 2:
155  kf_bfly2(Fout, fstride, m);
156  break;
157  case 3:
158  kf_bfly3(Fout, fstride, m);
159  break;
160  case 4:
161  kf_bfly4(Fout, fstride, m);
162  break;
163  case 5:
164  kf_bfly5(Fout, fstride, m);
165  break;
166  default:
167  kf_bfly_generic(Fout, fstride, m, p);
168  break;
169  }
170  }
171 
172  // these were #define macros in the original kiss_fft
173  void C_ADD(cpx_type& c, const cpx_type& a, const cpx_type& b)
174  {
175  c = a + b;
176  }
177  void C_MUL(cpx_type& c, const cpx_type& a, const cpx_type& b)
178  {
179  //c = a * b;
180  c = cpx_type(a.real() * b.real() - a.imag() * b.imag(), a.real() * b.imag() + a.imag() * b.real());
181  }
182  void C_SUB(cpx_type& c, const cpx_type& a, const cpx_type& b)
183  {
184  c = a - b;
185  }
186  void C_ADDTO(cpx_type& c, const cpx_type& a)
187  {
188  c += a;
189  }
190  void C_FIXDIV(cpx_type&, int)
191  {
192  } // NO-OP for float types
193  scalar_type S_MUL(const scalar_type& a, const scalar_type& b)
194  {
195  return a * b;
196  }
197  scalar_type HALF_OF(const scalar_type& a)
198  {
199  return a * .5;
200  }
201  void C_MULBYSCALAR(cpx_type& c, const scalar_type& a)
202  {
203  c *= a;
204  }
205 
206  void kf_bfly2(cpx_type* Fout, const size_t fstride, int m)
207  {
208  for(int k = 0; k < m; ++k) {
209  //cpx_type t = Fout[m + k] * _traits.twiddle(k * fstride);
210  cpx_type t;
211  C_MUL(t, Fout[m + k], _traits.twiddle(k * fstride));
212  Fout[m + k] = Fout[k] - t;
213  Fout[k] += t;
214  }
215  }
216 
217  void kf_bfly4(cpx_type* Fout, const size_t fstride, const size_t m)
218  {
219  cpx_type scratch[7];
220  int negative_if_inverse = _inverse * -2 + 1;
221  for(size_t k = 0; k < m; ++k) {
222  //scratch[0] = Fout[k + m] * _traits.twiddle(k * fstride);
223  C_MUL(scratch[0], Fout[k + m], _traits.twiddle(k * fstride));
224  C_MUL(scratch[1], Fout[k + 2 * m], _traits.twiddle(k * fstride * 2));
225  C_MUL(scratch[2], Fout[k + 3 * m], _traits.twiddle(k * fstride * 3));
226  scratch[5] = Fout[k] - scratch[1];
227 
228  Fout[k] += scratch[1];
229  scratch[3] = scratch[0] + scratch[2];
230  scratch[4] = scratch[0] - scratch[2];
231  scratch[4] = cpx_type(scratch[4].imag() * negative_if_inverse, -scratch[4].real() * negative_if_inverse);
232 
233  Fout[k + 2 * m] = Fout[k] - scratch[3];
234  Fout[k] += scratch[3];
235  Fout[k + m] = scratch[5] + scratch[4];
236  Fout[k + 3 * m] = scratch[5] - scratch[4];
237  }
238  }
239 
240  void kf_bfly3(cpx_type* Fout, const size_t fstride, const size_t m)
241  {
242  size_t k = m;
243  const size_t m2 = 2 * m;
244  cpx_type* tw1;
245  cpx_type* tw2;
246  cpx_type scratch[5];
247  cpx_type epi3;
248  epi3 = _twiddles[fstride * m];
249  tw1 = tw2 = &_twiddles[0];
250 
251  do {
252  C_FIXDIV(*Fout, 3);
253  C_FIXDIV(Fout[m], 3);
254  C_FIXDIV(Fout[m2], 3);
255 
256  C_MUL(scratch[1], Fout[m], *tw1);
257  C_MUL(scratch[2], Fout[m2], *tw2);
258 
259  C_ADD(scratch[3], scratch[1], scratch[2]);
260  C_SUB(scratch[0], scratch[1], scratch[2]);
261  tw1 += fstride;
262  tw2 += fstride * 2;
263 
264  Fout[m] = cpx_type(Fout->real() - HALF_OF(scratch[3].real()), Fout->imag() - HALF_OF(scratch[3].imag()));
265 
266  C_MULBYSCALAR(scratch[0], epi3.imag());
267 
268  C_ADDTO(*Fout, scratch[3]);
269 
270  Fout[m2] = cpx_type(Fout[m].real() + scratch[0].imag(), Fout[m].imag() - scratch[0].real());
271 
272  C_ADDTO(Fout[m], cpx_type(-scratch[0].imag(), scratch[0].real()));
273  ++Fout;
274  } while(--k);
275  }
276 
277  void kf_bfly5(cpx_type* Fout, const size_t fstride, const size_t m)
278  {
279  cpx_type* Fout0;
280  cpx_type* Fout1;
281  cpx_type* Fout2;
282  cpx_type* Fout3;
283  cpx_type* Fout4;
284  size_t u;
285  cpx_type scratch[13];
286  cpx_type* twiddles = &_twiddles[0];
287  cpx_type* tw;
288  cpx_type ya, yb;
289  ya = twiddles[fstride * m];
290  yb = twiddles[fstride * 2 * m];
291 
292  Fout0 = Fout;
293  Fout1 = Fout0 + m;
294  Fout2 = Fout0 + 2 * m;
295  Fout3 = Fout0 + 3 * m;
296  Fout4 = Fout0 + 4 * m;
297 
298  tw = twiddles;
299  for(u = 0; u < m; ++u) {
300  C_FIXDIV(*Fout0, 5);
301  C_FIXDIV(*Fout1, 5);
302  C_FIXDIV(*Fout2, 5);
303  C_FIXDIV(*Fout3, 5);
304  C_FIXDIV(*Fout4, 5);
305  scratch[0] = *Fout0;
306 
307  C_MUL(scratch[1], *Fout1, tw[u * fstride]);
308  C_MUL(scratch[2], *Fout2, tw[2 * u * fstride]);
309  C_MUL(scratch[3], *Fout3, tw[3 * u * fstride]);
310  C_MUL(scratch[4], *Fout4, tw[4 * u * fstride]);
311 
312  C_ADD(scratch[7], scratch[1], scratch[4]);
313  C_SUB(scratch[10], scratch[1], scratch[4]);
314  C_ADD(scratch[8], scratch[2], scratch[3]);
315  C_SUB(scratch[9], scratch[2], scratch[3]);
316 
317  C_ADDTO(*Fout0, scratch[7]);
318  C_ADDTO(*Fout0, scratch[8]);
319 
320  scratch[5] = scratch[0] + cpx_type(S_MUL(scratch[7].real(), ya.real()) + S_MUL(scratch[8].real(), yb.real()), S_MUL(scratch[7].imag(), ya.real())
321  + S_MUL(scratch[8].imag(), yb.real()));
322 
323  scratch[6] = cpx_type(S_MUL(scratch[10].imag(), ya.imag()) + S_MUL(scratch[9].imag(), yb.imag()), -S_MUL(scratch[10].real(), ya.imag()) - S_MUL(
324  scratch[9].real(), yb.imag()));
325 
326  C_SUB(*Fout1, scratch[5], scratch[6]);
327  C_ADD(*Fout4, scratch[5], scratch[6]);
328 
329  scratch[11] = scratch[0] + cpx_type(S_MUL(scratch[7].real(), yb.real()) + S_MUL(scratch[8].real(), ya.real()), S_MUL(scratch[7].imag(), yb.real())
330  + S_MUL(scratch[8].imag(), ya.real()));
331 
332  scratch[12] = cpx_type(-S_MUL(scratch[10].imag(), yb.imag()) + S_MUL(scratch[9].imag(), ya.imag()), S_MUL(scratch[10].real(), yb.imag()) - S_MUL(
333  scratch[9].real(), ya.imag()));
334 
335  C_ADD(*Fout2, scratch[11], scratch[12]);
336  C_SUB(*Fout3, scratch[11], scratch[12]);
337 
338  ++Fout0;
339  ++Fout1;
340  ++Fout2;
341  ++Fout3;
342  ++Fout4;
343  }
344  }
345 
346  /* perform the butterfly for one stage of a mixed radix FFT */
347  void kf_bfly_generic(cpx_type* Fout, const size_t fstride, int m, int p)
348  {
349  int u;
350  int k;
351  int q1;
352  int q;
353  cpx_type* twiddles = &_twiddles[0];
354  cpx_type t;
355  int Norig = _nfft;
356  cpx_type* scratchbuf = new cpx_type[p];
357 
358  for(u = 0; u < m; ++u) {
359  k = u;
360  for(q1 = 0; q1 < p; ++q1) {
361  scratchbuf[q1] = Fout[k];
362  C_FIXDIV(scratchbuf[q1], p);
363  k += m;
364  }
365 
366  k = u;
367  for(q1 = 0; q1 < p; ++q1) {
368  int twidx = 0;
369  Fout[k] = scratchbuf[0];
370  for(q = 1; q < p; ++q) {
371  twidx += fstride * k;
372  if(twidx >= Norig)
373  twidx -= Norig;
374  C_MUL(t, scratchbuf[q], twiddles[twidx]);
375  C_ADDTO(Fout[k], t);
376  }
377  k += m;
378  }
379  }
380 
381  delete[] scratchbuf;
382  }
383 
384  int _nfft;
385  bool _inverse;
386  std::vector<cpx_type> _twiddles;
387  std::vector<int> _stageRadix;
388  std::vector<int> _stageRemainder;
389  traits_type _traits;
390 };
391 #endif
void C_SUB(cpx_type &c, const cpx_type &a, const cpx_type &b)
Definition: kissfft.h:182
traits_type::cpx_type cpx_type
Definition: kissfft.h:97
void fill_twiddles(std::complex< T_scalar > *dst, int nfft, bool inverse)
Definition: kissfft.h:44
void kf_bfly3(cpx_type *Fout, const size_t fstride, const size_t m)
Definition: kissfft.h:240
scalar_type S_MUL(const scalar_type &a, const scalar_type &b)
Definition: kissfft.h:193
scalar_type HALF_OF(const scalar_type &a)
Definition: kissfft.h:197
T_scalar scalar_type
Definition: kissfft.h:42
int _nfft
Definition: kissfft.h:384
std::vector< int > _stageRadix
Definition: kissfft.h:387
void kf_bfly_generic(cpx_type *Fout, const size_t fstride, int m, int p)
Definition: kissfft.h:347
Fixed< IntType, IntBits > exp(Fixed< IntType, IntBits > const &x)
Definition: fixed.h:2289
kissfft(int nfft, bool inverse, const traits_type &traits=traits_type())
Definition: kissfft.h:103
std::vector< cpx_type > _twiddles
Definition: kissfft.h:386
T_traits traits_type
Definition: kissfft.h:95
void C_ADD(cpx_type &c, const cpx_type &a, const cpx_type &b)
Definition: kissfft.h:173
void prepare(std::vector< std::complex< T_scalar > > &dst, int nfft, bool inverse, std::vector< int > &stageRadix, std::vector< int > &stageRemainder)
Definition: kissfft.h:51
void C_MUL(cpx_type &c, const cpx_type &a, const cpx_type &b)
Definition: kissfft.h:177
void configure(int nfft, bool inverse, const traits_type &traits=traits_type())
Definition: kissfft.h:109
traits_type::scalar_type scalar_type
Definition: kissfft.h:96
std::vector< int > _stageRemainder
Definition: kissfft.h:388
T_complex cpx_type
Definition: kissfft.h:43
int32_t i
Definition: decimators.h:244
void C_FIXDIV(cpx_type &, int)
Definition: kissfft.h:190
void kf_bfly5(cpx_type *Fout, const size_t fstride, const size_t m)
Definition: kissfft.h:277
const cpx_type twiddle(int i)
Definition: kissfft.h:84
void kf_bfly4(cpx_type *Fout, const size_t fstride, const size_t m)
Definition: kissfft.h:217
void transform(const cpx_type *src, cpx_type *dst)
Definition: kissfft.h:121
kissfft()
Definition: kissfft.h:99
void C_MULBYSCALAR(cpx_type &c, const scalar_type &a)
Definition: kissfft.h:201
std::vector< cpx_type > _twiddles
Definition: kissfft.h:82
void kf_bfly2(cpx_type *Fout, const size_t fstride, int m)
Definition: kissfft.h:206
traits_type _traits
Definition: kissfft.h:389
void kf_work(int stage, cpx_type *Fout, const cpx_type *f, size_t fstride, size_t in_stride)
Definition: kissfft.h:127
bool _inverse
Definition: kissfft.h:385
void C_ADDTO(cpx_type &c, const cpx_type &a)
Definition: kissfft.h:186