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;
  44
  45     m_DFType = Config.DFType;
  46     m_stepSize = Config.stepSize;
  47
  48     m_whiten = Config.adaptiveWhitening;
  49     m_whitenRelaxCoeff = Config.whiteningRelaxCoeff;
  50     m_whitenFloor = Config.whiteningFloor;
  51     if (m_whitenRelaxCoeff < 0) m_whitenRelaxCoeff = 0.9997;
  52     if (m_whitenFloor < 0) m_whitenFloor = 0.01;
  53
  54     m_magHistory = new double[ m_halfLength ];
  55     memset(m_magHistory,0, m_halfLength*sizeof(double));
  56
  57     m_phaseHistory = new double[ m_halfLength ];
  58     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
  59
  60     m_phaseHistoryOld = new double[ m_halfLength ];
  61     memset(m_phaseHistoryOld,0, m_halfLength*sizeof(double));
  62
  63     m_magPeaks = new double[ m_halfLength ];
  64     memset(m_magPeaks,0, m_halfLength*sizeof(double));
  65
  66     // See note in process(const double *) below
  67     int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
  68     m_phaseVoc = new PhaseVocoder(actualLength);
  69
  70     m_DFWindowedFrame = new double[ m_dataLength ];
  71     m_magnitude = new double[ m_halfLength ];
  72     m_thetaAngle = new double[ m_halfLength ];
  73
  74     m_window = new Window<double>(HanningWindow, 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_DFWindowedFrame;
  87     delete [] m_magnitude;
  88     delete [] m_thetaAngle;
  89
  90     delete m_window;
  91 }
  92
  93 double DetectionFunction::process( const double *TDomain )
  94 {
  95     m_window->cut( TDomain, m_DFWindowedFrame );
  96
  97     // Our own FFT implementation supports power-of-two sizes only.
  98     // If we have to use this implementation (as opposed to the
  99     // version of process() below that operates on frequency domain
 100     // data directly), we will have to use the next smallest power of
 101     // two from the block size.  Results may vary accordingly!
 102
 103     unsigned int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
 104
 105     if (actualLength != m_dataLength) {
 106         // Pre-fill mag and phase vectors with zero, as the FFT output
 107         // will not fill the arrays
 108         for (unsigned int i = actualLength/2; i < m_dataLength/2; ++i) {
 109             m_magnitude[i] = 0;
 110             m_thetaAngle[0] = 0;
 111         }
 112     }
 113
 114     m_phaseVoc->process(m_DFWindowedFrame, m_magnitude, m_thetaAngle);
 115
 116     if (m_whiten) whiten();
 117
 118     return runDF();
 119 }
 120
 121 double DetectionFunction::process( const double *magnitudes, const double *phases )
 122 {
 123     for (size_t i = 0; i < m_halfLength; ++i) {
 124         m_magnitude[i] = magnitudes[i];
 125         m_thetaAngle[i] = phases[i];
 126     }
 127
 128     if (m_whiten) whiten();
 129
 130     return runDF();
 131 }
 132
 133 void DetectionFunction::whiten()
 134 {
 135     for (unsigned int i = 0; i < m_halfLength; ++i) {
 136         double m = m_magnitude[i];
 137         if (m < m_magPeaks[i]) {
 138             m = m + (m_magPeaks[i] - m) * m_whitenRelaxCoeff;
 139         }
 140         if (m < m_whitenFloor) m = m_whitenFloor;
 141         m_magPeaks[i] = m;
 142         m_magnitude[i] /= m;
 143     }
 144 }
 145
 146 double DetectionFunction::runDF()
 147 {
 148     double retVal = 0;
 149
 150     switch( m_DFType )
 151     {
 152     case DF_HFC:
 153         retVal = HFC( m_halfLength, m_magnitude);
 154         break;
 155
 156     case DF_SPECDIFF:
 157         retVal = specDiff( m_halfLength, m_magnitude);
 158         break;
 159
 160     case DF_PHASEDEV:
 161         retVal = phaseDev( m_halfLength, m_thetaAngle);
 162         break;
 163
 164     case DF_COMPLEXSD:
 165         retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
 166         break;
 167
 168     case DF_BROADBAND:
 169         retVal = broadband( m_halfLength, m_magnitude);
 170         break;
 171     }
 172
 173     return retVal;
 174 }
 175
 176 double DetectionFunction::HFC(unsigned int length, double *src)
 177 {
 178     unsigned int i;
 179     double val = 0;
 180
 181     for( i = 0; i < length; i++)
 182     {
 183         val += src[ i ] * ( i + 1);
 184     }
 185     return val;
 186 }
 187
 188 double DetectionFunction::specDiff(unsigned int length, double *src)
 189 {
 190     unsigned int i;
 191     double val = 0.0;
 192     double temp = 0.0;
 193     double diff = 0.0;
 194
 195     for( i = 0; i < length; i++)
 196     {
 197         temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
 198
 199         diff= sqrt(temp);
 200
 201         // (See note in phaseDev below.)
 202
 203         val += diff;
 204
 205         m_magHistory[ i ] = src[ i ];
 206     }
 207
 208     return val;
 209 }
 210
 211
 212 double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
 213 {
 214     unsigned int i;
 215     double tmpPhase = 0;
 216     double tmpVal = 0;
 217     double val = 0;
 218
 219     double dev = 0;
 220
 221     for( i = 0; i < length; i++)
 222     {
 223         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
 224         dev = MathUtilities::princarg( tmpPhase );
 225
 226         // A previous version of this code only counted the value here
 227         // if the magnitude exceeded 0.1.  My impression is that
 228         // doesn't greatly improve the results for "loud" music (so
 229         // long as the peak picker is reasonably sophisticated), but
 230         // does significantly damage its ability to work with quieter
 231         // music, so I'm removing it and counting the result always.
 232         // Same goes for the spectral difference measure above.
 233
 234         tmpVal  = fabs(dev);
 235         val += tmpVal ;
 236
 237         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
 238         m_phaseHistory[ i ] = srcPhase[ i ];
 239     }
 240
 241
 242     return val;
 243 }
 244
 245
 246 double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
 247 {
 248     unsigned int i;
 249     double val = 0;
 250     double tmpPhase = 0;
 251     double tmpReal = 0;
 252     double tmpImag = 0;
 253
 254     double dev = 0;
 255     ComplexData meas = ComplexData( 0, 0 );
 256     ComplexData j = ComplexData( 0, 1 );
 257
 258     for( i = 0; i < length; i++)
 259     {
 260         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
 261         dev= MathUtilities::princarg( tmpPhase );
 262
 263         meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
 264
 265         tmpReal = real( meas );
 266         tmpImag = imag( meas );
 267
 268         val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
 269
 270         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
 271         m_phaseHistory[ i ] = srcPhase[ i ];
 272         m_magHistory[ i ] = srcMagnitude[ i ];
 273     }
 274
 275     return val;
 276 }
 277
 278 double DetectionFunction::broadband(unsigned int length, double *src)
 279 {
 280     double val = 0;
 281     for (unsigned int i = 0; i < length; ++i) {
 282         double sqrmag = src[i] * src[i];
 283         if (m_magHistory[i] > 0.0) {
 284             double diff = 10.0 * log10(sqrmag / m_magHistory[i]);
 285             if (diff > m_dbRise) val = val + 1;
 286         }
 287         m_magHistory[i] = sqrmag;
 288     }
 289     return val;
 290 }
 291
 292 double* DetectionFunction::getSpectrumMagnitude()
 293 {
 294     return m_magnitude;
 295 }
 296