Krotos Modules 3
Loading...
Searching...
No Matches
AudioFFT.cpp
Go to the documentation of this file.
1// ==================================================================================
2// Copyright (c) 2017 HiFi-LoFi
3//
4// Permission is hereby granted, free of charge, to any person obtaining a copy
5// of this software and associated documentation files (the "Software"), to deal
6// in the Software without restriction, including without limitation the rights
7// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
8// copies of the Software, and to permit persons to whom the Software is furnished
9// to do so, subject to the following conditions:
10//
11// The above copyright notice and this permission notice shall be included in
12// all copies or substantial portions of the Software.
13//
14// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
16// FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
17// COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
18// IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
19// WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
20// ==================================================================================
21
22#include "AudioFFT.h"
23
24#include <cassert>
25#include <cmath>
26#include <cstring>
27
28#if defined(AUDIOFFT_APPLE_ACCELERATE)
29#define AUDIOFFT_APPLE_ACCELERATE_USED
30#include <Accelerate/Accelerate.h>
31#include <vector>
32#elif defined(AUDIOFFT_FFTW3)
33#define AUDIOFFT_FFTW3_USED
34#include <fftw3.h>
35#else
36#if !defined(AUDIOFFT_OOURA)
37#define AUDIOFFT_OOURA
38#endif
39#define AUDIOFFT_OOURA_USED
40#include <vector>
41#endif
42
43namespace audiofft
44{
45
46 namespace detail
47 {
48
49 class AudioFFTImpl
50 {
51 public:
52 AudioFFTImpl() = default;
53 AudioFFTImpl(const AudioFFTImpl&) = delete;
54 AudioFFTImpl& operator=(const AudioFFTImpl&) = delete;
55 virtual ~AudioFFTImpl() = default;
56 virtual void init(size_t size) = 0;
57 virtual void fft(const float* data, float* re, float* im) = 0;
58 virtual void ifft(float* data, const float* re, const float* im) = 0;
59 };
60
61 constexpr bool IsPowerOf2(size_t val) { return (val == 1 || (val & (val - 1)) == 0); }
62
63 template <typename TypeDest, typename TypeSrc>
64 void ConvertBuffer(TypeDest* dest, const TypeSrc* src, size_t len)
65 {
66 for (size_t i = 0; i < len; ++i)
67 {
68 dest[i] = static_cast<TypeDest>(src[i]);
69 }
70 }
71
72 template <typename TypeDest, typename TypeSrc, typename TypeFactor>
73 void ScaleBuffer(TypeDest* dest, const TypeSrc* src, const TypeFactor factor, size_t len)
74 {
75 for (size_t i = 0; i < len; ++i)
76 {
77 dest[i] = static_cast<TypeDest>(static_cast<TypeFactor>(src[i]) * factor);
78 }
79 }
80
81 } // End of namespace detail
82
83 // ================================================================
84
85#ifdef AUDIOFFT_OOURA_USED
86
92 class OouraFFT : public detail::AudioFFTImpl
93 {
94 public:
95 OouraFFT() : detail::AudioFFTImpl(), _size(0), _ip(), _w(), _buffer() {}
96
97 OouraFFT(const OouraFFT&) = delete;
98 OouraFFT& operator=(const OouraFFT&) = delete;
99
100 virtual void init(size_t size) override
101 {
102 if (_size != size)
103 {
104 _ip.resize(2 + static_cast<size_t>(std::sqrt(static_cast<double>(size))));
105 _w.resize(size / 2);
106 _buffer.resize(size);
107 _size = size;
108
109 const int size4 = static_cast<int>(_size) / 4;
110 makewt(size4, _ip.data(), _w.data());
111 makect(size4, _ip.data(), _w.data() + size4);
112 }
113 }
114
115 virtual void fft(const float* data, float* re, float* im) override
116 {
117 // Convert into the format as required by the Ooura FFT
118 detail::ConvertBuffer(_buffer.data(), data, _size);
119
120 rdft(static_cast<int>(_size), +1, _buffer.data(), _ip.data(), _w.data());
121
122 // Convert back to split-complex
123 {
124 double* b = _buffer.data();
125 double* bEnd = b + _size;
126 float* r = re;
127 float* i = im;
128 while (b != bEnd)
129 {
130 *(r++) = static_cast<float>(*(b++));
131 *(i++) = static_cast<float>(-(*(b++)));
132 }
133 }
134 const size_t size2 = _size / 2;
135 re[size2] = -im[0];
136 im[0] = 0.0;
137 im[size2] = 0.0;
138 }
139
140 virtual void ifft(float* data, const float* re, const float* im) override
141 {
142 // Convert into the format as required by the Ooura FFT
143 {
144 double* b = _buffer.data();
145 double* bEnd = b + _size;
146 const float* r = re;
147 const float* i = im;
148 while (b != bEnd)
149 {
150 *(b++) = static_cast<double>(*(r++));
151 *(b++) = -static_cast<double>(*(i++));
152 }
153 _buffer[1] = re[_size / 2];
154 }
155
156 rdft(static_cast<int>(_size), -1, _buffer.data(), _ip.data(), _w.data());
157
158 // Convert back to split-complex
159 detail::ScaleBuffer(data, _buffer.data(), 2.0 / static_cast<double>(_size), _size);
160 }
161
162 private:
163 size_t _size;
164 std::vector<int> _ip;
165 std::vector<double> _w;
166 std::vector<double> _buffer;
167
168 void rdft(int n, int isgn, double* a, int* ip, double* w)
169 {
170 int nw = ip[0];
171 int nc = ip[1];
172
173 if (isgn >= 0)
174 {
175 if (n > 4)
176 {
177 bitrv2(n, ip + 2, a);
178 cftfsub(n, a, w);
179 rftfsub(n, a, nc, w + nw);
180 }
181 else if (n == 4)
182 {
183 cftfsub(n, a, w);
184 }
185 double xi = a[0] - a[1];
186 a[0] += a[1];
187 a[1] = xi;
188 }
189 else
190 {
191 a[1] = 0.5 * (a[0] - a[1]);
192 a[0] -= a[1];
193 if (n > 4)
194 {
195 rftbsub(n, a, nc, w + nw);
196 bitrv2(n, ip + 2, a);
197 cftbsub(n, a, w);
198 }
199 else if (n == 4)
200 {
201 cftfsub(n, a, w);
202 }
203 }
204 }
205
206 /* -------- initializing routines -------- */
207
208 void makewt(int nw, int* ip, double* w)
209 {
210 int j, nwh;
211 double delta, x, y;
212
213 ip[0] = nw;
214 ip[1] = 1;
215 if (nw > 2)
216 {
217 nwh = nw >> 1;
218 delta = atan(1.0) / nwh;
219 w[0] = 1;
220 w[1] = 0;
221 w[nwh] = cos(delta * nwh);
222 w[nwh + 1] = w[nwh];
223 if (nwh > 2)
224 {
225 for (j = 2; j < nwh; j += 2)
226 {
227 x = cos(delta * j);
228 y = sin(delta * j);
229 w[j] = x;
230 w[j + 1] = y;
231 w[nw - j] = y;
232 w[nw - j + 1] = x;
233 }
234 bitrv2(nw, ip + 2, w);
235 }
236 }
237 }
238
239 void makect(int nc, int* ip, double* c)
240 {
241 int j, nch;
242 double delta;
243
244 ip[1] = nc;
245 if (nc > 1)
246 {
247 nch = nc >> 1;
248 delta = atan(1.0) / nch;
249 c[0] = cos(delta * nch);
250 c[nch] = 0.5 * c[0];
251 for (j = 1; j < nch; j++)
252 {
253 c[j] = 0.5 * cos(delta * j);
254 c[nc - j] = 0.5 * sin(delta * j);
255 }
256 }
257 }
258
259 /* -------- child routines -------- */
260
261 void bitrv2(int n, int* ip, double* a)
262 {
263 int j, j1, k, k1, l, m, m2;
264 double xr, xi, yr, yi;
265
266 ip[0] = 0;
267 l = n;
268 m = 1;
269 while ((m << 3) < l)
270 {
271 l >>= 1;
272 for (j = 0; j < m; j++)
273 {
274 ip[m + j] = ip[j] + l;
275 }
276 m <<= 1;
277 }
278 m2 = 2 * m;
279 if ((m << 3) == l)
280 {
281 for (k = 0; k < m; k++)
282 {
283 for (j = 0; j < k; j++)
284 {
285 j1 = 2 * j + ip[k];
286 k1 = 2 * k + ip[j];
287 xr = a[j1];
288 xi = a[j1 + 1];
289 yr = a[k1];
290 yi = a[k1 + 1];
291 a[j1] = yr;
292 a[j1 + 1] = yi;
293 a[k1] = xr;
294 a[k1 + 1] = xi;
295 j1 += m2;
296 k1 += 2 * m2;
297 xr = a[j1];
298 xi = a[j1 + 1];
299 yr = a[k1];
300 yi = a[k1 + 1];
301 a[j1] = yr;
302 a[j1 + 1] = yi;
303 a[k1] = xr;
304 a[k1 + 1] = xi;
305 j1 += m2;
306 k1 -= m2;
307 xr = a[j1];
308 xi = a[j1 + 1];
309 yr = a[k1];
310 yi = a[k1 + 1];
311 a[j1] = yr;
312 a[j1 + 1] = yi;
313 a[k1] = xr;
314 a[k1 + 1] = xi;
315 j1 += m2;
316 k1 += 2 * m2;
317 xr = a[j1];
318 xi = a[j1 + 1];
319 yr = a[k1];
320 yi = a[k1 + 1];
321 a[j1] = yr;
322 a[j1 + 1] = yi;
323 a[k1] = xr;
324 a[k1 + 1] = xi;
325 }
326 j1 = 2 * k + m2 + ip[k];
327 k1 = j1 + m2;
328 xr = a[j1];
329 xi = a[j1 + 1];
330 yr = a[k1];
331 yi = a[k1 + 1];
332 a[j1] = yr;
333 a[j1 + 1] = yi;
334 a[k1] = xr;
335 a[k1 + 1] = xi;
336 }
337 }
338 else
339 {
340 for (k = 1; k < m; k++)
341 {
342 for (j = 0; j < k; j++)
343 {
344 j1 = 2 * j + ip[k];
345 k1 = 2 * k + ip[j];
346 xr = a[j1];
347 xi = a[j1 + 1];
348 yr = a[k1];
349 yi = a[k1 + 1];
350 a[j1] = yr;
351 a[j1 + 1] = yi;
352 a[k1] = xr;
353 a[k1 + 1] = xi;
354 j1 += m2;
355 k1 += m2;
356 xr = a[j1];
357 xi = a[j1 + 1];
358 yr = a[k1];
359 yi = a[k1 + 1];
360 a[j1] = yr;
361 a[j1 + 1] = yi;
362 a[k1] = xr;
363 a[k1 + 1] = xi;
364 }
365 }
366 }
367 }
368
369 void cftfsub(int n, double* a, double* w)
370 {
371 int j, j1, j2, j3, l;
372 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
373
374 l = 2;
375 if (n > 8)
376 {
377 cft1st(n, a, w);
378 l = 8;
379 while ((l << 2) < n)
380 {
381 cftmdl(n, l, a, w);
382 l <<= 2;
383 }
384 }
385 if ((l << 2) == n)
386 {
387 for (j = 0; j < l; j += 2)
388 {
389 j1 = j + l;
390 j2 = j1 + l;
391 j3 = j2 + l;
392 x0r = a[j] + a[j1];
393 x0i = a[j + 1] + a[j1 + 1];
394 x1r = a[j] - a[j1];
395 x1i = a[j + 1] - a[j1 + 1];
396 x2r = a[j2] + a[j3];
397 x2i = a[j2 + 1] + a[j3 + 1];
398 x3r = a[j2] - a[j3];
399 x3i = a[j2 + 1] - a[j3 + 1];
400 a[j] = x0r + x2r;
401 a[j + 1] = x0i + x2i;
402 a[j2] = x0r - x2r;
403 a[j2 + 1] = x0i - x2i;
404 a[j1] = x1r - x3i;
405 a[j1 + 1] = x1i + x3r;
406 a[j3] = x1r + x3i;
407 a[j3 + 1] = x1i - x3r;
408 }
409 }
410 else
411 {
412 for (j = 0; j < l; j += 2)
413 {
414 j1 = j + l;
415 x0r = a[j] - a[j1];
416 x0i = a[j + 1] - a[j1 + 1];
417 a[j] += a[j1];
418 a[j + 1] += a[j1 + 1];
419 a[j1] = x0r;
420 a[j1 + 1] = x0i;
421 }
422 }
423 }
424
425 void cftbsub(int n, double* a, double* w)
426 {
427 int j, j1, j2, j3, l;
428 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
429
430 l = 2;
431 if (n > 8)
432 {
433 cft1st(n, a, w);
434 l = 8;
435 while ((l << 2) < n)
436 {
437 cftmdl(n, l, a, w);
438 l <<= 2;
439 }
440 }
441 if ((l << 2) == n)
442 {
443 for (j = 0; j < l; j += 2)
444 {
445 j1 = j + l;
446 j2 = j1 + l;
447 j3 = j2 + l;
448 x0r = a[j] + a[j1];
449 x0i = -a[j + 1] - a[j1 + 1];
450 x1r = a[j] - a[j1];
451 x1i = -a[j + 1] + a[j1 + 1];
452 x2r = a[j2] + a[j3];
453 x2i = a[j2 + 1] + a[j3 + 1];
454 x3r = a[j2] - a[j3];
455 x3i = a[j2 + 1] - a[j3 + 1];
456 a[j] = x0r + x2r;
457 a[j + 1] = x0i - x2i;
458 a[j2] = x0r - x2r;
459 a[j2 + 1] = x0i + x2i;
460 a[j1] = x1r - x3i;
461 a[j1 + 1] = x1i - x3r;
462 a[j3] = x1r + x3i;
463 a[j3 + 1] = x1i + x3r;
464 }
465 }
466 else
467 {
468 for (j = 0; j < l; j += 2)
469 {
470 j1 = j + l;
471 x0r = a[j] - a[j1];
472 x0i = -a[j + 1] + a[j1 + 1];
473 a[j] += a[j1];
474 a[j + 1] = -a[j + 1] - a[j1 + 1];
475 a[j1] = x0r;
476 a[j1 + 1] = x0i;
477 }
478 }
479 }
480
481 void cft1st(int n, double* a, double* w)
482 {
483 int j, k1, k2;
484 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
485 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
486
487 x0r = a[0] + a[2];
488 x0i = a[1] + a[3];
489 x1r = a[0] - a[2];
490 x1i = a[1] - a[3];
491 x2r = a[4] + a[6];
492 x2i = a[5] + a[7];
493 x3r = a[4] - a[6];
494 x3i = a[5] - a[7];
495 a[0] = x0r + x2r;
496 a[1] = x0i + x2i;
497 a[4] = x0r - x2r;
498 a[5] = x0i - x2i;
499 a[2] = x1r - x3i;
500 a[3] = x1i + x3r;
501 a[6] = x1r + x3i;
502 a[7] = x1i - x3r;
503 wk1r = w[2];
504 x0r = a[8] + a[10];
505 x0i = a[9] + a[11];
506 x1r = a[8] - a[10];
507 x1i = a[9] - a[11];
508 x2r = a[12] + a[14];
509 x2i = a[13] + a[15];
510 x3r = a[12] - a[14];
511 x3i = a[13] - a[15];
512 a[8] = x0r + x2r;
513 a[9] = x0i + x2i;
514 a[12] = x2i - x0i;
515 a[13] = x0r - x2r;
516 x0r = x1r - x3i;
517 x0i = x1i + x3r;
518 a[10] = wk1r * (x0r - x0i);
519 a[11] = wk1r * (x0r + x0i);
520 x0r = x3i + x1r;
521 x0i = x3r - x1i;
522 a[14] = wk1r * (x0i - x0r);
523 a[15] = wk1r * (x0i + x0r);
524 k1 = 0;
525 for (j = 16; j < n; j += 16)
526 {
527 k1 += 2;
528 k2 = 2 * k1;
529 wk2r = w[k1];
530 wk2i = w[k1 + 1];
531 wk1r = w[k2];
532 wk1i = w[k2 + 1];
533 wk3r = wk1r - 2 * wk2i * wk1i;
534 wk3i = 2 * wk2i * wk1r - wk1i;
535 x0r = a[j] + a[j + 2];
536 x0i = a[j + 1] + a[j + 3];
537 x1r = a[j] - a[j + 2];
538 x1i = a[j + 1] - a[j + 3];
539 x2r = a[j + 4] + a[j + 6];
540 x2i = a[j + 5] + a[j + 7];
541 x3r = a[j + 4] - a[j + 6];
542 x3i = a[j + 5] - a[j + 7];
543 a[j] = x0r + x2r;
544 a[j + 1] = x0i + x2i;
545 x0r -= x2r;
546 x0i -= x2i;
547 a[j + 4] = wk2r * x0r - wk2i * x0i;
548 a[j + 5] = wk2r * x0i + wk2i * x0r;
549 x0r = x1r - x3i;
550 x0i = x1i + x3r;
551 a[j + 2] = wk1r * x0r - wk1i * x0i;
552 a[j + 3] = wk1r * x0i + wk1i * x0r;
553 x0r = x1r + x3i;
554 x0i = x1i - x3r;
555 a[j + 6] = wk3r * x0r - wk3i * x0i;
556 a[j + 7] = wk3r * x0i + wk3i * x0r;
557 wk1r = w[k2 + 2];
558 wk1i = w[k2 + 3];
559 wk3r = wk1r - 2 * wk2r * wk1i;
560 wk3i = 2 * wk2r * wk1r - wk1i;
561 x0r = a[j + 8] + a[j + 10];
562 x0i = a[j + 9] + a[j + 11];
563 x1r = a[j + 8] - a[j + 10];
564 x1i = a[j + 9] - a[j + 11];
565 x2r = a[j + 12] + a[j + 14];
566 x2i = a[j + 13] + a[j + 15];
567 x3r = a[j + 12] - a[j + 14];
568 x3i = a[j + 13] - a[j + 15];
569 a[j + 8] = x0r + x2r;
570 a[j + 9] = x0i + x2i;
571 x0r -= x2r;
572 x0i -= x2i;
573 a[j + 12] = -wk2i * x0r - wk2r * x0i;
574 a[j + 13] = -wk2i * x0i + wk2r * x0r;
575 x0r = x1r - x3i;
576 x0i = x1i + x3r;
577 a[j + 10] = wk1r * x0r - wk1i * x0i;
578 a[j + 11] = wk1r * x0i + wk1i * x0r;
579 x0r = x1r + x3i;
580 x0i = x1i - x3r;
581 a[j + 14] = wk3r * x0r - wk3i * x0i;
582 a[j + 15] = wk3r * x0i + wk3i * x0r;
583 }
584 }
585
586 void cftmdl(int n, int l, double* a, double* w)
587 {
588 int j, j1, j2, j3, k, k1, k2, m, m2;
589 double wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
590 double x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
591
592 m = l << 2;
593 for (j = 0; j < l; j += 2)
594 {
595 j1 = j + l;
596 j2 = j1 + l;
597 j3 = j2 + l;
598 x0r = a[j] + a[j1];
599 x0i = a[j + 1] + a[j1 + 1];
600 x1r = a[j] - a[j1];
601 x1i = a[j + 1] - a[j1 + 1];
602 x2r = a[j2] + a[j3];
603 x2i = a[j2 + 1] + a[j3 + 1];
604 x3r = a[j2] - a[j3];
605 x3i = a[j2 + 1] - a[j3 + 1];
606 a[j] = x0r + x2r;
607 a[j + 1] = x0i + x2i;
608 a[j2] = x0r - x2r;
609 a[j2 + 1] = x0i - x2i;
610 a[j1] = x1r - x3i;
611 a[j1 + 1] = x1i + x3r;
612 a[j3] = x1r + x3i;
613 a[j3 + 1] = x1i - x3r;
614 }
615 wk1r = w[2];
616 for (j = m; j < l + m; j += 2)
617 {
618 j1 = j + l;
619 j2 = j1 + l;
620 j3 = j2 + l;
621 x0r = a[j] + a[j1];
622 x0i = a[j + 1] + a[j1 + 1];
623 x1r = a[j] - a[j1];
624 x1i = a[j + 1] - a[j1 + 1];
625 x2r = a[j2] + a[j3];
626 x2i = a[j2 + 1] + a[j3 + 1];
627 x3r = a[j2] - a[j3];
628 x3i = a[j2 + 1] - a[j3 + 1];
629 a[j] = x0r + x2r;
630 a[j + 1] = x0i + x2i;
631 a[j2] = x2i - x0i;
632 a[j2 + 1] = x0r - x2r;
633 x0r = x1r - x3i;
634 x0i = x1i + x3r;
635 a[j1] = wk1r * (x0r - x0i);
636 a[j1 + 1] = wk1r * (x0r + x0i);
637 x0r = x3i + x1r;
638 x0i = x3r - x1i;
639 a[j3] = wk1r * (x0i - x0r);
640 a[j3 + 1] = wk1r * (x0i + x0r);
641 }
642 k1 = 0;
643 m2 = 2 * m;
644 for (k = m2; k < n; k += m2)
645 {
646 k1 += 2;
647 k2 = 2 * k1;
648 wk2r = w[k1];
649 wk2i = w[k1 + 1];
650 wk1r = w[k2];
651 wk1i = w[k2 + 1];
652 wk3r = wk1r - 2 * wk2i * wk1i;
653 wk3i = 2 * wk2i * wk1r - wk1i;
654 for (j = k; j < l + k; j += 2)
655 {
656 j1 = j + l;
657 j2 = j1 + l;
658 j3 = j2 + l;
659 x0r = a[j] + a[j1];
660 x0i = a[j + 1] + a[j1 + 1];
661 x1r = a[j] - a[j1];
662 x1i = a[j + 1] - a[j1 + 1];
663 x2r = a[j2] + a[j3];
664 x2i = a[j2 + 1] + a[j3 + 1];
665 x3r = a[j2] - a[j3];
666 x3i = a[j2 + 1] - a[j3 + 1];
667 a[j] = x0r + x2r;
668 a[j + 1] = x0i + x2i;
669 x0r -= x2r;
670 x0i -= x2i;
671 a[j2] = wk2r * x0r - wk2i * x0i;
672 a[j2 + 1] = wk2r * x0i + wk2i * x0r;
673 x0r = x1r - x3i;
674 x0i = x1i + x3r;
675 a[j1] = wk1r * x0r - wk1i * x0i;
676 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
677 x0r = x1r + x3i;
678 x0i = x1i - x3r;
679 a[j3] = wk3r * x0r - wk3i * x0i;
680 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
681 }
682 wk1r = w[k2 + 2];
683 wk1i = w[k2 + 3];
684 wk3r = wk1r - 2 * wk2r * wk1i;
685 wk3i = 2 * wk2r * wk1r - wk1i;
686 for (j = k + m; j < l + (k + m); j += 2)
687 {
688 j1 = j + l;
689 j2 = j1 + l;
690 j3 = j2 + l;
691 x0r = a[j] + a[j1];
692 x0i = a[j + 1] + a[j1 + 1];
693 x1r = a[j] - a[j1];
694 x1i = a[j + 1] - a[j1 + 1];
695 x2r = a[j2] + a[j3];
696 x2i = a[j2 + 1] + a[j3 + 1];
697 x3r = a[j2] - a[j3];
698 x3i = a[j2 + 1] - a[j3 + 1];
699 a[j] = x0r + x2r;
700 a[j + 1] = x0i + x2i;
701 x0r -= x2r;
702 x0i -= x2i;
703 a[j2] = -wk2i * x0r - wk2r * x0i;
704 a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
705 x0r = x1r - x3i;
706 x0i = x1i + x3r;
707 a[j1] = wk1r * x0r - wk1i * x0i;
708 a[j1 + 1] = wk1r * x0i + wk1i * x0r;
709 x0r = x1r + x3i;
710 x0i = x1i - x3r;
711 a[j3] = wk3r * x0r - wk3i * x0i;
712 a[j3 + 1] = wk3r * x0i + wk3i * x0r;
713 }
714 }
715 }
716
717 void rftfsub(int n, double* a, int nc, double* c)
718 {
719 int j, k, kk, ks, m;
720 double wkr, wki, xr, xi, yr, yi;
721
722 m = n >> 1;
723 ks = 2 * nc / m;
724 kk = 0;
725 for (j = 2; j < m; j += 2)
726 {
727 k = n - j;
728 kk += ks;
729 wkr = 0.5 - c[nc - kk];
730 wki = c[kk];
731 xr = a[j] - a[k];
732 xi = a[j + 1] + a[k + 1];
733 yr = wkr * xr - wki * xi;
734 yi = wkr * xi + wki * xr;
735 a[j] -= yr;
736 a[j + 1] -= yi;
737 a[k] += yr;
738 a[k + 1] -= yi;
739 }
740 }
741
742 void rftbsub(int n, double* a, int nc, double* c)
743 {
744 int j, k, kk, ks, m;
745 double wkr, wki, xr, xi, yr, yi;
746
747 a[1] = -a[1];
748 m = n >> 1;
749 ks = 2 * nc / m;
750 kk = 0;
751 for (j = 2; j < m; j += 2)
752 {
753 k = n - j;
754 kk += ks;
755 wkr = 0.5 - c[nc - kk];
756 wki = c[kk];
757 xr = a[j] - a[k];
758 xi = a[j + 1] + a[k + 1];
759 yr = wkr * xr + wki * xi;
760 yi = wkr * xi - wki * xr;
761 a[j] -= yr;
762 a[j + 1] = yi - a[j + 1];
763 a[k] += yr;
764 a[k + 1] = yi - a[k + 1];
765 }
766 a[m + 1] = -a[m + 1];
767 }
768 };
769
774 typedef OouraFFT AudioFFTImplementation;
775
776#endif // AUDIOFFT_OOURA_USED
777
778 // ================================================================
779
780#ifdef AUDIOFFT_APPLE_ACCELERATE_USED
781
787 class AppleAccelerateFFT : public detail::AudioFFTImpl
788 {
789 public:
790 AppleAccelerateFFT() : detail::AudioFFTImpl(), _size(0), _powerOf2(0), _fftSetup(0), _re(), _im() {}
791
792 AppleAccelerateFFT(const AppleAccelerateFFT&) = delete;
793 AppleAccelerateFFT& operator=(const AppleAccelerateFFT&) = delete;
794
795 virtual ~AppleAccelerateFFT() { init(0); }
796
797 virtual void init(size_t size) override
798 {
799 if (_fftSetup)
800 {
801 vDSP_destroy_fftsetup(_fftSetup);
802 _size = 0;
803 _powerOf2 = 0;
804 _fftSetup = 0;
805 _re.clear();
806 _im.clear();
807 }
808
809 if (size > 0)
810 {
811 _size = size;
812 _powerOf2 = 0;
813 while ((1 << _powerOf2) < _size)
814 {
815 ++_powerOf2;
816 }
817 _fftSetup = vDSP_create_fftsetup(_powerOf2, FFT_RADIX2);
818 _re.resize(_size / 2);
819 _im.resize(_size / 2);
820 }
821 }
822
823 virtual void fft(const float* data, float* re, float* im) override
824 {
825 const size_t size2 = _size / 2;
826 DSPSplitComplex splitComplex;
827 splitComplex.realp = re;
828 splitComplex.imagp = im;
829 vDSP_ctoz(reinterpret_cast<const COMPLEX*>(data), 2, &splitComplex, 1, size2);
830 vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_FORWARD);
831 const float factor = 0.5f;
832 vDSP_vsmul(re, 1, &factor, re, 1, size2);
833 vDSP_vsmul(im, 1, &factor, im, 1, size2);
834 re[size2] = im[0];
835 im[0] = 0.0f;
836 im[size2] = 0.0f;
837 }
838
839 virtual void ifft(float* data, const float* re, const float* im) override
840 {
841 const size_t size2 = _size / 2;
842 ::memcpy(_re.data(), re, size2 * sizeof(float));
843 ::memcpy(_im.data(), im, size2 * sizeof(float));
844 _im[0] = re[size2];
845 DSPSplitComplex splitComplex;
846 splitComplex.realp = _re.data();
847 splitComplex.imagp = _im.data();
848 vDSP_fft_zrip(_fftSetup, &splitComplex, 1, _powerOf2, FFT_INVERSE);
849 vDSP_ztoc(&splitComplex, 1, reinterpret_cast<COMPLEX*>(data), 2, size2);
850 const float factor = 1.0f / static_cast<float>(_size);
851 vDSP_vsmul(data, 1, &factor, data, 1, _size);
852 }
853
854 private:
855 size_t _size;
856 size_t _powerOf2;
857 FFTSetup _fftSetup;
858 std::vector<float> _re;
859 std::vector<float> _im;
860 };
861
866 typedef AppleAccelerateFFT AudioFFTImplementation;
867
868#endif // AUDIOFFT_APPLE_ACCELERATE_USED
869
870 // ================================================================
871
872#ifdef AUDIOFFT_FFTW3_USED
873
879 class FFTW3FFT : public detail::AudioFFTImpl
880 {
881 public:
882 FFTW3FFT()
883 : detail::AudioFFTImpl(), _size(0), _complexSize(0), _planForward(0), _planBackward(0), _data(0), _re(0),
884 _im(0)
885 {
886 }
887
888 FFTW3FFT(const FFTW3FFT&) = delete;
889 FFTW3FFT& operator=(const FFTW3FFT&) = delete;
890
891 virtual ~FFTW3FFT() { init(0); }
892
893 virtual void init(size_t size) override
894 {
895 if (_size != size)
896 {
897 if (_size > 0)
898 {
899 fftwf_destroy_plan(_planForward);
900 fftwf_destroy_plan(_planBackward);
901 _planForward = 0;
902 _planBackward = 0;
903 _size = 0;
904 _complexSize = 0;
905
906 if (_data)
907 {
908 fftwf_free(_data);
909 _data = 0;
910 }
911
912 if (_re)
913 {
914 fftwf_free(_re);
915 _re = 0;
916 }
917
918 if (_im)
919 {
920 fftwf_free(_im);
921 _im = 0;
922 }
923 }
924
925 if (size > 0)
926 {
927 _size = size;
928 _complexSize = AudioFFT::ComplexSize(_size);
929 const size_t complexSize = AudioFFT::ComplexSize(_size);
930 _data = reinterpret_cast<float*>(fftwf_malloc(_size * sizeof(float)));
931 _re = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));
932 _im = reinterpret_cast<float*>(fftwf_malloc(complexSize * sizeof(float)));
933
934 fftw_iodim dim;
935 dim.n = static_cast<int>(size);
936 dim.is = 1;
937 dim.os = 1;
938 _planForward = fftwf_plan_guru_split_dft_r2c(1, &dim, 0, 0, _data, _re, _im, FFTW_MEASURE);
939 _planBackward = fftwf_plan_guru_split_dft_c2r(1, &dim, 0, 0, _re, _im, _data, FFTW_MEASURE);
940 }
941 }
942 }
943
944 virtual void fft(const float* data, float* re, float* im) override
945 {
946 ::memcpy(_data, data, _size * sizeof(float));
947 fftwf_execute_split_dft_r2c(_planForward, _data, _re, _im);
948 ::memcpy(re, _re, _complexSize * sizeof(float));
949 ::memcpy(im, _im, _complexSize * sizeof(float));
950 }
951
952 virtual void ifft(float* data, const float* re, const float* im) override
953 {
954 ::memcpy(_re, re, _complexSize * sizeof(float));
955 ::memcpy(_im, im, _complexSize * sizeof(float));
956 fftwf_execute_split_dft_c2r(_planBackward, _re, _im, _data);
957 detail::ScaleBuffer(data, _data, 1.0f / static_cast<float>(_size), _size);
958 }
959
960 private:
961 size_t _size;
962 size_t _complexSize;
963 fftwf_plan _planForward;
964 fftwf_plan _planBackward;
965 float* _data;
966 float* _re;
967 float* _im;
968 };
969
974 typedef FFTW3FFT AudioFFTImplementation;
975
976#endif // AUDIOFFT_FFTW3_USED
977
978 // =============================================================
979
980 AudioFFT::AudioFFT() : _impl(new AudioFFTImplementation()) {}
981
982 AudioFFT::~AudioFFT() {}
983
984 void AudioFFT::init(size_t size)
985 {
986 assert(detail::IsPowerOf2(size));
987 _impl->init(size);
988 }
989
990 void AudioFFT::fft(const float* data, float* re, float* im) { _impl->fft(data, re, im); }
991
992 void AudioFFT::ifft(float* data, const float* re, const float* im) { _impl->ifft(data, re, im); }
993
994 size_t AudioFFT::ComplexSize(size_t size) { return (size / 2) + 1; }
995
996} // namespace audiofft
audiofft::AudioFFT::_impl
std::unique_ptr< detail::AudioFFTImpl > _impl
Definition AudioFFT.h:151
audiofft::AudioFFT::init
void init(size_t size)
Initializes the FFT object.
Definition AudioFFT.cpp:984
audiofft::AudioFFT::fft
void fft(const float *data, float *re, float *im)
Performs the forward FFT.
Definition AudioFFT.cpp:990
audiofft::AudioFFT::ifft
void ifft(float *data, const float *re, const float *im)
Performs the inverse FFT.
Definition AudioFFT.cpp:992
audiofft::AudioFFT::~AudioFFT
~AudioFFT()
Destructor.
Definition AudioFFT.cpp:982
audiofft::AudioFFT::AudioFFT
AudioFFT()
Constructor.
Definition AudioFFT.cpp:980
audiofft::AudioFFT::ComplexSize
static size_t ComplexSize(size_t size)
Calculates the necessary size of the real/imaginary complex arrays.
Definition AudioFFT.cpp:994
audiofft::OouraFFT
Definition AudioFFT.cpp:93
audiofft::OouraFFT::makewt
void makewt(int nw, int *ip, double *w)
Definition AudioFFT.cpp:208
audiofft::OouraFFT::fft
virtual void fft(const float *data, float *re, float *im) override
Definition AudioFFT.cpp:115
audiofft::OouraFFT::makect
void makect(int nc, int *ip, double *c)
Definition AudioFFT.cpp:239
audiofft::OouraFFT::cft1st
void cft1st(int n, double *a, double *w)
Definition AudioFFT.cpp:481
audiofft::OouraFFT::rdft
void rdft(int n, int isgn, double *a, int *ip, double *w)
Definition AudioFFT.cpp:168
audiofft::OouraFFT::OouraFFT
OouraFFT(const OouraFFT &)=delete
audiofft::OouraFFT::init
virtual void init(size_t size) override
Definition AudioFFT.cpp:100
audiofft::OouraFFT::rftbsub
void rftbsub(int n, double *a, int nc, double *c)
Definition AudioFFT.cpp:742
audiofft::OouraFFT::_buffer
std::vector< double > _buffer
Definition AudioFFT.cpp:166
audiofft::OouraFFT::cftmdl
void cftmdl(int n, int l, double *a, double *w)
Definition AudioFFT.cpp:586
audiofft::OouraFFT::_w
std::vector< double > _w
Definition AudioFFT.cpp:165
audiofft::OouraFFT::ifft
virtual void ifft(float *data, const float *re, const float *im) override
Definition AudioFFT.cpp:140
audiofft::OouraFFT::bitrv2
void bitrv2(int n, int *ip, double *a)
Definition AudioFFT.cpp:261
audiofft::OouraFFT::rftfsub
void rftfsub(int n, double *a, int nc, double *c)
Definition AudioFFT.cpp:717
audiofft::OouraFFT::cftfsub
void cftfsub(int n, double *a, double *w)
Definition AudioFFT.cpp:369
audiofft::OouraFFT::OouraFFT
OouraFFT()
Definition AudioFFT.cpp:95
audiofft::OouraFFT::operator=
OouraFFT & operator=(const OouraFFT &)=delete
audiofft::OouraFFT::cftbsub
void cftbsub(int n, double *a, double *w)
Definition AudioFFT.cpp:425
audiofft::OouraFFT::_size
size_t _size
Definition AudioFFT.cpp:163
audiofft::OouraFFT::_ip
std::vector< int > _ip
Definition AudioFFT.cpp:164
audiofft::detail::AudioFFTImpl
Definition AudioFFT.cpp:50
audiofft::detail::AudioFFTImpl::~AudioFFTImpl
virtual ~AudioFFTImpl()=default
audiofft::detail::AudioFFTImpl::ifft
virtual void ifft(float *data, const float *re, const float *im)=0
audiofft::detail::AudioFFTImpl::init
virtual void init(size_t size)=0
audiofft::detail::AudioFFTImpl::operator=
AudioFFTImpl & operator=(const AudioFFTImpl &)=delete
audiofft::detail::AudioFFTImpl::AudioFFTImpl
AudioFFTImpl(const AudioFFTImpl &)=delete
audiofft::detail::AudioFFTImpl::fft
virtual void fft(const float *data, float *re, float *im)=0
audiofft::detail::IsPowerOf2
constexpr bool IsPowerOf2(size_t val)
Definition AudioFFT.cpp:61
audiofft::detail::ConvertBuffer
void ConvertBuffer(TypeDest *dest, const TypeSrc *src, size_t len)
Definition AudioFFT.cpp:64
audiofft::detail::ScaleBuffer
void ScaleBuffer(TypeDest *dest, const TypeSrc *src, const TypeFactor factor, size_t len)
Definition AudioFFT.cpp:73
audiofft
Definition AudioFFT.cpp:44
audiofft::AudioFFTImplementation
OouraFFT AudioFFTImplementation
Definition AudioFFT.cpp:774