# HG changeset patch
# User Chris Cannam
# Date 1434128937 -3600
# Node ID bdebff3265ae4ede7cad64b512d6fa43be12ea27
# Parent  420fc961c0c460ea259af1e2b80c3b16da95525f
Simplest naive FFTModel implementation (+ fill in tests)

diff -r 420fc961c0c4 -r bdebff3265ae data/fft/FFTapi.h
--- a/data/fft/FFTapi.h	Fri Jun 12 14:51:46 2015 +0100
+++ b/data/fft/FFTapi.h	Fri Jun 12 18:08:57 2015 +0100
@@ -50,5 +50,46 @@
 
 #endif
 
+#include <vector>
+#include <complex>
+
+class FFTForward // with fft shift but not window
+{
+public:
+    FFTForward(int size) :
+        m_size(size),
+        m_input((float *)fftf_malloc(size * sizeof(float))),
+        m_output((fftf_complex *)fftf_malloc((size/2 + 1) * sizeof(fftf_complex))),
+        m_plan(fftf_plan_dft_r2c_1d(size, m_input, m_output, FFTW_MEASURE))
+    { }
+
+    ~FFTForward() {
+        fftf_destroy_plan(m_plan);
+        fftf_free(m_input);
+        fftf_free(m_output);
+    }
+
+    std::vector<std::complex<float> > process(std::vector<float> in) const {
+        const int hs = m_size/2;
+        for (int i = 0; i < hs; ++i) {
+            m_input[i] = in[i + hs];
+            m_input[i + hs] = in[i];
+        }
+        fftf_execute(m_plan);
+        std::vector<std::complex<float> > result;
+        result.reserve(hs + 1);
+        for (int i = 0; i <= hs; ++i) {
+            result[i] = { m_output[i][0], m_output[i][1] };
+        }
+        return result;
+    }
+
+private:
+    int m_size;
+    float *m_input;
+    fftf_complex *m_output;
+    fftf_plan m_plan;
+};
+
 #endif
 
diff -r 420fc961c0c4 -r bdebff3265ae data/model/FFTModel.cpp
--- a/data/model/FFTModel.cpp	Fri Jun 12 14:51:46 2015 +0100
+++ b/data/model/FFTModel.cpp	Fri Jun 12 18:08:57 2015 +0100
@@ -42,8 +42,14 @@
     m_windowSize(windowSize),
     m_windowIncrement(windowIncrement),
     m_fftSize(fftSize),
-    m_windower(windowType, windowSize)
+    m_windower(windowType, windowSize),
+    m_fft(fftSize)
 {
+    if (m_windowSize > m_fftSize) {
+        cerr << "ERROR: FFTModel::FFTModel: window size (" << m_windowSize
+             << ") must be at least FFT size (" << m_fftSize << ")" << endl;
+        throw invalid_argument("FFTModel window size must be at least FFT size");
+    }
 }
 
 FFTModel::~FFTModel()
@@ -59,6 +65,20 @@
     }
 }
 
+int
+FFTModel::getWidth() const
+{
+    if (!m_model) return 0;
+    return int((m_model->getEndFrame() - m_model->getStartFrame())
+               / m_windowIncrement) + 1;
+}
+
+int
+FFTModel::getHeight() const
+{
+    return m_fftSize / 2 + 1;
+}
+
 QString
 FFTModel::getBinName(int n) const
 {
@@ -68,6 +88,131 @@
     return name;
 }
 
