8b5b1108142f — Chris Cannam 8 days ago
Merge
11 files changed, 907 insertions(+), 292 deletions(-)

M .hgignore
A => .hgtags
M COPYING
M Makefile
M README.md
M bqfft/FFT.h
M build/Makefile.inc
M build/run-platform-tests.sh
M src/FFT.cpp
M test/TestFFT.cpp
A => test/timings.cpp
M .hgignore +1 -0
@@ 8,6 8,7 @@ syntax: glob
 *.lib
 *.bak
 test-fft
+timings
 moc_*
 test/Makefile
 

          
A => .hgtags +1 -0
@@ 0,0 1,1 @@ 
+a766fe47501b185bc46cffc210735304e28f2189 v1.0.0

          
M COPYING +1 -1
@@ 1,5 1,5 @@ 
 
-    Copyright 2007-2015 Particular Programs Ltd.
+    Copyright 2007-2019 Particular Programs Ltd.
 
     Permission is hereby granted, free of charge, to any person
     obtaining a copy of this software and associated documentation

          
M Makefile +10 -0
@@ 36,4 36,14 @@ THIRD_PARTY_INCLUDES	:=
 THIRD_PARTY_LIBS	:=
 
 
+# If you are including a set of bq libraries into a project, you can
+# override variables for all of them (including all of the above) in
+# the following file, which all bq* Makefiles will include if found
+
+-include ../Makefile.inc-bq
+
+
+# This project-local Makefile describes the source files and contains
+# no routinely user-modifiable parts
+
 include build/Makefile.inc

          
M README.md +23 -4
@@ 4,22 4,37 @@ bqfft
 
 A small library wrapping various FFT implementations for some common
 audio processing use cases. Note this is not a general FFT interface,
-as it handles only power-of-two FFT sizes and real inputs.
+as it handles only real inputs.
+
+Transforms of any length are supported, but if you request a length
+that bqfft does not know how to calculate using any of the libraries
+that have been compiled in, a simple slow DFT will be used instead. A
+warning will be printed to stderr if this happens.
+
+Of the most commonly used libraries, Accelerate and IPP support
+power-of-two FFT lengths only, KissFFT supports any multiple of two
+(because we only use the real-input interface), and FFTW supports any
+length. You can compile in more than one library, so for example if
+you compile in Accelerate and KissFFT, the former will be used for
+powers of two and the latter for other even lengths.
 
 Requires the bqvec library.
 
 This code originated as part of the Rubber Band Library written by the
