libs/rubberband/src/FFT.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     Rubber Band
   5     An audio time-stretching and pitch-shifting library.
   6     Copyright 2007-2008 Chris Cannam.
   7
   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.
  13 */
  14
  15 #include "FFT.h"
  16 #include "Thread.h"
  17 #include "Profiler.h"
  18
  19 //#define FFT_MEASUREMENT 1
  20
  21 #define HAVE_FFTW3 // for Ardour
  22
  23 #ifdef HAVE_FFTW3
  24 #include <fftw3.h>
  25 #endif
  26
  27 #include <cstdlib>
  28
  29 #ifdef USE_KISSFFT
  30 #include "bsd-3rdparty/kissfft/kiss_fftr.h"
  31 #endif
  32
  33 #ifndef HAVE_FFTW3
  34 #ifndef USE_KISSFFT
  35 #ifndef USE_BUILTIN_FFT
  36 #error No FFT implementation selected!
  37 #endif
  38 #endif
  39 #endif
  40
  41 #include <cmath>
  42 #include <iostream>
  43 #include <map>
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <vector>
  47
  48 namespace RubberBand {
  49
  50 class FFTImpl
  51 {
  52 public:
  53     virtual ~FFTImpl() { }
  54
  55     virtual void initFloat() = 0;
  56     virtual void initDouble() = 0;
  57
  58     virtual void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) = 0;
  59     virtual void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) = 0;
  60     virtual void forwardMagnitude(const double *R__ realIn, double *R__ magOut) = 0;
  61
  62     virtual void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) = 0;
  63     virtual void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) = 0;
  64     virtual void forwardMagnitude(const float *R__ realIn, float *R__ magOut) = 0;
  65
  66     virtual void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) = 0;
  67     virtual void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) = 0;
  68     virtual void inverseCepstral(const double *R__ magIn, double *R__ cepOut) = 0;
  69
  70     virtual void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) = 0;
  71     virtual void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) = 0;
  72     virtual void inverseCepstral(const float *R__ magIn, float *R__ cepOut) = 0;
  73
  74     virtual float *getFloatTimeBuffer() = 0;
  75     virtual double *getDoubleTimeBuffer() = 0;
  76 };
  77
  78 namespace FFTs {
  79
  80
  81 #ifdef HAVE_FFTW3
  82
  83 // Define FFTW_DOUBLE_ONLY to make all uses of FFTW functions be
  84 // double-precision (so "float" FFTs are calculated by casting to
  85 // doubles and using the double-precision FFTW function).
  86 //
  87 // Define FFTW_FLOAT_ONLY to make all uses of FFTW functions be
  88 // single-precision (so "double" FFTs are calculated by casting to
  89 // floats and using the single-precision FFTW function).
  90 //
  91 // Neither of these flags is terribly desirable -- FFTW_FLOAT_ONLY
  92 // obviously loses you precision, and neither is handled in the most
  93 // efficient way so any performance improvement will be small at best.
  94 // The only real reason to define either flag would be to avoid
  95 // linking against both fftw3 and fftw3f libraries.
  96
  97 //#define FFTW_DOUBLE_ONLY 1
  98 //#define FFTW_FLOAT_ONLY 1
  99
 100 #if defined(FFTW_DOUBLE_ONLY) && defined(FFTW_FLOAT_ONLY)
 101 // Can't meaningfully define both
 102 #undef FFTW_DOUBLE_ONLY
 103 #undef FFTW_FLOAT_ONLY
 104 #endif
 105
 106 #ifdef FFTW_DOUBLE_ONLY
 107 #define fft_float_type double
 108 #define fftwf_complex fftw_complex
 109 #define fftwf_plan fftw_plan
 110 #define fftwf_plan_dft_r2c_1d fftw_plan_dft_r2c_1d
 111 #define fftwf_plan_dft_c2r_1d fftw_plan_dft_c2r_1d
 112 #define fftwf_destroy_plan fftw_destroy_plan
 113 #define fftwf_malloc fftw_malloc
 114 #define fftwf_free fftw_free
 115 #define fftwf_execute fftw_execute
 116 #define atan2f atan2
 117 #define sqrtf sqrt
 118 #define cosf cos
 119 #define sinf sin
 120 #else
 121 #define fft_float_type float
 122 #endif /* FFTW_DOUBLE_ONLY */
 123
 124 #ifdef FFTW_FLOAT_ONLY
 125 #define fft_double_type float
 126 #define fftw_complex fftwf_complex
 127 #define fftw_plan fftwf_plan
 128 #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
 129 #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
 130 #define fftw_destroy_plan fftwf_destroy_plan
 131 #define fftw_malloc fftwf_malloc
 132 #define fftw_free fftwf_free
 133 #define fftw_execute fftwf_execute
 134 #define atan2 atan2f
 135 #define sqrt sqrtf
 136 #define cos cosf
 137 #define sin sinf
 138 #else
 139 #define fft_double_type double
 140 #endif /* FFTW_FLOAT_ONLY */
 141
 142 class D_FFTW : public FFTImpl
 143 {
 144 public:
 145     D_FFTW(int size) : m_fplanf(0)
 146 #ifdef FFTW_DOUBLE_ONLY
 147                               , m_frb(0)
 148 #endif
 149                               , m_dplanf(0)
 150 #ifdef FFTW_FLOAT_ONLY
 151                               , m_drb(0)
 152 #endif
 153                               , m_size(size)
 154     {
 155     }
 156
 157     ~D_FFTW() {
 158         if (m_fplanf) {
 159             bool save = false;
 160             m_extantMutex.lock();
 161             if (m_extantf > 0 && --m_extantf == 0) save = true;
 162             m_extantMutex.unlock();
 163 #ifndef FFTW_DOUBLE_ONLY
 164             if (save) saveWisdom('f');
 165 #endif
 166             fftwf_destroy_plan(m_fplanf);
 167             fftwf_destroy_plan(m_fplani);
 168             fftwf_free(m_fbuf);
 169             fftwf_free(m_fpacked);
 170 #ifdef FFTW_DOUBLE_ONLY
 171             if (m_frb) fftw_free(m_frb);
 172 #endif
 173         }
 174         if (m_dplanf) {
 175             bool save = false;
 176             m_extantMutex.lock();
 177             if (m_extantd > 0 && --m_extantd == 0) save = true;
 178             m_extantMutex.unlock();
 179 #ifndef FFTW_FLOAT_ONLY
 180             if (save) saveWisdom('d');
 181 #endif
 182             fftw_destroy_plan(m_dplanf);
 183             fftw_destroy_plan(m_dplani);
 184             fftw_free(m_dbuf);
 185             fftw_free(m_dpacked);
 186 #ifdef FFTW_FLOAT_ONLY
 187             if (m_drb) fftwf_free(m_drb);
 188 #endif
 189         }
 190     }
 191
 192     void initFloat() {
 193         if (m_fplanf) return;
 194         bool load = false;
 195         m_extantMutex.lock();
 196         if (m_extantf++ == 0) load = true;
 197         m_extantMutex.unlock();
 198 #ifdef FFTW_DOUBLE_ONLY
 199         if (load) loadWisdom('d');
 200 #else
 201         if (load) loadWisdom('f');
 202 #endif
 203         m_fbuf = (fft_float_type *)fftw_malloc(m_size * sizeof(fft_float_type));
 204         m_fpacked = (fftwf_complex *)fftw_malloc
 205             ((m_size/2 + 1) * sizeof(fftwf_complex));
 206         m_fplanf = fftwf_plan_dft_r2c_1d
 207             (m_size, m_fbuf, m_fpacked, FFTW_MEASURE);
 208         m_fplani = fftwf_plan_dft_c2r_1d
 209             (m_size, m_fpacked, m_fbuf, FFTW_MEASURE);
 210     }
 211
 212     void initDouble() {
 213         if (m_dplanf) return;
 214         bool load = false;
 215         m_extantMutex.lock();
 216         if (m_extantd++ == 0) load = true;
 217         m_extantMutex.unlock();
 218 #ifdef FFTW_FLOAT_ONLY
 219         if (load) loadWisdom('f');
 220 #else
 221         if (load) loadWisdom('d');
 222 #endif
 223         m_dbuf = (fft_double_type *)fftw_malloc(m_size * sizeof(fft_double_type));
 224         m_dpacked = (fftw_complex *)fftw_malloc
 225             ((m_size/2 + 1) * sizeof(fftw_complex));
 226         m_dplanf = fftw_plan_dft_r2c_1d
 227             (m_size, m_dbuf, m_dpacked, FFTW_MEASURE);
 228         m_dplani = fftw_plan_dft_c2r_1d
 229             (m_size, m_dpacked, m_dbuf, FFTW_MEASURE);
 230     }
 231
 232     void loadWisdom(char type) { wisdom(false, type); }
 233     void saveWisdom(char type) { wisdom(true, type); }
 234
 235     void wisdom(bool save, char type) {
 236
 237 #ifdef FFTW_DOUBLE_ONLY
 238         if (type == 'f') return;
 239 #endif
 240 #ifdef FFTW_FLOAT_ONLY
 241         if (type == 'd') return;
 242 #endif
 243
 244         const char *home = getenv("HOME");
 245         if (!home) return;
 246
 247         char fn[256];
 248         snprintf(fn, 256, "%s/%s.%c", home, ".rubberband.wisdom", type);
 249
 250         FILE *f = fopen(fn, save ? "wb" : "rb");
 251         if (!f) return;
 252
 253         if (save) {
 254             switch (type) {
 255 #ifdef FFTW_DOUBLE_ONLY
 256             case 'f': break;
 257 #else
 258             case 'f': fftwf_export_wisdom_to_file(f); break;
 259 #endif
 260 #ifdef FFTW_FLOAT_ONLY
 261             case 'd': break;
 262 #else
 263             case 'd': fftw_export_wisdom_to_file(f); break;
 264 #endif
 265             default: break;
 266             }
 267         } else {
 268             switch (type) {
 269 #ifdef FFTW_DOUBLE_ONLY
 270             case 'f': break;
 271 #else
 272             case 'f': fftwf_import_wisdom_from_file(f); break;
 273 #endif
 274 #ifdef FFTW_FLOAT_ONLY
 275             case 'd': break;
 276 #else
 277             case 'd': fftw_import_wisdom_from_file(f); break;
 278 #endif
 279             default: break;
 280             }
 281         }
 282
 283         fclose(f);
 284     }
 285
 286     void packFloat(const float *R__ re, const float *R__ im) {
 287         const int hs = m_size/2;
 288         fftwf_complex *const R__ fpacked = m_fpacked;
 289         for (int i = 0; i <= hs; ++i) {
 290             fpacked[i][0] = re[i];
 291         }
 292         if (im) {
 293             for (int i = 0; i <= hs; ++i) {
 294                 fpacked[i][1] = im[i];
 295             }
 296         } else {
 297             for (int i = 0; i <= hs; ++i) {
 298                 fpacked[i][1] = 0.f;
 299             }
 300         }
 301     }
 302
 303     void packDouble(const double *R__ re, const double *R__ im) {
 304         const int hs = m_size/2;
 305         fftw_complex *const R__ dpacked = m_dpacked;
 306         for (int i = 0; i <= hs; ++i) {
 307             dpacked[i][0] = re[i];
 308         }
 309         if (im) {
 310             for (int i = 0; i <= hs; ++i) {
 311                 dpacked[i][1] = im[i];
 312             }
 313         } else {
 314             for (int i = 0; i <= hs; ++i) {
 315                 dpacked[i][1] = 0.0;
 316             }
 317         }
 318     }
 319
 320     void unpackFloat(float *R__ re, float *R__ im) {
 321         const int hs = m_size/2;
 322         for (int i = 0; i <= hs; ++i) {
 323             re[i] = m_fpacked[i][0];
 324         }
 325         if (im) {
 326             for (int i = 0; i <= hs; ++i) {
 327                 im[i] = m_fpacked[i][1];
 328             }
 329         }
 330     }
 331
 332     void unpackDouble(double *R__ re, double *R__ im) {
 333         const int hs = m_size/2;
 334         for (int i = 0; i <= hs; ++i) {
 335             re[i] = m_dpacked[i][0];
 336         }
 337         if (im) {
 338             for (int i = 0; i <= hs; ++i) {
 339                 im[i] = m_dpacked[i][1];
 340             }
 341         }
 342     }
 343
 344     void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
 345         if (!m_dplanf) initDouble();
 346         const int sz = m_size;
 347         fft_double_type *const R__ dbuf = m_dbuf;
 348 #ifndef FFTW_FLOAT_ONLY
 349         if (realIn != dbuf)
 350 #endif
 351             for (int i = 0; i < sz; ++i) {
 352                 dbuf[i] = realIn[i];
 353             }
 354         fftw_execute(m_dplanf);
 355         unpackDouble(realOut, imagOut);
 356     }
 357
 358     void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
 359         if (!m_dplanf) initDouble();
 360         fft_double_type *const R__ dbuf = m_dbuf;
 361         const int sz = m_size;
 362 #ifndef FFTW_FLOAT_ONLY
 363         if (realIn != dbuf)
 364 #endif
 365             for (int i = 0; i < sz; ++i) {
 366                 dbuf[i] = realIn[i];
 367             }
 368         fftw_execute(m_dplanf);
 369         const int hs = m_size/2;
 370         for (int i = 0; i <= hs; ++i) {
 371             magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
 372                              m_dpacked[i][1] * m_dpacked[i][1]);
 373         }
 374         for (int i = 0; i <= hs; ++i) {
 375             phaseOut[i] = atan2(m_dpacked[i][1], m_dpacked[i][0]);
 376         }
 377     }
 378
 379     void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
 380         if (!m_dplanf) initDouble();
 381         fft_double_type *const R__ dbuf = m_dbuf;
 382         const int sz = m_size;
 383 #ifndef FFTW_FLOAT_ONLY
 384         if (realIn != m_dbuf)
 385 #endif
 386             for (int i = 0; i < sz; ++i) {
 387                 dbuf[i] = realIn[i];
 388             }
 389         fftw_execute(m_dplanf);
 390         const int hs = m_size/2;
 391         for (int i = 0; i <= hs; ++i) {
 392             magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
 393                              m_dpacked[i][1] * m_dpacked[i][1]);
 394         }
 395     }
 396
 397     void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
 398         if (!m_fplanf) initFloat();
 399         fft_float_type *const R__ fbuf = m_fbuf;
 400         const int sz = m_size;
 401 #ifndef FFTW_DOUBLE_ONLY
 402         if (realIn != fbuf)
 403 #endif
 404             for (int i = 0; i < sz; ++i) {
 405                 fbuf[i] = realIn[i];
 406             }
 407         fftwf_execute(m_fplanf);
 408         unpackFloat(realOut, imagOut);
 409     }
 410
 411     void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
 412         if (!m_fplanf) initFloat();
 413         fft_float_type *const R__ fbuf = m_fbuf;
 414         const int sz = m_size;
 415 #ifndef FFTW_DOUBLE_ONLY
 416         if (realIn != fbuf)
 417 #endif
 418             for (int i = 0; i < sz; ++i) {
 419                 fbuf[i] = realIn[i];
 420             }
 421         fftwf_execute(m_fplanf);
 422         const int hs = m_size/2;
 423         for (int i = 0; i <= hs; ++i) {
 424             magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
 425                               m_fpacked[i][1] * m_fpacked[i][1]);
 426         }
 427         for (int i = 0; i <= hs; ++i) {
 428             phaseOut[i] = atan2f(m_fpacked[i][1], m_fpacked[i][0]) ;
 429         }
 430     }
 431
 432     void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
 433         if (!m_fplanf) initFloat();
 434         fft_float_type *const R__ fbuf = m_fbuf;
 435         const int sz = m_size;
 436 #ifndef FFTW_DOUBLE_ONLY
 437         if (realIn != fbuf)
 438 #endif
 439             for (int i = 0; i < sz; ++i) {
 440                 fbuf[i] = realIn[i];
 441             }
 442         fftwf_execute(m_fplanf);
 443         const int hs = m_size/2;
 444         for (int i = 0; i <= hs; ++i) {
 445             magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
 446                               m_fpacked[i][1] * m_fpacked[i][1]);
 447         }
 448     }
 449
 450     void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
 451         if (!m_dplanf) initDouble();
 452         packDouble(realIn, imagIn);
 453         fftw_execute(m_dplani);
 454         const int sz = m_size;
 455         fft_double_type *const R__ dbuf = m_dbuf;
 456 #ifndef FFTW_FLOAT_ONLY
 457         if (realOut != dbuf)
 458 #endif
 459             for (int i = 0; i < sz; ++i) {
 460                 realOut[i] = dbuf[i];
 461             }
 462     }
 463
 464     void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
 465         if (!m_dplanf) initDouble();
 466         const int hs = m_size/2;
 467         fftw_complex *const R__ dpacked = m_dpacked;
 468         for (int i = 0; i <= hs; ++i) {
 469             dpacked[i][0] = magIn[i] * cos(phaseIn[i]);
 470         }
 471         for (int i = 0; i <= hs; ++i) {
 472             dpacked[i][1] = magIn[i] * sin(phaseIn[i]);
 473         }
 474         fftw_execute(m_dplani);
 475         const int sz = m_size;
 476         fft_double_type *const R__ dbuf = m_dbuf;
 477 #ifndef FFTW_FLOAT_ONLY
 478         if (realOut != dbuf)
 479 #endif
 480             for (int i = 0; i < sz; ++i) {
 481                 realOut[i] = dbuf[i];
 482             }
 483     }
 484
 485     void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
 486         if (!m_dplanf) initDouble();
 487         fft_double_type *const R__ dbuf = m_dbuf;
 488         fftw_complex *const R__ dpacked = m_dpacked;
 489         const int hs = m_size/2;
 490         for (int i = 0; i <= hs; ++i) {
 491             dpacked[i][0] = log(magIn[i] + 0.000001);
 492         }
 493         for (int i = 0; i <= hs; ++i) {
 494             dpacked[i][1] = 0.0;
 495         }
 496         fftw_execute(m_dplani);
 497         const int sz = m_size;
 498 #ifndef FFTW_FLOAT_ONLY
 499         if (cepOut != dbuf)
 500 #endif
 501             for (int i = 0; i < sz; ++i) {
 502                 cepOut[i] = dbuf[i];
 503             }
 504     }
 505
 506     void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
 507         if (!m_fplanf) initFloat();
 508         packFloat(realIn, imagIn);
 509         fftwf_execute(m_fplani);
 510         const int sz = m_size;
 511         fft_float_type *const R__ fbuf = m_fbuf;
 512 #ifndef FFTW_DOUBLE_ONLY
 513         if (realOut != fbuf)
 514 #endif
 515             for (int i = 0; i < sz; ++i) {
 516                 realOut[i] = fbuf[i];
 517             }
 518     }
 519
 520     void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
 521         if (!m_fplanf) initFloat();
 522         const int hs = m_size/2;
 523         fftwf_complex *const R__ fpacked = m_fpacked;
 524         for (int i = 0; i <= hs; ++i) {
 525             fpacked[i][0] = magIn[i] * cosf(phaseIn[i]);
 526         }
 527         for (int i = 0; i <= hs; ++i) {
 528             fpacked[i][1] = magIn[i] * sinf(phaseIn[i]);
 529         }
 530         fftwf_execute(m_fplani);
 531         const int sz = m_size;
 532         fft_float_type *const R__ fbuf = m_fbuf;
 533 #ifndef FFTW_DOUBLE_ONLY
 534         if (realOut != fbuf)
 535 #endif
 536             for (int i = 0; i < sz; ++i) {
 537                 realOut[i] = fbuf[i];
 538             }
 539     }
 540
 541     void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
 542         if (!m_fplanf) initFloat();
 543         const int hs = m_size/2;
 544         fftwf_complex *const R__ fpacked = m_fpacked;
 545         for (int i = 0; i <= hs; ++i) {
 546             fpacked[i][0] = logf(magIn[i] + 0.000001f);
 547         }
 548         for (int i = 0; i <= hs; ++i) {
 549             fpacked[i][1] = 0.f;
 550         }
 551         fftwf_execute(m_fplani);
 552         const int sz = m_size;
 553         fft_float_type *const R__ fbuf = m_fbuf;
 554 #ifndef FFTW_DOUBLE_ONLY
 555         if (cepOut != fbuf)
 556 #endif
 557             for (int i = 0; i < sz; ++i) {
 558                 cepOut[i] = fbuf[i];
 559             }
 560     }
 561
 562     float *getFloatTimeBuffer() {
 563         initFloat();
 564 #ifdef FFTW_DOUBLE_ONLY
 565         if (!m_frb) m_frb = (float *)fftw_malloc(m_size * sizeof(float));
 566         return m_frb;
 567 #else
 568         return m_fbuf;
 569 #endif
 570     }
 571
 572     double *getDoubleTimeBuffer() {
 573         initDouble();
 574 #ifdef FFTW_FLOAT_ONLY
 575         if (!m_drb) m_drb = (double *)fftwf_malloc(m_size * sizeof(double));
 576         return m_drb;
 577 #else
 578         return m_dbuf;
 579 #endif
 580     }
 581
 582 private:
 583     fftwf_plan m_fplanf;
 584     fftwf_plan m_fplani;
 585 #ifdef FFTW_DOUBLE_ONLY
 586     float *m_frb;
 587     double *m_fbuf;
 588 #else
 589     float *m_fbuf;
 590 #endif
 591     fftwf_complex *m_fpacked;
 592     fftw_plan m_dplanf;
 593     fftw_plan m_dplani;
 594 #ifdef FFTW_FLOAT_ONLY
 595     float *m_dbuf;
 596     double *m_drb;
 597 #else
 598     double *m_dbuf;
 599 #endif
 600     fftw_complex * m_dpacked;
 601     const int m_size;
 602     static int m_extantf;
 603     static int m_extantd;
 604     static Mutex m_extantMutex;
 605 };
 606
 607 int
 608 D_FFTW::m_extantf = 0;
 609
 610 int
 611 D_FFTW::m_extantd = 0;
 612
 613 Mutex
 614 D_FFTW::m_extantMutex;
 615
 616 #endif /* HAVE_FFTW3 */
 617
 618 #ifdef USE_KISSFFT
 619
 620 class D_KISSFFT : public FFTImpl
 621 {
 622 public:
 623     D_KISSFFT(int size) :
 624         m_size(size),
 625         m_frb(0),
 626         m_drb(0),
 627         m_fplanf(0),
 628         m_fplani(0)
 629     {
 630 #ifdef FIXED_POINT
 631 #error KISSFFT is not configured for float values
 632 #endif
 633         if (sizeof(kiss_fft_scalar) != sizeof(float)) {
 634             std::cerr << "ERROR: KISSFFT is not configured for float values"
 635                       << std::endl;
 636         }
 637
 638         m_fbuf = new kiss_fft_scalar[m_size + 2];
 639         m_fpacked = new kiss_fft_cpx[m_size + 2];
 640         m_fplanf = kiss_fftr_alloc(m_size, 0, NULL, NULL);
 641         m_fplani = kiss_fftr_alloc(m_size, 1, NULL, NULL);
 642     }
 643
 644     ~D_KISSFFT() {
 645         kiss_fftr_free(m_fplanf);
 646         kiss_fftr_free(m_fplani);
 647         kiss_fft_cleanup();
 648
 649         delete[] m_fbuf;
 650         delete[] m_fpacked;
 651
 652         if (m_frb) delete[] m_frb;
 653         if (m_drb) delete[] m_drb;
 654     }
 655
 656     void initFloat() { }
 657     void initDouble() { }
 658
 659     void packFloat(const float *R__ re, const float *R__ im) {
 660         const int hs = m_size/2;
 661         for (int i = 0; i <= hs; ++i) {
 662             m_fpacked[i].r = re[i];
 663             m_fpacked[i].i = im[i];
 664         }
 665     }
 666
 667     void unpackFloat(float *R__ re, float *R__ im) {
 668         const int hs = m_size/2;
 669         for (int i = 0; i <= hs; ++i) {
 670             re[i] = m_fpacked[i].r;
 671             im[i] = m_fpacked[i].i;
 672         }
 673     }
 674
 675     void packDouble(const double *R__ re, const double *R__ im) {
 676         const int hs = m_size/2;
 677         for (int i = 0; i <= hs; ++i) {
 678             m_fpacked[i].r = float(re[i]);
 679             m_fpacked[i].i = float(im[i]);
 680         }
 681     }
 682
 683     void unpackDouble(double *R__ re, double *R__ im) {
 684         const int hs = m_size/2;
 685         for (int i = 0; i <= hs; ++i) {
 686             re[i] = double(m_fpacked[i].r);
 687             im[i] = double(m_fpacked[i].i);
 688         }
 689     }
 690
 691     void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
 692
 693         for (int i = 0; i < m_size; ++i) {
 694             m_fbuf[i] = float(realIn[i]);
 695         }
 696
 697         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
 698         unpackDouble(realOut, imagOut);
 699     }
 700
 701     void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
 702
 703         for (int i = 0; i < m_size; ++i) {
 704             m_fbuf[i] = float(realIn[i]);
 705         }
 706
 707         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
 708
 709         const int hs = m_size/2;
 710
 711         for (int i = 0; i <= hs; ++i) {
 712             magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
 713                              double(m_fpacked[i].i) * double(m_fpacked[i].i));
 714         }
 715
 716         for (int i = 0; i <= hs; ++i) {
 717             phaseOut[i] = atan2(double(m_fpacked[i].i), double(m_fpacked[i].r));
 718         }
 719     }
 720
 721     void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
 722
 723         for (int i = 0; i < m_size; ++i) {
 724             m_fbuf[i] = float(realIn[i]);
 725         }
 726
 727         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
 728
 729         const int hs = m_size/2;
 730
 731         for (int i = 0; i <= hs; ++i) {
 732             magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
 733                              double(m_fpacked[i].i) * double(m_fpacked[i].i));
 734         }
 735     }
 736
 737     void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
 738
 739         kiss_fftr(m_fplanf, realIn, m_fpacked);
 740         unpackFloat(realOut, imagOut);
 741     }
 742
 743     void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
 744
 745         kiss_fftr(m_fplanf, realIn, m_fpacked);
 746
 747         const int hs = m_size/2;
 748
 749         for (int i = 0; i <= hs; ++i) {
 750             magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
 751                               m_fpacked[i].i * m_fpacked[i].i);
 752         }
 753
 754         for (int i = 0; i <= hs; ++i) {
 755             phaseOut[i] = atan2f(m_fpacked[i].i, m_fpacked[i].r);
 756         }
 757     }
 758
 759     void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
 760
 761         kiss_fftr(m_fplanf, realIn, m_fpacked);
 762
 763         const int hs = m_size/2;
 764
 765         for (int i = 0; i <= hs; ++i) {
 766             magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
 767                               m_fpacked[i].i * m_fpacked[i].i);
 768         }
 769     }
 770
 771     void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
 772
 773         packDouble(realIn, imagIn);
 774
 775         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
 776
 777         for (int i = 0; i < m_size; ++i) {
 778             realOut[i] = m_fbuf[i];
 779         }
 780     }
 781
 782     void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
 783
 784         const int hs = m_size/2;
 785
 786         for (int i = 0; i <= hs; ++i) {
 787             m_fpacked[i].r = float(magIn[i] * cos(phaseIn[i]));
 788             m_fpacked[i].i = float(magIn[i] * sin(phaseIn[i]));
 789         }
 790
 791         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
 792
 793         for (int i = 0; i < m_size; ++i) {
 794             realOut[i] = m_fbuf[i];
 795         }
 796     }
 797
 798     void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
 799
 800         const int hs = m_size/2;
 801
 802         for (int i = 0; i <= hs; ++i) {
 803             m_fpacked[i].r = float(log(magIn[i] + 0.000001));
 804             m_fpacked[i].i = 0.0f;
 805         }
 806
 807         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
 808
 809         for (int i = 0; i < m_size; ++i) {
 810             cepOut[i] = m_fbuf[i];
 811         }
 812     }
 813
 814     void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
 815
 816         packFloat(realIn, imagIn);
 817         kiss_fftri(m_fplani, m_fpacked, realOut);
 818     }
 819
 820     void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
 821
 822         const int hs = m_size/2;
 823
 824         for (int i = 0; i <= hs; ++i) {
 825             m_fpacked[i].r = magIn[i] * cosf(phaseIn[i]);
 826             m_fpacked[i].i = magIn[i] * sinf(phaseIn[i]);
 827         }
 828
 829         kiss_fftri(m_fplani, m_fpacked, realOut);
 830     }
 831
 832     void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
 833
 834         const int hs = m_size/2;
 835
 836         for (int i = 0; i <= hs; ++i) {
 837             m_fpacked[i].r = logf(magIn[i] + 0.000001f);
 838             m_fpacked[i].i = 0.0f;
 839         }
 840
 841         kiss_fftri(m_fplani, m_fpacked, cepOut);
 842     }
 843
 844     float *getFloatTimeBuffer() {
 845         if (!m_frb) m_frb = new float[m_size];
 846         return m_frb;
 847     }
 848
 849     double *getDoubleTimeBuffer() {
 850         if (!m_drb) m_drb = new double[m_size];
 851         return m_drb;
 852     }
 853
 854 private:
 855     const int m_size;
 856     float* m_frb;
 857     double* m_drb;
 858     kiss_fftr_cfg m_fplanf;
 859     kiss_fftr_cfg m_fplani;
 860     kiss_fft_scalar *m_fbuf;
 861     kiss_fft_cpx *m_fpacked;
 862 };
 863
 864 #endif /* USE_KISSFFT */
 865
 866 #ifdef USE_BUILTIN_FFT
 867
 868 class D_Cross : public FFTImpl
 869 {
 870 public:
 871     D_Cross(int size) : m_size(size), m_table(0), m_frb(0), m_drb(0) {
 872
 873         m_a = new double[size];
 874         m_b = new double[size];
 875         m_c = new double[size];
 876         m_d = new double[size];
 877
 878         m_table = new int[m_size];
 879
 880         int bits;
 881         int i, j, k, m;
 882
 883         for (i = 0; ; ++i) {
 884             if (m_size & (1 << i)) {
 885                 bits = i;
 886                 break;
 887             }
 888         }
 889
 890         for (i = 0; i < m_size; ++i) {
 891
 892             m = i;
 893
 894             for (j = k = 0; j < bits; ++j) {
 895                 k = (k << 1) | (m & 1);
 896                 m >>= 1;
 897             }
 898
 899             m_table[i] = k;
 900         }
 901     }
 902
 903     ~D_Cross() {
 904         delete[] m_table;
 905         delete[] m_a;
 906         delete[] m_b;
 907         delete[] m_c;
 908         delete[] m_d;
 909         delete[] m_frb;
 910         delete[] m_drb;
 911     }
 912
 913     void initFloat() { }
 914     void initDouble() { }
 915
 916     void forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut) {
 917         basefft(false, realIn, 0, m_c, m_d);
 918         const int hs = m_size/2;
 919         for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
 920         if (imagOut) {
 921             for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
 922         }
 923     }
 924
 925     void forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut) {
 926         basefft(false, realIn, 0, m_c, m_d);
 927         const int hs = m_size/2;
 928         for (int i = 0; i <= hs; ++i) {
 929             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
 930             phaseOut[i] = atan2(m_d[i], m_c[i]) ;
 931         }
 932     }
 933
 934     void forwardMagnitude(const double *R__ realIn, double *R__ magOut) {
 935         basefft(false, realIn, 0, m_c, m_d);
 936         const int hs = m_size/2;
 937         for (int i = 0; i <= hs; ++i) {
 938             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
 939         }
 940     }
 941
 942     void forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut) {
 943         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
 944         basefft(false, m_a, 0, m_c, m_d);
 945         const int hs = m_size/2;
 946         for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
 947         if (imagOut) {
 948             for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
 949         }
 950     }
 951
 952     void forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut) {
 953         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
 954         basefft(false, m_a, 0, m_c, m_d);
 955         const int hs = m_size/2;
 956         for (int i = 0; i <= hs; ++i) {
 957             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
 958             phaseOut[i] = atan2(m_d[i], m_c[i]) ;
 959         }
 960     }
 961
 962     void forwardMagnitude(const float *R__ realIn, float *R__ magOut) {
 963         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
 964         basefft(false, m_a, 0, m_c, m_d);
 965         const int hs = m_size/2;
 966         for (int i = 0; i <= hs; ++i) {
 967             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
 968         }
 969     }
 970
 971     void inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut) {
 972         const int hs = m_size/2;
 973         for (int i = 0; i <= hs; ++i) {
 974             double real = realIn[i];
 975             double imag = imagIn[i];
 976             m_a[i] = real;
 977             m_b[i] = imag;
 978             if (i > 0) {
 979                 m_a[m_size-i] = real;
 980                 m_b[m_size-i] = -imag;
 981             }
 982         }
 983         basefft(true, m_a, m_b, realOut, m_d);
 984     }
 985
 986     void inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut) {
 987         const int hs = m_size/2;
 988         for (int i = 0; i <= hs; ++i) {
 989             double real = magIn[i] * cos(phaseIn[i]);
 990             double imag = magIn[i] * sin(phaseIn[i]);
 991             m_a[i] = real;
 992             m_b[i] = imag;
 993             if (i > 0) {
 994                 m_a[m_size-i] = real;
 995                 m_b[m_size-i] = -imag;
 996             }
 997         }
 998         basefft(true, m_a, m_b, realOut, m_d);
 999     }