+FFTModel::Column
+FFTModel::getColumn(int x) const
+{
+    auto cplx = getFFTColumn(x);
+    Column col;
+    col.reserve(int(cplx.size()));
+    for (auto c: cplx) col.push_back(abs(c));
+    return col;
+}
+
+float
+FFTModel::getMagnitudeAt(int x, int y) const
+{
+    //!!! 
+    return abs(getFFTColumn(x)[y]);
+}
+
+float
+FFTModel::getMaximumMagnitudeAt(int x) const
+{
+    Column col(getColumn(x));
+    auto itr = max_element(col.begin(), col.end());
+    if (itr == col.end()) return 0.f;
+    else return *itr;
+}
+
+float
+FFTModel::getPhaseAt(int x, int y) const
+{
+    //!!! 
+    return arg(getFFTColumn(x)[y]);
+}
+
+void
+FFTModel::getValuesAt(int x, int y, float &re, float &im) const
+{
+    auto col = getFFTColumn(x);
+    re = col[y].real();
+    im = col[y].imag();
+}
+
+bool
+FFTModel::isColumnAvailable(int ) const
+{
+    //!!!
+    return true;
+}
+
+bool
+FFTModel::getMagnitudesAt(int x, float *values, int minbin, int count) const
+{
+    if (count == 0) count = getHeight();
+    auto col = getFFTColumn(x);
+    for (int i = 0; i < count; ++i) {
+        values[i] = abs(col[minbin + i]);
+    }
+    return true;
+}
+
+bool
+FFTModel::getNormalizedMagnitudesAt(int x, float *values, int minbin, int count) const
+{
+    //!!! WRONG
+    return getMagnitudesAt(x, values, minbin, count);
+}
+
+bool
+FFTModel::getPhasesAt(int x, float *values, int minbin, int count) const
+{
+    if (count == 0) count = getHeight();
+    auto col = getFFTColumn(x);
+    for (int i = 0; i < count; ++i) {
+        values[i] = arg(col[minbin + i]);
+    }
+    return true;
+}
+
+bool
+FFTModel::getValuesAt(int x, float *reals, float *imags, int minbin, int count) const
+{
+    if (count == 0) count = getHeight();
+    auto col = getFFTColumn(x);
+    for (int i = 0; i < count; ++i) {
+        reals[i] = col[minbin + i].real();
+    }
+    for (int i = 0; i < count; ++i) {
+        imags[i] = col[minbin + i].imag();
+    }
+    return true;
+}
+
+vector<float>
+FFTModel::getSourceSamples(int column) const
+{
+    auto range = getSourceSampleRange(column);
+    vector<float> samples(m_fftSize, 0.f);
+    int off = (m_fftSize - m_windowSize) / 2;
+    decltype(range.first) pfx = 0;
+    if (range.first < 0) {
+        pfx = -range.first;
+        range = { 0, range.second };
+    }
+    (void) m_model->getData(m_channel,
+                            range.first,
+                            range.second - range.first,
+                            &samples[off + pfx]);
+    if (m_channel == -1) {
+	int channels = m_model->getChannelCount();
+	if (channels > 1) {
+	    for (int i = 0; i < range.second - range.first; ++i) {
+		samples[off + pfx + i] /= float(channels);
+	    }
+	}
+    }
+    return samples;
+}
+
+vector<complex<float>>
+FFTModel::getFFTColumn(int column) const
+{
+    auto samples = getSourceSamples(column);
+    m_windower.cut(&samples[0]);
+    return m_fft.process(samples);
+}
+
 bool
 FFTModel::estimateStableFrequency(int x, int y, double &frequency)
 {
@@ -239,14 +384,14 @@
     if (type == MajorPeaks) return 10;
     if (bin == 0) return 3;
 
-    double binfreq = (getSampleRate() * bin) / m_fftSize;
+    double binfreq = (sampleRate * bin) / m_fftSize;
     double hifreq = Pitch::getFrequencyForPitch(73, 0, binfreq);
 
-    int hibin = int(lrint((hifreq * m_fftSize) / getSampleRate()));
+    int hibin = int(lrint((hifreq * m_fftSize) / sampleRate));
     int medianWinSize = hibin - bin;
     if (medianWinSize < 3) medianWinSize = 3;
 
-    percentile = 0.5f + float(binfreq / getSampleRate());
+    percentile = 0.5f + float(binfreq / sampleRate);
 
     return medianWinSize;
 }
diff -r 420fc961c0c4 -r bdebff3265ae data/model/FFTModel.h
--- a/data/model/FFTModel.h	Fri Jun 12 14:51:46 2015 +0100
+++ b/data/model/FFTModel.h	Fri Jun 12 18:08:57 2015 +0100
@@ -21,8 +21,12 @@
 
 #include "base/Window.h"
 
+#include "data/fft/FFTapi.h"
+
 #include <set>
 #include <map>
