update qm-dsp library

[ardour.git] / libs / qm-dsp / dsp / onsets / DetectionFunction.cpp

diff --git a/libs/qm-dsp/dsp/onsets/DetectionFunction.cpp b/libs/qm-dsp/dsp/onsets/DetectionFunction.cpp
index af65af8cd28f1f7f23d21da03b5793405b35817b..dfcd4da921e3b2a250eda237c503dc41445cc02f 100644 (file)
--- a/libs/qm-dsp/dsp/onsets/DetectionFunction.cpp
+++ b/libs/qm-dsp/dsp/onsets/DetectionFunction.cpp
@@ -40,10 +40,11 @@ DetectionFunction::~DetectionFunction()
 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;
+    m_dbRise = Config.dbRise;
 
     m_whiten = Config.adaptiveWhitening;
     m_whitenRelaxCoeff = Config.whiteningRelaxCoeff;
@@ -53,7 +54,7 @@ void DetectionFunction::initialise( DFConfig Config )
 
     m_magHistory = new double[ m_halfLength ];
     memset(m_magHistory,0, m_halfLength*sizeof(double));
-
+               
     m_phaseHistory = new double[ m_halfLength ];
     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
 
@@ -63,15 +64,14 @@ void DetectionFunction::initialise( DFConfig Config )
     m_magPeaks = new double[ m_halfLength ];
     memset(m_magPeaks,0, m_halfLength*sizeof(double));
 
-    // See note in process(const double *) below
-    int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
-    m_phaseVoc = new PhaseVocoder(actualLength);
+    m_phaseVoc = new PhaseVocoder(m_dataLength, 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,47 +83,31 @@ void DetectionFunction::deInitialise()
 
     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 );
-
-    // 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!
-
-    unsigned int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
-
-    if (actualLength != m_dataLength) {
-        // Pre-fill mag and phase vectors with zero, as the FFT output
-        // will not fill the arrays
-        for (unsigned int i = actualLength/2; i < m_dataLength/2; ++i) {
-            m_magnitude[i] = 0;
-            m_thetaAngle[0] = 0;
-        }
-    }
+    m_window->cut(samples, m_windowed);
 
-    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();
 
@@ -152,15 +136,19 @@ double DetectionFunction::runDF()
     case DF_HFC:
        retVal = HFC( m_halfLength, m_magnitude);
        break;
-
+       
     case DF_SPECDIFF:
        retVal = specDiff( m_halfLength, m_magnitude);
        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;
-
+       
     case DF_COMPLEXSD:
        retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
        break;
@@ -169,7 +157,7 @@ double DetectionFunction::runDF()
         retVal = broadband( m_halfLength, m_magnitude);
         break;
     }
-
+       
     return retVal;
 }
 
@@ -195,7 +183,7 @@ double DetectionFunction::specDiff(unsigned int length, double *src)
     for( i = 0; i < length; i++)
     {
        temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
-
+               
        diff= sqrt(temp);
 
         // (See note in phaseDev below.)
@@ -230,15 +218,14 @@ double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
         // does significantly damage its ability to work with quieter
         // music, so I'm removing it and counting the result always.
         // Same goes for the spectral difference measure above.
-
+               
         tmpVal  = fabs(dev);
         val += tmpVal ;
 
        m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
        m_phaseHistory[ i ] = srcPhase[ i ];
     }
-
-
+       
     return val;
 }
 
@@ -250,7 +237,7 @@ double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, d
     double tmpPhase = 0;
     double tmpReal = 0;
     double tmpImag = 0;
-
+   
     double dev = 0;
     ComplexData meas = ComplexData( 0, 0 );
     ComplexData j = ComplexData( 0, 1 );
@@ -259,14 +246,14 @@ double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, d
     {
        tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
        dev= MathUtilities::princarg( tmpPhase );
-
+               
        meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
 
        tmpReal = real( meas );
        tmpImag = imag( meas );
 
        val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
-
+               
        m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
        m_phaseHistory[ i ] = srcPhase[ i ];
        m_magHistory[ i ] = srcMagnitude[ i ];
@@ -287,7 +274,7 @@ double DetectionFunction::broadband(unsigned int length, double *src)
         m_magHistory[i] = sqrmag;
     }
     return val;
-}
+}        
 
 double* DetectionFunction::getSpectrumMagnitude()
 {