qm-dsp  1.8
DetectionFunction.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 
00003 /*
00004     QM DSP Library
00005 
00006     Centre for Digital Music, Queen Mary, University of London.
00007     This file 2005-2006 Christian Landone.
00008 
00009     This program is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU General Public License as
00011     published by the Free Software Foundation; either version 2 of the
00012     License, or (at your option) any later version.  See the file
00013     COPYING included with this distribution for more information.
00014 */
00015 
00016 #include "DetectionFunction.h"
00017 #include <cstring>
00018 
00020 // Construction/Destruction
00022 
00023 DetectionFunction::DetectionFunction( DFConfig Config ) :
00024     m_window(0)
00025 {
00026     m_magHistory = NULL;
00027     m_phaseHistory = NULL;
00028     m_phaseHistoryOld = NULL;
00029     m_magPeaks = NULL;
00030 
00031     initialise( Config );
00032 }
00033 
00034 DetectionFunction::~DetectionFunction()
00035 {
00036     deInitialise();
00037 }
00038 
00039 
00040 void DetectionFunction::initialise( DFConfig Config )
00041 {
00042     m_dataLength = Config.frameLength;
00043     m_halfLength = m_dataLength/2 + 1;
00044 
00045     m_DFType = Config.DFType;
00046     m_stepSize = Config.stepSize;
00047 
00048     m_whiten = Config.adaptiveWhitening;
00049     m_whitenRelaxCoeff = Config.whiteningRelaxCoeff;
00050     m_whitenFloor = Config.whiteningFloor;
00051     if (m_whitenRelaxCoeff < 0) m_whitenRelaxCoeff = 0.9997;
00052     if (m_whitenFloor < 0) m_whitenFloor = 0.01;
00053 
00054     m_magHistory = new double[ m_halfLength ];
00055     memset(m_magHistory,0, m_halfLength*sizeof(double));
00056                 
00057     m_phaseHistory = new double[ m_halfLength ];
00058     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
00059 
00060     m_phaseHistoryOld = new double[ m_halfLength ];
00061     memset(m_phaseHistoryOld,0, m_halfLength*sizeof(double));
00062 
00063     m_magPeaks = new double[ m_halfLength ];
00064     memset(m_magPeaks,0, m_halfLength*sizeof(double));
00065 
00066     m_phaseVoc = new PhaseVocoder(m_dataLength, m_stepSize);
00067 
00068     m_magnitude = new double[ m_halfLength ];
00069     m_thetaAngle = new double[ m_halfLength ];
00070     m_unwrapped = new double[ m_halfLength ];
00071 
00072     m_window = new Window<double>(HanningWindow, m_dataLength);
00073     m_windowed = new double[ m_dataLength ];
00074 }
00075 
00076 void DetectionFunction::deInitialise()
00077 {
00078     delete [] m_magHistory ;
00079     delete [] m_phaseHistory ;
00080     delete [] m_phaseHistoryOld ;
00081     delete [] m_magPeaks ;
00082 
00083     delete m_phaseVoc;
00084 
00085     delete [] m_magnitude;
00086     delete [] m_thetaAngle;
00087     delete [] m_windowed;
00088     delete [] m_unwrapped;
00089 
00090     delete m_window;
00091 }
00092 
00093 double DetectionFunction::processTimeDomain(const double *samples)
00094 {
00095     m_window->cut(samples, m_windowed);
00096 
00097     m_phaseVoc->processTimeDomain(m_windowed, 
00098                                   m_magnitude, m_thetaAngle, m_unwrapped);
00099 
00100     if (m_whiten) whiten();
00101 
00102     return runDF();
00103 }
00104 
00105 double DetectionFunction::processFrequencyDomain(const double *reals,
00106                                                  const double *imags)
00107 {
00108     m_phaseVoc->processFrequencyDomain(reals, imags,
00109                                        m_magnitude, m_thetaAngle, m_unwrapped);
00110 
00111     if (m_whiten) whiten();
00112 
00113     return runDF();
00114 }
00115 
00116 void DetectionFunction::whiten()
00117 {
00118     for (unsigned int i = 0; i < m_halfLength; ++i) {
00119         double m = m_magnitude[i];
00120         if (m < m_magPeaks[i]) {
00121             m = m + (m_magPeaks[i] - m) * m_whitenRelaxCoeff;
00122         }
00123         if (m < m_whitenFloor) m = m_whitenFloor;
00124         m_magPeaks[i] = m;
00125         m_magnitude[i] /= m;
00126     }
00127 }
00128 
00129 double DetectionFunction::runDF()
00130 {
00131     double retVal = 0;
00132 
00133     switch( m_DFType )
00134     {
00135     case DF_HFC:
00136         retVal = HFC( m_halfLength, m_magnitude);
00137         break;
00138         
00139     case DF_SPECDIFF:
00140         retVal = specDiff( m_halfLength, m_magnitude);
00141         break;
00142         
00143     case DF_PHASEDEV:
00144         // Using the instantaneous phases here actually provides the
00145         // same results (for these calculations) as if we had used
00146         // unwrapped phases, but without the possible accumulation of
00147         // phase error over time
00148         retVal = phaseDev( m_halfLength, m_thetaAngle);
00149         break;
00150         
00151     case DF_COMPLEXSD:
00152         retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
00153         break;
00154 
00155     case DF_BROADBAND:
00156         retVal = broadband( m_halfLength, m_magnitude);
00157         break;
00158     }
00159         
00160     return retVal;
00161 }
00162 
00163 double DetectionFunction::HFC(unsigned int length, double *src)
00164 {
00165     unsigned int i;
00166     double val = 0;
00167 
00168     for( i = 0; i < length; i++)
00169     {
00170         val += src[ i ] * ( i + 1);
00171     }
00172     return val;
00173 }
00174 
00175 double DetectionFunction::specDiff(unsigned int length, double *src)
00176 {
00177     unsigned int i;
00178     double val = 0.0;
00179     double temp = 0.0;
00180     double diff = 0.0;
00181 
00182     for( i = 0; i < length; i++)
00183     {
00184         temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
00185                 
00186         diff= sqrt(temp);
00187 
00188         // (See note in phaseDev below.)
00189 
00190         val += diff;
00191 
00192         m_magHistory[ i ] = src[ i ];
00193     }
00194 
00195     return val;
00196 }
00197 
00198 
00199 double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
00200 {
00201     unsigned int i;
00202     double tmpPhase = 0;
00203     double tmpVal = 0;
00204     double val = 0;
00205 
00206     double dev = 0;
00207 
00208     for( i = 0; i < length; i++)
00209     {
00210         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
00211         dev = MathUtilities::princarg( tmpPhase );
00212 
00213         // A previous version of this code only counted the value here
00214         // if the magnitude exceeded 0.1.  My impression is that
00215         // doesn't greatly improve the results for "loud" music (so
00216         // long as the peak picker is reasonably sophisticated), but
00217         // does significantly damage its ability to work with quieter
00218         // music, so I'm removing it and counting the result always.
00219         // Same goes for the spectral difference measure above.
00220                 
00221         tmpVal  = fabs(dev);
00222         val += tmpVal ;
00223 
00224         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
00225         m_phaseHistory[ i ] = srcPhase[ i ];
00226     }
00227         
00228     return val;
00229 }
00230 
00231 
00232 double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
00233 {
00234     unsigned int i;
00235     double val = 0;
00236     double tmpPhase = 0;
00237     double tmpReal = 0;
00238     double tmpImag = 0;
00239    
00240     double dev = 0;
00241     ComplexData meas = ComplexData( 0, 0 );
00242     ComplexData j = ComplexData( 0, 1 );
00243 
00244     for( i = 0; i < length; i++)
00245     {
00246         tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
00247         dev= MathUtilities::princarg( tmpPhase );
00248                 
00249         meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
00250 
00251         tmpReal = real( meas );
00252         tmpImag = imag( meas );
00253 
00254         val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
00255                 
00256         m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
00257         m_phaseHistory[ i ] = srcPhase[ i ];
00258         m_magHistory[ i ] = srcMagnitude[ i ];
00259     }
00260 
00261     return val;
00262 }
00263 
00264 double DetectionFunction::broadband(unsigned int length, double *src)
00265 {
00266     double val = 0;
00267     for (unsigned int i = 0; i < length; ++i) {
00268         double sqrmag = src[i] * src[i];
00269         if (m_magHistory[i] > 0.0) {
00270             double diff = 10.0 * log10(sqrmag / m_magHistory[i]);
00271             if (diff > m_dbRise) val = val + 1;
00272         }
00273         m_magHistory[i] = sqrmag;
00274     }
00275     return val;
00276 }        
00277 
00278 double* DetectionFunction::getSpectrumMagnitude()
00279 {
00280     return m_magnitude;
00281 }
00282