Mercurial > hg > qm-dsp

changeset 114:f6ccde089491 pvoc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Tidy real-to-complex FFT -- forward and inverse have different arguments, so make them separate functions; document
author Chris Cannam
date Wed, 02 Oct 2013 15:04:38 +0100
parents 3cb359d043f0
children f3c69325cca2
files dsp/chromagram/Chromagram.cpp dsp/mfcc/MFCC.cpp dsp/segmentation/ClusterMeltSegmenter.cpp dsp/tempotracking/DownBeat.cpp dsp/transforms/FFT.cpp dsp/transforms/FFT.h tests/TestFFT.cpp
diffstat 7 files changed, 371 insertions(+), 128 deletions(-) [+]
line wrap: on
line diff
--- a/dsp/chromagram/Chromagram.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/chromagram/Chromagram.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -139,8 +139,7 @@
     }
     m_window->cut(m_windowbuf);
 
-    // FFT of current frame
-    m_FFT->process(false, m_windowbuf, m_FFTRe, m_FFTIm);
+    m_FFT->forward(m_windowbuf, m_FFTRe, m_FFTIm);
 
     return process(m_FFTRe, m_FFTIm);
 }
@@ -158,7 +157,6 @@
 
     double cmax = 0.0;
     double cval = 0;
-
     // Calculate ConstantQ frame
     m_ConstantQ->process( real, imag, m_CQRe, m_CQIm );
 	
--- a/dsp/mfcc/MFCC.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/mfcc/MFCC.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -210,7 +210,7 @@
     window->cut(inputData);
   
     /* Calculate the fft on the input frame */
-    fft->process(0, inputData, realOut, imagOut);
+    fft->forward(inputData, realOut, imagOut);
 
     free(inputData);
 
--- a/dsp/segmentation/ClusterMeltSegmenter.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/segmentation/ClusterMeltSegmenter.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -209,7 +209,7 @@
 
         window->cut(frame);
         
-        fft->process(false, frame, real, imag);
+        fft->forward(frame, real, imag);
         
         constq->process(real, imag, cqre, cqim);
 	
--- a/dsp/tempotracking/DownBeat.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/tempotracking/DownBeat.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -193,7 +193,7 @@
 
         // Now FFT beat frame
         
-        m_fft->process(false, m_beatframe, m_fftRealOut, m_fftImagOut);
+        m_fft->forward(m_beatframe, m_fftRealOut, m_fftImagOut);
         
         // Calculate magnitudes
 
--- a/dsp/transforms/FFT.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/transforms/FFT.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -15,9 +15,8 @@
 
 #include <iostream>
 
