Chris@29: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@29: 
Chris@29: /*
Chris@29:     bqfft
Chris@29: 
Chris@29:     A small library wrapping various FFT implementations for some
Chris@29:     common audio processing use cases.
Chris@29: 
Chris@29:     Copyright 2007-2015 Particular Programs Ltd.
Chris@29: 
Chris@29:     Permission is hereby granted, free of charge, to any person
Chris@29:     obtaining a copy of this software and associated documentation
Chris@29:     files (the "Software"), to deal in the Software without
Chris@29:     restriction, including without limitation the rights to use, copy,
Chris@29:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@29:     of the Software, and to permit persons to whom the Software is
Chris@29:     furnished to do so, subject to the following conditions:
Chris@29: 
Chris@29:     The above copyright notice and this permission notice shall be
Chris@29:     included in all copies or substantial portions of the Software.
Chris@29: 
Chris@29:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@29:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@29:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@29:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
Chris@29:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@29:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@29:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@29: 
Chris@29:     Except as contained in this notice, the names of Chris Cannam and
Chris@29:     Particular Programs Ltd shall not be used in advertising or
Chris@29:     otherwise to promote the sale, use or other dealings in this
Chris@29:     Software without prior written authorization.
Chris@29: */
Chris@29: 
Chris@29: #include "bqfft/FFT.h"
Chris@29: 
Chris@29: //#include "system/Thread.h"
Chris@29: 
Chris@29: #include <bqvec/Allocators.h>
Chris@29: #include <bqvec/VectorOps.h>
Chris@29: #include <bqvec/VectorOpsComplex.h>
Chris@29: 
Chris@29: //#define FFT_MEASUREMENT 1
Chris@29: 
Chris@29: #ifdef FFT_MEASUREMENT
Chris@29: #include <sstream>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_IPP
Chris@29: #include <ipps.h>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_FFTW3
Chris@29: #include <fftw3.h>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_VDSP
Chris@29: #include <Accelerate/Accelerate.h>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_MEDIALIB
Chris@29: #include <mlib_signal.h>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_OPENMAX
Chris@29: #include <omxSP.h>
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_SFFT
Chris@29: extern "C" {
Chris@29: #include <sfft.h>
Chris@29: }
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_KISSFFT
Chris@29: #include "kissfft/kiss_fftr.h"
Chris@29: #endif
Chris@29: 
Chris@29: #ifndef HAVE_IPP
Chris@29: #ifndef HAVE_FFTW3
Chris@29: #ifndef HAVE_KISSFFT
Chris@29: #ifndef USE_BUILTIN_FFT
Chris@29: #ifndef HAVE_VDSP
Chris@29: #ifndef HAVE_MEDIALIB
Chris@29: #ifndef HAVE_OPENMAX
Chris@29: #ifndef HAVE_SFFT
Chris@29: #error No FFT implementation selected!
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: #endif
Chris@29: 
Chris@29: #include <cmath>
Chris@29: #include <iostream>
Chris@29: #include <map>
Chris@29: #include <cstdio>
Chris@29: #include <cstdlib>
Chris@29: #include <vector>
Chris@29: 
Chris@29: #ifdef FFT_MEASUREMENT
Chris@29: #ifndef _WIN32
Chris@29: #include <unistd.h>
Chris@29: #endif
Chris@29: #endif
Chris@29: 
Chris@29: namespace breakfastquay {
Chris@29: 
Chris@29: class FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     virtual ~FFTImpl() { }
Chris@29: 
Chris@29:     virtual FFT::Precisions getSupportedPrecisions() const = 0;
Chris@29: 
Chris@29:     virtual void initFloat() = 0;
Chris@29:     virtual void initDouble() = 0;
Chris@29: 
Chris@29:     virtual void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) = 0;
Chris@29:     virtual void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) = 0;
Chris@29:     virtual void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) = 0;
Chris@29:     virtual void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) = 0;
Chris@29: 
Chris@29:     virtual void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) = 0;
Chris@29:     virtual void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) = 0;
Chris@29:     virtual void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) = 0;
Chris@29:     virtual void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) = 0;
Chris@29: 
Chris@29:     virtual void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) = 0;
Chris@29:     virtual void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) = 0;
Chris@29:     virtual void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) = 0;
Chris@29:     virtual void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) = 0;
Chris@29: 
Chris@29:     virtual void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) = 0;
Chris@29:     virtual void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) = 0;
Chris@29:     virtual void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) = 0;
Chris@29:     virtual void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) = 0;
Chris@29: };    
Chris@29: 
Chris@29: namespace FFTs {
Chris@29: 
Chris@29: #ifdef HAVE_IPP
Chris@29: 
Chris@29: class D_IPP : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_IPP(int size) :
Chris@29:         m_size(size), m_fspec(0), m_dspec(0)
Chris@29:     { 
Chris@29:         for (int i = 0; ; ++i) {
Chris@29:             if (m_size & (1 << i)) {
Chris@29:                 m_order = i;
Chris@29:                 break;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     ~D_IPP() {
Chris@29:         if (m_fspec) {
Chris@29:             ippsFFTFree_R_32f(m_fspec);
Chris@29:             ippsFree(m_fbuf);
Chris@29:             ippsFree(m_fpacked);
Chris@29:             ippsFree(m_fspare);
Chris@29:         }
Chris@29:         if (m_dspec) {
Chris@29:             ippsFFTFree_R_64f(m_dspec);
Chris@29:             ippsFree(m_dbuf);
Chris@29:             ippsFree(m_dpacked);
Chris@29:             ippsFree(m_dspare);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::SinglePrecision | FFT::DoublePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     //!!! rv check
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         if (m_fspec) return;
Chris@29:         int specSize, specBufferSize, bufferSize;
Chris@29:         ippsFFTGetSize_R_32f(m_order, IPP_FFT_NODIV_BY_ANY, ippAlgHintFast,
Chris@29:                              &specSize, &specBufferSize, &bufferSize);
Chris@29:         m_fbuf = ippsMalloc_8u(bufferSize);
Chris@29:         m_fpacked = ippsMalloc_32f(m_size + 2);
Chris@29:         m_fspare = ippsMalloc_32f(m_size / 2 + 1);
Chris@29:         ippsFFTInitAlloc_R_32f(&m_fspec, m_order, IPP_FFT_NODIV_BY_ANY, 
Chris@29:                                ippAlgHintFast);
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         if (m_dspec) return;
Chris@29:         int specSize, specBufferSize, bufferSize;
Chris@29:         ippsFFTGetSize_R_64f(m_order, IPP_FFT_NODIV_BY_ANY, ippAlgHintFast,
Chris@29:                              &specSize, &specBufferSize, &bufferSize);
Chris@29:         m_dbuf = ippsMalloc_8u(bufferSize);
Chris@29:         m_dpacked = ippsMalloc_64f(m_size + 2);
Chris@29:         m_dspare = ippsMalloc_64f(m_size / 2 + 1);
Chris@29:         ippsFFTInitAlloc_R_64f(&m_dspec, m_order, IPP_FFT_NODIV_BY_ANY, 
Chris@29:                                ippAlgHintFast);
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[index++] = re[i];
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_fpacked[index++] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_fpacked[index++] = 0.f;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im) {
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_dpacked[index++] = re[i];
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_dpacked[index++] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_dpacked[index++] = 0.0;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(float *re, float *BQ_R__ im) { // re may be equal to m_fpacked
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = m_fpacked[index++];
Chris@29:             }
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_fpacked[index++];
Chris@29:             index++;
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackDouble(double *re, double *BQ_R__ im) { // re may be equal to m_dpacked
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = m_dpacked[index++];
Chris@29:             }
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_dpacked[index++];
Chris@29:             index++;
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsFFTFwd_RToCCS_64f(realIn, m_dpacked, m_dspec, m_dbuf);
Chris@29:         unpackDouble(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsFFTFwd_RToCCS_64f(realIn, complexOut, m_dspec, m_dbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsFFTFwd_RToCCS_64f(realIn, m_dpacked, m_dspec, m_dbuf);
Chris@29:         unpackDouble(m_dpacked, m_dspare);
Chris@29:         ippsCartToPolar_64f(m_dpacked, m_dspare, magOut, phaseOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsFFTFwd_RToCCS_64f(realIn, m_dpacked, m_dspec, m_dbuf);
Chris@29:         unpackDouble(m_dpacked, m_dspare);
Chris@29:         ippsMagnitude_64f(m_dpacked, m_dspare, magOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsFFTFwd_RToCCS_32f(realIn, m_fpacked, m_fspec, m_fbuf);
Chris@29:         unpackFloat(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsFFTFwd_RToCCS_32f(realIn, complexOut, m_fspec, m_fbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsFFTFwd_RToCCS_32f(realIn, m_fpacked, m_fspec, m_fbuf);
Chris@29:         unpackFloat(m_fpacked, m_fspare);
Chris@29:         ippsCartToPolar_32f(m_fpacked, m_fspare, magOut, phaseOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsFFTFwd_RToCCS_32f(realIn, m_fpacked, m_fspec, m_fbuf);
Chris@29:         unpackFloat(m_fpacked, m_fspare);
Chris@29:         ippsMagnitude_32f(m_fpacked, m_fspare, magOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         packDouble(realIn, imagIn);
Chris@29:         ippsFFTInv_CCSToR_64f(m_dpacked, realOut, m_dspec, m_dbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsFFTInv_CCSToR_64f(complexIn, realOut, m_dspec, m_dbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         ippsPolarToCart_64f(magIn, phaseIn, realOut, m_dspare, m_size/2+1);
Chris@29:         packDouble(realOut, m_dspare); // to m_dpacked
Chris@29:         ippsFFTInv_CCSToR_64f(m_dpacked, realOut, m_dspec, m_dbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         ippsCopy_64f(magIn, m_dspare, hs1);
Chris@29:         ippsAddC_64f_I(0.000001, m_dspare, hs1);
Chris@29:         ippsLn_64f_I(m_dspare, hs1);
Chris@29:         packDouble(m_dspare, 0);
Chris@29:         ippsFFTInv_CCSToR_64f(m_dpacked, cepOut, m_dspec, m_dbuf);
Chris@29:     }
Chris@29:     
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         packFloat(realIn, imagIn);
Chris@29:         ippsFFTInv_CCSToR_32f(m_fpacked, realOut, m_fspec, m_fbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsFFTInv_CCSToR_32f(complexIn, realOut, m_fspec, m_fbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         ippsPolarToCart_32f(magIn, phaseIn, realOut, m_fspare, m_size/2+1);
Chris@29:         packFloat(realOut, m_fspare); // to m_fpacked
Chris@29:         ippsFFTInv_CCSToR_32f(m_fpacked, realOut, m_fspec, m_fbuf);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         ippsCopy_32f(magIn, m_fspare, hs1);
Chris@29:         ippsAddC_32f_I(0.000001f, m_fspare, hs1);
Chris@29:         ippsLn_32f_I(m_fspare, hs1);
Chris@29:         packFloat(m_fspare, 0);
Chris@29:         ippsFFTInv_CCSToR_32f(m_fpacked, cepOut, m_fspec, m_fbuf);
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     int m_order;
Chris@29:     IppsFFTSpec_R_32f *m_fspec;
Chris@29:     IppsFFTSpec_R_64f *m_dspec;
Chris@29:     Ipp8u *m_fbuf;
Chris@29:     Ipp8u *m_dbuf;
Chris@29:     float *m_fpacked;
Chris@29:     float *m_fspare;
Chris@29:     double *m_dpacked;
Chris@29:     double *m_dspare;
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_IPP */
Chris@29: 
Chris@29: #ifdef HAVE_VDSP
Chris@29: 
Chris@29: class D_VDSP : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_VDSP(int size) :
Chris@29:         m_size(size), m_fspec(0), m_dspec(0),
Chris@29:         m_fpacked(0), m_fspare(0),
Chris@29:         m_dpacked(0), m_dspare(0)
Chris@29:     { 
Chris@29:         for (int i = 0; ; ++i) {
Chris@29:             if (m_size & (1 << i)) {
Chris@29:                 m_order = i;
Chris@29:                 break;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     ~D_VDSP() {
Chris@29:         if (m_fspec) {
Chris@29:             vDSP_destroy_fftsetup(m_fspec);
Chris@29:             deallocate(m_fspare);
Chris@29:             deallocate(m_fspare2);
Chris@29:             deallocate(m_fbuf->realp);
Chris@29:             deallocate(m_fbuf->imagp);
Chris@29:             delete m_fbuf;
Chris@29:             deallocate(m_fpacked->realp);
Chris@29:             deallocate(m_fpacked->imagp);
Chris@29:             delete m_fpacked;
Chris@29:         }
Chris@29:         if (m_dspec) {
Chris@29:             vDSP_destroy_fftsetupD(m_dspec);
Chris@29:             deallocate(m_dspare);
Chris@29:             deallocate(m_dspare2);
Chris@29:             deallocate(m_dbuf->realp);
Chris@29:             deallocate(m_dbuf->imagp);
Chris@29:             delete m_dbuf;
Chris@29:             deallocate(m_dpacked->realp);
Chris@29:             deallocate(m_dpacked->imagp);
Chris@29:             delete m_dpacked;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::SinglePrecision | FFT::DoublePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     //!!! rv check
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         if (m_fspec) return;
Chris@29:         m_fspec = vDSP_create_fftsetup(m_order, FFT_RADIX2);
Chris@29:         m_fbuf = new DSPSplitComplex;
Chris@29:         //!!! "If possible, tempBuffer->realp and tempBuffer->imagp should be 32-byte aligned for best performance."
Chris@29:         m_fbuf->realp = allocate<float>(m_size);
Chris@29:         m_fbuf->imagp = allocate<float>(m_size);
Chris@29:         m_fpacked = new DSPSplitComplex;
Chris@29:         m_fpacked->realp = allocate<float>(m_size / 2 + 1);
Chris@29:         m_fpacked->imagp = allocate<float>(m_size / 2 + 1);
Chris@29:         m_fspare = allocate<float>(m_size + 2);
Chris@29:         m_fspare2 = allocate<float>(m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         if (m_dspec) return;
Chris@29:         m_dspec = vDSP_create_fftsetupD(m_order, FFT_RADIX2);
Chris@29:         m_dbuf = new DSPDoubleSplitComplex;
Chris@29:         //!!! "If possible, tempBuffer->realp and tempBuffer->imagp should be 32-byte aligned for best performance."
Chris@29:         m_dbuf->realp = allocate<double>(m_size);
Chris@29:         m_dbuf->imagp = allocate<double>(m_size);
Chris@29:         m_dpacked = new DSPDoubleSplitComplex;
Chris@29:         m_dpacked->realp = allocate<double>(m_size / 2 + 1);
Chris@29:         m_dpacked->imagp = allocate<double>(m_size / 2 + 1);
Chris@29:         m_dspare = allocate<double>(m_size + 2);
Chris@29:         m_dspare2 = allocate<double>(m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void packReal(const float *BQ_R__ const re) {
Chris@29:         // Pack input for forward transform 
Chris@29:         vDSP_ctoz((DSPComplex *)re, 2, m_fpacked, 1, m_size/2);
Chris@29:     }
Chris@29:     void packComplex(const float *BQ_R__ const re, const float *BQ_R__ const im) {
Chris@29:         // Pack input for inverse transform 
Chris@29:         if (re) v_copy(m_fpacked->realp, re, m_size/2 + 1);
Chris@29:         else v_zero(m_fpacked->realp, m_size/2 + 1);
Chris@29:         if (im) v_copy(m_fpacked->imagp, im, m_size/2 + 1);
Chris@29:         else v_zero(m_fpacked->imagp, m_size/2 + 1);
Chris@29:         fnyq();
Chris@29:     }
Chris@29: 
Chris@29:     void unpackReal(float *BQ_R__ const re) {
Chris@29:         // Unpack output for inverse transform
Chris@29:         vDSP_ztoc(m_fpacked, 1, (DSPComplex *)re, 2, m_size/2);
Chris@29:     }
Chris@29:     void unpackComplex(float *BQ_R__ const re, float *BQ_R__ const im) {
Chris@29:         // Unpack output for forward transform
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         float two = 2.f;
Chris@29:         vDSP_vsdiv(m_fpacked->realp, 1, &two, re, 1, m_size/2 + 1);
Chris@29:         vDSP_vsdiv(m_fpacked->imagp, 1, &two, im, 1, m_size/2 + 1);
Chris@29:     }
Chris@29:     void unpackComplex(float *BQ_R__ const cplx) {
Chris@29:         // Unpack output for forward transform
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         for (int i = 0; i < hs1; ++i) {
Chris@29:             cplx[i*2] = m_fpacked->realp[i] / 2.f;
Chris@29:             cplx[i*2+1] = m_fpacked->imagp[i] / 2.f;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packReal(const double *BQ_R__ const re) {
Chris@29:         // Pack input for forward transform
Chris@29:         vDSP_ctozD((DSPDoubleComplex *)re, 2, m_dpacked, 1, m_size/2);
Chris@29:     }
Chris@29:     void packComplex(const double *BQ_R__ const re, const double *BQ_R__ const im) {
Chris@29:         // Pack input for inverse transform
Chris@29:         if (re) v_copy(m_dpacked->realp, re, m_size/2 + 1);
Chris@29:         else v_zero(m_dpacked->realp, m_size/2 + 1);
Chris@29:         if (im) v_copy(m_dpacked->imagp, im, m_size/2 + 1);
Chris@29:         else v_zero(m_dpacked->imagp, m_size/2 + 1);
Chris@29:         dnyq();
Chris@29:     }
Chris@29: 
Chris@29:     void unpackReal(double *BQ_R__ const re) {
Chris@29:         // Unpack output for inverse transform
Chris@29:         vDSP_ztocD(m_dpacked, 1, (DSPDoubleComplex *)re, 2, m_size/2);
Chris@29:     }
Chris@29:     void unpackComplex(double *BQ_R__ const re, double *BQ_R__ const im) {
Chris@29:         // Unpack output for forward transform
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         double two = 2.0;
Chris@29:         vDSP_vsdivD(m_dpacked->realp, 1, &two, re, 1, m_size/2 + 1);
Chris@29:         vDSP_vsdivD(m_dpacked->imagp, 1, &two, im, 1, m_size/2 + 1);
Chris@29:     }
Chris@29:     void unpackComplex(double *BQ_R__ const cplx) {
Chris@29:         // Unpack output for forward transform
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         for (int i = 0; i < hs1; ++i) {
Chris@29:             cplx[i*2] = m_dpacked->realp[i] / 2.0;
Chris@29:             cplx[i*2+1] = m_dpacked->imagp[i] / 2.0;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void fdenyq() {
Chris@29:         // for fft result in packed form, unpack the DC and Nyquist bins
Chris@29:         const int hs = m_size/2;
Chris@29:         m_fpacked->realp[hs] = m_fpacked->imagp[0];
Chris@29:         m_fpacked->imagp[hs] = 0.f;
Chris@29:         m_fpacked->imagp[0] = 0.f;
Chris@29:     }
Chris@29:     void ddenyq() {
Chris@29:         // for fft result in packed form, unpack the DC and Nyquist bins
Chris@29:         const int hs = m_size/2;
Chris@29:         m_dpacked->realp[hs] = m_dpacked->imagp[0];
Chris@29:         m_dpacked->imagp[hs] = 0.;
Chris@29:         m_dpacked->imagp[0] = 0.;
Chris@29:     }
Chris@29: 
Chris@29:     void fnyq() {
Chris@29:         // for ifft input in packed form, pack the DC and Nyquist bins
Chris@29:         const int hs = m_size/2;
Chris@29:         m_fpacked->imagp[0] = m_fpacked->realp[hs];
Chris@29:         m_fpacked->realp[hs] = 0.f;
Chris@29:         m_fpacked->imagp[hs] = 0.f;
Chris@29:     }
Chris@29:     void dnyq() {
Chris@29:         // for ifft input in packed form, pack the DC and Nyquist bins
Chris@29:         const int hs = m_size/2;
Chris@29:         m_dpacked->imagp[0] = m_dpacked->realp[hs];
Chris@29:         m_dpacked->realp[hs] = 0.;
Chris@29:         m_dpacked->imagp[hs] = 0.;
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_FORWARD);
Chris@29:         ddenyq();
Chris@29:         unpackComplex(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_FORWARD);
Chris@29:         ddenyq();
Chris@29:         unpackComplex(complexOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         const int hs1 = m_size/2+1;
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_FORWARD);
Chris@29:         ddenyq();
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         for (int i = 0; i < hs1; ++i) m_dpacked->realp[i] /= 2.0;
Chris@29:         for (int i = 0; i < hs1; ++i) m_dpacked->imagp[i] /= 2.0;
Chris@29:         v_cartesian_to_polar(magOut, phaseOut,
Chris@29:                              m_dpacked->realp, m_dpacked->imagp, hs1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_FORWARD);
Chris@29:         ddenyq();
Chris@29:         const int hs1 = m_size/2+1;
Chris@29:         vDSP_zvmagsD(m_dpacked, 1, m_dspare, 1, hs1);
Chris@29:         vvsqrt(m_dspare2, m_dspare, &hs1);
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         double two = 2.0;
Chris@29:         vDSP_vsdivD(m_dspare2, 1, &two, magOut, 1, hs1);
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_FORWARD);
Chris@29:         fdenyq();
Chris@29:         unpackComplex(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_FORWARD);
Chris@29:         fdenyq();
Chris@29:         unpackComplex(complexOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         const int hs1 = m_size/2+1;
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_FORWARD);
Chris@29:         fdenyq();
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         for (int i = 0; i < hs1; ++i) m_fpacked->realp[i] /= 2.f;
Chris@29:         for (int i = 0; i < hs1; ++i) m_fpacked->imagp[i] /= 2.f;
Chris@29:         v_cartesian_to_polar(magOut, phaseOut,
Chris@29:                              m_fpacked->realp, m_fpacked->imagp, hs1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         packReal(realIn);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_FORWARD);
Chris@29:         fdenyq();
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         vDSP_zvmags(m_fpacked, 1, m_fspare, 1, hs1);
Chris@29:         vvsqrtf(m_fspare2, m_fspare, &hs1);
Chris@29:         // vDSP forward FFTs are scaled 2x (for some reason)
Chris@29:         float two = 2.f;
Chris@29:         vDSP_vsdiv(m_fspare2, 1, &two, magOut, 1, hs1);
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         packComplex(realIn, imagIn);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         double *d[2] = { m_dpacked->realp, m_dpacked->imagp };
Chris@29:         v_deinterleave(d, complexIn, 2, m_size/2 + 1);
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         const int hs1 = m_size/2+1;
Chris@29:         vvsincos(m_dpacked->imagp, m_dpacked->realp, phaseIn, &hs1);
Chris@29:         double *const rp = m_dpacked->realp;
Chris@29:         double *const ip = m_dpacked->imagp;
Chris@29:         for (int i = 0; i < hs1; ++i) rp[i] *= magIn[i];
Chris@29:         for (int i = 0; i < hs1; ++i) ip[i] *= magIn[i];
Chris@29:         dnyq();
Chris@29:         vDSP_fft_zriptD(m_dspec, m_dpacked, 1, m_dbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_dspec) initDouble();
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         v_copy(m_dspare, magIn, hs1);
Chris@29:         for (int i = 0; i < hs1; ++i) m_dspare[i] += 0.000001;
Chris@29:         vvlog(m_dspare2, m_dspare, &hs1);
Chris@29:         inverse(m_dspare2, 0, cepOut);
Chris@29:     }
Chris@29:     
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         packComplex(realIn, imagIn);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         float *f[2] = { m_fpacked->realp, m_fpacked->imagp };
Chris@29:         v_deinterleave(f, complexIn, 2, m_size/2 + 1);
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29: 
Chris@29:         const int hs1 = m_size/2+1;
Chris@29:         vvsincosf(m_fpacked->imagp, m_fpacked->realp, phaseIn, &hs1);
Chris@29:         float *const rp = m_fpacked->realp;
Chris@29:         float *const ip = m_fpacked->imagp;
Chris@29:         for (int i = 0; i < hs1; ++i) rp[i] *= magIn[i];
Chris@29:         for (int i = 0; i < hs1; ++i) ip[i] *= magIn[i];
Chris@29:         fnyq();
Chris@29:         vDSP_fft_zript(m_fspec, m_fpacked, 1, m_fbuf, m_order, FFT_INVERSE);
Chris@29:         unpackReal(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_fspec) initFloat();
Chris@29:         const int hs1 = m_size/2 + 1;
Chris@29:         v_copy(m_fspare, magIn, hs1);
Chris@29:         for (int i = 0; i < hs1; ++i) m_fspare[i] += 0.000001f;
Chris@29:         vvlogf(m_fspare2, m_fspare, &hs1);
Chris@29:         inverse(m_fspare2, 0, cepOut);
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     int m_order;
Chris@29:     FFTSetup m_fspec;
Chris@29:     FFTSetupD m_dspec;
Chris@29:     DSPSplitComplex *m_fbuf;
Chris@29:     DSPDoubleSplitComplex *m_dbuf;
Chris@29:     DSPSplitComplex *m_fpacked;
Chris@29:     float *m_fspare;
Chris@29:     float *m_fspare2;
Chris@29:     DSPDoubleSplitComplex *m_dpacked;
Chris@29:     double *m_dspare;
Chris@29:     double *m_dspare2;
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_VDSP */
Chris@29: 
Chris@29: #ifdef HAVE_MEDIALIB
Chris@29: 
Chris@29: class D_MEDIALIB : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_MEDIALIB(int size) :
Chris@29:         m_size(size),
Chris@29:         m_dpacked(0), m_fpacked(0)
Chris@29:     { 
Chris@29:         for (int i = 0; ; ++i) {
Chris@29:             if (m_size & (1 << i)) {
Chris@29:                 m_order = i;
Chris@29:                 break;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     ~D_MEDIALIB() {
Chris@29:         if (m_dpacked) {
Chris@29:             deallocate(m_dpacked);
Chris@29:         }
Chris@29:         if (m_fpacked) {
Chris@29:             deallocate(m_fpacked);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::SinglePrecision | FFT::DoublePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     //!!! rv check
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         m_fpacked = allocate<float>(m_size*2);
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         m_dpacked = allocate<double>(m_size*2);
Chris@29:     }
Chris@29: 
Chris@29:     void packFloatConjugates() {
Chris@29:         const int hs = m_size / 2;
Chris@29:         for (int i = 1; i <= hs; ++i) {
Chris@29:             m_fpacked[(m_size-i)*2] = m_fpacked[2*i];
Chris@29:             m_fpacked[(m_size-i)*2 + 1] = -m_fpacked[2*i + 1];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packDoubleConjugates() {
Chris@29:         const int hs = m_size / 2;
Chris@29:         for (int i = 1; i <= hs; ++i) {
Chris@29:             m_dpacked[(m_size-i)*2] = m_dpacked[2*i];
Chris@29:             m_dpacked[(m_size-i)*2 + 1] = -m_dpacked[2*i + 1];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[index++] = re[i];
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_fpacked[index++] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_fpacked[index++] = 0.f;
Chris@29:             }
Chris@29:         }
Chris@29:         packFloatConjugates();
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im) {
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_dpacked[index++] = re[i];
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_dpacked[index++] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_dpacked[index++] = 0.0;
Chris@29:             }
Chris@29:         }
Chris@29:         packDoubleConjugates();
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(float *re, float *BQ_R__ im) { // re may be equal to m_fpacked
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = m_fpacked[index++];
Chris@29:             }
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_fpacked[index++];
Chris@29:             index++;
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackDouble(double *re, double *BQ_R__ im) { // re may be equal to m_dpacked
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = m_dpacked[index++];
Chris@29:             }
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_dpacked[index++];
Chris@29:             index++;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         mlib_SignalFFT_1_D64C_D64(m_dpacked, realIn, m_order);
Chris@29:         unpackDouble(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         // mlib FFT gives the whole redundant complex result
Chris@29:         mlib_SignalFFT_1_D64C_D64(m_dpacked, realIn, m_order);
Chris@29:         v_copy(complexOut, m_dpacked, m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         mlib_SignalFFT_1_D64C_D64(m_dpacked, realIn, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         int index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = index;
Chris@29:             ++index;
Chris@29:             magOut[i] = sqrt(m_dpacked[reali] * m_dpacked[reali] +
Chris@29:                              m_dpacked[index] * m_dpacked[index]);
Chris@29:             phaseOut[i] = atan2(m_dpacked[index], m_dpacked[reali]) ;
Chris@29:             ++index;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         mlib_SignalFFT_1_D64C_D64(m_dpacked, realIn, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         int index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = index;
Chris@29:             ++index;
Chris@29:             magOut[i] = sqrt(m_dpacked[reali] * m_dpacked[reali] +
Chris@29:                              m_dpacked[index] * m_dpacked[index]);
Chris@29:             ++index;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         mlib_SignalFFT_1_F32C_F32(m_fpacked, realIn, m_order);
Chris@29:         unpackFloat(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         // mlib FFT gives the whole redundant complex result
Chris@29:         mlib_SignalFFT_1_F32C_F32(m_fpacked, realIn, m_order);
Chris@29:         v_copy(complexOut, m_fpacked, m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         mlib_SignalFFT_1_F32C_F32(m_fpacked, realIn, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         int index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = index;
Chris@29:             ++index;
Chris@29:             magOut[i] = sqrtf(m_fpacked[reali] * m_fpacked[reali] +
Chris@29:                               m_fpacked[index] * m_fpacked[index]);
Chris@29:             phaseOut[i] = atan2f(m_fpacked[index], m_fpacked[reali]);
Chris@29:             ++index;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         mlib_SignalFFT_1_F32C_F32(m_fpacked, realIn, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         int index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = index;
Chris@29:             ++index;
Chris@29:             magOut[i] = sqrtf(m_fpacked[reali] * m_fpacked[reali] +
Chris@29:                               m_fpacked[index] * m_fpacked[index]);
Chris@29:             ++index;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         packDouble(realIn, imagIn);
Chris@29:         mlib_SignalIFFT_2_D64_D64C(realOut, m_dpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         v_copy(m_dpacked, complexIn, m_size + 2);
Chris@29:         packDoubleConjugates();
Chris@29:         mlib_SignalIFFT_2_D64_D64C(realOut, m_dpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = magIn[i] * cos(phaseIn[i]);
Chris@29:             double imag = magIn[i] * sin(phaseIn[i]);
Chris@29:             m_dpacked[i*2] = real;
Chris@29:             m_dpacked[i*2 + 1] = imag;
Chris@29:         }
Chris@29:         packDoubleConjugates();
Chris@29:         mlib_SignalIFFT_2_D64_D64C(realOut, m_dpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_dpacked) initDouble();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_dpacked[i*2] = log(magIn[i] + 0.000001);
Chris@29:             m_dpacked[i*2 + 1] = 0.0;
Chris@29:         }
Chris@29:         packDoubleConjugates();
Chris@29:         mlib_SignalIFFT_2_D64_D64C(cepOut, m_dpacked, m_order);
Chris@29:     }
Chris@29:     
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         packFloat(realIn, imagIn);
Chris@29:         mlib_SignalIFFT_2_F32_F32C(realOut, m_fpacked, m_order);
Chris@29:     }
Chris@29:     
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         v_convert(m_fpacked, complexIn, m_size + 2);
Chris@29:         packFloatConjugates();
Chris@29:         mlib_SignalIFFT_2_F32_F32C(realOut, m_fpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = magIn[i] * cos(phaseIn[i]);
Chris@29:             double imag = magIn[i] * sin(phaseIn[i]);
Chris@29:             m_fpacked[i*2] = real;
Chris@29:             m_fpacked[i*2 + 1] = imag;
Chris@29:         }
Chris@29:         packFloatConjugates();
Chris@29:         mlib_SignalIFFT_2_F32_F32C(realOut, m_fpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_fpacked) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i*2] = logf(magIn[i] + 0.000001);
Chris@29:             m_fpacked[i*2 + 1] = 0.f;
Chris@29:         }
Chris@29:         packFloatConjugates();
Chris@29:         mlib_SignalIFFT_2_F32_F32C(cepOut, m_fpacked, m_order);
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     int m_order;
Chris@29:     double *m_dpacked;
Chris@29:     float *m_fpacked;
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_MEDIALIB */
Chris@29: 
Chris@29: #ifdef HAVE_OPENMAX
Chris@29: 
Chris@29: class D_OPENMAX : public FFTImpl
Chris@29: {
Chris@29:     // Convert a signed 32-bit integer to a float in the range [-1,1)
Chris@29:     static inline float i2f(OMX_S32 i)
Chris@29:     {
Chris@29:         return float(i) / float(OMX_MAX_S32);
Chris@29:     }
Chris@29: 
Chris@29:     // Convert a signed 32-bit integer to a double in the range [-1,1)
Chris@29:     static inline double i2d(OMX_S32 i)
Chris@29:     {
Chris@29:         return double(i) / double(OMX_MAX_S32);
Chris@29:     }
Chris@29: 
Chris@29:     // Convert a float in the range [-1,1) to a signed 32-bit integer
Chris@29:     static inline OMX_S32 f2i(float f)
Chris@29:     {
Chris@29:         return OMX_S32(f * OMX_MAX_S32);
Chris@29:     }
Chris@29: 
Chris@29:     // Convert a double in the range [-1,1) to a signed 32-bit integer
Chris@29:     static inline OMX_S32 d2i(double d)
Chris@29:     {
Chris@29:         return OMX_S32(d * OMX_MAX_S32);
Chris@29:     }
Chris@29: 
Chris@29: public:
Chris@29:     D_OPENMAX(int size) :
Chris@29:         m_size(size),
Chris@29:         m_packed(0)
Chris@29:     { 
Chris@29:         for (int i = 0; ; ++i) {
Chris@29:             if (m_size & (1 << i)) {
Chris@29:                 m_order = i;
Chris@29:                 break;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     ~D_OPENMAX() {
Chris@29:         if (m_packed) {
Chris@29:             deallocate(m_packed);
Chris@29:             deallocate(m_buf);
Chris@29:             deallocate(m_fbuf);
Chris@29:             deallocate(m_spec);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::SinglePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     //!!! rv check
Chris@29: 
Chris@29:     // The OpenMAX implementation uses a fixed-point representation in
Chris@29:     // 32-bit signed integers, with a downward scaling factor (0-32
Chris@29:     // bits) supplied as an argument to the FFT function.
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         initDouble();
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         if (!m_packed) {
Chris@29:             m_buf = allocate<OMX_S32>(m_size);
Chris@29:             m_packed = allocate<OMX_S32>(m_size*2 + 2);
Chris@29:             m_fbuf = allocate<float>(m_size*2 + 2);
Chris@29:             OMX_INT sz = 0;
Chris@29:             omxSP_FFTGetBufSize_R_S32(m_order, &sz);
Chris@29:             m_spec = (OMXFFTSpec_R_S32 *)allocate<char>(sz);
Chris@29:             omxSP_FFTInit_R_S32(m_spec, m_order);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re) {
Chris@29:         // prepare fixed point input for forward transform
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             m_buf[i] = f2i(re[i]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re) {
Chris@29:         // prepare fixed point input for forward transform
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             m_buf[i] = d2i(re[i]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(float *BQ_R__ re, float *BQ_R__ im) {
Chris@29:         // convert fixed point output for forward transform
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = i2f(m_packed[index++]);
Chris@29:             }
Chris@29:             v_scale(im, m_size, hs + 1);
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = i2f(m_packed[index++]);
Chris@29:             index++;
Chris@29:         }
Chris@29:         v_scale(re, m_size, hs + 1);
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackDouble(double *BQ_R__ re, double *BQ_R__ im) {
Chris@29:         // convert fixed point output for forward transform
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 im[i] = i2d(m_packed[index++]);
Chris@29:             }
Chris@29:             v_scale(im, m_size, hs + 1);
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = i2d(m_packed[index++]);
Chris@29:             index++;
Chris@29:         }
Chris@29:         v_scale(re, m_size, hs + 1);
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloatInterleaved(float *BQ_R__ cplx) {
Chris@29:         // convert fixed point output for forward transform
Chris@29:         for (int i = 0; i < m_size + 2; ++i) {
Chris@29:             cplx[i] = i2f(m_packed[i]);
Chris@29:         }            
Chris@29:         v_scale(cplx, m_size, m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void unpackDoubleInterleaved(double *BQ_R__ cplx) {
Chris@29:         // convert fixed point output for forward transform
Chris@29:         for (int i = 0; i < m_size + 2; ++i) {
Chris@29:             cplx[i] = i2d(m_packed[i]);
Chris@29:         }            
Chris@29:         v_scale(cplx, m_size, m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {
Chris@29:         // prepare fixed point input for inverse transform
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_packed[index++] = f2i(re[i]);
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_packed[index++] = f2i(im[i]);
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_packed[index++] = 0;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im) {
Chris@29:         // prepare fixed point input for inverse transform
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_packed[index++] = d2i(re[i]);
Chris@29:             index++;
Chris@29:         }
Chris@29:         index = 0;
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_packed[index++] = d2i(im[i]);
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 index++;
Chris@29:                 m_packed[index++] = 0;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void convertFloat(const float *BQ_R__ f) {
Chris@29:         // convert interleaved input for inverse interleaved transform
Chris@29:         const int n = m_size + 2;
Chris@29:         for (int i = 0; i < n; ++i) {
Chris@29:             m_packed[i] = f2i(f[i]);
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void convertDouble(const double *BQ_R__ d) {
Chris@29:         // convert interleaved input for inverse interleaved transform
Chris@29:         const int n = m_size + 2;
Chris@29:         for (int i = 0; i < n; ++i) {
Chris@29:             m_packed[i] = d2i(d[i]);
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackFloat(float *BQ_R__ re) {
Chris@29:         // convert fixed point output for inverse transform
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             re[i] = i2f(m_buf[i]) * m_size;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackDouble(double *BQ_R__ re) {
Chris@29:         // convert fixed point output for inverse transform
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             re[i] = i2d(m_buf[i]) * m_size;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         packDouble(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackDouble(realOut, imagOut);
Chris@29:     }
Chris@29:     
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         packDouble(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackDoubleInterleaved(complexOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         packDouble(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackDouble(magOut, phaseOut); // temporarily
Chris@29:         // at this point we actually have real/imag in the mag/phase arrays
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = magOut[i];
Chris@29:             double imag = phaseOut[i];
Chris@29:             c_magphase(magOut + i, phaseOut + i, real, imag);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         packDouble(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = i * 2;
Chris@29:             int imagi = reali + 1;
Chris@29:             double real = i2d(m_packed[reali]) * m_size;
Chris@29:             double imag = i2d(m_packed[imagi]) * m_size;
Chris@29:             magOut[i] = sqrt(real * real + imag * imag);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         packFloat(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackFloat(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         packFloat(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackFloatInterleaved(complexOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29: 
Chris@29:         packFloat(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         unpackFloat(magOut, phaseOut); // temporarily
Chris@29:         // at this point we actually have real/imag in the mag/phase arrays
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             float real = magOut[i];
Chris@29:             float imag = phaseOut[i];
Chris@29:             c_magphase(magOut + i, phaseOut + i, real, imag);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         packFloat(realIn);
Chris@29:         omxSP_FFTFwd_RToCCS_S32_Sfs(m_buf, m_packed, m_spec, m_order);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             int reali = i * 2;
Chris@29:             int imagi = reali + 1;
Chris@29:             float real = i2f(m_packed[reali]) * m_size;
Chris@29:             float imag = i2f(m_packed[imagi]) * m_size;
Chris@29:             magOut[i] = sqrtf(real * real + imag * imag);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         packDouble(realIn, imagIn);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackDouble(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         convertDouble(complexIn);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackDouble(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         int index = 0;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real, imag;
Chris@29:             c_phasor(&real, &imag, phaseIn[i]);
Chris@29:             m_fbuf[index++] = float(real);
Chris@29:             m_fbuf[index++] = float(imag);
Chris@29:         }
Chris@29:         convertFloat(m_fbuf);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackDouble(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_packed) initDouble();
Chris@29:         //!!! implement
Chris@29: #warning OpenMAX implementation lacks cepstral transforms
Chris@29:     }
Chris@29:     
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         packFloat(realIn, imagIn);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackFloat(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         convertFloat(complexIn);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackFloat(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         v_polar_to_cartesian_interleaved(m_fbuf, magIn, phaseIn, hs+1);
Chris@29:         convertFloat(m_fbuf);
Chris@29:         omxSP_FFTInv_CCSToR_S32_Sfs(m_packed, m_buf, m_spec, 0);
Chris@29:         unpackFloat(realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_packed) initFloat();
Chris@29:         //!!! implement
Chris@29: #warning OpenMAX implementation lacks cepstral transforms
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     int m_order;
Chris@29:     OMX_S32 *m_packed;
Chris@29:     OMX_S32 *m_buf;
Chris@29:     float *m_fbuf;
Chris@29:     OMXFFTSpec_R_S32 *m_spec;
Chris@29: 
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_OPENMAX */
Chris@29: 
Chris@29: #ifdef HAVE_FFTW3
Chris@29: 
Chris@29: /*
Chris@29:  Define FFTW_DOUBLE_ONLY to make all uses of FFTW functions be
Chris@29:  double-precision (so "float" FFTs are calculated by casting to
Chris@29:  doubles and using the double-precision FFTW function).
Chris@29: 
Chris@29:  Define FFTW_SINGLE_ONLY to make all uses of FFTW functions be
Chris@29:  single-precision (so "double" FFTs are calculated by casting to
Chris@29:  floats and using the single-precision FFTW function).
Chris@29: 
Chris@29:  Neither of these flags is desirable for either performance or
Chris@29:  precision. The main reason to define either flag is to avoid linking
Chris@29:  against both fftw3 and fftw3f libraries.
Chris@29: */
Chris@29: 
Chris@29: //#define FFTW_DOUBLE_ONLY 1
Chris@29: //#define FFTW_SINGLE_ONLY 1
Chris@29: 
Chris@29: #if defined(FFTW_DOUBLE_ONLY) && defined(FFTW_SINGLE_ONLY)
Chris@29: // Can't meaningfully define both
Chris@29: #error Can only define one of FFTW_DOUBLE_ONLY and FFTW_SINGLE_ONLY
Chris@29: #endif
Chris@29: 
Chris@29: #if defined(FFTW_FLOAT_ONLY)
Chris@29: #warning FFTW_FLOAT_ONLY is deprecated, use FFTW_SINGLE_ONLY instead
Chris@29: #define FFTW_SINGLE_ONLY 1
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29: #define fft_float_type double
Chris@29: #define fftwf_complex fftw_complex
Chris@29: #define fftwf_plan fftw_plan
Chris@29: #define fftwf_plan_dft_r2c_1d fftw_plan_dft_r2c_1d
Chris@29: #define fftwf_plan_dft_c2r_1d fftw_plan_dft_c2r_1d
Chris@29: #define fftwf_destroy_plan fftw_destroy_plan
Chris@29: #define fftwf_malloc fftw_malloc
Chris@29: #define fftwf_free fftw_free
Chris@29: #define fftwf_execute fftw_execute
Chris@29: #define atan2f atan2
Chris@29: #define sqrtf sqrt
Chris@29: #define cosf cos
Chris@29: #define sinf sin
Chris@29: #else
Chris@29: #define fft_float_type float
Chris@29: #endif /* FFTW_DOUBLE_ONLY */
Chris@29: 
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29: #define fft_double_type float
Chris@29: #define fftw_complex fftwf_complex
Chris@29: #define fftw_plan fftwf_plan
Chris@29: #define fftw_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
Chris@29: #define fftw_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
Chris@29: #define fftw_destroy_plan fftwf_destroy_plan
Chris@29: #define fftw_malloc fftwf_malloc
Chris@29: #define fftw_free fftwf_free
Chris@29: #define fftw_execute fftwf_execute
Chris@29: #define atan2 atan2f
Chris@29: #define sqrt sqrtf
Chris@29: #define cos cosf
Chris@29: #define sin sinf
Chris@29: #else
Chris@29: #define fft_double_type double
Chris@29: #endif /* FFTW_SINGLE_ONLY */
Chris@29: 
Chris@29: class D_FFTW : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_FFTW(int size) :
Chris@29:         m_fplanf(0), m_dplanf(0), m_size(size)
Chris@29:     {
Chris@29:         initMutex();
Chris@29:     }
Chris@29: 
Chris@29:     ~D_FFTW() {
Chris@29:         if (m_fplanf) {
Chris@29:             lock();
Chris@29:             bool save = false;
Chris@29:             if (m_extantf > 0 && --m_extantf == 0) save = true;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:             if (save) saveWisdom('f');
Chris@29: #endif
Chris@29:             fftwf_destroy_plan(m_fplanf);
Chris@29:             fftwf_destroy_plan(m_fplani);
Chris@29:             fftwf_free(m_fbuf);
Chris@29:             fftwf_free(m_fpacked);
Chris@29:             unlock();
Chris@29:         }
Chris@29:         if (m_dplanf) {
Chris@29:             lock();
Chris@29:             bool save = false;
Chris@29:             if (m_extantd > 0 && --m_extantd == 0) save = true;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:             if (save) saveWisdom('d');
Chris@29: #endif
Chris@29:             fftw_destroy_plan(m_dplanf);
Chris@29:             fftw_destroy_plan(m_dplani);
Chris@29:             fftw_free(m_dbuf);
Chris@29:             fftw_free(m_dpacked);
Chris@29:             unlock();
Chris@29:         }
Chris@29:         destroyMutex();
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:         return FFT::SinglePrecision;
Chris@29: #else
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:         return FFT::DoublePrecision;
Chris@29: #else
Chris@29:         return FFT::SinglePrecision | FFT::DoublePrecision;
Chris@29: #endif
Chris@29: #endif
Chris@29:     }
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         if (m_fplanf) return;
Chris@29:         bool load = false;
Chris@29:         lock();
Chris@29:         if (m_extantf++ == 0) load = true;
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:         if (load) loadWisdom('d');
Chris@29: #else
Chris@29:         if (load) loadWisdom('f');
Chris@29: #endif
Chris@29:         m_fbuf = (fft_float_type *)fftw_malloc(m_size * sizeof(fft_float_type));
Chris@29:         m_fpacked = (fftwf_complex *)fftw_malloc
Chris@29:             ((m_size/2 + 1) * sizeof(fftwf_complex));
Chris@29:         m_fplanf = fftwf_plan_dft_r2c_1d
Chris@29:             (m_size, m_fbuf, m_fpacked, FFTW_MEASURE);
Chris@29:         m_fplani = fftwf_plan_dft_c2r_1d
Chris@29:             (m_size, m_fpacked, m_fbuf, FFTW_MEASURE);
Chris@29:         unlock();
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         if (m_dplanf) return;
Chris@29:         bool load = false;
Chris@29:         lock();
Chris@29:         if (m_extantd++ == 0) load = true;
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:         if (load) loadWisdom('f');
Chris@29: #else
Chris@29:         if (load) loadWisdom('d');
Chris@29: #endif
Chris@29:         m_dbuf = (fft_double_type *)fftw_malloc(m_size * sizeof(fft_double_type));
Chris@29:         m_dpacked = (fftw_complex *)fftw_malloc
Chris@29:             ((m_size/2 + 1) * sizeof(fftw_complex));
Chris@29:         m_dplanf = fftw_plan_dft_r2c_1d
Chris@29:             (m_size, m_dbuf, m_dpacked, FFTW_MEASURE);
Chris@29:         m_dplani = fftw_plan_dft_c2r_1d
Chris@29:             (m_size, m_dpacked, m_dbuf, FFTW_MEASURE);
Chris@29:         unlock();
Chris@29:     }
Chris@29: 
Chris@29:     void loadWisdom(char type) { wisdom(false, type); }
Chris@29:     void saveWisdom(char type) { wisdom(true, type); }
Chris@29: 
Chris@29:     void wisdom(bool save, char type) {
Chris@29: 
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:         if (type == 'f') return;
Chris@29: #endif
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:         if (type == 'd') return;
Chris@29: #endif
Chris@29: 
Chris@29:         const char *home = getenv("HOME");
Chris@29:         if (!home) return;
Chris@29: 
Chris@29:         char fn[256];
Chris@29:         snprintf(fn, 256, "%s/%s.%c", home, ".turbot.wisdom", type);
Chris@29: 
Chris@29:         FILE *f = fopen(fn, save ? "wb" : "rb");
Chris@29:         if (!f) return;
Chris@29: 
Chris@29:         if (save) {
Chris@29:             switch (type) {
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:             case 'f': break;
Chris@29: #else
Chris@29:             case 'f': fftwf_export_wisdom_to_file(f); break;
Chris@29: #endif
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:             case 'd': break;
Chris@29: #else
Chris@29:             case 'd': fftw_export_wisdom_to_file(f); break;
Chris@29: #endif
Chris@29:             default: break;
Chris@29:             }
Chris@29:         } else {
Chris@29:             switch (type) {
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:             case 'f': break;
Chris@29: #else
Chris@29:             case 'f': fftwf_import_wisdom_from_file(f); break;
Chris@29: #endif
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:             case 'd': break;
Chris@29: #else
Chris@29:             case 'd': fftw_import_wisdom_from_file(f); break;
Chris@29: #endif
Chris@29:             default: break;
Chris@29:             }
Chris@29:         }
Chris@29: 
Chris@29:         fclose(f);
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         fftwf_complex *const BQ_R__ fpacked = m_fpacked; 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             fpacked[i][0] = re[i];
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 fpacked[i][1] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 fpacked[i][1] = 0.f;
Chris@29:             }
Chris@29:         }                
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         fftw_complex *const BQ_R__ dpacked = m_dpacked; 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             dpacked[i][0] = re[i];
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 dpacked[i][1] = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 dpacked[i][1] = 0.0;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(float *BQ_R__ re, float *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_fpacked[i][0];
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 im[i] = m_fpacked[i][1];
Chris@29:             }
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackDouble(double *BQ_R__ re, double *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_dpacked[i][0];
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 im[i] = m_dpacked[i][1];
Chris@29:             }
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         const int sz = m_size;
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realIn != dbuf) 
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 dbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftw_execute(m_dplanf);
Chris@29:         unpackDouble(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         const int sz = m_size;
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realIn != dbuf) 
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 dbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftw_execute(m_dplanf);
Chris@29:         v_convert(complexOut, (fft_double_type *)m_dpacked, sz + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realIn != dbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 dbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftw_execute(m_dplanf);
Chris@29:         v_cartesian_interleaved_to_polar(magOut, phaseOut,
Chris@29:                                          (double *)m_dpacked, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realIn != m_dbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 dbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftw_execute(m_dplanf);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_dpacked[i][0] * m_dpacked[i][0] +
Chris@29:                              m_dpacked[i][1] * m_dpacked[i][1]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realIn != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 fbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftwf_execute(m_fplanf);
Chris@29:         unpackFloat(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realIn != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 fbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftwf_execute(m_fplanf);
Chris@29:         v_convert(complexOut, (fft_float_type *)m_fpacked, sz + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realIn != fbuf) 
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 fbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftwf_execute(m_fplanf);
Chris@29:         v_cartesian_interleaved_to_polar(magOut, phaseOut,
Chris@29:                                          (float *)m_fpacked, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realIn != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 fbuf[i] = realIn[i];
Chris@29:             }
Chris@29:         fftwf_execute(m_fplanf);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrtf(m_fpacked[i][0] * m_fpacked[i][0] +
Chris@29:                               m_fpacked[i][1] * m_fpacked[i][1]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, imagIn);
Chris@29:         fftw_execute(m_dplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realOut != dbuf) 
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = dbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         v_convert((double *)m_dpacked, complexIn, m_size + 2);
Chris@29:         fftw_execute(m_dplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realOut != dbuf) 
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = dbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         const int hs = m_size/2;
Chris@29:         fftw_complex *const BQ_R__ dpacked = m_dpacked;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             dpacked[i][0] = magIn[i] * cos(phaseIn[i]);
Chris@29:         }
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             dpacked[i][1] = magIn[i] * sin(phaseIn[i]);
Chris@29:         }
Chris@29:         fftw_execute(m_dplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (realOut != dbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = dbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         fft_double_type *const BQ_R__ dbuf = m_dbuf;
Chris@29:         fftw_complex *const BQ_R__ dpacked = m_dpacked;
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             dpacked[i][0] = log(magIn[i] + 0.000001);
Chris@29:         }
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             dpacked[i][1] = 0.0;
Chris@29:         }
Chris@29:         fftw_execute(m_dplani);
Chris@29:         const int sz = m_size;
Chris@29: #ifndef FFTW_SINGLE_ONLY
Chris@29:         if (cepOut != dbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 cepOut[i] = dbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, imagIn);
Chris@29:         fftwf_execute(m_fplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realOut != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = fbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         v_copy((float *)m_fpacked, complexIn, m_size + 2);
Chris@29:         fftwf_execute(m_fplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realOut != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = fbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         fftwf_complex *const BQ_R__ fpacked = m_fpacked;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             fpacked[i][0] = magIn[i] * cosf(phaseIn[i]);
Chris@29:         }
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             fpacked[i][1] = magIn[i] * sinf(phaseIn[i]);
Chris@29:         }
Chris@29:         fftwf_execute(m_fplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (realOut != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 realOut[i] = fbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         fftwf_complex *const BQ_R__ fpacked = m_fpacked;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             fpacked[i][0] = logf(magIn[i] + 0.000001f);
Chris@29:         }
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             fpacked[i][1] = 0.f;
Chris@29:         }
Chris@29:         fftwf_execute(m_fplani);
Chris@29:         const int sz = m_size;
Chris@29:         fft_float_type *const BQ_R__ fbuf = m_fbuf;
Chris@29: #ifndef FFTW_DOUBLE_ONLY
Chris@29:         if (cepOut != fbuf)
Chris@29: #endif
Chris@29:             for (int i = 0; i < sz; ++i) {
Chris@29:                 cepOut[i] = fbuf[i];
Chris@29:             }
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     fftwf_plan m_fplanf;
Chris@29:     fftwf_plan m_fplani;
Chris@29: #ifdef FFTW_DOUBLE_ONLY
Chris@29:     double *m_fbuf;
Chris@29: #else
Chris@29:     float *m_fbuf;
Chris@29: #endif
Chris@29:     fftwf_complex *m_fpacked;
Chris@29:     fftw_plan m_dplanf;
Chris@29:     fftw_plan m_dplani;
Chris@29: #ifdef FFTW_SINGLE_ONLY
Chris@29:     float *m_dbuf;
Chris@29: #else
Chris@29:     double *m_dbuf;
Chris@29: #endif
Chris@29:     fftw_complex *m_dpacked;
Chris@29:     const int m_size;
Chris@29:     static int m_extantf;
Chris@29:     static int m_extantd;
Chris@29: #ifdef NO_THREADING
Chris@29:     void initMutex() {}
Chris@29:     void destroyMutex() {}
Chris@29:     void lock() {}
Chris@29:     void unlock() {}
Chris@29: #else
Chris@29: #ifdef _WIN32
Chris@29:     static HANDLE m_commonMutex;
Chris@29:     void initMutex() { m_commonMutex = CreateMutex(NULL, FALSE, NULL); }
Chris@29:     void destroyMutex() { CloseHandle(m_commonMutex); }
Chris@29:     void lock() { WaitForSingleObject(m_commonMutex, INFINITE); }
Chris@29:     void unlock() { ReleaseMutex(m_commonMutex); }
Chris@29: #else
Chris@29:     static pthread_mutex_t m_commonMutex;
Chris@29:     void initMutex() { pthread_mutex_init(&m_commonMutex, 0); }
Chris@29:     void destroyMutex() { pthread_mutex_destroy(&m_commonMutex); }
Chris@29:     void lock() { pthread_mutex_lock(&m_commonMutex); }
Chris@29:     void unlock() { pthread_mutex_unlock(&m_commonMutex); }
Chris@29: #endif
Chris@29: #endif
Chris@29: };
Chris@29: 
Chris@29: int
Chris@29: D_FFTW::m_extantf = 0;
Chris@29: 
Chris@29: int
Chris@29: D_FFTW::m_extantd = 0;
Chris@29: 
Chris@29: #ifndef NO_THREADING
Chris@29: #ifdef _WIN32
Chris@29: HANDLE D_FFTW::m_commonMutex;
Chris@29: #else
Chris@29: pthread_mutex_t D_FFTW::m_commonMutex;
Chris@29: #endif
Chris@29: #endif
Chris@29: 
Chris@29: #endif /* HAVE_FFTW3 */
Chris@29: 
Chris@29: #ifdef HAVE_SFFT
Chris@29: 
Chris@29: /*
Chris@29:  Define SFFT_DOUBLE_ONLY to make all uses of SFFT functions be
Chris@29:  double-precision (so "float" FFTs are calculated by casting to
Chris@29:  doubles and using the double-precision SFFT function).
Chris@29: 
Chris@29:  Define SFFT_SINGLE_ONLY to make all uses of SFFT functions be
Chris@29:  single-precision (so "double" FFTs are calculated by casting to
Chris@29:  floats and using the single-precision SFFT function).
Chris@29: 
Chris@29:  Neither of these flags is desirable for either performance or
Chris@29:  precision.
Chris@29: */
Chris@29: 
Chris@29: //#define SFFT_DOUBLE_ONLY 1
Chris@29: //#define SFFT_SINGLE_ONLY 1
Chris@29: 
Chris@29: #if defined(SFFT_DOUBLE_ONLY) && defined(SFFT_SINGLE_ONLY)
Chris@29: // Can't meaningfully define both
Chris@29: #error Can only define one of SFFT_DOUBLE_ONLY and SFFT_SINGLE_ONLY
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef SFFT_DOUBLE_ONLY
Chris@29: #define fft_float_type double
Chris@29: #define FLAG_SFFT_FLOAT SFFT_DOUBLE
Chris@29: #define atan2f atan2
Chris@29: #define sqrtf sqrt
Chris@29: #define cosf cos
Chris@29: #define sinf sin
Chris@29: #define logf log
Chris@29: #else
Chris@29: #define FLAG_SFFT_FLOAT SFFT_FLOAT
Chris@29: #define fft_float_type float
Chris@29: #endif /* SFFT_DOUBLE_ONLY */
Chris@29: 
Chris@29: #ifdef SFFT_SINGLE_ONLY
Chris@29: #define fft_double_type float
Chris@29: #define FLAG_SFFT_DOUBLE SFFT_FLOAT
Chris@29: #define atan2 atan2f
Chris@29: #define sqrt sqrtf
Chris@29: #define cos cosf
Chris@29: #define sin sinf
Chris@29: #define log logf
Chris@29: #else
Chris@29: #define FLAG_SFFT_DOUBLE SFFT_DOUBLE
Chris@29: #define fft_double_type double
Chris@29: #endif /* SFFT_SINGLE_ONLY */
Chris@29: 
Chris@29: class D_SFFT : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_SFFT(int size) :
Chris@29:         m_fplanf(0), m_fplani(0), m_dplanf(0), m_dplani(0), m_size(size)
Chris@29:     {
Chris@29:     }
Chris@29: 
Chris@29:     ~D_SFFT() {
Chris@29:         if (m_fplanf) {
Chris@29:             sfft_free(m_fplanf);
Chris@29:             sfft_free(m_fplani);
Chris@29:             deallocate(m_fbuf);
Chris@29:             deallocate(m_fresult);
Chris@29:         }
Chris@29:         if (m_dplanf) {
Chris@29:             sfft_free(m_dplanf);
Chris@29:             sfft_free(m_dplani);
Chris@29:             deallocate(m_dbuf);
Chris@29:             deallocate(m_dresult);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29: #ifdef SFFT_SINGLE_ONLY
Chris@29:         return FFT::SinglePrecision;
Chris@29: #else
Chris@29: #ifdef SFFT_DOUBLE_ONLY
Chris@29:         return FFT::DoublePrecision;
Chris@29: #else
Chris@29:         return FFT::SinglePrecision | FFT::DoublePrecision;
Chris@29: #endif
Chris@29: #endif
Chris@29:     }
Chris@29: 
Chris@29:     void initFloat() {
Chris@29:         if (m_fplanf) return;
Chris@29:         m_fbuf = allocate<fft_float_type>(2 * m_size);
Chris@29:         m_fresult = allocate<fft_float_type>(2 * m_size);
Chris@29:         m_fplanf = sfft_init(m_size, SFFT_FORWARD | FLAG_SFFT_FLOAT);
Chris@29:         m_fplani = sfft_init(m_size, SFFT_BACKWARD | FLAG_SFFT_FLOAT);
Chris@29:         if (!m_fplanf || !m_fplani) {
Chris@29:             if (!m_fplanf) {
Chris@29:                 std::cerr << "D_SFFT: Failed to construct forward float transform for size " << m_size << " (check SFFT library's target configuration)" << std::endl;
Chris@29:             } else {
Chris@29:                 std::cerr << "D_SFFT: Failed to construct inverse float transform for size " << m_size << " (check SFFT library's target configuration)" << std::endl;
Chris@29:             }
Chris@29: #ifndef NO_EXCEPTIONS
Chris@29:             throw FFT::InternalError;
Chris@29: #else
Chris@29:             abort();
Chris@29: #endif
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void initDouble() {
Chris@29:         if (m_dplanf) return;
Chris@29:         m_dbuf = allocate<fft_double_type>(2 * m_size);
Chris@29:         m_dresult = allocate<fft_double_type>(2 * m_size);
Chris@29:         m_dplanf = sfft_init(m_size, SFFT_FORWARD | FLAG_SFFT_DOUBLE);
Chris@29:         m_dplani = sfft_init(m_size, SFFT_BACKWARD | FLAG_SFFT_DOUBLE);
Chris@29:         if (!m_dplanf || !m_dplani) {
Chris@29:             if (!m_dplanf) {
Chris@29:                 std::cerr << "D_SFFT: Failed to construct forward double transform for size " << m_size << " (check SFFT library's target configuration)" << std::endl;
Chris@29:             } else {
Chris@29:                 std::cerr << "D_SFFT: Failed to construct inverse double transform for size " << m_size << " (check SFFT library's target configuration)" << std::endl;
Chris@29:             }
Chris@29: #ifndef NO_EXCEPTIONS
Chris@29:             throw FFT::InternalError;
Chris@29: #else
Chris@29:             abort();
Chris@29: #endif
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im, fft_float_type *target, int n) {
Chris@29:         for (int i = 0; i < n; ++i) target[i*2] = re[i];
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i < n; ++i) target[i*2+1] = im[i]; 
Chris@29:         } else {
Chris@29:             for (int i = 0; i < n; ++i) target[i*2+1] = 0.f;
Chris@29:         }                
Chris@29:     }
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im, fft_double_type *target, int n) {
Chris@29:         for (int i = 0; i < n; ++i) target[i*2] = re[i];
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i < n; ++i) target[i*2+1] = im[i];
Chris@29:         } else {
Chris@29:             for (int i = 0; i < n; ++i) target[i*2+1] = 0.0;
Chris@29:         }                
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(const fft_float_type *source, float *BQ_R__ re, float *BQ_R__ im, int n) {
Chris@29:         for (int i = 0; i < n; ++i) re[i] = source[i*2];
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i < n; ++i) im[i] = source[i*2+1];
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void unpackDouble(const fft_double_type *source, double *BQ_R__ re, double *BQ_R__ im, int n) {
Chris@29:         for (int i = 0; i < n; ++i) re[i] = source[i*2];
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i < n; ++i) im[i] = source[i*2+1];
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     template<typename T>
Chris@29:     void mirror(T *BQ_R__ cplx, int n) {
Chris@29:         for (int i = 1; i <= n/2; ++i) {
Chris@29:             int j = n-i;
Chris@29:             cplx[j*2] = cplx[i*2];
Chris@29:             cplx[j*2+1] = -cplx[i*2+1];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, 0, m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplanf, m_dbuf, m_dresult);
Chris@29:         unpackDouble(m_dresult, realOut, imagOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, 0, m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplanf, m_dbuf, m_dresult);
Chris@29:         v_convert(complexOut, m_dresult, m_size+2); // i.e. m_size/2+1 complex
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, 0, m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplanf, m_dbuf, m_dresult);
Chris@29:         v_cartesian_interleaved_to_polar(magOut, phaseOut,
Chris@29:                                          m_dresult, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, 0, m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplanf, m_dbuf, m_dresult);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_dresult[i*2] * m_dresult[i*2] +
Chris@29:                              m_dresult[i*2+1] * m_dresult[i*2+1]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, 0, m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplanf, m_fbuf, m_fresult);
Chris@29:         unpackFloat(m_fresult, realOut, imagOut, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, 0, m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplanf, m_fbuf, m_fresult);
Chris@29:         v_convert(complexOut, m_fresult, m_size+2); // i.e. m_size/2+1 complex
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, 0, m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplanf, m_fbuf, m_fresult);
Chris@29:         v_cartesian_interleaved_to_polar(magOut, phaseOut,
Chris@29:                                          m_fresult, m_size/2+1);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, 0, m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplanf, m_fbuf, m_fresult);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrtf(m_fresult[i*2] * m_fresult[i*2] +
Chris@29:                               m_fresult[i*2+1] * m_fresult[i*2+1]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         packDouble(realIn, imagIn, m_dbuf, m_size/2+1);
Chris@29:         mirror(m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplani, m_dbuf, m_dresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_dresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         v_convert((double *)m_dbuf, complexIn, m_size + 2);
Chris@29:         mirror(m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplani, m_dbuf, m_dresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_dresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_dbuf[i*2] = magIn[i] * cos(phaseIn[i]);
Chris@29:             m_dbuf[i*2+1] = magIn[i] * sin(phaseIn[i]);
Chris@29:         }
Chris@29:         mirror(m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplani, m_dbuf, m_dresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_dresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         if (!m_dplanf) initDouble();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_dbuf[i*2] = log(magIn[i] + 0.000001);
Chris@29:             m_dbuf[i*2+1] = 0.0;
Chris@29:         }
Chris@29:         mirror(m_dbuf, m_size);
Chris@29:         sfft_execute(m_dplani, m_dbuf, m_dresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             cepOut[i] = m_dresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         packFloat(realIn, imagIn, m_fbuf, m_size/2+1);
Chris@29:         mirror(m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplani, m_fbuf, m_fresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         v_convert((float *)m_fbuf, complexIn, m_size + 2);
Chris@29:         mirror(m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplani, m_fbuf, m_fresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fbuf[i*2] = magIn[i] * cosf(phaseIn[i]);
Chris@29:             m_fbuf[i*2+1] = magIn[i] * sinf(phaseIn[i]);
Chris@29:         }
Chris@29:         mirror(m_fbuf, m_size);
Chris@29:         sfft_execute(m_fplani, m_fbuf, m_fresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         if (!m_fplanf) initFloat();
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fbuf[i*2] = logf(magIn[i] + 0.00001);
Chris@29:             m_fbuf[i*2+1] = 0.0f;
Chris@29:         }
Chris@29:         sfft_execute(m_fplani, m_fbuf, m_fresult);
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             cepOut[i] = m_fresult[i*2];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     sfft_plan_t *m_fplanf;
Chris@29:     sfft_plan_t *m_fplani;
Chris@29:     fft_float_type *m_fbuf;
Chris@29:     fft_float_type *m_fresult;
Chris@29: 
Chris@29:     sfft_plan_t *m_dplanf;
Chris@29:     sfft_plan_t *m_dplani;
Chris@29:     fft_double_type *m_dbuf;
Chris@29:     fft_double_type *m_dresult;
Chris@29: 
Chris@29:     const int m_size;
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_SFFT */
Chris@29: 
Chris@29: #ifdef HAVE_KISSFFT
Chris@29: 
Chris@29: class D_KISSFFT : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_KISSFFT(int size) :
Chris@29:         m_size(size),
Chris@29:         m_fplanf(0),  
Chris@29:         m_fplani(0)
Chris@29:     {
Chris@29: #ifdef FIXED_POINT
Chris@29: #error KISSFFT is not configured for float values
Chris@29: #endif
Chris@29:         if (sizeof(kiss_fft_scalar) != sizeof(float)) {
Chris@29:             std::cerr << "ERROR: KISSFFT is not configured for float values"
Chris@29:                       << std::endl;
Chris@29:         }
Chris@29: 
Chris@29:         m_fbuf = new kiss_fft_scalar[m_size + 2];
Chris@29:         m_fpacked = new kiss_fft_cpx[m_size + 2];
Chris@29:         m_fplanf = kiss_fftr_alloc(m_size, 0, NULL, NULL);
Chris@29:         m_fplani = kiss_fftr_alloc(m_size, 1, NULL, NULL);
Chris@29:     }
Chris@29: 
Chris@29:     ~D_KISSFFT() {
Chris@29:         kiss_fftr_free(m_fplanf);
Chris@29:         kiss_fftr_free(m_fplani);
Chris@29:         kiss_fft_cleanup();
Chris@29: 
Chris@29:         delete[] m_fbuf;
Chris@29:         delete[] m_fpacked;
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::SinglePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     void initFloat() { }
Chris@29:     void initDouble() { }
Chris@29: 
Chris@29:     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = re[i];
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 m_fpacked[i].i = im[i];
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 m_fpacked[i].i = 0.f;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackFloat(float *BQ_R__ re, float *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = m_fpacked[i].r;
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 im[i] = m_fpacked[i].i;
Chris@29:             }
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void packDouble(const double *BQ_R__ re, const double *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = float(re[i]);
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 m_fpacked[i].i = float(im[i]);
Chris@29:             }
Chris@29:         } else {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 m_fpacked[i].i = 0.f;
Chris@29:             }
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void unpackDouble(double *BQ_R__ re, double *BQ_R__ im) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             re[i] = double(m_fpacked[i].r);
Chris@29:         }
Chris@29:         if (im) {
Chris@29:             for (int i = 0; i <= hs; ++i) {
Chris@29:                 im[i] = double(m_fpacked[i].i);
Chris@29:             }
Chris@29:         }
Chris@29:     }        
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29: 
Chris@29:         v_convert(m_fbuf, realIn, m_size);
Chris@29:         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
Chris@29:         unpackDouble(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29: 
Chris@29:         v_convert(m_fbuf, realIn, m_size);
Chris@29:         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
Chris@29:         v_convert(complexOut, (float *)m_fpacked, m_size + 2);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             m_fbuf[i] = float(realIn[i]);
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
Chris@29:                              double(m_fpacked[i].i) * double(m_fpacked[i].i));
Chris@29:         }
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             phaseOut[i] = atan2(double(m_fpacked[i].i), double(m_fpacked[i].r));
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             m_fbuf[i] = float(realIn[i]);
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, m_fbuf, m_fpacked);
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(double(m_fpacked[i].r) * double(m_fpacked[i].r) +
Chris@29:                              double(m_fpacked[i].i) * double(m_fpacked[i].i));
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, realIn, m_fpacked);
Chris@29:         unpackFloat(realOut, imagOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, realIn, (kiss_fft_cpx *)complexOut);
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, realIn, m_fpacked);
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
Chris@29:                               m_fpacked[i].i * m_fpacked[i].i);
Chris@29:         }
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             phaseOut[i] = atan2f(m_fpacked[i].i, m_fpacked[i].r);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29: 
Chris@29:         kiss_fftr(m_fplanf, realIn, m_fpacked);
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrtf(m_fpacked[i].r * m_fpacked[i].r +
Chris@29:                               m_fpacked[i].i * m_fpacked[i].i);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29: 
Chris@29:         packDouble(realIn, imagIn);
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fbuf[i];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29: 
Chris@29:         v_convert((float *)m_fpacked, complexIn, m_size + 2);
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fbuf[i];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = float(magIn[i] * cos(phaseIn[i]));
Chris@29:             m_fpacked[i].i = float(magIn[i] * sin(phaseIn[i]));
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             realOut[i] = m_fbuf[i];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = float(log(magIn[i] + 0.000001));
Chris@29:             m_fpacked[i].i = 0.0f;
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, m_fbuf);
Chris@29: 
Chris@29:         for (int i = 0; i < m_size; ++i) {
Chris@29:             cepOut[i] = m_fbuf[i];
Chris@29:         }
Chris@29:     }
Chris@29:     
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29: 
Chris@29:         packFloat(realIn, imagIn);
Chris@29:         kiss_fftri(m_fplani, m_fpacked, realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29: 
Chris@29:         v_copy((float *)m_fpacked, complexIn, m_size + 2);
Chris@29:         kiss_fftri(m_fplani, m_fpacked, realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = magIn[i] * cosf(phaseIn[i]);
Chris@29:             m_fpacked[i].i = magIn[i] * sinf(phaseIn[i]);
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, realOut);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29: 
Chris@29:         const int hs = m_size/2;
Chris@29: 
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             m_fpacked[i].r = logf(magIn[i] + 0.000001f);
Chris@29:             m_fpacked[i].i = 0.0f;
Chris@29:         }
Chris@29: 
Chris@29:         kiss_fftri(m_fplani, m_fpacked, cepOut);
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     kiss_fftr_cfg m_fplanf;
Chris@29:     kiss_fftr_cfg m_fplani;
Chris@29:     kiss_fft_scalar *m_fbuf;
Chris@29:     kiss_fft_cpx *m_fpacked;
Chris@29: };
Chris@29: 
Chris@29: #endif /* HAVE_KISSFFT */
Chris@29: 
Chris@29: #ifdef USE_BUILTIN_FFT
Chris@29: 
Chris@29: class D_Cross : public FFTImpl
Chris@29: {
Chris@29: public:
Chris@29:     D_Cross(int size) : m_size(size), m_table(0) {
Chris@29:         
Chris@29:         m_a = new double[size];
Chris@29:         m_b = new double[size];
Chris@29:         m_c = new double[size];
Chris@29:         m_d = new double[size];
Chris@29: 
Chris@29:         m_table = new int[m_size];
Chris@29:     
Chris@29:         int bits;
Chris@29:         int i, j, k, m;
Chris@29: 
Chris@29:         for (i = 0; ; ++i) {
Chris@29:             if (m_size & (1 << i)) {
Chris@29:                 bits = i;
Chris@29:                 break;
Chris@29:             }
Chris@29:         }
Chris@29:         
Chris@29:         for (i = 0; i < m_size; ++i) {
Chris@29:             
Chris@29:             m = i;
Chris@29:             
Chris@29:             for (j = k = 0; j < bits; ++j) {
Chris@29:                 k = (k << 1) | (m & 1);
Chris@29:                 m >>= 1;
Chris@29:             }
Chris@29:             
Chris@29:             m_table[i] = k;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     ~D_Cross() {
Chris@29:         delete[] m_table;
Chris@29:         delete[] m_a;
Chris@29:         delete[] m_b;
Chris@29:         delete[] m_c;
Chris@29:         delete[] m_d;
Chris@29:     }
Chris@29: 
Chris@29:     FFT::Precisions
Chris@29:     getSupportedPrecisions() const {
Chris@29:         return FFT::DoublePrecision;
Chris@29:     }
Chris@29: 
Chris@29:     void initFloat() { }
Chris@29:     void initDouble() { }
Chris@29: 
Chris@29:     void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
Chris@29:         basefft(false, realIn, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
Chris@29:         if (imagOut) {
Chris@29:             for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
Chris@29:         basefft(false, realIn, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) complexOut[i*2] = m_c[i];
Chris@29:         for (int i = 0; i <= hs; ++i) complexOut[i*2+1] = m_d[i];
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
Chris@29:         basefft(false, realIn, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
Chris@29:             phaseOut[i] = atan2(m_d[i], m_c[i]) ;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
Chris@29:         basefft(false, realIn, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
Chris@29:         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
Chris@29:         basefft(false, m_a, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) realOut[i] = m_c[i];
Chris@29:         if (imagOut) {
Chris@29:             for (int i = 0; i <= hs; ++i) imagOut[i] = m_d[i];
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
Chris@29:         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
Chris@29:         basefft(false, m_a, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) complexOut[i*2] = m_c[i];
Chris@29:         for (int i = 0; i <= hs; ++i) complexOut[i*2+1] = m_d[i];
Chris@29:     }
Chris@29: 
Chris@29:     void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
Chris@29:         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
Chris@29:         basefft(false, m_a, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
Chris@29:             phaseOut[i] = atan2(m_d[i], m_c[i]) ;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
Chris@29:         for (int i = 0; i < m_size; ++i) m_a[i] = realIn[i];
Chris@29:         basefft(false, m_a, 0, m_c, m_d);
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             magOut[i] = sqrt(m_c[i] * m_c[i] + m_d[i] * m_d[i]);
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = realIn[i];
Chris@29:             double imag = imagIn[i];
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, realOut, m_d);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = complexIn[i*2];
Chris@29:             double imag = complexIn[i*2+1];
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, realOut, m_d);
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = magIn[i] * cos(phaseIn[i]);
Chris@29:             double imag = magIn[i] * sin(phaseIn[i]);
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, realOut, m_d);
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             double real = log(magIn[i] + 0.000001);
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = 0.0;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = 0.0;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, cepOut, m_d);
Chris@29:     }
Chris@29: 
Chris@29:     void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             float real = realIn[i];
Chris@29:             float imag = imagIn[i];
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, m_c, m_d);
Chris@29:         for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
Chris@29:     }
Chris@29: 
Chris@29:     void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             float real = complexIn[i*2];
Chris@29:             float imag = complexIn[i*2+1];
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, m_c, m_d);
Chris@29:         for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
Chris@29:     }
Chris@29: 
Chris@29:     void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             float real = magIn[i] * cosf(phaseIn[i]);
Chris@29:             float imag = magIn[i] * sinf(phaseIn[i]);
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = imag;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = -imag;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, m_c, m_d);
Chris@29:         for (int i = 0; i < m_size; ++i) realOut[i] = m_c[i];
Chris@29:     }
Chris@29: 
Chris@29:     void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
Chris@29:         const int hs = m_size/2;
Chris@29:         for (int i = 0; i <= hs; ++i) {
Chris@29:             float real = logf(magIn[i] + 0.000001);
Chris@29:             m_a[i] = real;
Chris@29:             m_b[i] = 0.0;
Chris@29:             if (i > 0) {
Chris@29:                 m_a[m_size-i] = real;
Chris@29:                 m_b[m_size-i] = 0.0;
Chris@29:             }
Chris@29:         }
Chris@29:         basefft(true, m_a, m_b, m_c, m_d);
Chris@29:         for (int i = 0; i < m_size; ++i) cepOut[i] = m_c[i];
Chris@29:     }
Chris@29: 
Chris@29: private:
Chris@29:     const int m_size;
Chris@29:     int *m_table;
Chris@29:     double *m_a;
Chris@29:     double *m_b;
Chris@29:     double *m_c;
Chris@29:     double *m_d;
Chris@29:     void basefft(bool inverse, const double *BQ_R__ ri, const double *BQ_R__ ii, double *BQ_R__ ro, double *BQ_R__ io);
Chris@29: };
Chris@29: 
Chris@29: void
Chris@29: D_Cross::basefft(bool inverse, const double *BQ_R__ ri, const double *BQ_R__ ii, double *BQ_R__ ro, double *BQ_R__ io)
Chris@29: {
Chris@29:     if (!ri || !ro || !io) return;
Chris@29: 
Chris@29:     int i, j, k, m;
Chris@29:     int blockSize, blockEnd;
Chris@29: 
Chris@29:     double tr, ti;
Chris@29: 
Chris@29:     double angle = 2.0 * M_PI;
Chris@29:     if (inverse) angle = -angle;
Chris@29: 
Chris@29:     const int n = m_size;
Chris@29: 
Chris@29:     if (ii) {
Chris@29: 	for (i = 0; i < n; ++i) {
Chris@29: 	    ro[m_table[i]] = ri[i];
Chris@29:         }
Chris@29: 	for (i = 0; i < n; ++i) {
Chris@29: 	    io[m_table[i]] = ii[i];
Chris@29: 	}
Chris@29:     } else {
Chris@29: 	for (i = 0; i < n; ++i) {
Chris@29: 	    ro[m_table[i]] = ri[i];
Chris@29:         }
Chris@29: 	for (i = 0; i < n; ++i) {
Chris@29: 	    io[m_table[i]] = 0.0;
Chris@29: 	}
Chris@29:     }
Chris@29: 
Chris@29:     blockEnd = 1;
Chris@29: 
Chris@29:     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
Chris@29: 
Chris@29: 	double delta = angle / (double)blockSize;
Chris@29: 	double sm2 = -sin(-2 * delta);
Chris@29: 	double sm1 = -sin(-delta);
Chris@29: 	double cm2 = cos(-2 * delta);
Chris@29: 	double cm1 = cos(-delta);
Chris@29: 	double w = 2 * cm1;
Chris@29: 	double ar[3], ai[3];
Chris@29: 
Chris@29: 	for (i = 0; i < n; i += blockSize) {
Chris@29: 
Chris@29: 	    ar[2] = cm2;
Chris@29: 	    ar[1] = cm1;
Chris@29: 
Chris@29: 	    ai[2] = sm2;
Chris@29: 	    ai[1] = sm1;
Chris@29: 
Chris@29: 	    for (j = i, m = 0; m < blockEnd; j++, m++) {
Chris@29: 
Chris@29: 		ar[0] = w * ar[1] - ar[2];
Chris@29: 		ar[2] = ar[1];
Chris@29: 		ar[1] = ar[0];
Chris@29: 
Chris@29: 		ai[0] = w * ai[1] - ai[2];
Chris@29: 		ai[2] = ai[1];
Chris@29: 		ai[1] = ai[0];
Chris@29: 
Chris@29: 		k = j + blockEnd;
Chris@29: 		tr = ar[0] * ro[k] - ai[0] * io[k];
Chris@29: 		ti = ar[0] * io[k] + ai[0] * ro[k];
Chris@29: 
Chris@29: 		ro[k] = ro[j] - tr;
Chris@29: 		io[k] = io[j] - ti;
Chris@29: 
Chris@29: 		ro[j] += tr;
Chris@29: 		io[j] += ti;
Chris@29: 	    }
Chris@29: 	}
Chris@29: 
Chris@29: 	blockEnd = blockSize;
Chris@29:     }
Chris@29: 
Chris@29: /* fftw doesn't rescale, so nor will we
Chris@29: 
Chris@29:     if (inverse) {
Chris@29: 
Chris@29: 	double denom = (double)n;
Chris@29: 
Chris@29: 	for (i = 0; i < n; i++) {
Chris@29: 	    ro[i] /= denom;
Chris@29: 	    io[i] /= denom;
Chris@29: 	}
Chris@29:     }
Chris@29: */
Chris@29: }
Chris@29: 
Chris@29: #endif /* USE_BUILTIN_FFT */
Chris@29: 
Chris@29: } /* end namespace FFTs */
Chris@29: 
Chris@29: std::string
Chris@29: FFT::m_implementation;
Chris@29: 
Chris@29: std::set<std::string>
Chris@29: FFT::getImplementations()
Chris@29: {
Chris@29:     std::set<std::string> impls;
Chris@29: #ifdef HAVE_IPP
Chris@29:     impls.insert("ipp");
Chris@29: #endif
Chris@29: #ifdef HAVE_FFTW3
Chris@29:     impls.insert("fftw");
Chris@29: #endif
Chris@29: #ifdef HAVE_KISSFFT
Chris@29:     impls.insert("kissfft");
Chris@29: #endif
Chris@29: #ifdef HAVE_VDSP
Chris@29:     impls.insert("vdsp");
Chris@29: #endif
Chris@29: #ifdef HAVE_MEDIALIB
Chris@29:     impls.insert("medialib");
Chris@29: #endif
Chris@29: #ifdef HAVE_OPENMAX
Chris@29:     impls.insert("openmax");
Chris@29: #endif
Chris@29: #ifdef HAVE_SFFT
Chris@29:     impls.insert("sfft");
Chris@29: #endif
Chris@29: #ifdef USE_BUILTIN_FFT
Chris@29:     impls.insert("cross");
Chris@29: #endif
Chris@29:     return impls;
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::pickDefaultImplementation()
Chris@29: {
Chris@29:     if (m_implementation != "") return;
Chris@29: 
Chris@29:     std::set<std::string> impls = getImplementations();
Chris@29: 
Chris@29:     std::string best = "cross";
Chris@29:     if (impls.find("kissfft") != impls.end()) best = "kissfft";
Chris@29:     if (impls.find("medialib") != impls.end()) best = "medialib";
Chris@29:     if (impls.find("openmax") != impls.end()) best = "openmax";
Chris@29:     if (impls.find("sfft") != impls.end()) best = "sfft";
Chris@29:     if (impls.find("fftw") != impls.end()) best = "fftw";
Chris@29:     if (impls.find("vdsp") != impls.end()) best = "vdsp";
Chris@29:     if (impls.find("ipp") != impls.end()) best = "ipp";
Chris@29:     
Chris@29:     m_implementation = best;
Chris@29: }
Chris@29: 
Chris@29: std::string
Chris@29: FFT::getDefaultImplementation()
Chris@29: {
Chris@29:     return m_implementation;
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::setDefaultImplementation(std::string i)
Chris@29: {
Chris@29:     m_implementation = i;
Chris@29: }
Chris@29: 
Chris@29: FFT::FFT(int size, int debugLevel) :
Chris@29:     d(0)
Chris@29: {
Chris@29:     if ((size < 2) ||
Chris@29:         (size & (size-1))) {
Chris@29:         std::cerr << "FFT::FFT(" << size << "): power-of-two sizes only supported, minimum size 2" << std::endl;
Chris@29: #ifndef NO_EXCEPTIONS
Chris@29:         throw InvalidSize;
Chris@29: #else
Chris@29:         abort();
Chris@29: #endif
Chris@29:     }
Chris@29: 
Chris@29:     if (m_implementation == "") pickDefaultImplementation();
Chris@29:     std::string impl = m_implementation;
Chris@29: 
Chris@29:     if (debugLevel > 0) {
Chris@29:         std::cerr << "FFT::FFT(" << size << "): using implementation: "
Chris@29:                   << impl << std::endl;
Chris@29:     }
Chris@29: 
Chris@29:     if (impl == "ipp") {
Chris@29: #ifdef HAVE_IPP
Chris@29:         d = new FFTs::D_IPP(size);
Chris@29: #endif
Chris@29:     } else if (impl == "fftw") {
Chris@29: #ifdef HAVE_FFTW3
Chris@29:         d = new FFTs::D_FFTW(size);
Chris@29: #endif
Chris@29:     } else if (impl == "kissfft") {        
Chris@29: #ifdef HAVE_KISSFFT
Chris@29:         d = new FFTs::D_KISSFFT(size);
Chris@29: #endif
Chris@29:     } else if (impl == "vdsp") {
Chris@29: #ifdef HAVE_VDSP
Chris@29:         d = new FFTs::D_VDSP(size);
Chris@29: #endif
Chris@29:     } else if (impl == "medialib") {
Chris@29: #ifdef HAVE_MEDIALIB
Chris@29:         d = new FFTs::D_MEDIALIB(size);
Chris@29: #endif
Chris@29:     } else if (impl == "openmax") {
Chris@29: #ifdef HAVE_OPENMAX
Chris@29:         d = new FFTs::D_OPENMAX(size);
Chris@29: #endif
Chris@29:     } else if (impl == "sfft") {
Chris@29: #ifdef HAVE_SFFT
Chris@29:         d = new FFTs::D_SFFT(size);
Chris@29: #endif
Chris@29:     } else if (impl == "cross") {
Chris@29: #ifdef USE_BUILTIN_FFT
Chris@29:         d = new FFTs::D_Cross(size);
Chris@29: #endif
Chris@29:     }
Chris@29: 
Chris@29:     if (!d) {
Chris@29:         std::cerr << "FFT::FFT(" << size << "): ERROR: implementation "
Chris@29:                   << impl << " is not compiled in" << std::endl;
Chris@29: #ifndef NO_EXCEPTIONS
Chris@29:         throw InvalidImplementation;
Chris@29: #else
Chris@29:         abort();
Chris@29: #endif
Chris@29:     }
Chris@29: }
Chris@29: 
Chris@29: FFT::~FFT()
Chris@29: {
Chris@29:     delete d;
Chris@29: }
Chris@29: 
Chris@29: #ifndef NO_EXCEPTIONS
Chris@29: #define CHECK_NOT_NULL(x) \
Chris@29:     if (!(x)) { \
Chris@29:         std::cerr << "FFT: ERROR: Null argument " #x << std::endl;  \
Chris@29:         throw NullArgument; \
Chris@29:     }
Chris@29: #else
Chris@29: #define CHECK_NOT_NULL(x) \
Chris@29:     if (!(x)) { \
Chris@29:         std::cerr << "FFT: ERROR: Null argument " #x << std::endl;  \
Chris@29:         std::cerr << "FFT: Would be throwing NullArgument here, if exceptions were not disabled" << std::endl;  \
Chris@29:         return; \
Chris@29:     }
Chris@29: #endif
Chris@29: 
Chris@29: void
Chris@29: FFT::forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     CHECK_NOT_NULL(imagOut);
Chris@29:     d->forward(realIn, realOut, imagOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(complexOut);
Chris@29:     d->forwardInterleaved(realIn, complexOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(magOut);
Chris@29:     CHECK_NOT_NULL(phaseOut);
Chris@29:     d->forwardPolar(realIn, magOut, phaseOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(magOut);
Chris@29:     d->forwardMagnitude(realIn, magOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     CHECK_NOT_NULL(imagOut);
Chris@29:     d->forward(realIn, realOut, imagOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(complexOut);
Chris@29:     d->forwardInterleaved(realIn, complexOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(magOut);
Chris@29:     CHECK_NOT_NULL(phaseOut);
Chris@29:     d->forwardPolar(realIn, magOut, phaseOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(magOut);
Chris@29:     d->forwardMagnitude(realIn, magOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(imagIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inverse(realIn, imagIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(complexIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inverseInterleaved(complexIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(magIn);
Chris@29:     CHECK_NOT_NULL(phaseIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inversePolar(magIn, phaseIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(magIn);
Chris@29:     CHECK_NOT_NULL(cepOut);
Chris@29:     d->inverseCepstral(magIn, cepOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(realIn);
Chris@29:     CHECK_NOT_NULL(imagIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inverse(realIn, imagIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(complexIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inverseInterleaved(complexIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(magIn);
Chris@29:     CHECK_NOT_NULL(phaseIn);
Chris@29:     CHECK_NOT_NULL(realOut);
Chris@29:     d->inversePolar(magIn, phaseIn, realOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut)
Chris@29: {
Chris@29:     CHECK_NOT_NULL(magIn);
Chris@29:     CHECK_NOT_NULL(cepOut);
Chris@29:     d->inverseCepstral(magIn, cepOut);
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::initFloat() 
Chris@29: {
Chris@29:     d->initFloat();
Chris@29: }
Chris@29: 
Chris@29: void
Chris@29: FFT::initDouble() 
Chris@29: {
Chris@29:     d->initDouble();
Chris@29: }
Chris@29: 
Chris@29: FFT::Precisions
Chris@29: FFT::getSupportedPrecisions() const
Chris@29: {
Chris@29:     return d->getSupportedPrecisions();
Chris@29: }
Chris@29: 
Chris@29: #ifdef FFT_MEASUREMENT
Chris@29: 
Chris@29: std::string
Chris@29: FFT::tune()
Chris@29: {
Chris@29:     std::ostringstream os;
Chris@29:     os << "FFT::tune()..." << std::endl;
Chris@29: 
Chris@29:     std::vector<int> sizes;
Chris@29:     std::map<FFTImpl *, int> candidates;
Chris@29:     std::map<int, int> wins;
Chris@29: 
Chris@29:     sizes.push_back(512);
Chris@29:     sizes.push_back(1024);
Chris@29:     sizes.push_back(4096);
Chris@29:     
Chris@29:     for (unsigned int si = 0; si < sizes.size(); ++si) {
Chris@29: 
Chris@29:         int size = sizes[si];
Chris@29: 
Chris@29:         while (!candidates.empty()) {
Chris@29:             delete candidates.begin()->first;
Chris@29:             candidates.erase(candidates.begin());
Chris@29:         }
Chris@29: 
Chris@29:         FFTImpl *d;
Chris@29:         
Chris@29: #ifdef HAVE_IPP
Chris@29:         std::cerr << "Constructing new IPP FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_IPP(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 0;
Chris@29: #endif
Chris@29:         
Chris@29: #ifdef HAVE_FFTW3
Chris@29:         os << "Constructing new FFTW3 FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_FFTW(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 1;
Chris@29: #endif
Chris@29: 
Chris@29: #ifdef HAVE_KISSFFT
Chris@29:         os << "Constructing new KISSFFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_KISSFFT(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 2;
Chris@29: #endif        
Chris@29: 
Chris@29: #ifdef USE_BUILTIN_FFT
Chris@29:         os << "Constructing new Cross FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_Cross(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 3;
Chris@29: #endif
Chris@29:         
Chris@29: #ifdef HAVE_VDSP
Chris@29:         os << "Constructing new vDSP FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_VDSP(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 4;
Chris@29: #endif
Chris@29:         
Chris@29: #ifdef HAVE_MEDIALIB
Chris@29:         std::cerr << "Constructing new MediaLib FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_MEDIALIB(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 5;
Chris@29: #endif
Chris@29:         
Chris@29: #ifdef HAVE_OPENMAX
Chris@29:         os << "Constructing new OpenMAX FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_OPENMAX(size);
Chris@29:         d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 6;
Chris@29: #endif
Chris@29:         
Chris@29: #ifdef HAVE_SFFT
Chris@29:         os << "Constructing new SFFT FFT object for size " << size << "..." << std::endl;
Chris@29:         d = new FFTs::D_SFFT(size);
Chris@29: //        d->initFloat();
Chris@29:         d->initDouble();
Chris@29:         candidates[d] = 6;
Chris@29: #endif
Chris@29: 
Chris@29:         os << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << std::endl;
Chris@29:         float divisor = float(CLOCKS_PER_SEC) / 1000.f;
Chris@29:         
Chris@29:         os << "Timing order is: ";
Chris@29:         for (std::map<FFTImpl *, int>::iterator ci = candidates.begin();
Chris@29:              ci != candidates.end(); ++ci) {
Chris@29:             os << ci->second << " ";
Chris@29:         }
Chris@29:         os << std::endl;
Chris@29: 
Chris@29:         int iterations = 500;
Chris@29:         os << "Iterations: " << iterations << std::endl;
Chris@29: 
Chris@29:         double *da = new double[size];
Chris@29:         double *db = new double[size];
Chris@29:         double *dc = new double[size];
Chris@29:         double *dd = new double[size];
Chris@29:         double *di = new double[size + 2];
Chris@29:         double *dj = new double[size + 2];
Chris@29: 
Chris@29:         float *fa = new float[size];
Chris@29:         float *fb = new float[size];
Chris@29:         float *fc = new float[size];
Chris@29:         float *fd = new float[size];
Chris@29:         float *fi = new float[size + 2];
Chris@29:         float *fj = new float[size + 2];
Chris@29: 
Chris@29:         for (int type = 0; type < 16; ++type) {
Chris@29:     
Chris@29:             //!!!
Chris@29:             if ((type > 3 && type < 8) ||
Chris@29:                 (type > 11)) {
Chris@29:                 continue;
Chris@29:             }
Chris@29: 
Chris@29:             if (type > 7) {
Chris@29:                 // inverse transform: bigger inputs, to simulate the
Chris@29:                 // fact that the forward transform is unscaled
Chris@29:                 for (int i = 0; i < size; ++i) {
Chris@29:                     da[i] = drand48() * size;
Chris@29:                     fa[i] = da[i];
Chris@29:                     db[i] = drand48() * size;
Chris@29:                     fb[i] = db[i];
Chris@29:                 }
Chris@29:             } else {    
Chris@29:                 for (int i = 0; i < size; ++i) {
Chris@29:                     da[i] = drand48();
Chris@29:                     fa[i] = da[i];
Chris@29:                     db[i] = drand48();
Chris@29:                     fb[i] = db[i];
Chris@29:                 }
Chris@29:             }
Chris@29:                 
Chris@29:             for (int i = 0; i < size + 2; ++i) {
Chris@29:                 di[i] = drand48();
Chris@29:                 fi[i] = di[i];
Chris@29:             }
Chris@29: 
Chris@29:             int low = -1;
Chris@29:             int lowscore = 0;
Chris@29: 
Chris@29:             const char *names[] = {
Chris@29: 
Chris@29:                 "Forward Cartesian Double",
Chris@29:                 "Forward Interleaved Double",
Chris@29:                 "Forward Polar Double",
Chris@29:                 "Forward Magnitude Double",
Chris@29:                 "Forward Cartesian Float",
Chris@29:                 "Forward Interleaved Float",
Chris@29:                 "Forward Polar Float",
Chris@29:                 "Forward Magnitude Float",
Chris@29: 
Chris@29:                 "Inverse Cartesian Double",
Chris@29:                 "Inverse Interleaved Double",
Chris@29:                 "Inverse Polar Double",
Chris@29:                 "Inverse Cepstral Double",
Chris@29:                 "Inverse Cartesian Float",
Chris@29:                 "Inverse Interleaved Float",
Chris@29:                 "Inverse Polar Float",
Chris@29:                 "Inverse Cepstral Float"
Chris@29:             };
Chris@29:             os << names[type] << " :: ";
Chris@29: 
Chris@29:             for (std::map<FFTImpl *, int>::iterator ci = candidates.begin();
Chris@29:                  ci != candidates.end(); ++ci) {
Chris@29: 
Chris@29:                 FFTImpl *d = ci->first;
Chris@29: 
Chris@29:                 double mean = 0;
Chris@29: 
Chris@29:                 clock_t start = clock();
Chris@29:                 
Chris@29:                 for (int i = 0; i < iterations; ++i) {
Chris@29: 
Chris@29:                     if (i == 0) {
Chris@29:                         for (int j = 0; j < size; ++j) {
Chris@29:                             dc[j] = 0;
Chris@29:                             dd[j] = 0;
Chris@29:                             fc[j] = 0;
Chris@29:                             fd[j] = 0;
Chris@29:                             fj[j] = 0;
Chris@29:                             dj[j] = 0;
Chris@29:                         }
Chris@29:                     }
Chris@29: 
Chris@29:                     switch (type) {
Chris@29:                     case 0: d->forward(da, dc, dd); break;
Chris@29:                     case 1: d->forwardInterleaved(da, dj); break;
Chris@29:                     case 2: d->forwardPolar(da, dc, dd); break;
Chris@29:                     case 3: d->forwardMagnitude(da, dc); break;
Chris@29:                     case 4: d->forward(fa, fc, fd); break;
Chris@29:                     case 5: d->forwardInterleaved(fa, fj); break;
Chris@29:                     case 6: d->forwardPolar(fa, fc, fd); break;
Chris@29:                     case 7: d->forwardMagnitude(fa, fc); break;
Chris@29:                     case 8: d->inverse(da, db, dc); break;
Chris@29:                     case 9: d->inverseInterleaved(di, dc); break;
Chris@29:                     case 10: d->inversePolar(da, db, dc); break;
Chris@29:                     case 11: d->inverseCepstral(da, dc); break;
Chris@29:                     case 12: d->inverse(fa, fb, fc); break;
Chris@29:                     case 13: d->inverseInterleaved(fi, fc); break;
Chris@29:                     case 14: d->inversePolar(fa, fb, fc); break;
Chris@29:                     case 15: d->inverseCepstral(fa, fc); break;
Chris@29:                     }
Chris@29: 
Chris@29:                     if (i == 0) {
Chris@29:                         mean = 0;
Chris@29:                         for (int j = 0; j < size; ++j) {
Chris@29:                             mean += dc[j];
Chris@29:                             mean += dd[j];
Chris@29:                             mean += fc[j];
Chris@29:                             mean += fd[j];
Chris@29:                             mean += fj[j];
Chris@29:                             mean += dj[j];
Chris@29:                         }
Chris@29:                         mean /= size * 6;
Chris@29:                     }
Chris@29:                 }
Chris@29: 
Chris@29:                 clock_t end = clock();
Chris@29: 
Chris@29:                 os << float(end - start)/divisor << " (" << mean << ") ";
Chris@29: 
Chris@29:                 if (low == -1 || (end - start) < lowscore) {
Chris@29:                     low = ci->second;
Chris@29:                     lowscore = end - start;
Chris@29:                 }
Chris@29:             }
Chris@29: 
Chris@29:             os << std::endl;
Chris@29: 
Chris@29:             os << "  size " << size << ", type " << type << ": fastest is " << low << " (time " << float(lowscore)/divisor << ")" << std::endl;
Chris@29: 
Chris@29:             wins[low]++;
Chris@29:         }
Chris@29:         
Chris@29:         delete[] fa;
Chris@29:         delete[] fb;
Chris@29:         delete[] fc;
Chris@29:         delete[] fd;
Chris@29:         delete[] da;
Chris@29:         delete[] db;
Chris@29:         delete[] dc;
Chris@29:         delete[] dd;
Chris@29:     }
Chris@29: 
Chris@29:     while (!candidates.empty()) {
Chris@29:         delete candidates.begin()->first;
Chris@29:         candidates.erase(candidates.begin());
Chris@29:     }
Chris@29: 
Chris@29:     int bestscore = 0;
Chris@29:     int best = -1;
Chris@29: 
Chris@29:     for (std::map<int, int>::iterator wi = wins.begin(); wi != wins.end(); ++wi) {
Chris@29:         if (best == -1 || wi->second > bestscore) {
Chris@29:             best = wi->first;
Chris@29:             bestscore = wi->second;
Chris@29:         }
Chris@29:     }
Chris@29: 
Chris@29:     os << "overall winner is " << best << " with " << bestscore << " wins" << std::endl;
Chris@29: 
Chris@29:     return os.str();
Chris@29: }
Chris@29: 
Chris@29: #endif
Chris@29: 
Chris@29: }