1000
1001     void inverseCepstral(const double *R__ magIn, double *R__ cepOut) {
1002         const int hs = m_size/2;
1003         for (int i = 0; i <= hs; ++i) {
1004             double real = log(magIn[i] + 0.000001);
1005             m_a[i] = real;
1006             m_b[i] = 0.0;
1007             if (i > 0) {
1008                 m_a[m_size-i] = real;
1009                 m_b[m_size-i] = 0.0;
1010             }
1011         }
1012         basefft(true, m_a, m_b, cepOut, m_d);
1013     }
1014
1015     void inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut) {
1016         const int hs = m_size/2;
1017         for (int i = 0; i <= hs; ++i) {
1018             float real = realIn[i];
1019             float imag = imagIn[i];
1020             m_a[i] = real;
1021             m_b[i] = imag;
1022             if (i > 0) {
1023                 m_a[m_size-i] = real;
1024                 m_b[m_size-i] = -imag;
1025             }
1026         }
1027         basefft(true, m_a, m_b, m_c, m_d);
1028         for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
1029     }
1030
1031     void inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut) {
1032         const int hs = m_size/2;
1033         for (int i = 0; i <= hs; ++i) {
1034             float real = magIn[i] * cosf(phaseIn[i]);
1035             float imag = magIn[i] * sinf(phaseIn[i]);
1036             m_a[i] = real;
1037             m_b[i] = imag;
1038             if (i > 0) {
1039                 m_a[m_size-i] = real;
1040                 m_b[m_size-i] = -imag;
1041             }
1042         }
1043         basefft(true, m_a, m_b, m_c, m_d);
1044         for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
1045     }
1046
1047     void inverseCepstral(const float *R__ magIn, float *R__ cepOut) {
1048         const int hs = m_size/2;
1049         for (int i = 0; i <= hs; ++i) {
1050             float real = logf(magIn[i] + 0.000001);
1051             m_a[i] = real;
1052             m_b[i] = 0.0;
1053             if (i > 0) {
1054                 m_a[m_size-i] = real;
1055                 m_b[m_size-i] = 0.0;
1056             }
1057         }
1058         basefft(true, m_a, m_b, m_c, m_d);
1059         for (int i = 0; i < m_size; ++i) cepOut[i] = m_c[i];
1060     }
1061
1062     float *getFloatTimeBuffer() {
1063         if (!m_frb) m_frb = new float[m_size];
1064         return m_frb;
1065     }
1066
1067     double *getDoubleTimeBuffer() {
1068         if (!m_drb) m_drb = new double[m_size];
1069         return m_drb;
1070     }
1071
1072 private:
1073     const int m_size;
1074     int *m_table;
1075     float *m_frb;
1076     double *m_drb;
1077     double *m_a;
1078     double *m_b;
1079     double *m_c;
1080     double *m_d;
1081     void basefft(bool inverse, const double *R__ ri, const double *R__ ii, double *R__ ro, double *R__ io);
1082 };
1083
1084 void
1085 D_Cross::basefft(bool inverse, const double *R__ ri, const double *R__ ii, double *R__ ro, double *R__ io)
1086 {
1087     if (!ri || !ro || !io) return;
1088
1089     int i, j, k, m;
1090     int blockSize, blockEnd;
1091
1092     double tr, ti;
1093
1094     double angle = 2.0 * M_PI;
1095     if (inverse) angle = -angle;
1096
1097     const int n = m_size;
1098
1099     if (ii) {
1100         for (i = 0; i < n; ++i) {
1101             ro[m_table[i]] = ri[i];
1102         }
1103         for (i = 0; i < n; ++i) {
1104             io[m_table[i]] = ii[i];
1105         }
1106     } else {
1107         for (i = 0; i < n; ++i) {
1108             ro[m_table[i]] = ri[i];
1109         }
1110         for (i = 0; i < n; ++i) {
1111             io[m_table[i]] = 0.0;
1112         }
1113     }
1114
1115     blockEnd = 1;
1116
1117     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
1118
1119         double delta = angle / (double)blockSize;
1120         double sm2 = -sin(-2 * delta);
1121         double sm1 = -sin(-delta);
1122         double cm2 = cos(-2 * delta);
1123         double cm1 = cos(-delta);
1124         double w = 2 * cm1;
1125         double ar[3], ai[3];
1126
1127         for (i = 0; i < n; i += blockSize) {
1128
1129             ar[2] = cm2;
1130             ar[1] = cm1;
1131
1132             ai[2] = sm2;
1133             ai[1] = sm1;
1134
1135             for (j = i, m = 0; m < blockEnd; j++, m++) {
1136
1137                 ar[0] = w * ar[1] - ar[2];
1138                 ar[2] = ar[1];
1139                 ar[1] = ar[0];
1140
1141                 ai[0] = w * ai[1] - ai[2];
1142                 ai[2] = ai[1];
1143                 ai[1] = ai[0];
1144
1145                 k = j + blockEnd;
1146                 tr = ar[0] * ro[k] - ai[0] * io[k];
1147                 ti = ar[0] * io[k] + ai[0] * ro[k];
1148
1149                 ro[k] = ro[j] - tr;
1150                 io[k] = io[j] - ti;
1151
1152                 ro[j] += tr;
1153                 io[j] += ti;
1154             }
1155         }
1156
1157         blockEnd = blockSize;
1158     }
1159
1160 /* fftw doesn't rescale, so nor will we
1161
1162     if (inverse) {
1163
1164         double denom = (double)n;
1165
1166         for (i = 0; i < n; i++) {
1167             ro[i] /= denom;
1168             io[i] /= denom;
1169         }
1170     }
1171 */
1172 }
1173
1174 #endif /* USE_BUILTIN_FFT */
1175
1176 } /* end namespace FFTs */
1177
1178 int
1179 FFT::m_method = -1;
1180
1181 FFT::FFT(int size, int debugLevel)
1182 {
1183     if ((size < 2) ||
1184         (size & (size-1))) {
1185         std::cerr << "FFT::FFT(" << size << "): power-of-two sizes only supported, minimum size 2" << std::endl;
1186         throw InvalidSize;
1187     }
1188
1189     if (m_method == -1) {
1190         m_method = 3;
1191 #ifdef USE_KISSFFT
1192         m_method = 2;
1193 #endif
1194 #ifdef HAVE_FFTW3
1195         m_method = 1;
1196 #endif
1197     }
1198
1199     switch (m_method) {
1200
1201     case 0:
1202         std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1203 #ifdef USE_BUILTIN_FFT
1204         d = new FFTs::D_Cross(size);
1205 #else
1206         std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1207         abort();
1208 #endif
1209         break;
1210
1211     case 1:
1212 #ifdef HAVE_FFTW3
1213         if (debugLevel > 0) {
1214             std::cerr << "FFT::FFT(" << size << "): using FFTW3 implementation"
1215                       << std::endl;
1216         }
1217         d = new FFTs::D_FFTW(size);
1218 #else
1219         std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1220 #ifdef USE_BUILTIN_FFT
1221         d = new FFTs::D_Cross(size);
1222 #else
1223         std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1224         abort();
1225 #endif
1226 #endif
1227         break;
1228
1229     case 2:
1230 #ifdef USE_KISSFFT
1231         if (debugLevel > 0) {
1232             std::cerr << "FFT::FFT(" << size << "): using KISSFFT implementation"
1233                       << std::endl;
1234         }
1235         d = new FFTs::D_KISSFFT(size);
1236 #else
1237         std::cerr << "FFT::FFT(" << size << "): WARNING: Selected implemention not available" << std::endl;
1238 #ifdef USE_BUILTIN_FFT
1239         d = new FFTs::D_Cross(size);
1240 #else
1241         std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1242         abort();
1243 #endif
1244 #endif
1245         break;
1246
1247     default:
1248 #ifdef USE_BUILTIN_FFT
1249         std::cerr << "FFT::FFT(" << size << "): WARNING: using slow built-in implementation" << std::endl;
1250         d = new FFTs::D_Cross(size);
1251 #else
1252         std::cerr << "FFT::FFT(" << size << "): ERROR: Fallback implementation not available!" << std::endl;
1253         abort();
1254 #endif
1255         break;
1256     }
1257 }
1258
1259 FFT::~FFT()
1260 {
1261     delete d;
1262 }
1263
1264 void
1265 FFT::forward(const double *R__ realIn, double *R__ realOut, double *R__ imagOut)
1266 {
1267     d->forward(realIn, realOut, imagOut);
1268 }
1269
1270 void
1271 FFT::forwardPolar(const double *R__ realIn, double *R__ magOut, double *R__ phaseOut)
1272 {
1273     d->forwardPolar(realIn, magOut, phaseOut);
1274 }
1275
1276 void
1277 FFT::forwardMagnitude(const double *R__ realIn, double *R__ magOut)
1278 {
1279     d->forwardMagnitude(realIn, magOut);
1280 }
1281
1282 void
1283 FFT::forward(const float *R__ realIn, float *R__ realOut, float *R__ imagOut)
1284 {
1285     d->forward(realIn, realOut, imagOut);
1286 }
1287
1288 void
1289 FFT::forwardPolar(const float *R__ realIn, float *R__ magOut, float *R__ phaseOut)
1290 {
1291     d->forwardPolar(realIn, magOut, phaseOut);
1292 }
1293
1294 void
1295 FFT::forwardMagnitude(const float *R__ realIn, float *R__ magOut)
1296 {
1297     d->forwardMagnitude(realIn, magOut);
1298 }
1299
1300 void
1301 FFT::inverse(const double *R__ realIn, const double *R__ imagIn, double *R__ realOut)
1302 {
1303     d->inverse(realIn, imagIn, realOut);
1304 }
1305
1306 void
1307 FFT::inversePolar(const double *R__ magIn, const double *R__ phaseIn, double *R__ realOut)
1308 {
1309     d->inversePolar(magIn, phaseIn, realOut);
1310 }
1311
1312 void
1313 FFT::inverseCepstral(const double *R__ magIn, double *R__ cepOut)
1314 {
1315     d->inverseCepstral(magIn, cepOut);
1316 }
1317
1318 void
1319 FFT::inverse(const float *R__ realIn, const float *R__ imagIn, float *R__ realOut)
1320 {
1321     d->inverse(realIn, imagIn, realOut);
1322 }
1323
1324 void
1325 FFT::inversePolar(const float *R__ magIn, const float *R__ phaseIn, float *R__ realOut)
1326 {
1327     d->inversePolar(magIn, phaseIn, realOut);
1328 }
1329
1330 void
1331 FFT::inverseCepstral(const float *R__ magIn, float *R__ cepOut)
1332 {
1333     d->inverseCepstral(magIn, cepOut);
1334 }
1335
1336 void
1337 FFT::initFloat()
1338 {
1339     d->initFloat();
1340 }
1341
1342 void
1343 FFT::initDouble()
1344 {
1345     d->initDouble();
1346 }
1347
1348 float *
1349 FFT::getFloatTimeBuffer()
1350 {
1351     return d->getFloatTimeBuffer();
1352 }
1353
1354 double *
1355 FFT::getDoubleTimeBuffer()
1356 {
1357     return d->getDoubleTimeBuffer();
1358 }
1359
1360
1361 void
1362 FFT::tune()
1363 {
1364 }
1365
1366
1367 }