libs/qm-dsp/dsp/onsets/DetectionFunction.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file 2005-2006 Christian Landone.
   8
   9     This program is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU General Public License as
  11     published by the Free Software Foundation; either version 2 of the
  12     License, or (at your option) any later version.  See the file
  13     COPYING included with this distribution for more information.
  14 */
  15
  16 #include "DetectionFunction.h"
  17 #include <cstring>
  18
  19 //////////////////////////////////////////////////////////////////////
  20 // Construction/Destruction
  21 //////////////////////////////////////////////////////////////////////
  22
  23 DetectionFunction::DetectionFunction( DFConfig Config ) :
  24     m_window(0)
  25 {
  26     m_magHistory = NULL;
  27     m_phaseHistory = NULL;
  28     m_phaseHistoryOld = NULL;
  29     m_magPeaks = NULL;
  30
  31     initialise( Config );
  32 }
  33
  34 DetectionFunction::~DetectionFunction()
  35 {
  36     deInitialise();
  37 }
  38
  39
  40 void DetectionFunction::initialise( DFConfig Config )
  41 {
  42     m_dataLength = Config.frameLength;
  43     m_halfLength = m_dataLength/2 + 1;
  44
  45     m_DFType = Config.DFType;
  46     m_stepSize = Config.stepSize;
  47     m_dbRise = Config.dbRise;
  48
  49     m_whiten = Config.adaptiveWhitening;
  50     m_whitenRelaxCoeff = Config.whiteningRelaxCoeff;
  51     m_whitenFloor = Config.whiteningFloor;
  52     if (m_whitenRelaxCoeff < 0) m_whitenRelaxCoeff = 0.9997;
  53     if (m_whitenFloor < 0) m_whitenFloor = 0.01;
  54
  55     m_magHistory = new double[ m_halfLength ];
  56     memset(m_magHistory,0, m_halfLength*sizeof(double));
  57
  58     m_phaseHistory = new double[ m_halfLength ];
  59     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
  60
  61     m_phaseHistoryOld = new double[ m_halfLength ];
  62     memset(m_phaseHistoryOld,0, m_halfLength*sizeof(double));
  63
  64     m_magPeaks = new double[ m_halfLength ];
  65     memset(m_magPeaks,0, m_halfLength*sizeof(double));
  66
  67     m_phaseVoc = new PhaseVocoder(m_dataLength, m_stepSize);
  68
  69     m_magnitude = new double[ m_halfLength ];
  70     m_thetaAngle = new double[ m_halfLength ];
  71     m_unwrapped = new double[ m_halfLength ];
  72
  73     m_window = new Window<double>(HanningWindow, m_dataLength);
  74     m_windowed = new double[ m_dataLength ];
  75 }
  76
  77 void DetectionFunction::deInitialise()
  78 {
  79     delete [] m_magHistory ;
  80     delete [] m_phaseHistory ;
  81     delete [] m_phaseHistoryOld ;
  82     delete [] m_magPeaks ;
  83
  84     delete m_phaseVoc;
  85
  86     delete [] m_magnitude;
  87     delete [] m_thetaAngle;
  88     delete [] m_windowed;
  89     delete [] m_unwrapped;
  90
  91     delete m_window;
  92 }
  93
  94 double DetectionFunction::processTimeDomain(const double *samples)
  95 {
  96     m_window->cut(samples, m_windowed);
  97
  98     m_phaseVoc->processTimeDomain(m_windowed,
  99                                   m_magnitude, m_thetaAngle, m_unwrapped);
 100
 101     if (m_whiten) whiten();
 102
 103     return runDF();
 104 }
 105
 106 double DetectionFunction::processFrequencyDomain(const double *reals,
 107                                                  const double *imags)
 108 {
 109     m_phaseVoc->processFrequencyDomain(reals, imags,
 110                                        m_magnitude, m_thetaAngle, m_unwrapped);
 111
 112     if (m_whiten) whiten();
 113
 114     return runDF();
 115 }
 116
 117 void DetectionFunction::whiten()
 118 {
 119     for (unsigned int i = 0; i < m_halfLength; ++i) {
 120         double m = m_magnitude[i];
 121         if (m < m_magPeaks[i]) {
 122             m = m + (m_magPeaks[i] - m) * m_whitenRelaxCoeff;
 123         }
 124         if (m < m_whitenFloor) m = m_whitenFloor;
 125         m_magPeaks[i] = m;
 126         m_magnitude[i] /= m;
 127     }
 128 }
 129
 130 double DetectionFunction::runDF()
 131 {
 132     double retVal = 0;
 133
 134     switch( m_DFType )
 135     {
 136     case DF_HFC:
 137         retVal = HFC( m_halfLength, m_magnitude);
 138         break;
 139
 140     case DF_SPECDIFF:
 141         retVal = specDiff( m_halfLength, m_magnitude);
 142         break;
 143
 144     case DF_PHASEDEV:
 145         // Using the instantaneous phases here actually provides the
 146         // same results (for these calculations) as if we had used
 147         // unwrapped phases, but without the possible accumulation of
 148         // phase error over time
 149         retVal = phaseDev( m_halfLength, m_thetaAngle);
 150         break;
 151
 152     case DF_COMPLEXSD:
 153         retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
 154         break;
 155
 156     case DF_BROADBAND:
 157         retVal = broadband( m_halfLength, m_magnitude);
 158         break;
 159     }
 160
 161     return retVal;
 162 }
 163
 164 double DetectionFunction::HFC(unsigned int length, double *src)
 165 {
 166     unsigned int i;
 167     double val = 0;
 168
 169     for( i = 0; i < length; i++)
 170     {
 171         val += src[ i ] * ( i + 1);
 172     }
 173     return val;
 174 }
 175
 176 double DetectionFunction::specDiff(unsigned int length, double *src)
 177 {
 178     unsigned int i;
 179     double val = 0.0;
 180     double temp = 0.0;
 181     double diff = 0.0;
 182
 183     for( i = 0; i < length; i++)
 184     {
 185         temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
 186
 187         diff= sqrt(temp);
 188
 189         // (See note in phaseDev below.)
 190
 191         val += diff;
 192
 193         m_magHistory[ i ] = src[ i ];
 194     }
 195
 196     return val;
 197 }
 198
 199
 200 double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
 201 {
 202     unsigned int i;
 203     double tmpPhase = 0;
 204     double tmpVal = 0;
 205     double val = 0;
 206
 207     double dev = 0;
 208
 209     for( i = 0; i < length; i++)
 210     {
 211         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
 212         dev = MathUtilities::princarg( tmpPhase );
 213
 214         // A previous version of this code only counted the value here
 215         // if the magnitude exceeded 0.1.  My impression is that
 216         // doesn't greatly improve the results for "loud" music (so
 217         // long as the peak picker is reasonably sophisticated), but
 218         // does significantly damage its ability to work with quieter
 219         // music, so I'm removing it and counting the result always.
 220         // Same goes for the spectral difference measure above.
 221
 222         tmpVal  = fabs(dev);
 223         val += tmpVal ;
 224
 225         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
 226         m_phaseHistory[ i ] = srcPhase[ i ];
 227     }
 228
 229     return val;
 230 }
 231
 232
 233 double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
 234 {
 235     unsigned int i;
 236     double val = 0;
 237     double tmpPhase = 0;
 238     double tmpReal = 0;
 239     double tmpImag = 0;
 240
 241     double dev = 0;
 242     ComplexData meas = ComplexData( 0, 0 );
 243     ComplexData j = ComplexData( 0, 1 );
 244
 245     for( i = 0; i < length; i++)
 246     {
 247         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
 248         dev= MathUtilities::princarg( tmpPhase );
 249
 250         meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
 251
 252         tmpReal = real( meas );
 253         tmpImag = imag( meas );
 254
 255         val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
 256
 257         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
 258         m_phaseHistory[ i ] = srcPhase[ i ];
 259         m_magHistory[ i ] = srcMagnitude[ i ];
 260     }
 261
 262     return val;
 263 }
 264
 265 double DetectionFunction::broadband(unsigned int length, double *src)
 266 {
 267     double val = 0;
 268     for (unsigned int i = 0; i < length; ++i) {
 269         double sqrmag = src[i] * src[i];
 270         if (m_magHistory[i] > 0.0) {
 271             double diff = 10.0 * log10(sqrmag / m_magHistory[i]);
 272             if (diff > m_dbRise) val = val + 1;
 273         }
 274         m_magHistory[i] = sqrmag;
 275     }
 276     return val;
 277 }
 278
 279 double* DetectionFunction::getSpectrumMagnitude()
 280 {
 281     return m_magnitude;
 282 }
 283