Mercurial > hg > qm-dsp
changeset 358:7a225d665ed2
Merge from branch kissfft
| author | Chris Cannam <c.cannam@qmul.ac.uk> |
|---|---|
| date | Wed, 16 Oct 2013 12:52:44 +0100 |
| parents | 5ff562f5b521 (current diff) 650bbacf8288 (diff) |
| children | 0ea56b09577e |
| files | .hgsubstate |
| diffstat | 11 files changed, 247 insertions(+), 199 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.hgsub Wed Oct 16 12:52:44 2013 +0100 @@ -0,0 +1,1 @@ +ext/kissfft = http://hg.code.sf.net/p/kissfft/code
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.hgsubstate Wed Oct 16 12:52:44 2013 +0100 @@ -0,0 +1,1 @@ +fbe1bb0bc7b94ec252842b8b7e3f3347ec75d92f ext/kissfft
--- a/dsp/onsets/DetectionFunction.cpp Wed Oct 16 12:52:37 2013 +0100 +++ b/dsp/onsets/DetectionFunction.cpp Wed Oct 16 12:52:44 2013 +0100 @@ -63,9 +63,7 @@ m_magPeaks = new double[ m_halfLength ]; memset(m_magPeaks,0, m_halfLength*sizeof(double)); - // See note in processTimeDomain below - int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength); - m_phaseVoc = new PhaseVocoder(actualLength, m_stepSize); + m_phaseVoc = new PhaseVocoder(m_dataLength, m_stepSize); m_magnitude = new double[ m_halfLength ]; m_thetaAngle = new double[ m_halfLength ]; @@ -96,23 +94,6 @@ { m_window->cut(samples, m_windowed); - // Our own FFT implementation supports power-of-two sizes only. - // If we have to use this implementation (as opposed to the - // version of process() below that operates on frequency domain - // data directly), we will have to use the next smallest power of - // two from the block size. Results may vary accordingly! - - int actualLength = MathUtilities::previousPowerOfTwo((int)m_dataLength); - - if (actualLength != (int)m_dataLength) { - // Pre-fill mag and phase vectors with zero, as the FFT output - // will not fill the arrays - for (int i = actualLength/2; i < (int)m_dataLength/2; ++i) { - m_magnitude[i] = 0; - m_thetaAngle[0] = 0; - } - } - m_phaseVoc->processTimeDomain(m_windowed, m_magnitude, m_thetaAngle, m_unwrapped);
--- a/dsp/onsets/DetectionFunction.h Wed Oct 16 12:52:37 2013 +0100 +++ b/dsp/onsets/DetectionFunction.h Wed Oct 16 12:52:44 2013 +0100 @@ -29,7 +29,7 @@ struct DFConfig{ unsigned int stepSize; // DF step in samples - unsigned int frameLength; // DF analysis window - usually 2*step + unsigned int frameLength; // DF analysis window - usually 2*step. Must be even! int DFType; // type of detection function ( see defines ) double dbRise; // only used for broadband df (and required for it) bool adaptiveWhitening; // perform adaptive whitening
--- a/dsp/phasevocoder/PhaseVocoder.h Wed Oct 16 12:52:37 2013 +0100 +++ b/dsp/phasevocoder/PhaseVocoder.h Wed Oct 16 12:52:44 2013 +0100 @@ -21,7 +21,6 @@ class PhaseVocoder { public: - //!!! review: size must be a power of two, or not? PhaseVocoder(int size, int hop); virtual ~PhaseVocoder();
--- a/dsp/transforms/FFT.cpp Wed Oct 16 12:52:37 2013 +0100 +++ b/dsp/transforms/FFT.cpp Wed Oct 16 12:52:44 2013 +0100 @@ -11,195 +11,190 @@ #include "maths/MathUtilities.h" +#include "kiss_fft.h" +#include "kiss_fftr.h" + #include <cmath> #include <iostream> +#include <stdexcept> + +class FFT::D +{ +public: + D(int n) : m_n(n) { + m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL); + m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL); + m_kin = new kiss_fft_cpx[m_n]; + m_kout = new kiss_fft_cpx[m_n]; + } + + ~D() { + kiss_fft_free(m_planf); + kiss_fft_free(m_plani); + delete[] m_kin; + delete[] m_kout; + } + + void process(bool inverse, + const double *ri, + const double *ii, + double *ro, + double *io) { + + for (int i = 0; i < m_n; ++i) { + m_kin[i].r = ri[i]; + m_kin[i].i = (ii ? ii[i] : 0.0); + } + + if (!inverse) { + + kiss_fft(m_planf, m_kin, m_kout); + + for (int i = 0; i < m_n; ++i) { + ro[i] = m_kout[i].r; + io[i] = m_kout[i].i; + } + + } else { + + kiss_fft(m_plani, m_kin, m_kout); + + double scale = 1.0 / m_n; + + for (int i = 0; i < m_n; ++i) { + ro[i] = m_kout[i].r * scale; + io[i] = m_kout[i].i * scale; + } + } + } + +private: + int m_n; + kiss_fft_cfg m_planf; + kiss_fft_cfg m_plani; + kiss_fft_cpx *m_kin; + kiss_fft_cpx *m_kout; +}; + FFT::FFT(int n) : - m_n(n) + m_d(new D(n)) { - if( !MathUtilities::isPowerOfTwo(m_n) ) - { - std::cerr << "ERROR: FFT: Non-power-of-two FFT size " - << m_n << " not supported in this implementation" - << std::endl; - return; - } } FFT::~FFT() { - + delete m_d; } +void +FFT::process(bool inverse, + const double *p_lpRealIn, const double *p_lpImagIn, + double *p_lpRealOut, double *p_lpImagOut) +{ + m_d->process(inverse, + p_lpRealIn, p_lpImagIn, + p_lpRealOut, p_lpImagOut); +} + +class FFTReal::D +{ +public: + D(int n) : m_n(n) { + if (n % 2) { + throw std::invalid_argument + ("nsamples must be even in FFTReal constructor"); + } + m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL); + m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL); + m_c = new kiss_fft_cpx[m_n]; + } + + ~D() { + kiss_fftr_free(m_planf); + kiss_fftr_free(m_plani); + delete[] m_c; + } + + void forward(const double *ri, double *ro, double *io) { + + kiss_fftr(m_planf, ri, m_c); + + for (int i = 0; i <= m_n/2; ++i) { + ro[i] = m_c[i].r; + io[i] = m_c[i].i; + } + + for (int i = 0; i + 1 < m_n/2; ++i) { + ro[m_n - i - 1] = ro[i + 1]; + io[m_n - i - 1] = -io[i + 1]; + } + } + + void forwardMagnitude(const double *ri, double *mo) { + + double *io = new double[m_n]; + + forward(ri, mo, io); + + for (int i = 0; i < m_n; ++i) { + mo[i] = sqrt(mo[i] * mo[i] + io[i] * io[i]); + } + + delete[] io; + } + + void inverse(const double *ri, const double *ii, double *ro) { + + for (int i = 0; i < m_n; ++i) { + m_c[i].r = ri[i]; + m_c[i].i = ii[i]; + } + + kiss_fftri(m_plani, m_c, ro); + + double scale = 1.0 / m_n; + + for (int i = 0; i < m_n; ++i) { + ro[i] *= scale; + } + } + +private: + int m_n; + kiss_fftr_cfg m_planf; + kiss_fftr_cfg m_plani; + kiss_fft_cpx *m_c; +}; + FFTReal::FFTReal(int n) : - m_n(n), - m_fft(new FFT(n)), - m_r(new double[n]), - m_i(new double[n]), - m_discard(new double[n]) + m_d(new D(n)) { - for (int i = 0; i < n; ++i) { - m_r[i] = 0; - m_i[i] = 0; - m_discard[i] = 0; - } } FFTReal::~FFTReal() { - delete m_fft; - delete[] m_discard; - delete[] m_r; - delete[] m_i; + delete m_d; } void -FFTReal::forward(const double *realIn, - double *realOut, double *imagOut) +FFTReal::forward(const double *ri, double *ro, double *io) { - m_fft->process(false, realIn, 0, realOut, imagOut); + m_d->forward(ri, ro, io); } void -FFTReal::inverse(const double *realIn, const double *imagIn, - double *realOut) +FFTReal::forwardMagnitude(const double *ri, double *mo) { - 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; - - if( p_nSamples < 2 ) - { - return 0; - } - - for ( i=0; ; i++ ) - { - if( p_nSamples & (1 << i) ) return i; - } -} - -static int reverseBits(int p_nIndex, int p_nBits) -{ - int i, rev; - - for(i=rev=0; i < p_nBits; i++) - { - rev = (rev << 1) | (p_nIndex & 1); - p_nIndex >>= 1; - } - - return rev; + m_d->forwardMagnitude(ri, mo); } void -FFT::process(bool p_bInverseTransform, - const double *p_lpRealIn, const double *p_lpImagIn, - double *p_lpRealOut, double *p_lpImagOut) +FFTReal::inverse(const double *ri, const double *ii, double *ro) { - if (!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return; - -// std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl; - - int NumBits; - int i, j, k, n; - int BlockSize, BlockEnd; - - double angle_numerator = 2.0 * M_PI; - double tr, ti; - - if( !MathUtilities::isPowerOfTwo(m_n) ) - { - std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size " - << m_n << " not supported in this implementation" - << std::endl; - return; - } - - if( p_bInverseTransform ) angle_numerator = -angle_numerator; - - NumBits = numberOfBitsNeeded ( m_n ); - - - - for( i=0; i < m_n; i++ ) - { - j = reverseBits ( i, NumBits ); - p_lpRealOut[j] = p_lpRealIn[i]; - p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i]; - } - - - BlockEnd = 1; - for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 ) - { - double delta_angle = angle_numerator / (double)BlockSize; - double sm2 = -sin ( -2 * delta_angle ); - double sm1 = -sin ( -delta_angle ); - double cm2 = cos ( -2 * delta_angle ); - double cm1 = cos ( -delta_angle ); - double w = 2 * cm1; - double ar[3], ai[3]; - - for( i=0; i < m_n; i += BlockSize ) - { - - ar[2] = cm2; - ar[1] = cm1; - - ai[2] = sm2; - ai[1] = sm1; - - for ( j=i, n=0; n < BlockEnd; j++, n++ ) - { - - ar[0] = w*ar[1] - ar[2]; - ar[2] = ar[1]; - ar[1] = ar[0]; - - ai[0] = w*ai[1] - ai[2]; - ai[2] = ai[1]; - ai[1] = ai[0]; - - k = j + BlockEnd; - tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k]; - ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k]; - - p_lpRealOut[k] = p_lpRealOut[j] - tr; - p_lpImagOut[k] = p_lpImagOut[j] - ti; - - p_lpRealOut[j] += tr; - p_lpImagOut[j] += ti; - - } - } - - BlockEnd = BlockSize; - - } - - - if( p_bInverseTransform ) - { - double denom = (double)m_n; - - for ( i=0; i < m_n; i++ ) - { - p_lpRealOut[i] /= denom; - p_lpImagOut[i] /= denom; - } - } + m_d->inverse(ri, ii, ro); } + +
--- a/dsp/transforms/FFT.h Wed Oct 16 12:52:37 2013 +0100 +++ b/dsp/transforms/FFT.h Wed Oct 16 12:52:44 2013 +0100 @@ -14,8 +14,8 @@ public: /** * Construct an FFT object to carry out complex-to-complex - * transforms of size nsamples. nsamples must be a power of two in - * this implementation. + * transforms of size nsamples. nsamples does not have to be a + * power of two. */ FFT(int nsamples); ~FFT(); @@ -39,7 +39,8 @@ double *realOut, double *imagOut); private: - int m_n; + class D; + D *m_d; }; class FFTReal @@ -47,8 +48,10 @@ public: /** * Construct an FFT object to carry out real-to-complex transforms - * of size nsamples. nsamples must be a power of two in this - * implementation. + * of size nsamples. nsamples does not have to be a power of two, + * but it does have to be even. (Use the complex-complex FFT above + * if you need an odd FFT size. This constructor will throw + * std::invalid_argument if nsamples is odd.) */ FFTReal(int nsamples); ~FFTReal(); @@ -66,6 +69,18 @@ double *realOut, double *imagOut); /** + * Carry out a forward real-to-complex transform of size nsamples, + * where nsamples is the value provided to the constructor + * above. Return only the magnitudes of the complex output values. + * + * realIn and magOut 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 forwardMagnitude(const double *realIn, double *magOut); + + /** * Carry out an inverse real transform (i.e. complex-to-real) of * size nsamples, where nsamples is the value provided to the * constructor above. @@ -83,11 +98,8 @@ double *realOut); private: - int m_n; - FFT *m_fft; - double *m_r; - double *m_i; - double *m_discard; + class D; + D *m_d; }; #endif
--- a/qm-dsp.pro Wed Oct 16 12:52:37 2013 +0100 +++ b/qm-dsp.pro Wed Oct 16 12:52:44 2013 +0100 @@ -5,6 +5,12 @@ OBJECTS_DIR = tmp_obj MOC_DIR = tmp_moc +# This one is vital +QMAKE_CFLAGS += -Dkiss_fft_scalar=double +QMAKE_CXXFLAGS += -Dkiss_fft_scalar=double + +INCLUDEPATH += ./ext/kissfft/tools ./ext/kissfft + linux-g++* { QMAKE_CXXFLAGS_RELEASE += -DNDEBUG -O3 -fno-exceptions -fPIC -ffast-math -msse -mfpmath=sse -ftree-vectorize -fomit-frame-pointer DEFINES += USE_PTHREADS @@ -82,7 +88,9 @@ maths/pca/pca.h \ thread/AsynchronousTask.h \ thread/BlockAllocator.h \ - thread/Thread.h + thread/Thread.h \ + ext/kissfft/kiss_fft.h \ + ext/kissfft/tools/kiss_fftr.h SOURCES += base/Pitch.cpp \ base/KaiserWindow.cpp \ base/SincWindow.cpp \ @@ -117,4 +125,6 @@ maths/KLDivergence.cpp \ maths/MathUtilities.cpp \ maths/pca/pca.c \ - thread/Thread.cpp + thread/Thread.cpp \ + ext/kissfft/kiss_fft.c \ + ext/kissfft/tools/kiss_fftr.c
--- a/tests/Makefile Wed Oct 16 12:52:37 2013 +0100 +++ b/tests/Makefile Wed Oct 16 12:52:44 2013 +0100 @@ -7,8 +7,10 @@ TESTS := test-mathutilities test-window test-fft test-pvoc +#VG := valgrind + all: $(TESTS) - for t in $(TESTS); do echo "Running $$t"; ./"$$t" || exit 1; done + for t in $(TESTS); do echo "Running $$t"; $(VG) ./"$$t" || exit 1; done test-mathutilities: TestMathUtilities.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS)
--- a/tests/TestFFT.cpp Wed Oct 16 12:52:37 2013 +0100 +++ b/tests/TestFFT.cpp Wed Oct 16 12:52:44 2013 +0100 @@ -7,6 +7,8 @@ #include <boost/test/unit_test.hpp> +#include <stdexcept> + BOOST_AUTO_TEST_SUITE(TestFFT) #define COMPARE_CONST(a, n) \ @@ -316,5 +318,48 @@ COMPARE_ARRAY(back, in); } +BOOST_AUTO_TEST_CASE(sizes) +{ + // Complex supports any size. A single test with an odd size + // will do here, without getting too much into our expectations + // about supported butterflies etc + + double in[] = { 1, 1, 1 }; + double re[] = { 999, 999, 999 }; + double im[] = { 999, 999, 999 }; + FFT(3).process(false, in, 0, re, im); + BOOST_CHECK_EQUAL(re[0], 3.0); + BOOST_CHECK_EQUAL(re[1], 0.0); + BOOST_CHECK_EQUAL(re[2], 0.0); + COMPARE_CONST(im, 0.0); + double back[3]; + double backim[3]; + FFT(3).process(true, re, im, back, backim); + COMPARE_ARRAY(back, in); + COMPARE_CONST(backim, 0.0); +} + +BOOST_AUTO_TEST_CASE(r_sizes) +{ + // Real supports any even size, but not odd ones + + BOOST_CHECK_THROW(FFTReal r(3), std::invalid_argument); + + double in[] = { 1, 1, 1, 1, 1, 1 }; + double re[] = { 999, 999, 999, 999, 999, 999 }; + double im[] = { 999, 999, 999, 999, 999, 999 }; + FFTReal(6).forward(in, re, im); + BOOST_CHECK_EQUAL(re[0], 6.0); + BOOST_CHECK_EQUAL(re[1], 0.0); + BOOST_CHECK_EQUAL(re[2], 0.0); + BOOST_CHECK_EQUAL(re[3], 0.0); + BOOST_CHECK_EQUAL(re[4], 0.0); + BOOST_CHECK_EQUAL(re[5], 0.0); + COMPARE_CONST(im, 0.0); + double back[6]; + FFTReal(6).inverse(re, im, back); + COMPARE_ARRAY(back, in); +} + BOOST_AUTO_TEST_SUITE_END()
--- a/tests/TestPhaseVocoder.cpp Wed Oct 16 12:52:37 2013 +0100 +++ b/tests/TestPhaseVocoder.cpp Wed Oct 16 12:52:44 2013 +0100 @@ -164,14 +164,16 @@ // phase of -pi/2 plus advance of 2pi). It has the same measured // phase as last time around, -pi/2, which is consistent with // the expected phase, so the unwrapped phase is 3pi/2. - //!!! - // * Bin 3 is a bit of a puzzle -- it has an effectively zero - // magnitude but a non-zero measured phase. Spectral leakage? + // + // * Bin 3 I don't really know about -- the magnitude here is 0, + // but we get non-zero measured phase whose sign is + // implementation-dependent // // * Bin 4 (Nyquist) has expected phase 4*pi, measured phase 0, // hence error 0 and unwrapped phase 4*pi. - double phaseExpected1[] = { 999, 0, -M_PI, -M_PI/2, M_PI, 0, 999 }; + phase[1+3] = 0.0; // Because we aren't testing for this one + double phaseExpected1[] = { 999, 0, -M_PI, -M_PI/2, 0, 0, 999 }; COMPARE_ARRAY(phase, phaseExpected1); double unwExpected1[] = { 999, 0, M_PI, 3*M_PI/2, 3*M_PI, 4*M_PI, 999 };