qm-dsp
1.8
|
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