Mercurial > hg > qm-dsp

--- 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<double>(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();
--- 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<double> *m_window;
     PhaseVocoder* m_phaseVoc;	// Phase Vocoder
--- 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 <math.h>

-//////////////////////////////////////////////////////////////////////
-// Construction/Destruction
-//////////////////////////////////////////////////////////////////////
+#include <iostream>
+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];
+    }
+}
+
--- 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
--- 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()