-same authors (see https://bitbucket.org/breakfastquay/rubberband/).
+same authors (see https://hg.sr.ht/~breakfastquay/rubberband/).
 It has been pulled out into a separate library and relicensed under a
 more permissive licence.
 
 C++ standard required: C++98 (does not use C++11 or newer features)
 
- * To compile on Linux: Edit Makefile to select implementation, then make test.
+ * To compile on Linux: Edit Makefile to select implementation, then make.
    Do read the notes in the Makefile, and don't attempt to use the default
    implementation, which is very slow
+   
+ * To compile on macOS: make -f build/Makefile.osx
 
- * To compile on macOS: make -f build/Makefile.osx test
+ * To build and run tests: as above, but add the "test" target -
+   requires Boost test headers installed
 
  * Depends on: [bqvec](https://hg.sr.ht/~breakfastquay/bqvec)
 

          
@@ 27,5 42,9 @@ C++ standard required: C++98 (does not u
 
 [![Build Status](https://travis-ci.org/breakfastquay/bqfft.svg?branch=master)](https://travis-ci.org/breakfastquay/bqfft)
 
+<<<<<<< working copy
 Copyright 2007-2020 Particular Programs Ltd. See the file COPYING for
+=======
+Copyright 2007-2019 Particular Programs Ltd. See the file COPYING for
+>>>>>>> merge rev
 (BSD/MIT-style) licence terms.

          
M bqfft/FFT.h +7 -3
@@ 128,13 128,17 @@ public:
     static void setDefaultImplementation(std::string);
 
 #ifdef FFT_MEASUREMENT
-    static std::string tune();
+    static
+#ifdef FFT_MEASUREMENT_RETURN_RESULT_TEXT
+    std::string
+#else
+    void
+#endif
+    tune();
 #endif
 
 protected:
     FFTImpl *d;
-    static std::string m_implementation;
-    static void pickDefaultImplementation();
 
 private:
     FFT(const FFT &); // not provided

          
M build/Makefile.inc +4 -0
@@ 32,6 32,9 @@ valgrind:	$(LIBRARY) test-fft
 test-fft:	test/TestFFT.o $(LIBRARY)
 	$(CXX) $(CXXFLAGS) -o $@ $^ -lboost_unit_test_framework -L. -lbqfft -L../bqvec -lbqvec $(THIRD_PARTY_LIBS)
 
+timings:       test/timings.o
+	$(CXX) $(CXXFLAGS) -o $@ $^
+
 clean:		
 	rm -f $(OBJECTS) $(TEST_OBJECTS)
 

          
@@ 46,3 49,4 @@ depend:
 
 src/FFT.o: bqfft/FFT.h
 test/TestFFT.o: bqfft/FFT.h
+test/timings.o: src/FFT.cpp bqfft/FFT.h

          
M build/run-platform-tests.sh +2 -2
@@ 81,8 81,8 @@ for mf in Makefile build/Makefile.$platf
 
     if [ "$have_valgrind" = "yes" ]; then
 	for t in test-* ; do
-	    if [ -x "$t" ]; then
-		run "no leaks are possible" valgrind --leak-check=full ./"$t"
+	    if [ -f "$t" -a -x "$t" ]; then
+		run "0 errors from 0 contexts" valgrind --leak-check=full ./"$t"
 	    fi
 	done
     fi

          
M src/FFT.cpp +385 -63
@@ 6,7 6,7 @@ 
     A small library wrapping various FFT implementations for some
     common audio processing use cases.
 
-    Copyright 2007-2016 Particular Programs Ltd.
+    Copyright 2007-2019 Particular Programs Ltd.
 
     Permission is hereby granted, free of charge, to any person
     obtaining a copy of this software and associated documentation

          
@@ 1703,7 1703,7 @@ public:
     void saveWisdom(char type) { wisdom(true, type); }
 
     void wisdom(bool save, char type) {
-
+#ifdef USE_FFTW_WISDOM
 #ifdef FFTW_DOUBLE_ONLY
         if (type == 'f') return;
 #endif

          
@@ 1751,6 1751,10 @@ public:
         }
 
         fclose(f);
+#else
+        (void)save;
+        (void)type;
+#endif
     }
 
     void packFloat(const float *BQ_R__ re, const float *BQ_R__ im) {

          
@@ 2999,88 3003,388 @@ D_Cross::basefft(bool inverse, const dou
 
 #endif /* USE_BUILTIN_FFT */
 
+class D_DFT : public FFTImpl
+{
+private:
+    template <typename T>
+    class DFT
+    {
+    public:
+        DFT(int size) : m_size(size), m_bins(size/2 + 1) {
+
+            m_sin = allocate_channels<double>(m_size, m_size);
+            m_cos = allocate_channels<double>(m_size, m_size);
+
+            for (int i = 0; i < m_size; ++i) {
+                for (int j = 0; j < m_size; ++j) {
+                    double arg = (double(i) * double(j) * M_PI * 2.0) / m_size;
+                    m_sin[i][j] = sin(arg);
+                    m_cos[i][j] = cos(arg);
+                }
+            }
+
+            m_tmp = allocate_channels<double>(2, m_size);
+        }
+
+        ~DFT() {
+            deallocate_channels(m_tmp, 2);
+            deallocate_channels(m_sin, m_size);
+            deallocate_channels(m_cos, m_size);
+        }
+
+        void forward(const T *BQ_R__ realIn, T *BQ_R__ realOut, T *BQ_R__ imagOut) {
+            for (int i = 0; i < m_bins; ++i) {
+                double re = 0.0, im = 0.0;
+                for (int j = 0; j < m_size; ++j) re += realIn[j] * m_cos[i][j];
+                for (int j = 0; j < m_size; ++j) im -= realIn[j] * m_sin[i][j];
+                realOut[i] = T(re);
+                imagOut[i] = T(im);
+            }
+        }
+
+        void forwardInterleaved(const T *BQ_R__ realIn, T *BQ_R__ complexOut) {
+            for (int i = 0; i < m_bins; ++i) {
+                double re = 0.0, im = 0.0;
+                for (int j = 0; j < m_size; ++j) re += realIn[j] * m_cos[i][j];
+                for (int j = 0; j < m_size; ++j) im -= realIn[j] * m_sin[i][j];
+                complexOut[i*2] = T(re);
+                complexOut[i*2 + 1] = T(im);
+            }
+        }
+
+        void forwardPolar(const T *BQ_R__ realIn, T *BQ_R__ magOut, T *BQ_R__ phaseOut) {
+            forward(realIn, magOut, phaseOut); // temporarily
+            for (int i = 0; i < m_bins; ++i) {
+                T re = magOut[i], im = phaseOut[i];
+                c_magphase(magOut + i, phaseOut + i, re, im);
+            }
+        }
+
+        void forwardMagnitude(const T *BQ_R__ realIn, T *BQ_R__ magOut) {
+            for (int i = 0; i < m_bins; ++i) {
+                double re = 0.0, im = 0.0;
+                for (int j = 0; j < m_size; ++j) re += realIn[j] * m_cos[i][j];
+                for (int j = 0; j < m_size; ++j) im -= realIn[j] * m_sin[i][j];
+                magOut[i] = T(sqrt(re * re + im * im));
+            }
+        }
+
+        void inverse(const T *BQ_R__ realIn, const T *BQ_R__ imagIn, T *BQ_R__ realOut) {
+            for (int i = 0; i < m_bins; ++i) {
+                m_tmp[0][i] = realIn[i];
+                m_tmp[1][i] = imagIn[i];
+            }
+            for (int i = m_bins; i < m_size; ++i) {
+                m_tmp[0][i] = realIn[m_size - i];
+                m_tmp[1][i] = -imagIn[m_size - i];
+            }
+            for (int i = 0; i < m_size; ++i) {
+                double re = 0.0;
+                const double *const cos = m_cos[i];
+                const double *const sin = m_sin[i];
+                for (int j = 0; j < m_size; ++j) re += m_tmp[0][j] * cos[j];
+                for (int j = 0; j < m_size; ++j) re -= m_tmp[1][j] * sin[j];
+                realOut[i] = T(re);
+            }
+        }
+
+        void inverseInterleaved(const T *BQ_R__ complexIn, T *BQ_R__ realOut) {
+            for (int i = 0; i < m_bins; ++i) {
+                m_tmp[0][i] = complexIn[i*2];
+                m_tmp[1][i] = complexIn[i*2+1];
+            }
+            for (int i = m_bins; i < m_size; ++i) {
+                m_tmp[0][i] = complexIn[(m_size - i) * 2];
+                m_tmp[1][i] = -complexIn[(m_size - i) * 2 + 1];
+            }
+            for (int i = 0; i < m_size; ++i) {
+                double re = 0.0;
+                const double *const cos = m_cos[i];
+                const double *const sin = m_sin[i];
+                for (int j = 0; j < m_size; ++j) re += m_tmp[0][j] * cos[j];
+                for (int j = 0; j < m_size; ++j) re -= m_tmp[1][j] * sin[j];
+                realOut[i] = T(re);
+            }
+        }
+
+        void inversePolar(const T *BQ_R__ magIn, const T *BQ_R__ phaseIn, T *BQ_R__ realOut) {
+            T *complexIn = allocate<T>(m_bins * 2);
+            v_polar_to_cartesian_interleaved(complexIn, magIn, phaseIn, m_bins);
+            inverseInterleaved(complexIn, realOut);
+            deallocate(complexIn);
+        }
+
+        void inverseCepstral(const T *BQ_R__ magIn, T *BQ_R__ cepOut) {
+            T *complexIn = allocate_and_zero<T>(m_bins * 2);
+            for (int i = 0; i < m_bins; ++i) {
+                complexIn[i*2] = T(log(magIn[i] + 0.000001));
+            }
+            inverseInterleaved(complexIn, cepOut);
+            deallocate(complexIn);
+        }
+
+    private:
+        const int m_size;
+        const int m_bins;
+        double **m_sin;
+        double **m_cos;
+        double **m_tmp;
+    };
+    
+public:
+    D_DFT(int size) : m_size(size), m_double(0), m_float(0) { }
+
+    ~D_DFT() {
+        delete m_double;
+        delete m_float;
+    }
+
+    int getSize() const {
+        return m_size;
+    }
+
+    FFT::Precisions
+    getSupportedPrecisions() const {
+        return FFT::DoublePrecision;
+    }
+
+    void initFloat() {
+        if (!m_float) {
+            m_float = new DFT<float>(m_size);
+        }
+    }
+        
+    void initDouble() {
+        if (!m_double) {
+            m_double = new DFT<double>(m_size);
+        }
+    }
+
+    void forward(const double *BQ_R__ realIn, double *BQ_R__ realOut, double *BQ_R__ imagOut) {
+        initDouble();
+        m_double->forward(realIn, realOut, imagOut);
+    }
+
+    void forwardInterleaved(const double *BQ_R__ realIn, double *BQ_R__ complexOut) {
+        initDouble();
+        m_double->forwardInterleaved(realIn, complexOut);
+    }
+
+    void forwardPolar(const double *BQ_R__ realIn, double *BQ_R__ magOut, double *BQ_R__ phaseOut) {
+        initDouble();
+        m_double->forwardPolar(realIn, magOut, phaseOut);
+    }
+
+    void forwardMagnitude(const double *BQ_R__ realIn, double *BQ_R__ magOut) {
+        initDouble();
+        m_double->forwardMagnitude(realIn, magOut);
+    }
+
+    void forward(const float *BQ_R__ realIn, float *BQ_R__ realOut, float *BQ_R__ imagOut) {
+        initFloat();
+        m_float->forward(realIn, realOut, imagOut);
+    }
+
+    void forwardInterleaved(const float *BQ_R__ realIn, float *BQ_R__ complexOut) {
+        initFloat();
+        m_float->forwardInterleaved(realIn, complexOut);
+    }
+
+    void forwardPolar(const float *BQ_R__ realIn, float *BQ_R__ magOut, float *BQ_R__ phaseOut) {
+        initFloat();
+        m_float->forwardPolar(realIn, magOut, phaseOut);
+    }
+
+    void forwardMagnitude(const float *BQ_R__ realIn, float *BQ_R__ magOut) {
+        initFloat();
+        m_float->forwardMagnitude(realIn, magOut);
+    }
+
+    void inverse(const double *BQ_R__ realIn, const double *BQ_R__ imagIn, double *BQ_R__ realOut) {
+        initDouble();
+        m_double->inverse(realIn, imagIn, realOut);
+    }
+
+    void inverseInterleaved(const double *BQ_R__ complexIn, double *BQ_R__ realOut) {
+        initDouble();
+        m_double->inverseInterleaved(complexIn, realOut);
+    }
+
+    void inversePolar(const double *BQ_R__ magIn, const double *BQ_R__ phaseIn, double *BQ_R__ realOut) {
+        initDouble();
+        m_double->inversePolar(magIn, phaseIn, realOut);
+    }
+
+    void inverseCepstral(const double *BQ_R__ magIn, double *BQ_R__ cepOut) {
+        initDouble();
+        m_double->inverseCepstral(magIn, cepOut);
+    }
+
+    void inverse(const float *BQ_R__ realIn, const float *BQ_R__ imagIn, float *BQ_R__ realOut) {
+        initFloat();
+        m_float->inverse(realIn, imagIn, realOut);
+    }
+
+    void inverseInterleaved(const float *BQ_R__ complexIn, float *BQ_R__ realOut) {
+        initFloat();
+        m_float->inverseInterleaved(complexIn, realOut);
+    }
+
+    void inversePolar(const float *BQ_R__ magIn, const float *BQ_R__ phaseIn, float *BQ_R__ realOut) {
+        initFloat();
+        m_float->inversePolar(magIn, phaseIn, realOut);
+    }
+
+    void inverseCepstral(const float *BQ_R__ magIn, float *BQ_R__ cepOut) {
+        initFloat();
+        m_float->inverseCepstral(magIn, cepOut);
+    }
+
+private:
+    int m_size;
+    DFT<double> *m_double;
+    DFT<float> *m_float;
+};
+
 } /* end namespace FFTs */
 
-std::string
-FFT::m_implementation;
+enum SizeConstraint {
+    SizeConstraintNone           = 0x0,
+    SizeConstraintEven           = 0x1,
+    SizeConstraintPowerOfTwo     = 0x2,
+    SizeConstraintEvenPowerOfTwo = 0x3 // i.e. 0x1 | 0x2. Excludes size 1 obvs
+};
+
+typedef std::map<std::string, SizeConstraint> ImplMap;
+
+static std::string defaultImplementation;
+
+static ImplMap
+getImplementationDetails()
+{
+    ImplMap impls;
+    
+#ifdef HAVE_IPP
+    impls["ipp"] = SizeConstraintEvenPowerOfTwo;
+#endif
+#ifdef HAVE_FFTW3
+    impls["fftw"] = SizeConstraintNone;
+#endif
+#ifdef HAVE_KISSFFT
+    impls["kissfft"] = SizeConstraintEven;
+#endif
+#ifdef HAVE_VDSP
+    impls["vdsp"] = SizeConstraintEvenPowerOfTwo;
+#endif
+#ifdef HAVE_MEDIALIB
+    impls["medialib"] = SizeConstraintEvenPowerOfTwo;
+#endif
+#ifdef HAVE_OPENMAX
+    impls["openmax"] = SizeConstraintEvenPowerOfTwo;
+#endif
+#ifdef HAVE_SFFT
+    impls["sfft"] = SizeConstraintEvenPowerOfTwo;
+#endif
+#ifdef USE_BUILTIN_FFT
+    impls["cross"] = SizeConstraintPowerOfTwo;
+#endif
+
+    impls["dft"] = SizeConstraintNone;
+
+    return impls;
+}
+
+static std::string
+pickImplementation(int size)
+{
+    ImplMap impls = getImplementationDetails();
+
+    bool isPowerOfTwo = !(size & (size-1));
+    bool isEven = !(size & 1);
+
+    if (defaultImplementation != "") {
+        ImplMap::const_iterator itr = impls.find(defaultImplementation);
+        if (itr != impls.end()) {
+            if (((itr->second & SizeConstraintPowerOfTwo) && !isPowerOfTwo) ||
+                ((itr->second & SizeConstraintEven) && !isEven)) {
+//                std::cerr << "NOTE: bqfft: Explicitly-set default "
+//                          << "implementation \"" << defaultImplementation
+//                          << "\" does not support size " << size
+//                          << ", trying other compiled-in implementations"
+//                          << std::endl;
+            } else {
+                return defaultImplementation;
+            }
+        } else {
+            std::cerr << "WARNING: bqfft: Default implementation \""
+                      << defaultImplementation << "\" is not compiled in"
+                      << std::endl;
+        }
+    } 
+    
+    std::string preference[] = {
+        "ipp", "vdsp", "fftw", "sfft", "openmax",
+        "medialib", "kissfft", "cross"
+    };
+
+    for (int i = 0; i < int(sizeof(preference)/sizeof(preference[0])); ++i) {
+        ImplMap::const_iterator itr = impls.find(preference[i]);
+        if (itr != impls.end()) {
+            if ((itr->second & SizeConstraintPowerOfTwo) && !isPowerOfTwo) {
+                continue;
+            }
+            if ((itr->second & SizeConstraintEven) && !isEven) {
+                continue;
+            }
+            return preference[i];
+        }
+    }
+
+    std::cerr << "WARNING: bqfft: No compiled-in implementation supports size "
+              << size << ", falling back to slow DFT" << std::endl;
+    
+    return "dft";
+}
 
 std::set<std::string>
 FFT::getImplementations()
 {
-    std::set<std::string> impls;
-#ifdef HAVE_IPP
-    impls.insert("ipp");
-#endif
-#ifdef HAVE_FFTW3
-    impls.insert("fftw");
-#endif
-#ifdef HAVE_KISSFFT
-    impls.insert("kissfft");
-#endif
-#ifdef HAVE_VDSP
-    impls.insert("vdsp");
-#endif
-#ifdef HAVE_MEDIALIB
-    impls.insert("medialib");
-#endif
-#ifdef HAVE_OPENMAX
-    impls.insert("openmax");
-#endif
-#ifdef HAVE_SFFT
-    impls.insert("sfft");
-#endif
-#ifdef USE_BUILTIN_FFT
-    impls.insert("cross");
-#endif
-    return impls;
-}
-
-void
-FFT::pickDefaultImplementation()
-{
-    if (m_implementation != "") return;
-
-    std::set<std::string> impls = getImplementations();
-
-    std::string best = "cross";
-    if (impls.find("kissfft") != impls.end()) best = "kissfft";
-    if (impls.find("medialib") != impls.end()) best = "medialib";
-    if (impls.find("openmax") != impls.end()) best = "openmax";
-    if (impls.find("sfft") != impls.end()) best = "sfft";
-    if (impls.find("fftw") != impls.end()) best = "fftw";
-    if (impls.find("vdsp") != impls.end()) best = "vdsp";
-    if (impls.find("ipp") != impls.end()) best = "ipp";
-    
-    m_implementation = best;
+    ImplMap impls = getImplementationDetails();
+    std::set<std::string> toReturn;
+    for (ImplMap::const_iterator i = impls.begin(); i != impls.end(); ++i) {
+        toReturn.insert(i->first);
+    }
+    return toReturn;
 }
 
 std::string
 FFT::getDefaultImplementation()
 {
-    return m_implementation;
+    return defaultImplementation;
 }
 
 void
 FFT::setDefaultImplementation(std::string i)
 {
-    m_implementation = i;
+    if (i == "") {
+        defaultImplementation = i;
+        return;
+    } 
+    ImplMap impls = getImplementationDetails();
+    ImplMap::const_iterator itr = impls.find(i);
+    if (itr == impls.end()) {
+        std::cerr << "WARNING: bqfft: setDefaultImplementation: "
+                  << "requested implementation \"" << i
+                  << "\" is not compiled in" << std::endl;
+    } else {
+        defaultImplementation = i;
+    }
 }
 
 FFT::FFT(int size, int debugLevel) :
     d(0)
 {
-    if ((size < 2) ||
-        (size & (size-1))) {
-        std::cerr << "FFT::FFT(" << size << "): power-of-two sizes only supported, minimum size 2" << std::endl;
-#ifndef NO_EXCEPTIONS
-        throw InvalidSize;
-#else
-        abort();
-#endif
-    }
-
-    if (m_implementation == "") pickDefaultImplementation();
-    std::string impl = m_implementation;
+    std::string impl = pickImplementation(size);
 
     if (debugLevel > 0) {
         std::cerr << "FFT::FFT(" << size << "): using implementation: "

          
@@ 3119,6 3423,8 @@ FFT::FFT(int size, int debugLevel) :
 #ifdef USE_BUILTIN_FFT
         d = new FFTs::D_Cross(size);
 #endif
+    } else if (impl == "dft") {
+        d = new FFTs::D_DFT(size);
     }
 
     if (!d) {

          
@@ 3314,10 3620,18 @@ FFT::getSupportedPrecisions() const
 
 #ifdef FFT_MEASUREMENT
 
+#ifdef FFT_MEASUREMENT_RETURN_RESULT_TEXT
 std::string
+#else
+void
+#endif
 FFT::tune()
 {
+#ifdef FFT_MEASUREMENT_RETURN_RESULT_TEXT
     std::ostringstream os;
+#else
+#define os std::cerr
+#endif
     os << "FFT::tune()..." << std::endl;
 
     std::vector<int> sizes;

          
@@ 3404,6 3718,12 @@ FFT::tune()
         candidates[d] = 6;
 #endif
 
+        os << "Constructing new DFT object for size " << size << "..." << std::endl;
+        d = new FFTs::D_DFT(size);
+        d->initFloat();
+        d->initDouble();
+        candidates[d] = 7;
+
         os << "CLOCKS_PER_SEC = " << CLOCKS_PER_SEC << std::endl;
         float divisor = float(CLOCKS_PER_SEC) / 1000.f;
         

          
@@ 3586,7 3906,9 @@ FFT::tune()
 
     os << "overall winner is " << best << " with " << bestscore << " wins" << std::endl;
 
+#ifdef FFT_MEASUREMENT_RETURN_RESULT_TEXT
     return os.str();
+#endif
 }
 
 #endif

          
M test/TestFFT.cpp +462 -219
@@ 49,13 49,19 @@ using namespace breakfastquay;
 
 BOOST_AUTO_TEST_SUITE(TestFFT)
 
-#ifdef TEST_FFT_IS_SINGLE_PRECISION_ONLY
-static const double eps = 1e-7;
-static const float epsf = 1e-6f;
-#else
-static const double eps = 1e-14;
-static const float epsf = 1e-7f;
-#endif
+#define DEFINE_EPS(fft) \
+    float epsf = 1e-6f; \
+    double eps; \
+    if (fft.getSupportedPrecisions() & FFT::DoublePrecision) { \
+        eps = 1e-14; \
+    } else { \
+        eps = epsf; \
+    } \
+    (void)epsf; (void)eps;
+
+#define USING_FFT(n) \
+    FFT fft(n); \
+    DEFINE_EPS(fft);
 
 #define COMPARE(a, b) BOOST_CHECK_SMALL(a-b, eps) 
 #define COMPARE_F(a, b) BOOST_CHECK_SMALL(a-b, epsf) 

          
@@ 81,65 87,117 @@ static const float epsf = 1e-7f;
         BOOST_CHECK_SMALL(a[cmp_i]/s - b[cmp_i], epsf); \
     }
 
-BOOST_AUTO_TEST_CASE(dc)
+#define ONE_IMPL_AUTO_TEST_CASE(name, impl) \
+    BOOST_AUTO_TEST_CASE(name##_##impl) \
+    { \
+        std::set<std::string> impls = FFT::getImplementations(); \
+        if (impls.find(#impl) == impls.end()) return; \
+        FFT::setDefaultImplementation(#impl); \
+        performTest_##name(); \
+        FFT::setDefaultImplementation(""); \
+    }
+
+// If you add an implementation in FFT.cpp, add it also to
+// ALL_IMPL_AUTO_TEST_CASE and all_implementations[] below
+
+#define ALL_IMPL_AUTO_TEST_CASE(name) \
+    void performTest_##name (); \
+    ONE_IMPL_AUTO_TEST_CASE(name, ipp); \
+    ONE_IMPL_AUTO_TEST_CASE(name, fftw); \
+    ONE_IMPL_AUTO_TEST_CASE(name, kissfft); \
+    ONE_IMPL_AUTO_TEST_CASE(name, vdsp); \
+    ONE_IMPL_AUTO_TEST_CASE(name, medialib); \
+    ONE_IMPL_AUTO_TEST_CASE(name, openmax); \
+    ONE_IMPL_AUTO_TEST_CASE(name, sfft); \
+    ONE_IMPL_AUTO_TEST_CASE(name, cross); \
+    ONE_IMPL_AUTO_TEST_CASE(name, dft); \
+    void performTest_##name ()
+
+std::string all_implementations[] = {
+    "ipp", "vdsp", "fftw", "sfft", "openmax",
+    "medialib", "kissfft", "cross", "dft"
+};
+
+BOOST_AUTO_TEST_CASE(showImplementations)
+{
+    std::set<std::string> impls = FFT::getImplementations();
+    std::cerr << "\nThe following implementations are compiled in and will be tested:" << std::endl;
+    for (int i = 0; i < int(sizeof(all_implementations)/sizeof(all_implementations[0])); ++i) {
+        if (impls.find(all_implementations[i]) != impls.end()) {
+            std::cerr << " +" << all_implementations[i];
+        }
+    }
+    std::cerr << std::endl << std::endl;
+    std::cerr << "The following implementations are NOT compiled in and will not be tested:" << std::endl;
+    for (int i = 0; i < int(sizeof(all_implementations)/sizeof(all_implementations[0])); ++i) {
+        if (impls.find(all_implementations[i]) == impls.end()) {
+            std::cerr << " -" << all_implementations[i];
+        }
+    }
+    std::cerr << std::endl << std::endl;
+}
+
+
+/*
+ * 1a. Simple synthetic signals, transforms to separate real/imag arrays,
+ *     double-precision
+ */
+
+ALL_IMPL_AUTO_TEST_CASE(dc)
 {
     // DC-only signal. The DC bin is purely real
     double in[] = { 1, 1, 1, 1 };
     double re[3], im[3];
-    FFT(4).forward(in, re, im);
+    USING_FFT(4);
+    fft.forward(in, re, im);
     COMPARE(re[0], 4.0);
     COMPARE_ZERO(re[1]);
     COMPARE_ZERO(re[2]);
     COMPARE_ALL(im, 0.0);
     double back[4];
-    FFT(4).inverse(re, im, back);
+    fft.inverse(re, im, back);
     COMPARE_SCALED(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(showImplementation)
-{
-    FFT f(4); // to ensure an implementation has been picked
-    (void)f;  // avoid compiler warning
-    std::cerr << "Using implementation: " << FFT::getDefaultImplementation()
-              << std::endl;
-}
-
-BOOST_AUTO_TEST_CASE(sine)
+ALL_IMPL_AUTO_TEST_CASE(sine)
 {
     // Sine. Output is purely imaginary
     double in[] = { 0, 1, 0, -1 };
     double re[3], im[3];
-    FFT(4).forward(in, re, im);
+    USING_FFT(4);
+    fft.forward(in, re, im);
     COMPARE_ALL(re, 0.0);
     COMPARE_ZERO(im[0]);
     COMPARE(im[1], -2.0);
     COMPARE_ZERO(im[2]);
     double back[4];
-    FFT(4).inverse(re, im, back);
+    fft.inverse(re, im, back);
     COMPARE_SCALED(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(cosine)
+ALL_IMPL_AUTO_TEST_CASE(cosine)
 {
     // Cosine. Output is purely real
     double in[] = { 1, 0, -1, 0 };
     double re[3], im[3];
-    FFT(4).forward(in, re, im);
+    USING_FFT(4);
+    fft.forward(in, re, im);
     COMPARE_ZERO(re[0]);
     COMPARE(re[1], 2.0);
     COMPARE_ZERO(re[2]);
     COMPARE_ALL(im, 0.0);
     double back[4];
-    FFT(4).inverse(re, im, back);
+    fft.inverse(re, im, back);
     COMPARE_SCALED(back, in, 4);
 }
 	
-BOOST_AUTO_TEST_CASE(sineCosine)
+ALL_IMPL_AUTO_TEST_CASE(sineCosine)
 {
     // Sine and cosine mixed
     double in[] = { 0.5, 1, -0.5, -1 };
     double re[3], im[3];
-    FFT(4).forward(in, re, im);
+    USING_FFT(4);
+    fft.forward(in, re, im);
     COMPARE_ZERO(re[0]);
     COMPARE(re[1], 1.0);
     COMPARE_ZERO(re[2]);

          
@@ 147,30 205,156 @@ BOOST_AUTO_TEST_CASE(sineCosine)
     COMPARE(im[1], -2.0);
     COMPARE_ZERO(im[2]);
     double back[4];
-    FFT(4).inverse(re, im, back);
+    fft.inverse(re, im, back);
     COMPARE_SCALED(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(nyquist)
+ALL_IMPL_AUTO_TEST_CASE(nyquist)
 {
     double in[] = { 1, -1, 1, -1 };
     double re[3], im[3];
-    FFT(4).forward(in, re, im);
+    USING_FFT(4);
+    fft.forward(in, re, im);
     COMPARE_ZERO(re[0]);
     COMPARE_ZERO(re[1]);
     COMPARE(re[2], 4.0);
     COMPARE_ALL(im, 0.0);
     double back[4];
-    FFT(4).inverse(re, im, back);
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 1.0);
+    COMPARE(re[1], 1.0);
+    COMPARE(re[2], 1.0);
+    COMPARE_ALL(im, 0.0);
+    double back[4];
+    fft.inverse(re, im, back);
     COMPARE_SCALED(back, in, 4);
 }
+
+
+/*
+ * 1b. Simple synthetic signals, transforms to separate real/imag arrays,
+ *     single-precision (i.e. single-precision version of 1a)
+ */
+
+ALL_IMPL_AUTO_TEST_CASE(dcF)
+{
+    // DC-only signal. The DC bin is purely real
+    float in[] = { 1, 1, 1, 1 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_F(re[0], 4.0f);
+    COMPARE_ZERO_F(re[1]);
+    COMPARE_ZERO_F(re[2]);
+    COMPARE_ALL_F(im, 0.0f);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(sineF)
+{
+    // Sine. Output is purely imaginary
+    float in[] = { 0, 1, 0, -1 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_ALL_F(re, 0.0f);
+    COMPARE_ZERO_F(im[0]);
+    COMPARE_F(im[1], -2.0f);
+    COMPARE_ZERO_F(im[2]);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(cosineF)
+{
+    // Cosine. Output is purely real
+    float in[] = { 1, 0, -1, 0 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_ZERO_F(re[0]);
+    COMPARE_F(re[1], 2.0f);
+    COMPARE_ZERO_F(re[2]);
+    COMPARE_ALL_F(im, 0.0f);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
 	
-BOOST_AUTO_TEST_CASE(interleaved)
+ALL_IMPL_AUTO_TEST_CASE(sineCosineF)
+{
+    // Sine and cosine mixed
+    float in[] = { 0.5, 1, -0.5, -1 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_ZERO_F(re[0]);
+    COMPARE_F(re[1], 1.0f);
+    COMPARE_ZERO_F(re[2]);
+    COMPARE_ZERO_F(im[0]);
+    COMPARE_F(im[1], -2.0f);
+    COMPARE_ZERO_F(im[2]);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(nyquistF)
+{
+    float in[] = { 1, -1, 1, -1 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_ZERO_F(re[0]);
+    COMPARE_ZERO_F(re[1]);
+    COMPARE_F(re[2], 4.0f);
+    COMPARE_ALL_F(im, 0.0f);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(diracF)
+{
+    float in[] = { 1, 0, 0, 0 };
+    float re[3], im[3];
+    USING_FFT(4);
+    fft.forward(in, re, im);
+    COMPARE_F(re[0], 1.0f);
+    COMPARE_F(re[1], 1.0f);
+    COMPARE_F(re[2], 1.0f);
+    COMPARE_ALL_F(im, 0.0f);
+    float back[4];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED_F(back, in, 4);
+}
+
+
+/*
+ * 2a. Subset of synthetic signals, testing different output formats
+ *     (interleaved complex, polar, magnitude-only, and our weird
+ *     cepstral thing), double-precision
+ */ 
+
+ALL_IMPL_AUTO_TEST_CASE(interleaved)
 {
     // Sine and cosine mixed, test output format
     double in[] = { 0.5, 1, -0.5, -1 };
     double out[6];
-    FFT(4).forwardInterleaved(in, out);
+    USING_FFT(4);
+    fft.forwardInterleaved(in, out);
     COMPARE_ZERO(out[0]);
     COMPARE_ZERO(out[1]);
     COMPARE(out[2], 1.0);

          
@@ 178,31 362,16 @@ BOOST_AUTO_TEST_CASE(interleaved)
     COMPARE_ZERO(out[4]);
     COMPARE_ZERO(out[5]);
     double back[4];
-    FFT(4).inverseInterleaved(out, back);
+    fft.inverseInterleaved(out, back);
     COMPARE_SCALED(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(cosinePolar)
-{
-    double in[] = { 1, 0, -1, 0 };
-    double mag[3], phase[3];
-    FFT(4).forwardPolar(in, mag, phase);
-    COMPARE_ZERO(mag[0]);
-    COMPARE(mag[1], 2.0);
-    COMPARE_ZERO(mag[2]);
-    // No meaningful tests for phase[i] where mag[i]==0 (phase
-    // could legitimately be anything)
-    COMPARE_ZERO(phase[1]);
-    double back[4];
-    FFT(4).inversePolar(mag, phase, back);
-    COMPARE_SCALED(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(sinePolar)
+ALL_IMPL_AUTO_TEST_CASE(sinePolar)
 {
     double in[] = { 0, 1, 0, -1 };
     double mag[3], phase[3];
-    FFT(4).forwardPolar(in, mag, phase);
+    USING_FFT(4);
+    fft.forwardPolar(in, mag, phase);
     COMPARE_ZERO(mag[0]);
     COMPARE(mag[1], 2.0);
     COMPARE_ZERO(mag[2]);

          
@@ 210,162 379,72 @@ BOOST_AUTO_TEST_CASE(sinePolar)
     // could legitimately be anything)
     COMPARE(phase[1], -M_PI/2.0);
     double back[4];
-    FFT(4).inversePolar(mag, phase, back);
+    fft.inversePolar(mag, phase, back);
     COMPARE_SCALED(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(magnitude)
+ALL_IMPL_AUTO_TEST_CASE(cosinePolar)
+{
+    double in[] = { 1, 0, -1, 0 };
+    double mag[3], phase[3];
+    USING_FFT(4);
+    fft.forwardPolar(in, mag, phase);
+    COMPARE_ZERO(mag[0]);
+    COMPARE(mag[1], 2.0);
+    COMPARE_ZERO(mag[2]);
+    // No meaningful tests for phase[i] where mag[i]==0 (phase
+    // could legitimately be anything)
+    COMPARE_ZERO(phase[1]);
+    double back[4];
+    fft.inversePolar(mag, phase, back);
+    COMPARE_SCALED(back, in, 4);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(magnitude)
 {
     // Sine and cosine mixed
     double in[] = { 0.5, 1, -0.5, -1 };
     double out[3];
-    FFT(4).forwardMagnitude(in, out);
+    USING_FFT(4);
+    fft.forwardMagnitude(in, out);
     COMPARE_ZERO(out[0]);
     COMPARE_F(float(out[1]), sqrtf(5.0));
     COMPARE_ZERO(out[2]);
 }
 
-BOOST_AUTO_TEST_CASE(dirac)
-{
-    double in[] = { 1, 0, 0, 0 };
-    double re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE(re[0], 1.0);
-    COMPARE(re[1], 1.0);
-    COMPARE(re[2], 1.0);
-    COMPARE_ALL(im, 0.0);
-    double back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(cepstrum)
+ALL_IMPL_AUTO_TEST_CASE(cepstrum)
 {
     double in[] = { 1, 0, 0, 0, 1, 0, 0, 0 };
     double mag[5];
-    FFT(8).forwardMagnitude(in, mag);
+    USING_FFT(8);
+    fft.forwardMagnitude(in, mag);
     double cep[8];
-    FFT(8).inverseCepstral(mag, cep);
-    COMPARE_ZERO(cep[1]);
-    COMPARE_ZERO(cep[2]);
-    COMPARE_ZERO(cep[3]);
-    COMPARE_ZERO(cep[5]);
-    COMPARE_ZERO(cep[6]);
-    COMPARE_ZERO(cep[7]);
+    fft.inverseCepstral(mag, cep);
+    BOOST_CHECK_SMALL(cep[1], 1e-9);
+    BOOST_CHECK_SMALL(cep[2], 1e-9);
+    BOOST_CHECK_SMALL(cep[3], 1e-9);
+    BOOST_CHECK_SMALL(cep[5], 1e-9);
+    BOOST_CHECK_SMALL(cep[6], 1e-9);
+    BOOST_CHECK_SMALL(cep[7], 1e-9);
     BOOST_CHECK_SMALL(-6.561181 - cep[0]/8, 0.000001);
     BOOST_CHECK_SMALL( 7.254329 - cep[4]/8, 0.000001);
 }
 
-BOOST_AUTO_TEST_CASE(forwardArrayBounds)
-{
-    // initialise bins to something recognisable, so we can tell
-    // if they haven't been written
-    double in[] = { 1, 1, -1, -1 };
-    double re[] = { 999, 999, 999, 999, 999 };
-    double im[] = { 999, 999, 999, 999, 999 };
-    FFT(4).forward(in, re+1, im+1);
-    // And check we haven't overrun the arrays
-    COMPARE(re[0], 999.0);
-    COMPARE(im[0], 999.0);
-    COMPARE(re[4], 999.0);
-    COMPARE(im[4], 999.0);
-}
 
-BOOST_AUTO_TEST_CASE(inverseArrayBounds)
-{
-    // initialise bins to something recognisable, so we can tell
-    // if they haven't been written
-    double re[] = { 0, 1, 0 };
-    double im[] = { 0, -2, 0 };
-    double out[] = { 999, 999, 999, 999, 999, 999 };
-    FFT(4).inverse(re, im, out+1);
-    // And check we haven't overrun the arrays
-    COMPARE(out[0], 999.0);
-    COMPARE(out[5], 999.0);
-}
-
-BOOST_AUTO_TEST_CASE(dcF)
-{
-    // DC-only signal. The DC bin is purely real
-    float in[] = { 1, 1, 1, 1 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_F(re[0], 4.0f);
-    COMPARE_ZERO_F(re[1]);
-    COMPARE_ZERO_F(re[2]);
-    COMPARE_ALL_F(im, 0.0f);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(sineF)
-{
-    // Sine. Output is purely imaginary
-    float in[] = { 0, 1, 0, -1 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_ALL_F(re, 0.0f);
-    COMPARE_ZERO_F(im[0]);
-    COMPARE_F(im[1], -2.0f);
-    COMPARE_ZERO_F(im[2]);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(cosineF)
-{
-    // Cosine. Output is purely real
-    float in[] = { 1, 0, -1, 0 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_ZERO_F(re[0]);
-    COMPARE_F(re[1], 2.0f);
-    COMPARE_ZERO_F(re[2]);
-    COMPARE_ALL_F(im, 0.0f);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
+/*
+ * 2b. Subset of synthetic signals, testing different output formats
+ *     (interleaved complex, polar, magnitude-only, and our weird
+ *     cepstral thing), single-precision (i.e. single-precision
+ *     version of 2a)
+ */ 
 	
-BOOST_AUTO_TEST_CASE(sineCosineF)
-{
-    // Sine and cosine mixed
-    float in[] = { 0.5, 1, -0.5, -1 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_ZERO_F(re[0]);
-    COMPARE_F(re[1], 1.0f);
-    COMPARE_ZERO_F(re[2]);
-    COMPARE_ZERO_F(im[0]);
-    COMPARE_F(im[1], -2.0f);
-    COMPARE_ZERO_F(im[2]);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(nyquistF)
-{
-    float in[] = { 1, -1, 1, -1 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_ZERO_F(re[0]);
-    COMPARE_ZERO_F(re[1]);
-    COMPARE_F(re[2], 4.0f);
-    COMPARE_ALL_F(im, 0.0f);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
-	
-BOOST_AUTO_TEST_CASE(interleavedF)
+ALL_IMPL_AUTO_TEST_CASE(interleavedF)
 {
     // Sine and cosine mixed, test output format
     float in[] = { 0.5, 1, -0.5, -1 };
     float out[6];
-    FFT(4).forwardInterleaved(in, out);
+    USING_FFT(4);
+    fft.forwardInterleaved(in, out);
     COMPARE_ZERO_F(out[0]);
     COMPARE_ZERO_F(out[1]);
     COMPARE_F(out[2], 1.0f);

          
@@ 373,15 452,16 @@ BOOST_AUTO_TEST_CASE(interleavedF)
     COMPARE_ZERO_F(out[4]);
     COMPARE_ZERO_F(out[5]);
     float back[4];
-    FFT(4).inverseInterleaved(out, back);
+    fft.inverseInterleaved(out, back);
     COMPARE_SCALED_F(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(cosinePolarF)
+ALL_IMPL_AUTO_TEST_CASE(cosinePolarF)
 {
     float in[] = { 1, 0, -1, 0 };
     float mag[3], phase[3];
-    FFT(4).forwardPolar(in, mag, phase);
+    USING_FFT(4);
+    fft.forwardPolar(in, mag, phase);
     COMPARE_ZERO_F(mag[0]);
     COMPARE_F(mag[1], 2.0f);
     COMPARE_ZERO_F(mag[2]);

          
@@ 389,15 469,16 @@ BOOST_AUTO_TEST_CASE(cosinePolarF)
     // could legitimately be anything)
     COMPARE_ZERO_F(phase[1]);
     float back[4];
-    FFT(4).inversePolar(mag, phase, back);
+    fft.inversePolar(mag, phase, back);
     COMPARE_SCALED_F(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(sinePolarF)
+ALL_IMPL_AUTO_TEST_CASE(sinePolarF)
 {
     float in[] = { 0, 1, 0, -1 };
     float mag[3], phase[3];
-    FFT(4).forwardPolar(in, mag, phase);
+    USING_FFT(4);
+    fft.forwardPolar(in, mag, phase);
     COMPARE_ZERO_F(mag[0]);
     COMPARE_F(mag[1], 2.0f);
     COMPARE_ZERO_F(mag[2]);

          
@@ 405,42 486,30 @@ BOOST_AUTO_TEST_CASE(sinePolarF)
     // could legitimately be anything)
     COMPARE_F(phase[1], -float(M_PI)/2.0f);
     float back[4];
-    FFT(4).inversePolar(mag, phase, back);
+    fft.inversePolar(mag, phase, back);
     COMPARE_SCALED_F(back, in, 4);
 }
 
-BOOST_AUTO_TEST_CASE(magnitudeF)
+ALL_IMPL_AUTO_TEST_CASE(magnitudeF)
 {
     // Sine and cosine mixed
     float in[] = { 0.5, 1, -0.5, -1 };
     float out[3];
-    FFT(4).forwardMagnitude(in, out);
+    USING_FFT(4);
+    fft.forwardMagnitude(in, out);
     COMPARE_ZERO_F(out[0]);
     COMPARE_F(float(out[1]), sqrtf(5.0f));
     COMPARE_ZERO_F(out[2]);
 }
 
-BOOST_AUTO_TEST_CASE(diracF)
-{
-    float in[] = { 1, 0, 0, 0 };
-    float re[3], im[3];
-    FFT(4).forward(in, re, im);
-    COMPARE_F(re[0], 1.0f);
-    COMPARE_F(re[1], 1.0f);
-    COMPARE_F(re[2], 1.0f);
-    COMPARE_ALL_F(im, 0.0f);
-    float back[4];
-    FFT(4).inverse(re, im, back);
-    COMPARE_SCALED_F(back, in, 4);
-}
-
-BOOST_AUTO_TEST_CASE(cepstrumF)
+ALL_IMPL_AUTO_TEST_CASE(cepstrumF)
 {
     float in[] = { 1, 0, 0, 0, 1, 0, 0, 0 };
     float mag[5];
-    FFT(8).forwardMagnitude(in, mag);
+    USING_FFT(8);
+    fft.forwardMagnitude(in, mag);
     float cep[8];
-    FFT(8).inverseCepstral(mag, cep);
+    fft.inverseCepstral(mag, cep);
     COMPARE_ZERO_F(cep[1]);
     COMPARE_ZERO_F(cep[2]);
     COMPARE_ZERO_F(cep[3]);

          
@@ 451,32 520,206 @@ BOOST_AUTO_TEST_CASE(cepstrumF)
     BOOST_CHECK_SMALL( 7.254329 - cep[4]/8, 0.000001);
 }
 
-BOOST_AUTO_TEST_CASE(forwardArrayBoundsF)
+
+/*
+ * 4. Bounds checking, double-precision and single-precision
+ */
+
+ALL_IMPL_AUTO_TEST_CASE(forwardArrayBounds)
+{
+    double in[] = { 1, 1, -1, -1 };
+
+    // Initialise output bins to something recognisable, so we can
+    // tell if they haven't been written
+    double re[] = { 999, 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999, 999 };
+    
+    USING_FFT(4);
+    fft.forward(in, re+1, im+1);
+    
+    // Check we haven't overrun the output arrays
+    COMPARE(re[0], 999.0);
+    COMPARE(im[0], 999.0);
+    COMPARE(re[4], 999.0);
+    COMPARE(im[4], 999.0);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(inverseArrayBounds)
 {
-    // initialise bins to something recognisable, so we can tell
-    // if they haven't been written
+    // The inverse transform is only supposed to refer to the first
+    // N/2+1 bins and synthesise the rest rather than read them - so
+    // initialise the next one to some value that would mess up the
+    // results if it were used
+    double re[] = { 0, 1, 0, 456 };
+    double im[] = { 0, -2, 0, 456 };
+
+    // Initialise output bins to something recognisable, so we can
+    // tell if they haven't been written
+    double out[] = { 999, 999, 999, 999, 999, 999 };
+    
+    USING_FFT(4);
+    fft.inverse(re, im, out+1);
+
+    // Check we haven't overrun the output arrays
+    COMPARE(out[0], 999.0);
+    COMPARE(out[5], 999.0);
+
+    // And check the results are as we expect, i.e. that we haven't
+    // used the bogus final bin
+    COMPARE(out[1] / 4, 0.5);
+    COMPARE(out[2] / 4, 1.0);
+    COMPARE(out[3] / 4, -0.5);
+    COMPARE(out[4] / 4, -1.0);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(forwardArrayBoundsF)
+{
     float in[] = { 1, 1, -1, -1 };
+    
+    // Initialise output bins to something recognisable, so we can
+    // tell if they haven't been written
     float re[] = { 999, 999, 999, 999, 999 };
     float im[] = { 999, 999, 999, 999, 999 };
-    FFT(4).forward(in, re+1, im+1);
-    // And check we haven't overrun the arrays
+    
+    USING_FFT(4);
+    fft.forward(in, re+1, im+1);
+    
+    // Check we haven't overrun the output arrays
     COMPARE_F(re[0], 999.0f);
     COMPARE_F(im[0], 999.0f);
     COMPARE_F(re[4], 999.0f);
     COMPARE_F(im[4], 999.0f);
 }
 
-BOOST_AUTO_TEST_CASE(inverseArrayBoundsF)
+ALL_IMPL_AUTO_TEST_CASE(inverseArrayBoundsF)
 {
-    // initialise bins to something recognisable, so we can tell
-    // if they haven't been written
-    float re[] = { 0, 1, 0 };
-    float im[] = { 0, -2, 0 };
+    // The inverse transform is only supposed to refer to the first
+    // N/2+1 bins and synthesise the rest rather than read them - so
+    // initialise the next one to some value that would mess up the
+    // results if it were used
+    float re[] = { 0, 1, 0, 456 };
+    float im[] = { 0, -2, 0, 456 };
+
+    // Initialise output bins to something recognisable, so we can
+    // tell if they haven't been written
     float out[] = { 999, 999, 999, 999, 999, 999 };
-    FFT(4).inverse(re, im, out+1);
-    // And check we haven't overrun the arrays
+    
+    USING_FFT(4);
+    fft.inverse(re, im, out+1);
+    
+    // Check we haven't overrun the output arrays
     COMPARE_F(out[0], 999.0f);
     COMPARE_F(out[5], 999.0f);
+
+    // And check the results are as we expect, i.e. that we haven't
+    // used the bogus final bin
+    COMPARE_F(out[1] / 4.0f, 0.5f);
+    COMPARE_F(out[2] / 4.0f, 1.0f);
+    COMPARE_F(out[3] / 4.0f, -0.5f);
+    COMPARE_F(out[4] / 4.0f, -1.0f);
+}
+
+
+/*
+ * 5. Less common transform lengths - we should always fall back on
+ *    some implementation that can handle these, even if the requested
+ *    one doesn't. Note that the dirac tests we do first are "easier"
+ *    in that they don't vary with length
+ */
+
+ALL_IMPL_AUTO_TEST_CASE(dirac_1)
+{
+    double in[] = { 1 };
+    double re[1], im[1];
+    USING_FFT(1);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 1.0);
+    COMPARE_ALL(im, 0.0);
+    double back[1];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 1);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(dirac_6)
+{
+    double in[] = { 1, 0, 0, 0, 0, 0 };
+    double re[4], im[4];
+    USING_FFT(6);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 1.0);
+    COMPARE(re[1], 1.0);
+    COMPARE(re[2], 1.0);
+    COMPARE(re[3], 1.0);
+    COMPARE_ALL(im, 0.0);
+    double back[6];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 6);
 }
 
+ALL_IMPL_AUTO_TEST_CASE(dirac_7)
+{
+    double in[] = { 1, 0, 0, 0, 0, 0, 0 };
+    double re[4], im[4];
+    USING_FFT(7);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 1.0);
+    COMPARE(re[1], 1.0);
+    COMPARE(re[2], 1.0);
+    COMPARE(re[3], 1.0);
+    COMPARE_ALL(im, 0.0);
+    double back[7];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 7);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(sineCosine_6)
+{
+    // Sine and cosine mixed, i.e. f(x) = 0.5 * cos(2*x*pi/6) + sin(2*x*pi/6)
+    double r = sqrt(3.0)/2.0;
+    double in[] = { 0.5, r + 0.25, r - 0.25, -0.5, -r - 0.25, -r + 0.25 };
+    double re[4], im[4];
+    USING_FFT(6);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 0.0);
+    COMPARE(re[1], 1.5);
+    COMPARE(re[2], 0.0);
+    COMPARE(re[3], 0.0);
+    COMPARE(im[0], 0.0);
+    COMPARE(im[1], -3.0);
+    COMPARE(im[2], 0.0);
+    COMPARE(im[3], 0.0);
+    double back[6];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 6);
+}
+
+ALL_IMPL_AUTO_TEST_CASE(sineCosine_7)
+{
+    // Sine and cosine mixed, i.e. f(x) = 0.5 * cos(2*x*pi/6) + sin(2*x*pi/6)
+    double in[] = {
+        0.5,
+        1.0935763833973966,
+        0.8636674452036665,
+        -0.016600694833651286,
+        -0.8843681730687676,
+        -1.086188379159981,
+        -0.47008658153866323
+    };
+    double re[4], im[4];
+    USING_FFT(7);
+    fft.forward(in, re, im);
+    COMPARE(re[0], 0.0);
+    COMPARE(re[1], 1.75);
+    COMPARE(re[2], 0.0);
+    COMPARE(re[3], 0.0);
+    COMPARE(im[0], 0.0);
+    COMPARE(im[1], -3.5);
+    COMPARE(im[2], 0.0);
+    COMPARE(im[3], 0.0);
+    double back[7];
+    fft.inverse(re, im, back);
+    COMPARE_SCALED(back, in, 7);
+}
+
+
 BOOST_AUTO_TEST_SUITE_END()

          
A => test/timings.cpp +11 -0
@@ 0,0 1,11 @@ 
+
+#define FFT_MEASUREMENT 1
+
+#include "../src/FFT.cpp"
+
+#include <iostream>
+
+int main(int, char **)
+{
+    breakfastquay::FFT::tune();
+}