-FFT::FFT(unsigned int n) :
-    m_n(n),
-    m_private(0)
+FFT::FFT(int n) :
+    m_n(n)
 {
     if( !MathUtilities::isPowerOfTwo(m_n) )
     {
@@ -33,27 +32,51 @@
 
 }
 
-FFTReal::FFTReal(unsigned int n) :
+FFTReal::FFTReal(int n) :
     m_n(n),
-    m_private(0)
+    m_fft(new FFT(n)),
+    m_r(new double[n]),
+    m_i(new double[n]),
+    m_discard(new double[n])
 {
-    m_private = new FFT(m_n);
+    for (int i = 0; i < n; ++i) {
+        m_r[i] = 0;
+        m_i[i] = 0;
+        m_discard[i] = 0;
+    }
 }
 
 FFTReal::~FFTReal()
 {
-    delete (FFT *)m_private;
+    delete m_fft;
+    delete[] m_discard;
+    delete[] m_r;
+    delete[] m_i;
 }
 
 void
-FFTReal::process(bool inverse,
-                 const double *realIn,
+FFTReal::forward(const double *realIn,
                  double *realOut, double *imagOut)
 {
-    ((FFT *)m_private)->process(inverse, realIn, 0, realOut, imagOut);
+    m_fft->process(false, realIn, 0, realOut, imagOut);
 }
 
-static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
+void
+FFTReal::inverse(const double *realIn, const double *imagIn,
+                 double *realOut)
+{
+    for (int i = 0; i < m_n/2 + 1; ++i) {
+        m_r[i] = realIn[i];
+        m_i[i] = imagIn[i];
+        if (i > 0 && i < m_n/2) {
+            m_r[m_n - i] = realIn[i];
+            m_i[m_n - i] = -imagIn[i];
+        }
+    }
+    m_fft->process(true, m_r, m_i, realOut, m_discard);
+}
+
+static int numberOfBitsNeeded(int p_nSamples)
 {	
     int i;
 
@@ -68,9 +91,9 @@
     }
 }
 
-static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
+static int reverseBits(int p_nIndex, int p_nBits)
 {
-    unsigned int i, rev;
+    int i, rev;
 
     for(i=rev=0; i < p_nBits; i++)
     {
@@ -90,9 +113,9 @@
 
 //    std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
 
-    unsigned int NumBits;
-    unsigned int i, j, k, n;
-    unsigned int BlockSize, BlockEnd;
+    int NumBits;
+    int i, j, k, n;
+    int BlockSize, BlockEnd;
 
     double angle_numerator = 2.0 * M_PI;
     double tr, ti;
@@ -110,6 +133,7 @@
     NumBits = numberOfBitsNeeded ( m_n );
 
 
+
     for( i=0; i < m_n; i++ )
     {
 	j = reverseBits ( i, NumBits );
--- a/dsp/transforms/FFT.h	Tue Oct 01 15:38:56 2013 +0100
+++ b/dsp/transforms/FFT.h	Wed Oct 02 15:04:38 2013 +0100
@@ -12,31 +12,82 @@
 class FFT  
 {
 public:
-    FFT(unsigned int nsamples);
+    /**
+     * Construct an FFT object to carry out complex-to-complex
+     * transforms of size nsamples. nsamples must be a power of two in
+     * this implementation.
+     */
+    FFT(int nsamples);
     ~FFT();
 
+    /**
+     * Carry out a forward or inverse transform (depending on the
+     * value of inverse) of size nsamples, where nsamples is the value
+     * provided to the constructor above.
+     *
+     * realIn and (where present) imagIn should contain nsamples each,
+     * and realOut and imagOut should point to enough space to receive
+     * nsamples each.
+     *
+     * imagIn may be NULL if the signal is real, but the other
+     * pointers must be valid.
+     *
+     * The inverse transform is scaled by 1/nsamples.
+     */
     void process(bool inverse,
                  const double *realIn, const double *imagIn,
                  double *realOut, double *imagOut);
     
 private:
-    unsigned int m_n;
-    void *m_private;
+    int m_n;
 };
 
 class FFTReal
 {
 public:
-    FFTReal(unsigned int nsamples);
+    /**
+     * Construct an FFT object to carry out real-to-complex transforms
+     * of size nsamples. nsamples must be a power of two in this
+     * implementation.
+     */
+    FFTReal(int nsamples);
     ~FFTReal();
 
-    void process(bool inverse,
-                 const double *realIn,
+    /**
+     * Carry out a forward real-to-complex transform of size nsamples,
+     * where nsamples is the value provided to the constructor above.
+     *
+     * realIn, realOut, and imagOut must point to (enough space for)
+     * nsamples values. For consistency with the FFT class above, and
+     * compatibility with existing code, the conjugate half of the
+     * output is returned even though it is redundant.
+     */
+    void forward(const double *realIn,
                  double *realOut, double *imagOut);
 
+    /**
+     * Carry out an inverse real transform (i.e. complex-to-real) of
+     * size nsamples, where nsamples is the value provided to the
+     * constructor above.
+     *
+     * realIn and imagIn should point to at least nsamples/2+1 values;
+     * if more are provided, only the first nsamples/2+1 values of
+     * each will be used (the conjugate half will always be deduced
+     * from the first nsamples/2+1 rather than being read from the
+     * input data).  realOut should point to enough space to receive
+     * nsamples values.
+     *
+     * The inverse transform is scaled by 1/nsamples.
+     */
+    void inverse(const double *realIn, const double *imagIn,
+                 double *realOut);
+
 private:
-    unsigned int m_n;
-    void *m_private;
+    int m_n;
+    FFT *m_fft;
+    double *m_r;
+    double *m_i;
+    double *m_discard;
 };    
 
 #endif
--- a/tests/TestFFT.cpp	Tue Oct 01 15:38:56 2013 +0100
+++ b/tests/TestFFT.cpp	Wed Oct 02 15:04:38 2013 +0100
@@ -18,101 +18,7 @@
         BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-14);			\
     }
 
-BOOST_AUTO_TEST_CASE(dc)
-{
-    // DC-only signal. The DC bin is purely real
-    double in[] = { 1, 1, 1, 1 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 4.0);
-    BOOST_CHECK_EQUAL(re[1], 0.0);
-    BOOST_CHECK_EQUAL(re[2], 0.0);
-    COMPARE_CONST(im, 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-
-BOOST_AUTO_TEST_CASE(sine)
-{
-    // Sine. Output is purely imaginary
-    double in[] = { 0, 1, 0, -1 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    COMPARE_CONST(re, 0.0);
-    BOOST_CHECK_EQUAL(im[0], 0.0);
-    BOOST_CHECK_EQUAL(im[1], -2.0);
-    BOOST_CHECK_EQUAL(im[2], 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-
-BOOST_AUTO_TEST_CASE(cosine)
-{
-    // Cosine. Output is purely real
-    double in[] = { 1, 0, -1, 0 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 0.0);
-    BOOST_CHECK_EQUAL(re[1], 2.0);
-    BOOST_CHECK_EQUAL(re[2], 0.0);
-    COMPARE_CONST(im, 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-	
-BOOST_AUTO_TEST_CASE(sineCosine)
-{
-    // Sine and cosine mixed
-    double in[] = { 0.5, 1, -0.5, -1 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 0.0);
-    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
-    BOOST_CHECK_EQUAL(re[2], 0.0);
-    BOOST_CHECK_EQUAL(im[0], 0.0);
-    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
-    BOOST_CHECK_EQUAL(im[2], 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-
-BOOST_AUTO_TEST_CASE(nyquist)
-{
-    double in[] = { 1, -1, 1, -1 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 0.0);
-    BOOST_CHECK_EQUAL(re[1], 0.0);
-    BOOST_CHECK_EQUAL(re[2], 4.0);
-    COMPARE_CONST(im, 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
-
-BOOST_AUTO_TEST_CASE(dirac)
-{
-    double in[] = { 1, 0, 0, 0 };
-    double re[4], im[4];
-    FFT(4).process(false, in, 0, re, im);
-    BOOST_CHECK_EQUAL(re[0], 1.0);
-    BOOST_CHECK_EQUAL(re[1], 1.0);
-    BOOST_CHECK_EQUAL(re[2], 1.0);
-    COMPARE_CONST(im, 0.0);
-    double back[4];
-    double backim[4];
-    FFT(4).process(true, re, im, back, backim);
-    COMPARE_ARRAY(back, in);
-}
+//!!! need at least one test with complex time-domain signal
 
 BOOST_AUTO_TEST_CASE(forwardArrayBounds)
 {
@@ -129,15 +35,30 @@
     BOOST_CHECK_EQUAL(im[5], 999.0);
 }
 
+BOOST_AUTO_TEST_CASE(r_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, 999 };
+    double im[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re+1, im+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(re[0], 999.0);
+    BOOST_CHECK_EQUAL(im[0], 999.0);
+    BOOST_CHECK_EQUAL(re[5], 999.0);
+    BOOST_CHECK_EQUAL(im[5], 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 re[] = { 0, 1, 0, 1 };
+    double im[] = { 0, -2, 0, 2 };
     double outre[] = { 999, 999, 999, 999, 999, 999 };
     double outim[] = { 999, 999, 999, 999, 999, 999 };
-    FFT(4).process(false, re, im, outre+1, outim+1);
+    FFT(4).process(true, re, im, outre+1, outim+1);
     // And check we haven't overrun the arrays
     BOOST_CHECK_EQUAL(outre[0], 999.0);
     BOOST_CHECK_EQUAL(outim[0], 999.0);
@@ -145,5 +66,254 @@
     BOOST_CHECK_EQUAL(outim[5], 999.0);
 }
 
+BOOST_AUTO_TEST_CASE(r_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 outre[] = { 999, 999, 999, 999, 999, 999 };
+    FFTReal(4).inverse(re, im, outre+1);
+    // And check we haven't overrun the arrays
+    BOOST_CHECK_EQUAL(outre[0], 999.0);
+    BOOST_CHECK_EQUAL(outre[5], 999.0);
+}
+
+BOOST_AUTO_TEST_CASE(dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double in[] = { 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_dc)
+{
+    // DC-only signal. The DC bin is purely real
+    double in[] = { 1, 1, 1, 1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 4.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(sine)
+{
+    // Sine. Output is purely imaginary
+    double in[] = { 0, 1, 0, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    COMPARE_CONST(re, 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], -2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 2.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_sine)
+{
+    // Sine. Output is purely imaginary
+    double in[] = { 0, 1, 0, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    COMPARE_CONST(re, 0.0);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_EQUAL(im[1], -2.0);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_EQUAL(im[3], 2.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(cosine)
+{
+    // Cosine. Output is purely real
+    double in[] = { 1, 0, -1, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 2.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_cosine)
+{
+    // Cosine. Output is purely real
+    double in[] = { 1, 0, -1, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 2.0);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_EQUAL(re[3], 2.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+	
+BOOST_AUTO_TEST_CASE(sineCosine)
+{
+    // Sine and cosine mixed
+    double in[] = { 0.5, 1, -0.5, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+	
+BOOST_AUTO_TEST_CASE(r_sineCosine)
+{
+    // Sine and cosine mixed
+    double in[] = { 0.5, 1, -0.5, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_CLOSE(re[1], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(re[2], 0.0);
+    BOOST_CHECK_CLOSE(re[3], 1.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[0], 0.0);
+    BOOST_CHECK_CLOSE(im[1], -2.0, 1e-12);
+    BOOST_CHECK_EQUAL(im[2], 0.0);
+    BOOST_CHECK_CLOSE(im[3], 2.0, 1e-12);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(nyquist)
+{
+    double in[] = { 1, -1, 1, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 4.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_nyquist)
+{
+    double in[] = { 1, -1, 1, -1 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 0.0);
+    BOOST_CHECK_EQUAL(re[1], 0.0);
+    BOOST_CHECK_EQUAL(re[2], 4.0);
+    BOOST_CHECK_EQUAL(re[3], 0.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
+BOOST_AUTO_TEST_CASE(dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFT(4).process(false, in, 0, re, im);
+    BOOST_CHECK_EQUAL(re[0], 1.0);
+    BOOST_CHECK_EQUAL(re[1], 1.0);
+    BOOST_CHECK_EQUAL(re[2], 1.0);
+    BOOST_CHECK_EQUAL(re[3], 1.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    double backim[4];
+    FFT(4).process(true, re, im, back, backim);
+    COMPARE_ARRAY(back, in);
+    COMPARE_CONST(backim, 0.0);
+}
+
+BOOST_AUTO_TEST_CASE(r_dirac)
+{
+    double in[] = { 1, 0, 0, 0 };
+    double re[] = { 999, 999, 999, 999 };
+    double im[] = { 999, 999, 999, 999 };
+    FFTReal(4).forward(in, re, im);
+    BOOST_CHECK_EQUAL(re[0], 1.0);
+    BOOST_CHECK_EQUAL(re[1], 1.0);
+    BOOST_CHECK_EQUAL(re[2], 1.0);
+    BOOST_CHECK_EQUAL(re[3], 1.0);
+    COMPARE_CONST(im, 0.0);
+    double back[4];
+    // check conjugates are reconstructed
+    re[3] = 999;
+    im[3] = 999;
+    FFTReal(4).inverse(re, im, back);
+    COMPARE_ARRAY(back, in);
+}
+
 BOOST_AUTO_TEST_SUITE_END()