+#include <vector>
+#include <complex>
 
 /**
  * An implementation of DenseThreeDimensionalModel that makes FFT data
@@ -88,7 +92,6 @@
     int getFFTSize() const { return m_fftSize; }
     
     float getMagnitudeAt(int x, int y) const;
-    float getNormalizedMagnitudeAt(int x, int y) const;
     float getMaximumMagnitudeAt(int x) const;
     float getPhaseAt(int x, int y) const;
     void getValuesAt(int x, int y, float &real, float &imaginary) const;
@@ -144,9 +147,22 @@
     int m_windowIncrement;
     int m_fftSize;
     Window<float> m_windower;
+    FFTForward m_fft;
     
     int getPeakPickWindowSize(PeakPickType type, sv_samplerate_t sampleRate,
                               int bin, float &percentile) const;
+
+    std::pair<sv_frame_t, sv_frame_t> getSourceSampleRange(int column) const {
+        sv_frame_t startFrame = m_windowIncrement * sv_frame_t(column);
+        sv_frame_t endFrame = startFrame + m_windowSize;
+        // Cols are centred on the audio sample (e.g. col 0 is centred at sample 0)
+        startFrame -= m_windowSize / 2;
+        endFrame -= m_windowSize / 2;
+        return { startFrame, endFrame };
+    }
+
+    std::vector<std::complex<float> > getFFTColumn(int column) const;
+    std::vector<float> getSourceSamples(int column) const;
 };
 
 #endif
diff -r 420fc961c0c4 -r bdebff3265ae data/model/test/MockWaveModel.cpp
--- a/data/model/test/MockWaveModel.cpp	Fri Jun 12 14:51:46 2015 +0100
+++ b/data/model/test/MockWaveModel.cpp	Fri Jun 12 18:08:57 2015 +0100
@@ -30,17 +30,17 @@
 {
     sv_frame_t i = 0;
 
-    cerr << "MockWaveModel::getData(" << channel << "," << start << "," << count << "): ";
+//    cerr << "MockWaveModel::getData(" << channel << "," << start << "," << count << "): ";
 
     while (i < count) {
 	sv_frame_t idx = start + i;
 	if (!in_range_for(m_data[channel], idx)) break;
 	buffer[i] = m_data[channel][idx];
-	cerr << buffer[i] << " ";
+//	cerr << buffer[i] << " ";
 	++i;
     }
 
-    cerr << endl;
+//    cerr << endl;
     
     return i;
 }
diff -r 420fc961c0c4 -r bdebff3265ae data/model/test/TestFFTModel.h
--- a/data/model/test/TestFFTModel.h	Fri Jun 12 14:51:46 2015 +0100
+++ b/data/model/test/TestFFTModel.h	Fri Jun 12 18:08:57 2015 +0100
@@ -40,40 +40,38 @@
               int columnNo, vector<vector<complex<float>>> expectedValues,
               int expectedWidth) {
         for (int ch = 0; in_range_for(expectedValues, ch); ++ch) {
-            for (int polar = 0; polar <= 1; ++polar) {
-                FFTModel fftm(model, ch, window, windowSize, windowIncrement,
-                              fftSize, bool(polar));
-                QCOMPARE(fftm.getWidth(), expectedWidth);
-                int hs1 = fftSize/2 + 1;
-                QCOMPARE(fftm.getHeight(), hs1);
-                vector<float> reals(hs1 + 1, 0.f);
-                vector<float> imags(hs1 + 1, 0.f);
-                reals[hs1] = 999.f; // overrun guards
-                imags[hs1] = 999.f;
-                fftm.getValuesAt(columnNo, &reals[0], &imags[0]);
-                for (int i = 0; i < hs1; ++i) {
-                    float eRe = expectedValues[ch][i].real();
-                    float eIm = expectedValues[ch][i].imag();
-                    if (reals[i] != eRe || imags[i] != eIm) {
-                        cerr << "ERROR: output is not as expected for column "
-                             << i << " in channel " << ch << " (polar store = "
-                             << polar << ")" << endl;
-                        cerr << "expected : ";
-                        for (int j = 0; j < hs1; ++j) {
-                            cerr << expectedValues[ch][j] << " ";
-                        }
-                        cerr << "\nactual   : ";
-                        for (int j = 0; j < hs1; ++j) {
-                            cerr << complex<float>(reals[j], imags[j]) << " ";
-                        }
-                        cerr << endl;
+            FFTModel fftm(model, ch, window, windowSize, windowIncrement, fftSize);
+            QCOMPARE(fftm.getWidth(), expectedWidth);
+            int hs1 = fftSize/2 + 1;
+            QCOMPARE(fftm.getHeight(), hs1);
+            vector<float> reals(hs1 + 1, 0.f);
+            vector<float> imags(hs1 + 1, 0.f);
+            reals[hs1] = 999.f; // overrun guards
+            imags[hs1] = 999.f;
+            fftm.getValuesAt(columnNo, &reals[0], &imags[0]);
+            for (int i = 0; i < hs1; ++i) {
+                float eRe = expectedValues[ch][i].real();
+                float eIm = expectedValues[ch][i].imag();
+                float thresh = 1e-5f;
+                if (abs(reals[i] - eRe) > thresh ||
+                    abs(imags[i] - eIm) > thresh) {
+                    cerr << "ERROR: output is not as expected for column "
+                         << i << " in channel " << ch << endl;
+                    cerr << "expected : ";
+                    for (int j = 0; j < hs1; ++j) {
+                        cerr << expectedValues[ch][j] << " ";
                     }
-                    QCOMPARE(reals[i], eRe);
-                    QCOMPARE(imags[i], eIm);
+                    cerr << "\nactual   : ";
+                    for (int j = 0; j < hs1; ++j) {
+                        cerr << complex<float>(reals[j], imags[j]) << " ";
+                    }
+                    cerr << endl;
                 }
-                QCOMPARE(reals[hs1], 999.f);
-                QCOMPARE(imags[hs1], 999.f);
+                COMPARE_FUZZIER_F(reals[i], eRe);
+                COMPARE_FUZZIER_F(imags[i], eIm);
             }
+            QCOMPARE(reals[hs1], 999.f);
+            QCOMPARE(imags[hs1], 999.f);
         }
     }
 
@@ -117,6 +115,10 @@
     
     void sine_simple_rect() {
 	MockWaveModel mwm({ Sine }, 16, 4);
+        // Sine: output is purely imaginary. Note the sign is flipped
+        // (normally the first half of the output would have negative
+        // sign for a sine starting at 0) because the model does an
+        // FFT shift to centre the phase
         test(&mwm, RectangularWindow, 8, 8, 8, 0,
              { { {}, {}, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 1,
@@ -129,39 +131,62 @@
     
     void cosine_simple_rect() {
 	MockWaveModel mwm({ Cosine }, 16, 4);
+        // Cosine: output is purely real. Note the sign is flipped
+        // because the model does an FFT shift to centre the phase
         test(&mwm, RectangularWindow, 8, 8, 8, 0,
              { { {}, {}, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 1,
-             { { {}, { 2.f, 0.f }, {}, {}, {} } }, 4);
+             { { {}, { -2.f, 0.f }, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 2,
-             { { {}, { 2.f, 0.f }, {}, {}, {} } }, 4);
+             { { {}, { -2.f, 0.f }, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 3,
              { { {}, {}, {}, {}, {} } }, 4);
     }
     
     void nyquist_simple_rect() {
 	MockWaveModel mwm({ Nyquist }, 16, 4);
+        // Again, the sign is flipped. This has the same amount of
+        // energy as the DC example
         test(&mwm, RectangularWindow, 8, 8, 8, 0,
              { { {}, {}, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 1,
-             { { {}, {}, {}, {}, { 2.f, 0.f } } }, 4);
+             { { {}, {}, {}, {}, { -4.f, 0.f } } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 2,
-             { { {}, {}, {}, {}, { 2.f, 0.f } } }, 4);
+             { { {}, {}, {}, {}, { -4.f, 0.f } } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 3,
              { { {}, {}, {}, {}, {} } }, 4);
     }
     
     void dirac_simple_rect() {
 	MockWaveModel mwm({ Dirac }, 16, 4);
+        // The window scales by 0.5 and some signs are flipped. Only
+        // column 1 has any data (the single impulse).
         test(&mwm, RectangularWindow, 8, 8, 8, 0,
              { { {}, {}, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 1,
-             { { { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f } } }, 4);
+             { { { 0.5f, 0.f }, { -0.5f, 0.f }, { 0.5f, 0.f }, { -0.5f, 0.f }, { 0.5f, 0.f } } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 2,
-             { { { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f }, { 1.f, 0.f } } }, 4);
+             { { {}, {}, {}, {}, {} } }, 4);
         test(&mwm, RectangularWindow, 8, 8, 8, 3,
              { { {}, {}, {}, {}, {} } }, 4);
     }
+    
+    void dirac_simple_rect_2() {
+	MockWaveModel mwm({ Dirac }, 16, 8);
+        // With 8 samples padding, the FFT shift places the first
+        // Dirac impulse at the start of column 1, thus giving all
+        // positive values
+        test(&mwm, RectangularWindow, 8, 8, 8, 0,
+             { { {}, {}, {}, {}, {} } }, 5);
+        test(&mwm, RectangularWindow, 8, 8, 8, 1,
+             { { { 0.5f, 0.f }, { 0.5f, 0.f }, { 0.5f, 0.f }, { 0.5f, 0.f }, { 0.5f, 0.f } } }, 5);
+        test(&mwm, RectangularWindow, 8, 8, 8, 2,
+             { { {}, {}, {}, {}, {} } }, 5);
+        test(&mwm, RectangularWindow, 8, 8, 8, 3,
+             { { {}, {}, {}, {}, {} } }, 5);
+        test(&mwm, RectangularWindow, 8, 8, 8, 4,
+             { { {}, {}, {}, {}, {} } }, 5);
+    }
 
 };