Mercurial > hg > svcore
changeset 1091:bdebff3265ae simple-fft-model
Simplest naive FFTModel implementation (+ fill in tests)
| author | Chris Cannam |
|---|---|
| date | Fri, 12 Jun 2015 18:08:57 +0100 |
| parents | 420fc961c0c4 |
| children | 70f18770b72d |
| files | data/fft/FFTapi.h data/model/FFTModel.cpp data/model/FFTModel.h data/model/test/MockWaveModel.cpp data/model/test/TestFFTModel.h |
| diffstat | 5 files changed, 272 insertions(+), 45 deletions(-) [+] |
line wrap: on
line diff
--- 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
--- 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; }
--- 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
--- 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; }
--- 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); + } };