Mercurial > hg > qm-dsp
changeset 347:e3dedded9c4d
Merge from pvoc branch
| author | Chris Cannam <c.cannam@qmul.ac.uk> |
|---|---|
| date | Fri, 04 Oct 2013 16:43:44 +0100 |
| parents | f665f9ce2fd1 (current diff) 58ba20857a5e (diff) |
| children | 50fae18236ee |
| files | |
| diffstat | 15 files changed, 750 insertions(+), 194 deletions(-) [+] |
line wrap: on
line diff
--- a/dsp/chromagram/Chromagram.cpp Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/chromagram/Chromagram.cpp Fri Oct 04 16:43:44 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:15:59 2013 +0100 +++ b/dsp/mfcc/MFCC.cpp Fri Oct 04 16:43:44 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/onsets/DetectionFunction.cpp Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/onsets/DetectionFunction.cpp Fri Oct 04 16:43:44 2013 +0100 @@ -40,7 +40,7 @@ void DetectionFunction::initialise( DFConfig Config ) { m_dataLength = Config.frameLength; - m_halfLength = m_dataLength/2; + m_halfLength = m_dataLength/2 + 1; m_DFType = Config.DFType; m_stepSize = Config.stepSize; @@ -63,15 +63,16 @@ m_magPeaks = new double[ m_halfLength ]; memset(m_magPeaks,0, m_halfLength*sizeof(double)); - // See note in process(const double *) below + // See note in processTimeDomain below int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength); - m_phaseVoc = new PhaseVocoder(actualLength); + m_phaseVoc = new PhaseVocoder(actualLength, m_stepSize); - m_DFWindowedFrame = new double[ m_dataLength ]; m_magnitude = new double[ m_halfLength ]; m_thetaAngle = new double[ m_halfLength ]; + m_unwrapped = new double[ m_halfLength ]; m_window = new Window<double>(HanningWindow, m_dataLength); + m_windowed = new double[ m_dataLength ]; } void DetectionFunction::deInitialise() @@ -83,16 +84,17 @@ delete m_phaseVoc; - delete [] m_DFWindowedFrame; delete [] m_magnitude; delete [] m_thetaAngle; + delete [] m_windowed; + delete [] m_unwrapped; delete m_window; } -double DetectionFunction::process( const double *TDomain ) +double DetectionFunction::processTimeDomain(const double *samples) { - m_window->cut( TDomain, m_DFWindowedFrame ); + 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 @@ -111,19 +113,19 @@ } } - m_phaseVoc->process(m_DFWindowedFrame, m_magnitude, m_thetaAngle); + m_phaseVoc->processTimeDomain(m_windowed, + m_magnitude, m_thetaAngle, m_unwrapped); if (m_whiten) whiten(); return runDF(); } -double DetectionFunction::process( const double *magnitudes, const double *phases ) +double DetectionFunction::processFrequencyDomain(const double *reals, + const double *imags) { - for (size_t i = 0; i < m_halfLength; ++i) { - m_magnitude[i] = magnitudes[i]; - m_thetaAngle[i] = phases[i]; - } + m_phaseVoc->processFrequencyDomain(reals, imags, + m_magnitude, m_thetaAngle, m_unwrapped); if (m_whiten) whiten(); @@ -158,6 +160,10 @@ break; case DF_PHASEDEV: + // Using the instantaneous phases here actually provides the + // same results (for these calculations) as if we had used + // unwrapped phases, but without the possible accumulation of + // phase error over time retVal = phaseDev( m_halfLength, m_thetaAngle); break; @@ -238,7 +244,6 @@ m_phaseHistory[ i ] = srcPhase[ i ]; } - return val; }
--- a/dsp/onsets/DetectionFunction.h Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/onsets/DetectionFunction.h Fri Oct 04 16:43:44 2013 +0100 @@ -43,8 +43,18 @@ double* getSpectrumMagnitude(); DetectionFunction( DFConfig Config ); virtual ~DetectionFunction(); - double process( const double* TDomain ); - double process( const double* magnitudes, const double* phases ); + + /** + * Process a single time-domain frame of audio, provided as + * frameLength samples. + */ + double processTimeDomain(const double* samples); + + /** + * Process a single frequency-domain frame, provided as + * frameLength/2+1 real and imaginary component values. + */ + double processFrequencyDomain(const double* reals, const double* imags); private: void whiten(); @@ -74,9 +84,10 @@ double* m_phaseHistoryOld; double* m_magPeaks; - double* m_DFWindowedFrame; // Array for windowed analysis frame + double* m_windowed; // Array for windowed analysis frame double* m_magnitude; // Magnitude of analysis frame ( frequency domain ) double* m_thetaAngle;// Phase of analysis frame ( frequency domain ) + double* m_unwrapped; // Unwrapped phase of analysis frame Window<double> *m_window; PhaseVocoder* m_phaseVoc; // Phase Vocoder
--- a/dsp/phasevocoder/PhaseVocoder.cpp Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/phasevocoder/PhaseVocoder.cpp Fri Oct 04 16:43:44 2013 +0100 @@ -4,7 +4,7 @@ QM DSP Library Centre for Digital Music, Queen Mary, University of London. - This file 2005-2006 Christian Landone. + This file 2005-2006 Christian Landone, copyright 2013 QMUL. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as @@ -15,30 +15,47 @@ #include "PhaseVocoder.h" #include "dsp/transforms/FFT.h" +#include "maths/MathUtilities.h" #include <math.h> -////////////////////////////////////////////////////////////////////// -// Construction/Destruction -////////////////////////////////////////////////////////////////////// +#include <cassert> -PhaseVocoder::PhaseVocoder(unsigned int n) : - m_n(n) +#include <iostream> +using std::cerr; +using std::endl; + +PhaseVocoder::PhaseVocoder(int n, int hop) : + m_n(n), + m_hop(hop) { m_fft = new FFTReal(m_n); - m_realOut = new double[m_n]; - m_imagOut = new double[m_n]; + m_time = new double[m_n]; + m_real = new double[m_n]; + m_imag = new double[m_n]; + m_phase = new double[m_n/2 + 1]; + m_unwrapped = new double[m_n/2 + 1]; + + for (int i = 0; i < m_n/2 + 1; ++i) { + m_phase[i] = 0.0; + m_unwrapped[i] = 0.0; + } + + reset(); } PhaseVocoder::~PhaseVocoder() { - delete [] m_realOut; - delete [] m_imagOut; + delete[] m_unwrapped; + delete[] m_phase; + delete[] m_real; + delete[] m_imag; + delete[] m_time; delete m_fft; } -void PhaseVocoder::FFTShift(unsigned int size, double *src) +void PhaseVocoder::FFTShift(double *src) { - const int hs = size/2; + const int hs = m_n/2; for (int i = 0; i < hs; ++i) { double tmp = src[i]; src[i] = src[i + hs]; @@ -46,34 +63,73 @@ } } -void PhaseVocoder::process(double *src, double *mag, double *theta) +void PhaseVocoder::processTimeDomain(const double *src, + double *mag, double *theta, + double *unwrapped) { - FFTShift( m_n, src); - - m_fft->process(0, src, m_realOut, m_imagOut); - - getMagnitude( m_n/2, mag, m_realOut, m_imagOut); - getPhase( m_n/2, theta, m_realOut, m_imagOut); + for (int i = 0; i < m_n; ++i) { + m_time[i] = src[i]; + } + FFTShift(m_time); + m_fft->forward(m_time, m_real, m_imag); + getMagnitudes(mag); + getPhases(theta); + unwrapPhases(theta, unwrapped); } -void PhaseVocoder::getMagnitude(unsigned int size, double *mag, double *real, double *imag) -{ - unsigned int j; +void PhaseVocoder::processFrequencyDomain(const double *reals, + const double *imags, + double *mag, double *theta, + double *unwrapped) +{ + for (int i = 0; i < m_n/2 + 1; ++i) { + m_real[i] = reals[i]; + m_imag[i] = imags[i]; + } + getMagnitudes(mag); + getPhases(theta); + unwrapPhases(theta, unwrapped); +} - for( j = 0; j < size; j++) - { - mag[ j ] = sqrt( real[ j ] * real[ j ] + imag[ j ] * imag[ j ]); +void PhaseVocoder::reset() +{ + for (int i = 0; i < m_n/2 + 1; ++i) { + // m_phase stores the "previous" phase, so set to one step + // behind so that a signal with initial phase at zero matches + // the expected values. This is completely unnecessary for any + // analytical purpose, it's just tidier. + double omega = (2 * M_PI * m_hop * i) / m_n; + m_phase[i] = -omega; + m_unwrapped[i] = -omega; } } -void PhaseVocoder::getPhase(unsigned int size, double *theta, double *real, double *imag) +void PhaseVocoder::getMagnitudes(double *mag) +{ + for (int i = 0; i < m_n/2 + 1; i++) { + mag[i] = sqrt(m_real[i] * m_real[i] + m_imag[i] * m_imag[i]); + } +} + +void PhaseVocoder::getPhases(double *theta) { - unsigned int k; - - // Phase Angle "matlab" style - //Watch out for quadrant mapping !!! - for( k = 0; k < size; k++) - { - theta[ k ] = atan2( -imag[ k ], real[ k ]); + for (int i = 0; i < m_n/2 + 1; i++) { + theta[i] = atan2(m_imag[i], m_real[i]); } } + +void PhaseVocoder::unwrapPhases(double *theta, double *unwrapped) +{ + for (int i = 0; i < m_n/2 + 1; ++i) { + + double omega = (2 * M_PI * m_hop * i) / m_n; + double expected = m_phase[i] + omega; + double error = MathUtilities::princarg(theta[i] - expected); + + unwrapped[i] = m_unwrapped[i] + omega + error; + + m_phase[i] = theta[i]; + m_unwrapped[i] = unwrapped[i]; + } +} +
--- a/dsp/phasevocoder/PhaseVocoder.h Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/phasevocoder/PhaseVocoder.h Fri Oct 04 16:43:44 2013 +0100 @@ -4,7 +4,7 @@ QM DSP Library Centre for Digital Music, Queen Mary, University of London. - This file 2005-2006 Christian Landone. + This file 2005-2006 Christian Landone, copyright 2013 QMUL. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as @@ -21,22 +21,61 @@ class PhaseVocoder { public: - PhaseVocoder( unsigned int size ); + //!!! review: size must be a power of two, or not? + PhaseVocoder(int size, int hop); virtual ~PhaseVocoder(); - void process( double* src, double* mag, double* theta); + /** + * Given one frame of time-domain samples, FFT and return the + * magnitudes, instantaneous phases, and unwrapped phases. + * + * src must have size values (where size is the frame size value + * as passed to the PhaseVocoder constructor), and should have + * been windowed as necessary by the caller (but not fft-shifted). + * + * mag, phase, and unwrapped must each be non-NULL and point to + * enough space for size/2 + 1 values. The redundant conjugate + * half of the output is not returned. + */ + void processTimeDomain(const double *src, + double *mag, double *phase, double *unwrapped); + + /** + * Given one frame of frequency-domain samples, return the + * magnitudes, instantaneous phases, and unwrapped phases. + * + * reals and imags must each contain size/2+1 values (where size + * is the frame size value as passed to the PhaseVocoder + * constructor). + * + * mag, phase, and unwrapped must each be non-NULL and point to + * enough space for size/2+1 values. + */ + void processFrequencyDomain(const double *reals, const double *imags, + double *mag, double *phase, double *unwrapped); + + /** + * Reset the stored phases to zero. Note that this may be + * necessary occasionally (depending on the application) to avoid + * loss of floating-point precision in the accumulated unwrapped + * phase values as they grow. + */ + void reset(); protected: - void getPhase(unsigned int size, double *theta, double *real, double *imag); -// void coreFFT( unsigned int NumSamples, double *RealIn, double* ImagIn, double *RealOut, double *ImagOut); - void getMagnitude( unsigned int size, double* mag, double* real, double* imag); - void FFTShift( unsigned int size, double* src); + void FFTShift(double *src); + void getMagnitudes(double *mag); + void getPhases(double *theta); + void unwrapPhases(double *theta, double *unwrapped); - unsigned int m_n; + int m_n; + int m_hop; FFTReal *m_fft; - double *m_imagOut; - double *m_realOut; - + double *m_time; + double *m_imag; + double *m_real; + double *m_phase; + double *m_unwrapped; }; #endif
--- a/dsp/segmentation/ClusterMeltSegmenter.cpp Tue Oct 01 15:15:59 2013 +0100 +++ b/dsp/segmentation/ClusterMeltSegmenter.cpp Fri Oct 04 16:43:44 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:15:59 2013 +0100 +++ b/dsp/tempotracking/DownBeat.cpp Fri Oct 04 16:43:44 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:15:59 2013 +0100 +++ b/dsp/transforms/FFT.cpp Fri Oct 04 16:43:44 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:15:59 2013 +0100 +++ b/dsp/transforms/FFT.h Fri Oct 04 16:43:44 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/qm-dsp.pro Tue Oct 01 15:15:59 2013 +0100 +++ b/qm-dsp.pro Fri Oct 04 16:43:44 2013 +0100 @@ -1,6 +1,6 @@ TEMPLATE = lib -CONFIG += staticlib warn_on release +CONFIG += staticlib warn_on debug CONFIG -= qt OBJECTS_DIR = tmp_obj MOC_DIR = tmp_moc
--- a/tests/Makefile Tue Oct 01 15:15:59 2013 +0100 +++ b/tests/Makefile Fri Oct 04 16:43:44 2013 +0100 @@ -1,14 +1,14 @@ CFLAGS := -I.. $(CFLAGS) -CXXFLAGS := -I.. $(CXXFLAGS) +CXXFLAGS := -I.. -Wall -g $(CXXFLAGS) LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework LIBS := ../libqm-dsp.a -TESTS := test-window test-fft +TESTS := test-window test-fft test-pvoc all: $(TESTS) - for t in $(TESTS); do echo "Running $$t"; ./"$$t" || exit 1; done + for t in $(TESTS); do echo "Running $$t"; valgrind ./"$$t" || exit 1; done test-window: TestWindow.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) @@ -16,6 +16,13 @@ test-fft: TestFFT.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) +test-pvoc: TestPhaseVocoder.o $(LIBS) + $(CXX) -o $@ $^ $(LDFLAGS) + +TestWindow.o: $(LIBS) +TestFFT.o: $(LIBS) +TestPhaseVocoder.o: $(LIBS) + clean: - rm *.o $(TESTS) + rm -f *.o $(TESTS)
--- a/tests/TestFFT.cpp Tue Oct 01 15:15:59 2013 +0100 +++ b/tests/TestFFT.cpp Fri Oct 04 16:43:44 2013 +0100 @@ -1,3 +1,4 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ #include "dsp/transforms/FFT.h" @@ -18,101 +19,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 +36,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 +67,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()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/TestPhaseVocoder.cpp Fri Oct 04 16:43:44 2013 +0100 @@ -0,0 +1,193 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +#include "dsp/phasevocoder/PhaseVocoder.h" + +#include "base/Window.h" + +#include <iostream> + +using std::cerr; +using std::endl; + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MAIN + +#include <boost/test/unit_test.hpp> + +BOOST_AUTO_TEST_SUITE(TestFFT) + +#define COMPARE_CONST(a, n) \ + for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \ + BOOST_CHECK_SMALL(a[cmp_i] - n, 1e-7); \ + } + +#define COMPARE_ARRAY(a, b) \ + for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \ + BOOST_CHECK_SMALL(a[cmp_i] - b[cmp_i], 1e-7); \ + } + +#define COMPARE_ARRAY_EXACT(a, b) \ + for (int cmp_i = 0; cmp_i < (int)(sizeof(a)/sizeof(a[0])); ++cmp_i) { \ + BOOST_CHECK_EQUAL(a[cmp_i], b[cmp_i]); \ + } + +BOOST_AUTO_TEST_CASE(fullcycle) +{ + // Cosine with one cycle exactly equal to pvoc hopsize. This is + // pretty much the most trivial case -- in fact it's + // indistinguishable from totally silent input (in the phase + // values) because the measured phases are zero throughout. + + // We aren't windowing the input frame because (for once) it + // actually *is* just a short part of a continuous infinite + // sinusoid. + + double frame[] = { 1, 0, -1, 0, 1, 0, -1, 0 }; + + PhaseVocoder pvoc(8, 4); + + // Make these arrays one element too long at each end, so as to + // test for overruns. For frame size 8, we expect 8/2+1 = 5 + // mag/phase pairs. + double mag[] = { 999, 999, 999, 999, 999, 999, 999 }; + double phase[] = { 999, 999, 999, 999, 999, 999, 999 }; + double unw[] = { 999, 999, 999, 999, 999, 999, 999 }; + + pvoc.processTimeDomain(frame, mag + 1, phase + 1, unw + 1); + + double magExpected0[] = { 999, 0, 0, 4, 0, 0, 999 }; + COMPARE_ARRAY_EXACT(mag, magExpected0); + + double phaseExpected0[] = { 999, 0, 0, 0, 0, 0, 999 }; + COMPARE_ARRAY(phase, phaseExpected0); + + double unwExpected0[] = { 999, 0, 0, 0, 0, 0, 999 }; + COMPARE_ARRAY(unw, unwExpected0); + + pvoc.processTimeDomain(frame, mag + 1, phase + 1, unw + 1); + + double magExpected1[] = { 999, 0, 0, 4, 0, 0, 999 }; + COMPARE_ARRAY_EXACT(mag, magExpected1); + + double phaseExpected1[] = { 999, 0, 0, 0, 0, 0, 999 }; + COMPARE_ARRAY(phase, phaseExpected1); + + // Derivation of unwrapped values: + // + // * Bin 0 (DC) always has phase 0 and expected phase 0 + // + // * Bin 1 has expected phase pi (the hop size is half a cycle at + // its frequency), but measured phase 0 (because there is no + // signal in that bin). So it has phase error -pi, which is + // mapped into (-pi,pi] range as pi, giving an unwrapped phase + // of 2*pi. + // + // * Bin 2 has expected phase 2*pi, measured phase 0, hence error + // 0 and unwrapped phase 2*pi. + // + // * Bin 3 is like bin 1: it has expected phase 3*pi, measured + // phase 0, so phase error -pi and unwrapped phase 4*pi. + // + // * Bin 4 (Nyquist) has expected phase 4*pi, measured phase 0, + // hence error 0 and unwrapped phase 4*pi. + + double unwExpected1[] = { 999, 0, 2*M_PI, 2*M_PI, 4*M_PI, 4*M_PI, 999 }; + COMPARE_ARRAY(unw, unwExpected1); + + pvoc.processTimeDomain(frame, mag + 1, phase + 1, unw + 1); + + double magExpected2[] = { 999, 0, 0, 4, 0, 0, 999 }; + COMPARE_ARRAY_EXACT(mag, magExpected2); + + double phaseExpected2[] = { 999, 0, 0, 0, 0, 0, 999 }; + COMPARE_ARRAY(phase, phaseExpected2); + + double unwExpected2[] = { 999, 0, 4*M_PI, 4*M_PI, 8*M_PI, 8*M_PI, 999 }; + COMPARE_ARRAY(unw, unwExpected2); +} + +BOOST_AUTO_TEST_CASE(overlapping) +{ + // Sine (i.e. cosine starting at phase -pi/2) starting with the + // first sample, introducing a cosine of half the frequency + // starting at the fourth sample, i.e. the second hop. The cosine + // is introduced "by magic", i.e. it doesn't appear in the second + // half of the first frame (it would have quite strange effects on + // the first frame if it did). + + double data[32] = { // 3 x 8-sample frames which we pretend are overlapping + 0, 1, 0, -1, 0, 1, 0, -1, + 1, 1.70710678, 0, -1.70710678, -1, 0.29289322, 0, -0.29289322, + -1, 0.29289322, 0, -0.29289322, 1, 1.70710678, 0, -1.70710678, + }; + + PhaseVocoder pvoc(8, 4); + + // Make these arrays one element too long at each end, so as to + // test for overruns. For frame size 8, we expect 8/2+1 = 5 + // mag/phase pairs. + double mag[] = { 999, 999, 999, 999, 999, 999, 999 }; + double phase[] = { 999, 999, 999, 999, 999, 999, 999 }; + double unw[] = { 999, 999, 999, 999, 999, 999, 999 }; + + pvoc.processTimeDomain(data, mag + 1, phase + 1, unw + 1); + + double magExpected0[] = { 999, 0, 0, 4, 0, 0, 999 }; + COMPARE_ARRAY(mag, magExpected0); + + double phaseExpected0[] = { 999, 0, 0, -M_PI/2 , 0, 0, 999 }; + COMPARE_ARRAY(phase, phaseExpected0); + + double unwExpected0[] = { 999, 0, 0, -M_PI/2, 0, 0, 999 }; + COMPARE_ARRAY(unw, unwExpected0); + + pvoc.processTimeDomain(data + 8, mag + 1, phase + 1, unw + 1); + + double magExpected1[] = { 999, 0, 4, 4, 0, 0, 999 }; + COMPARE_ARRAY(mag, magExpected1); + + // Derivation of unwrapped values: + // + // * Bin 0 (DC) always has phase 0 and expected phase 0 + // + // * Bin 1 has a new signal, a cosine starting with phase 0. But + // because of the "FFT shift" which the phase vocoder carries + // out to place zero phase in the centre of the (usually + // windowed) frame, and because a single cycle at this frequency + // spans the whole frame, this bin actually has measured phase + // of either pi or -pi. (The shift doesn't affect those + // higher-frequency bins whose signals fit exact multiples of a + // cycle into a frame.) This maps into (-pi,pi] as pi, which + // matches the expected phase, hence unwrapped phase is also pi. + // + // * Bin 2 has expected phase 3pi/2 (being the previous measured + // 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 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 }; + COMPARE_ARRAY(phase, phaseExpected1); + + double unwExpected1[] = { 999, 0, M_PI, 3*M_PI/2, 3*M_PI, 4*M_PI, 999 }; + COMPARE_ARRAY(unw, unwExpected1); + + pvoc.processTimeDomain(data + 16, mag + 1, phase + 1, unw + 1); + + double magExpected2[] = { 999, 0, 4, 4, 0, 0, 999 }; + COMPARE_ARRAY(mag, magExpected2); + + double phaseExpected2[] = { 999, 0, 0, -M_PI/2, 0, 0, 999 }; + COMPARE_ARRAY(phase, phaseExpected2); + + double unwExpected2[] = { 999, 0, 2*M_PI, 7*M_PI/2, 6*M_PI, 8*M_PI, 999 }; + COMPARE_ARRAY(unw, unwExpected2); +} + +BOOST_AUTO_TEST_SUITE_END() +