# HG changeset patch # User Chris Cannam # Date 1380722734 -3600 # Node ID f3c69325cca2bdebd375b988d55de181b32e3dc9 # Parent f6ccde089491e445141a866bacf59c997a5f0b7f Do actual phase unwrapping in the phase vocoder! diff -r f6ccde089491 -r f3c69325cca2 dsp/onsets/DetectionFunction.cpp --- a/dsp/onsets/DetectionFunction.cpp Wed Oct 02 15:04:38 2013 +0100 +++ b/dsp/onsets/DetectionFunction.cpp Wed Oct 02 15:05:34 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; @@ -65,11 +65,12 @@ // See note in process(const double *) 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(HanningWindow, m_dataLength); } @@ -111,7 +112,8 @@ } } - m_phaseVoc->process(m_DFWindowedFrame, m_magnitude, m_thetaAngle); + m_phaseVoc->process(m_DFWindowedFrame, m_magnitude, + m_thetaAngle, m_unwrapped); if (m_whiten) whiten(); diff -r f6ccde089491 -r f3c69325cca2 dsp/onsets/DetectionFunction.h --- a/dsp/onsets/DetectionFunction.h Wed Oct 02 15:04:38 2013 +0100 +++ b/dsp/onsets/DetectionFunction.h Wed Oct 02 15:05:34 2013 +0100 @@ -77,6 +77,7 @@ double* m_DFWindowedFrame; // 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 *m_window; PhaseVocoder* m_phaseVoc; // Phase Vocoder diff -r f6ccde089491 -r f3c69325cca2 dsp/phasevocoder/PhaseVocoder.cpp --- a/dsp/phasevocoder/PhaseVocoder.cpp Wed Oct 02 15:04:38 2013 +0100 +++ b/dsp/phasevocoder/PhaseVocoder.cpp Wed Oct 02 15:05:34 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,43 @@ #include "PhaseVocoder.h" #include "dsp/transforms/FFT.h" +#include "maths/MathUtilities.h" #include -////////////////////////////////////////////////////////////////////// -// Construction/Destruction -////////////////////////////////////////////////////////////////////// +#include +using std::cerr; +using std::endl; -PhaseVocoder::PhaseVocoder(unsigned int n) : - m_n(n) +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_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_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 +59,61 @@ } } -void PhaseVocoder::process(double *src, double *mag, double *theta) +void PhaseVocoder::process(double *src, double *mag, double *theta, + double *unwrapped) { - FFTShift( m_n, src); + FFTShift(src); - m_fft->process(0, src, m_realOut, m_imagOut); + m_fft->forward(src, m_real, m_imag); - getMagnitude( m_n/2, mag, m_realOut, m_imagOut); - getPhase( m_n/2, theta, m_realOut, m_imagOut); + getMagnitudes(mag); + getPhases(theta); + unwrapPhases(theta, unwrapped); } -void PhaseVocoder::getMagnitude(unsigned int size, double *mag, double *real, double *imag) -{ - unsigned int j; - - 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) +{ + cerr << "PhaseVocoder::unwrapPhases" << endl; + + 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; + + cerr << "i = " << i << ", instantaneous phase = " << theta[i] << ", prev phase = " << m_phase[i] << ", omega = " << omega << ", expected = " << expected << ", error = " << error << ", unwrapped = " << unwrapped[i] << endl; + + m_phase[i] = theta[i]; + m_unwrapped[i] = unwrapped[i]; + } +} + diff -r f6ccde089491 -r f3c69325cca2 dsp/phasevocoder/PhaseVocoder.h --- a/dsp/phasevocoder/PhaseVocoder.h Wed Oct 02 15:04:38 2013 +0100 +++ b/dsp/phasevocoder/PhaseVocoder.h Wed Oct 02 15:05:34 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,44 @@ class PhaseVocoder { public: - PhaseVocoder( unsigned int size ); + 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 process(double *src, 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_imag; + double *m_real; + double *m_phase; + double *m_unwrapped; }; #endif diff -r f6ccde089491 -r f3c69325cca2 tests/TestPhaseVocoder.cpp --- a/tests/TestPhaseVocoder.cpp Wed Oct 02 15:04:38 2013 +0100 +++ b/tests/TestPhaseVocoder.cpp Wed Oct 02 15:05:34 2013 +0100 @@ -27,21 +27,27 @@ BOOST_AUTO_TEST_CASE(fullcycle) { - // Cosine with one cycle exactly equal to pvoc hopsize. We aren't - // windowing the input frame because (for once) it actually *is* - // just a short part of a continuous infinite sinusoid. + // 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); + 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.process(frame, mag + 1, phase + 1); + pvoc.process(frame, mag + 1, phase + 1, unw + 1); double magExpected0[] = { 999, 0, 0, 4, 0, 0, 999 }; COMPARE_ARRAY_EXACT(mag, magExpected0); @@ -49,22 +55,53 @@ double phaseExpected0[] = { 999, 0, 0, 0, 0, 0, 999 }; COMPARE_ARRAY_EXACT(phase, phaseExpected0); - pvoc.process(frame, mag + 1, phase + 1); + double unwExpected0[] = { 999, 0, 0, 0, 0, 0, 999 }; + COMPARE_ARRAY(unw, unwExpected0); + + pvoc.process(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, 2 * M_PI, 0, 0, 999 }; + double phaseExpected1[] = { 999, 0, 0, 0, 0, 0, 999 }; COMPARE_ARRAY(phase, phaseExpected1); - pvoc.process(frame, mag + 1, phase + 1); + // Derivation of 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 unwrapped 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) is like bin 2: 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.process(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, 4 * M_PI, 0, 0, 999 }; + 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); } +//!!! signal that starts mid-phase + + BOOST_AUTO_TEST_SUITE_END()