1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
5 An audio time-stretching and pitch-shifting library.
6 Copyright 2007-2008 Chris Cannam.
8 This program is free software; you can redistribute it and/or
9 modify it under the terms of the GNU General Public License as
10 published by the Free Software Foundation; either version 2 of the
11 License, or (at your option) any later version. See the file
12 COPYING included with this distribution for more information.
19 //#define FFT_MEASUREMENT 1
22 #define HAVE_FFTW3 // for Ardour
32 #include "bsd-3rdparty/kissfft/kiss_fftr.h"
37 #ifndef USE_BUILTIN_FFT
38 #error No FFT implementation selected!
50 namespace RubberBand {
55 virtual ~FFTImpl() { }
57 virtual void initFloat() = 0;
58 virtual void initDouble() = 0;
60 virtual void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) = 0;
61 virtual void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) = 0;
62 virtual void forwardMagnitude(const double *R__ realIn, double *R__ magOut) = 0;
64 virtual void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) = 0;
65 virtual void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) = 0;
66 virtual void forwardMagnitude(const float *R__ realIn, float *R__ magOut) = 0;
68 virtual void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) = 0;
69 virtual void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) = 0;
70 virtual void inverseCepstral(const double *R__ magIn, double *R__ cepOut) = 0;
72 virtual void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) = 0;
73 virtual void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) = 0;
74 virtual void inverseCepstral(const float *R__ magIn, float *R__ cepOut) = 0;
76 virtual float *getFloatTimeBuffer() = 0;
77 virtual double *getDoubleTimeBuffer() = 0;
85 // Define FFTW_DOUBLE_ONLY to make all uses of FFTW functions be
86 // double-precision (so "float" FFTs are calculated by casting to
87 // doubles and using the double-precision FFTW function).
89 // Define FFTW_FLOAT_ONLY to make all uses of FFTW functions be
90 // single-precision (so "double" FFTs are calculated by casting to
91 // floats and using the single-precision FFTW function).
93 // Neither of these flags is terribly desirable -- FFTW_FLOAT_ONLY
94 // obviously loses you precision, and neither is handled in the most
95 // efficient way so any performance improvement will be small at best.
96 // The only real reason to define either flag would be to avoid
97 // linking against both fftw3 and fftw3f libraries.
99 //#define FFTW_DOUBLE_ONLY 1
100 //#define FFTW_FLOAT_ONLY 1
102 #if defined(FFTW_DOUBLE_ONLY) && defined(FFTW_FLOAT_ONLY)
103 // Can't meaningfully define both
104 #undef FFTW_DOUBLE_ONLY
105 #undef FFTW_FLOAT_ONLY
108 #ifdef FFTW_DOUBLE_ONLY
109 #define fft_float_type double
110 #define fftwf_complex fftw_complex
111 #define fftwf_plan fftw_plan
112 #define fftwf_plan_dft_r2c_1d fftw_plan_dft_r2c_1d
113 #define fftwf_plan_dft_c2r_1d fftw_plan_dft_c2r_1d
114 #define fftwf_destroy_plan fftw_destroy_plan
115 #define fftwf_malloc fftw_malloc
116 #define fftwf_free fftw_free
117 #define fftwf_execute fftw_execute
123 #define fft_float_type float
124 #endif /* FFTW_DOUBLE_ONLY */
126 #ifdef FFTW_FLOAT_ONLY
127 #define fft_double_type float
128 #define fftw_complex fftwf_complex
129 #define fftw_plan fftwf_plan
130 #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
131 #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
132 #define fftw_destroy_plan fftwf_destroy_plan
133 #define fftw_malloc fftwf_malloc
134 #define fftw_free fftwf_free
135 #define fftw_execute fftwf_execute
141 #define fft_double_type double
142 #endif /* FFTW_FLOAT_ONLY */
144 class D_FFTW : public FFTImpl
147 D_FFTW(int size) : m_fplanf(0)
148 #ifdef FFTW_DOUBLE_ONLY
152 #ifdef FFTW_FLOAT_ONLY
162 m_extantMutex.lock();
163 if (m_extantf > 0 && --m_extantf == 0) save = true;
164 m_extantMutex.unlock();
165 #ifndef FFTW_DOUBLE_ONLY
166 if (save) saveWisdom('f');
168 fftwf_destroy_plan(m_fplanf);
169 fftwf_destroy_plan(m_fplani);
171 fftwf_free(m_fpacked);
172 #ifdef FFTW_DOUBLE_ONLY
173 if (m_frb) fftw_free(m_frb);
178 m_extantMutex.lock();
179 if (m_extantd > 0 && --m_extantd == 0) save = true;
180 m_extantMutex.unlock();
181 #ifndef FFTW_FLOAT_ONLY
182 if (save) saveWisdom('d');
184 fftw_destroy_plan(m_dplanf);
185 fftw_destroy_plan(m_dplani);
187 fftw_free(m_dpacked);
188 #ifdef FFTW_FLOAT_ONLY
189 if (m_drb) fftwf_free(m_drb);
195 if (m_fplanf) return;
197 m_extantMutex.lock();
198 if (m_extantf++ == 0) load = true;
199 m_extantMutex.unlock();
200 #ifdef FFTW_DOUBLE_ONLY
201 if (load) loadWisdom('d');
203 if (load) loadWisdom('f');
205 m_fbuf = (fft_float_type *)fftw_malloc(m_size * sizeof(fft_float_type));
206 m_fpacked = (fftwf_complex *)fftw_malloc
207 ((m_size/2 + 1) * sizeof(fftwf_complex));
208 m_fplanf = fftwf_plan_dft_r2c_1d
209 (m_size, m_fbuf, m_fpacked, FFTW_MEASURE);
210 m_fplani = fftwf_plan_dft_c2r_1d
211 (m_size, m_fpacked, m_fbuf, FFTW_MEASURE);
215 if (m_dplanf) return;
217 m_extantMutex.lock();
218 if (m_extantd++ == 0) load = true;
219 m_extantMutex.unlock();
220 #ifdef FFTW_FLOAT_ONLY
221 if (load) loadWisdom('f');
223 if (load) loadWisdom('d');
225 m_dbuf = (fft_double_type *)fftw_malloc(m_size * sizeof(fft_double_type));
226 m_dpacked = (fftw_complex *)fftw_malloc
227 ((m_size/2 + 1) * sizeof(fftw_complex));
228 m_dplanf = fftw_plan_dft_r2c_1d
229 (m_size, m_dbuf, m_dpacked, FFTW_MEASURE);
230 m_dplani = fftw_plan_dft_c2r_1d
231 (m_size, m_dpacked, m_dbuf, FFTW_MEASURE);
234 void loadWisdom(char type) { wisdom(false, type); }
235 void saveWisdom(char type) { wisdom(true, type); }
237 void wisdom(bool save, char type) {
239 #ifdef FFTW_DOUBLE_ONLY
240 if (type == 'f') return;
242 #ifdef FFTW_FLOAT_ONLY
243 if (type == 'd') return;
246 const char *home = getenv("HOME");
250 snprintf(fn, 256, "%s/%s.%c", home, ".rubberband.wisdom", type);
252 FILE *f = fopen(fn, save ? "wb" : "rb");
257 #ifdef FFTW_DOUBLE_ONLY
260 case 'f': fftwf_export_wisdom_to_file(f); break;
262 #ifdef FFTW_FLOAT_ONLY
265 case 'd': fftw_export_wisdom_to_file(f); break;
271 #ifdef FFTW_DOUBLE_ONLY
274 case 'f': fftwf_import_wisdom_from_file(f); break;
276 #ifdef FFTW_FLOAT_ONLY
279 case 'd': fftw_import_wisdom_from_file(f); break;
288 void packFloat(const float *R__ re, const float *R__ im) {
289 const int hs = m_size/2;
290 fftwf_complex *const R__ fpacked = m_fpacked;
291 for (int i = 0; i <= hs; ++i) {
292 fpacked[i][0] = re[i];
295 for (int i = 0; i <= hs; ++i) {
296 fpacked[i][1] = im[i];
299 for (int i = 0; i <= hs; ++i) {
305 void packDouble(const double *R__ re, const double *R__ im) {
306 const int hs = m_size/2;
307 fftw_complex *const R__ dpacked = m_dpacked;
308 for (int i = 0; i <= hs; ++i) {
309 dpacked[i][0] = re[i];
312 for (int i = 0; i <= hs; ++i) {
313 dpacked[i][1] = im[i];
316 for (int i = 0; i <= hs; ++i) {
322 void unpackFloat(float *R__ re, float *R__ im) {
323 const int hs = m_size/2;
324 for (int i = 0; i <= hs; ++i) {
325 re[i] = m_fpacked[i][0];
328 for (int i = 0; i <= hs; ++i) {
329 im[i] = m_fpacked[i][1];
334 void unpackDouble(double *R__ re, double *R__ im) {
335 const int hs = m_size/2;
336 for (int i = 0; i <= hs; ++i) {
337 re[i] = m_dpacked[i][0];
340 for (int i = 0; i <= hs; ++i) {
341 im[i] = m_dpacked[i][1];
346 void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
347 if (!m_dplanf) initDouble();
348 const int sz = m_size;
349 fft_double_type *const R__ dbuf = m_dbuf;
350 #ifndef FFTW_FLOAT_ONLY
353 for (int i = 0; i < sz; ++i) {
356 fftw_execute(m_dplanf);
357 unpackDouble(realOut, imagOut);
360 void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
361 if (!m_dplanf) initDouble();
362 fft_double_type *const R__ dbuf = m_dbuf;
363 const int sz = m_size;
364 #ifndef FFTW_FLOAT_ONLY
367 for (int i = 0; i < sz; ++i) {
370 fftw_execute(m_dplanf);
371 const int hs = m_size/2;
372 for (int i = 0; i <= hs; ++i) {
373 magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
374 m_dpacked[i][1] * m_dpacked[i][1]);
376 for (int i = 0; i <= hs; ++i) {
377 phaseOut[i] = atan2(m_dpacked[i][1], m_dpacked[i][0]);
381 void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
382 if (!m_dplanf) initDouble();
383 fft_double_type *const R__ dbuf = m_dbuf;
384 const int sz = m_size;
385 #ifndef FFTW_FLOAT_ONLY
386 if (realIn != m_dbuf)
388 for (int i = 0; i < sz; ++i) {
391 fftw_execute(m_dplanf);
392 const int hs = m_size/2;
393 for (int i = 0; i <= hs; ++i) {
394 magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
395 m_dpacked[i][1] * m_dpacked[i][1]);
399 void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
400 if (!m_fplanf) initFloat();
401 fft_float_type *const R__ fbuf = m_fbuf;
402 const int sz = m_size;
403 #ifndef FFTW_DOUBLE_ONLY
406 for (int i = 0; i < sz; ++i) {
409 fftwf_execute(m_fplanf);
410 unpackFloat(realOut, imagOut);
413 void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
414 if (!m_fplanf) initFloat();
415 fft_float_type *const R__ fbuf = m_fbuf;
416 const int sz = m_size;
417 #ifndef FFTW_DOUBLE_ONLY
420 for (int i = 0; i < sz; ++i) {
423 fftwf_execute(m_fplanf);
424 const int hs = m_size/2;
425 for (int i = 0; i <= hs; ++i) {
426 magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
427 m_fpacked[i][1] * m_fpacked[i][1]);
429 for (int i = 0; i <= hs; ++i) {
430 phaseOut[i] = atan2f(m_fpacked[i][1], m_fpacked[i][0]) ;
434 void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
435 if (!m_fplanf) initFloat();
436 fft_float_type *const R__ fbuf = m_fbuf;
437 const int sz = m_size;
438 #ifndef FFTW_DOUBLE_ONLY
441 for (int i = 0; i < sz; ++i) {
444 fftwf_execute(m_fplanf);
445 const int hs = m_size/2;
446 for (int i = 0; i <= hs; ++i) {
447 magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
448 m_fpacked[i][1] * m_fpacked[i][1]);
452 void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
453 if (!m_dplanf) initDouble();
454 packDouble(realIn, imagIn);
455 fftw_execute(m_dplani);
456 const int sz = m_size;
457 fft_double_type *const R__ dbuf = m_dbuf;
458 #ifndef FFTW_FLOAT_ONLY
461 for (int i = 0; i < sz; ++i) {
462 realOut[i] = dbuf[i];
466 void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
467 if (!m_dplanf) initDouble();
468 const int hs = m_size/2;
469 fftw_complex *const R__ dpacked = m_dpacked;
470 for (int i = 0; i <= hs; ++i) {
471 dpacked[i][0] = magIn[i] * cos(phaseIn[i]);
473 for (int i = 0; i <= hs; ++i) {
474 dpacked[i][1] = magIn[i] * sin(phaseIn[i]);
476 fftw_execute(m_dplani);
477 const int sz = m_size;
478 fft_double_type *const R__ dbuf = m_dbuf;
479 #ifndef FFTW_FLOAT_ONLY
482 for (int i = 0; i < sz; ++i) {
483 realOut[i] = dbuf[i];
487 void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
488 if (!m_dplanf) initDouble();
489 fft_double_type *const R__ dbuf = m_dbuf;
490 fftw_complex *const R__ dpacked = m_dpacked;
491 const int hs = m_size/2;
492 for (int i = 0; i <= hs; ++i) {
493 dpacked[i][0] = log(magIn[i] + 0.000001);
495 for (int i = 0; i <= hs; ++i) {
498 fftw_execute(m_dplani);
499 const int sz = m_size;
500 #ifndef FFTW_FLOAT_ONLY
503 for (int i = 0; i < sz; ++i) {
508 void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
509 if (!m_fplanf) initFloat();
510 packFloat(realIn, imagIn);
511 fftwf_execute(m_fplani);
512 const int sz = m_size;
513 fft_float_type *const R__ fbuf = m_fbuf;
514 #ifndef FFTW_DOUBLE_ONLY
517 for (int i = 0; i < sz; ++i) {
518 realOut[i] = fbuf[i];
522 void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
523 if (!m_fplanf) initFloat();
524 const int hs = m_size/2;
525 fftwf_complex *const R__ fpacked = m_fpacked;
526 for (int i = 0; i <= hs; ++i) {
527 fpacked[i][0] = magIn[i] * cosf(phaseIn[i]);
529 for (int i = 0; i <= hs; ++i) {
530 fpacked[i][1] = magIn[i] * sinf(phaseIn[i]);
532 fftwf_execute(m_fplani);
533 const int sz = m_size;
534 fft_float_type *const R__ fbuf = m_fbuf;
535 #ifndef FFTW_DOUBLE_ONLY
538 for (int i = 0; i < sz; ++i) {
539 realOut[i] = fbuf[i];
543 void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
544 if (!m_fplanf) initFloat();
545 const int hs = m_size/2;
546 fftwf_complex *const R__ fpacked = m_fpacked;
547 for (int i = 0; i <= hs; ++i) {
548 fpacked[i][0] = logf(magIn[i] + 0.000001f);
550 for (int i = 0; i <= hs; ++i) {
553 fftwf_execute(m_fplani);
554 const int sz = m_size;
555 fft_float_type *const R__ fbuf = m_fbuf;
556 #ifndef FFTW_DOUBLE_ONLY
559 for (int i = 0; i < sz; ++i) {
564 float *getFloatTimeBuffer() {
566 #ifdef FFTW_DOUBLE_ONLY
567 if (!m_frb) m_frb = (float *)fftw_malloc(m_size * sizeof(float));
574 double *getDoubleTimeBuffer() {
576 #ifdef FFTW_FLOAT_ONLY
577 if (!m_drb) m_drb = (double *)fftwf_malloc(m_size * sizeof(double));
587 #ifdef FFTW_DOUBLE_ONLY
593 fftwf_complex *m_fpacked;
596 #ifdef FFTW_FLOAT_ONLY
602 fftw_complex * m_dpacked;
604 static int m_extantf;
605 static int m_extantd;
606 static Mutex m_extantMutex;
610 D_FFTW::m_extantf = 0;
613 D_FFTW::m_extantd = 0;
616 D_FFTW::m_extantMutex;
618 #endif /* HAVE_FFTW3 */
622 class D_KISSFFT : public FFTImpl
625 D_KISSFFT(int size) :
633 #error KISSFFT is not configured for float values
635 if (sizeof(kiss_fft_scalar) != sizeof(float)) {
636 std::cerr << "ERROR: KISSFFT is not configured for float values"
640 m_fbuf = new kiss_fft_scalar[m_size + 2];
641 m_fpacked = new kiss_fft_cpx[m_size + 2];
642 m_fplanf = kiss_fftr_alloc(m_size, 0, NULL, NULL);
643 m_fplani = kiss_fftr_alloc(m_size, 1, NULL, NULL);
647 kiss_fftr_free(m_fplanf);
648 kiss_fftr_free(m_fplani);
654 if (m_frb) delete[] m_frb;
655 if (m_drb) delete[] m_drb;
659 void initDouble() { }
661 void packFloat(const float *R__ re, const float *R__ im) {
662 const int hs = m_size/2;
663 for (int i = 0; i <= hs; ++i) {
664 m_fpacked[i].r = re[i];
665 m_fpacked[i].i = im[i];
669 void unpackFloat(float *R__ re, float *R__ im) {
670 const int hs = m_size/2;
671 for (int i = 0; i <= hs; ++i) {
672 re[i] = m_fpacked[i].r;
673 im[i] = m_fpacked[i].i;
677 void packDouble(const double *R__ re, const double *R__ im) {
678 const int hs = m_size/2;
679 for (int i = 0; i <= hs; ++i) {
680 m_fpacked[i].r = float(re[i]);
681 m_fpacked[i].i = float(im[i]);
685 void unpackDouble(double *R__ re, double *R__ im) {
686 const int hs = m_size/2;
687 for (int i = 0; i <= hs; ++i) {
688 re[i] = double(m_fpacked[i].r);
689 im[i] = double(m_fpacked[i].i);
693 void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
695 for (int i = 0; i < m_size; ++i) {
696 m_fbuf[i] = float(realIn[i]);
699 kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
700 unpackDouble(realOut, imagOut);
703 void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
705 for (int i = 0; i < m_size; ++i) {
706 m_fbuf[i] = float(realIn[i]);
709 kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
711 const int hs = m_size/2;
713 for (int i = 0; i <= hs; ++i) {
714 magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
715 double(m_fpacked[i].i) * double(m_fpacked[i].i));
718 for (int i = 0; i <= hs; ++i) {
719 phaseOut[i] = atan2(double(m_fpacked[i].i), double(m_fpacked[i].r));
723 void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
725 for (int i = 0; i < m_size; ++i) {
726 m_fbuf[i] = float(realIn[i]);
729 kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
731 const int hs = m_size/2;
733 for (int i = 0; i <= hs; ++i) {
734 magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
735 double(m_fpacked[i].i) * double(m_fpacked[i].i));
739 void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
741 kiss_fftr(m_fplanf, realIn, m_fpacked);
742 unpackFloat(realOut, imagOut);
745 void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
747 kiss_fftr(m_fplanf, realIn, m_fpacked);
749 const int hs = m_size/2;
751 for (int i = 0; i <= hs; ++i) {
752 magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
753 m_fpacked[i].i * m_fpacked[i].i);
756 for (int i = 0; i <= hs; ++i) {
757 phaseOut[i] = atan2f(m_fpacked[i].i, m_fpacked[i].r);
761 void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
763 kiss_fftr(m_fplanf, realIn, m_fpacked);
765 const int hs = m_size/2;
767 for (int i = 0; i <= hs; ++i) {
768 magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
769 m_fpacked[i].i * m_fpacked[i].i);
773 void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
775 packDouble(realIn, imagIn);
777 kiss_fftri(m_fplani, m_fpacked, m_fbuf);
779 for (int i = 0; i < m_size; ++i) {
780 realOut[i] = m_fbuf[i];
784 void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
786 const int hs = m_size/2;
788 for (int i = 0; i <= hs; ++i) {
789 m_fpacked[i].r = float(magIn[i] * cos(phaseIn[i]));
790 m_fpacked[i].i = float(magIn[i] * sin(phaseIn[i]));
793 kiss_fftri(m_fplani, m_fpacked, m_fbuf);
795 for (int i = 0; i < m_size; ++i) {
796 realOut[i] = m_fbuf[i];
800 void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
802 const int hs = m_size/2;
804 for (int i = 0; i <= hs; ++i) {
805 m_fpacked[i].r = float(log(magIn[i] + 0.000001));
806 m_fpacked[i].i = 0.0f;
809 kiss_fftri(m_fplani, m_fpacked, m_fbuf);
811 for (int i = 0; i < m_size; ++i) {
812 cepOut[i] = m_fbuf[i];
816 void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
818 packFloat(realIn, imagIn);
819 kiss_fftri(m_fplani, m_fpacked, realOut);
822 void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
824 const int hs = m_size/2;
826 for (int i = 0; i <= hs; ++i) {
827 m_fpacked[i].r = magIn[i] * cosf(phaseIn[i]);
828 m_fpacked[i].i = magIn[i] * sinf(phaseIn[i]);
831 kiss_fftri(m_fplani, m_fpacked, realOut);
834 void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
836 const int hs = m_size/2;
838 for (int i = 0; i <= hs; ++i) {
839 m_fpacked[i].r = logf(magIn[i] + 0.000001f);
840 m_fpacked[i].i = 0.0f;
843 kiss_fftri(m_fplani, m_fpacked, cepOut);
846 float *getFloatTimeBuffer() {
847 if (!m_frb) m_frb = new float[m_size];
851 double *getDoubleTimeBuffer() {
852 if (!m_drb) m_drb = new double[m_size];
860 kiss_fftr_cfg m_fplanf;
861 kiss_fftr_cfg m_fplani;
862 kiss_fft_scalar *m_fbuf;
863 kiss_fft_cpx *m_fpacked;
866 #endif /* USE_KISSFFT */
868 #ifdef USE_BUILTIN_FFT
870 class D_Cross : public FFTImpl
873 D_Cross(int size) : m_size(size), m_table(0), m_frb(0), m_drb(0) {
875 m_a = new double[size];
876 m_b = new double[size];
877 m_c = new double[size];
878 m_d = new double[size];
880 m_table = new int[m_size];
886 if (m_size & (1 << i)) {
892 for (i = 0; i < m_size; ++i) {
896 for (j = k = 0; j < bits; ++j) {
897 k = (k << 1) | (m & 1);
916 void initDouble() { }
918 void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
919 basefft(false, realIn, 0, m_c, m_d);
920 const int hs = m_size/2;
921 for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
923 for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
927 void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
928 basefft(false, realIn, 0, m_c, m_d);
929 const int hs = m_size/2;
930 for (int i = 0; i <= hs; ++i) {
931 magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
932 phaseOut[i] = atan2(m_d[i], m_c[i]) ;
936 void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
937 basefft(false, realIn, 0, m_c, m_d);
938 const int hs = m_size/2;
939 for (int i = 0; i <= hs; ++i) {
940 magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
944 void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
945 for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
946 basefft(false, m_a, 0, m_c, m_d);
947 const int hs = m_size/2;
948 for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
950 for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
954 void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
955 for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
956 basefft(false, m_a, 0, m_c, m_d);
957 const int hs = m_size/2;
958 for (int i = 0; i <= hs; ++i) {
959 magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
960 phaseOut[i] = atan2(m_d[i], m_c[i]) ;
964 void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
965 for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
966 basefft(false, m_a, 0, m_c, m_d);
967 const int hs = m_size/2;
968 for (int i = 0; i <= hs; ++i) {
969 magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
973 void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
974 const int hs = m_size/2;
975 for (int i = 0; i <= hs; ++i) {
976 double real = realIn[i];
977 double imag = imagIn[i];
981 m_a[m_size-i] = real;
982 m_b[m_size-i] = -imag;
985 basefft(true, m_a, m_b, realOut, m_d);
988 void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
989 const int hs = m_size/2;
990 for (int i = 0; i <= hs; ++i) {
991 double real = magIn[i] * cos(phaseIn[i]);
992 double imag = magIn[i] * sin(phaseIn[i]);
996 m_a[m_size-i] = real;
997 m_b[m_size-i] = -imag;
1000 basefft(true, m_a, m_b, realOut, m_d);
1003 void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
1004 const int hs = m_size/2;
1005 for (int i = 0; i <= hs; ++i) {
1006 double real = log(magIn[i] + 0.000001);
1010 m_a[m_size-i] = real;
1011 m_b[m_size-i] = 0.0;
1014 basefft(true, m_a, m_b, cepOut, m_d);
1017 void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
1018 const int hs = m_size/2;
1019 for (int i = 0; i <= hs; ++i) {
1020 float real = realIn[i];
1021 float imag = imagIn[i];
1025 m_a[m_size-i] = real;
1026 m_b[m_size-i] = -imag;
1029 basefft(true, m_a, m_b, m_c, m_d);
1030 for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
1033 void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
1034 const int hs = m_size/2;
1035 for (int i = 0; i <= hs; ++i) {
1036 float real = magIn[i] * cosf(phaseIn[i]);
1037 float imag = magIn[i] * sinf(phaseIn[i]);
1041 m_a[m_size-i] = real;
1042 m_b[m_size-i] = -imag;
1045 basefft(true, m_a, m_b, m_c, m_d);
1046 for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
1049 void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
1050 const int hs = m_size/2;
1051 for (int i = 0; i <= hs; ++i) {
1052 float real = logf(magIn[i] + 0.000001);
1056 m_a[m_size-i] = real;
1057 m_b[m_size-i] = 0.0;
1060 basefft(true, m_a, m_b, m_c, m_d);
1061 for (int i = 0; i < m_size; ++i) cepOut[i] = m_c[i];
1064 float *getFloatTimeBuffer() {
1065 if (!m_frb) m_frb = new float[m_size];
1069 double *getDoubleTimeBuffer() {
1070 if (!m_drb) m_drb = new double[m_size];
1083 void basefft(bool inverse, const double *R__ ri, const double *R__ ii, double *R__ ro, double *R__ io);
1087 D_Cross::basefft(bool inverse, const double *R__ ri, const double *R__ ii, double *R__ ro, double *R__ io)
1089 if (!ri || !ro || !io) return;
1092 int blockSize, blockEnd;
1096 double angle = 2.0 * M_PI;
1097 if (inverse) angle = -angle;
1099 const int n = m_size;
1102 for (i = 0; i < n; ++i) {
1103 ro[m_table[i]] = ri[i];
1105 for (i = 0; i < n; ++i) {
1106 io[m_table[i]] = ii[i];
1109 for (i = 0; i < n; ++i) {
1110 ro[m_table[i]] = ri[i];
1112 for (i = 0; i < n; ++i) {
1113 io[m_table[i]] = 0.0;
1119 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
1121 double delta = angle / (double)blockSize;
1122 double sm2 = -sin(-2 * delta);
1123 double sm1 = -sin(-delta);
1124 double cm2 = cos(-2 * delta);
1125 double cm1 = cos(-delta);
1127 double ar[3], ai[3];
1129 for (i = 0; i < n; i += blockSize) {
1137 for (j = i, m = 0; m < blockEnd; j++, m++) {
1139 ar[0] = w * ar[1] - ar[2];
1143 ai[0] = w * ai[1] - ai[2];
1148 tr = ar[0] * ro[k] - ai[0] * io[k];
1149 ti = ar[0] * io[k] + ai[0] * ro[k];
1159 blockEnd = blockSize;
1162 /* fftw doesn't rescale, so nor will we
1166 double denom = (double)n;
1168 for (i = 0; i < n; i++) {
1176 #endif /* USE_BUILTIN_FFT */
1178 } /* end namespace FFTs */
1183 FFT::FFT(int size, int debugLevel)
1186 (size & (size-1))) {
1187 std::cerr << "FFT::FFT(" << size << "): power-of-two sizes only supported, minimum size 2" << std::endl;
1191 if (m_method == -1) {
1204 std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1205 #ifdef USE_BUILTIN_FFT
1206 d = new FFTs::D_Cross(size);
1208 std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1215 if (debugLevel > 0) {
1216 std::cerr << "FFT::FFT(" << size << "): using FFTW3 implementation"
1219 d = new FFTs::D_FFTW(size);
1221 std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1222 #ifdef USE_BUILTIN_FFT
1223 d = new FFTs::D_Cross(size);
1225 std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1233 if (debugLevel > 0) {
1234 std::cerr << "FFT::FFT(" << size << "): using KISSFFT implementation"
1237 d = new FFTs::D_KISSFFT(size);
1239 std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1240 #ifdef USE_BUILTIN_FFT
1241 d = new FFTs::D_Cross(size);
1243 std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1250 #ifdef USE_BUILTIN_FFT
1251 std::cerr << "FFT::FFT(" << size << "): WARNING: using slow built-in implementation" << std::endl;
1252 d = new FFTs::D_Cross(size);
1254 std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1267 FFT::forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut)
1269 d->forward(realIn, realOut, imagOut);
1273 FFT::forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut)
1275 d->forwardPolar(realIn, magOut, phaseOut);
1279 FFT::forwardMagnitude(const double *R__ realIn, double *R__ magOut)
1281 d->forwardMagnitude(realIn, magOut);
1285 FFT::forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut)
1287 d->forward(realIn, realOut, imagOut);
1291 FFT::forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut)
1293 d->forwardPolar(realIn, magOut, phaseOut);
1297 FFT::forwardMagnitude(const float *R__ realIn, float *R__ magOut)
1299 d->forwardMagnitude(realIn, magOut);
1303 FFT::inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut)
1305 d->inverse(realIn, imagIn, realOut);
1309 FFT::inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut)
1311 d->inversePolar(magIn, phaseIn, realOut);
1315 FFT::inverseCepstral(const double *R__ magIn, double *R__ cepOut)
1317 d->inverseCepstral(magIn, cepOut);
1321 FFT::inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut)
1323 d->inverse(realIn, imagIn, realOut);
1327 FFT::inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut)
1329 d->inversePolar(magIn, phaseIn, realOut);
1333 FFT::inverseCepstral(const float *R__ magIn, float *R__ cepOut)
1335 d->inverseCepstral(magIn, cepOut);
1351 FFT::getFloatTimeBuffer()
1353 return d->getFloatTimeBuffer();
1357 FFT::getDoubleTimeBuffer()
1359 return d->getDoubleTimeBuffer();