c@225: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@225: 
c@225: /*
c@225:     QM DSP Library
c@225: 
c@225:     Centre for Digital Music, Queen Mary, University of London.
c@309:     This file 2005-2006 Christian Landone.
c@309: 
c@309:     This program is free software; you can redistribute it and/or
c@309:     modify it under the terms of the GNU General Public License as
c@309:     published by the Free Software Foundation; either version 2 of the
c@309:     License, or (at your option) any later version.  See the file
c@309:     COPYING included with this distribution for more information.
c@225: */
c@225: 
c@225: #include "DetectionFunction.h"
c@272: #include <cstring>
c@225: 
c@225: //////////////////////////////////////////////////////////////////////
c@225: // Construction/Destruction
c@225: //////////////////////////////////////////////////////////////////////
c@225: 
c@225: DetectionFunction::DetectionFunction( DFConfig Config ) :
c@225:     m_window(0)
c@225: {
c@227:     m_magHistory = NULL;
c@227:     m_phaseHistory = NULL;
c@227:     m_phaseHistoryOld = NULL;
c@239:     m_magPeaks = NULL;
c@225: 
c@225:     initialise( Config );
c@225: }
c@225: 
c@225: DetectionFunction::~DetectionFunction()
c@225: {
c@225:     deInitialise();
c@225: }
c@225: 
c@225: 
c@225: void DetectionFunction::initialise( DFConfig Config )
c@225: {
c@225:     m_dataLength = Config.frameLength;
c@225:     m_halfLength = m_dataLength/2;
c@239: 
c@225:     m_DFType = Config.DFType;
c@238:     m_stepSize = Config.stepSize;
c@225: 
c@239:     m_whiten = Config.adaptiveWhitening;
c@239:     m_whitenRelaxCoeff = Config.whiteningRelaxCoeff;
c@239:     m_whitenFloor = Config.whiteningFloor;
c@239:     if (m_whitenRelaxCoeff < 0) m_whitenRelaxCoeff = 0.9997;
c@239:     if (m_whitenFloor < 0) m_whitenFloor = 0.01;
c@239: 
c@227:     m_magHistory = new double[ m_halfLength ];
c@227:     memset(m_magHistory,0, m_halfLength*sizeof(double));
c@225: 		
c@227:     m_phaseHistory = new double[ m_halfLength ];
c@227:     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
c@225: 
c@227:     m_phaseHistoryOld = new double[ m_halfLength ];
c@227:     memset(m_phaseHistoryOld,0, m_halfLength*sizeof(double));
c@225: 
c@239:     m_magPeaks = new double[ m_halfLength ];
c@239:     memset(m_magPeaks,0, m_halfLength*sizeof(double));
c@239: 
c@289:     // See note in process(const double *) below
c@289:     int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
c@289:     m_phaseVoc = new PhaseVocoder(actualLength);
c@225: 
c@225:     m_DFWindowedFrame = new double[ m_dataLength ];
c@225:     m_magnitude = new double[ m_halfLength ];
c@225:     m_thetaAngle = new double[ m_halfLength ];
c@225: 
c@225:     m_window = new Window<double>(HanningWindow, m_dataLength);
c@225: }
c@225: 
c@225: void DetectionFunction::deInitialise()
c@225: {
c@227:     delete [] m_magHistory ;
c@227:     delete [] m_phaseHistory ;
c@227:     delete [] m_phaseHistoryOld ;
c@239:     delete [] m_magPeaks ;
c@225: 
c@225:     delete m_phaseVoc;
c@225: 
c@225:     delete [] m_DFWindowedFrame;
c@225:     delete [] m_magnitude;
c@225:     delete [] m_thetaAngle;
c@225: 
c@225:     delete m_window;
c@225: }
c@225: 
c@280: double DetectionFunction::process( const double *TDomain )
c@225: {
c@225:     m_window->cut( TDomain, m_DFWindowedFrame );
c@280: 
c@280:     // Our own FFT implementation supports power-of-two sizes only.
c@280:     // If we have to use this implementation (as opposed to the
c@280:     // version of process() below that operates on frequency domain
c@280:     // data directly), we will have to use the next smallest power of
c@280:     // two from the block size.  Results may vary accordingly!
c@280: 
c@280:     int actualLength = MathUtilities::previousPowerOfTwo(m_dataLength);
c@280: 
c@280:     if (actualLength != m_dataLength) {
c@280:         // Pre-fill mag and phase vectors with zero, as the FFT output
c@280:         // will not fill the arrays
c@280:         for (int i = actualLength/2; i < m_dataLength/2; ++i) {
c@280:             m_magnitude[i] = 0;
c@280:             m_thetaAngle[0] = 0;
c@280:         }
c@280:     }
c@280: 
c@289:     m_phaseVoc->process(m_DFWindowedFrame, m_magnitude, m_thetaAngle);
c@225: 
c@239:     if (m_whiten) whiten();
c@239: 
c@227:     return runDF();
c@227: }
c@227: 
c@280: double DetectionFunction::process( const double *magnitudes, const double *phases )
c@227: {
c@227:     for (size_t i = 0; i < m_halfLength; ++i) {
c@227:         m_magnitude[i] = magnitudes[i];
c@227:         m_thetaAngle[i] = phases[i];
c@227:     }
c@227: 
c@239:     if (m_whiten) whiten();
c@239: 
c@227:     return runDF();
c@227: }
c@227: 
c@239: void DetectionFunction::whiten()
c@239: {
c@239:     for (unsigned int i = 0; i < m_halfLength; ++i) {
c@239:         double m = m_magnitude[i];
c@239:         if (m < m_magPeaks[i]) {
c@239:             m = m + (m_magPeaks[i] - m) * m_whitenRelaxCoeff;
c@239:         }
c@239:         if (m < m_whitenFloor) m = m_whitenFloor;
c@239:         m_magPeaks[i] = m;
c@239:         m_magnitude[i] /= m;
c@239:     }
c@239: }
c@239: 
c@227: double DetectionFunction::runDF()
c@227: {
c@227:     double retVal = 0;
c@227: 
c@225:     switch( m_DFType )
c@225:     {
c@225:     case DF_HFC:
c@225: 	retVal = HFC( m_halfLength, m_magnitude);
c@225: 	break;
c@225: 	
c@238:     case DF_SPECDIFF:
c@225: 	retVal = specDiff( m_halfLength, m_magnitude);
c@225: 	break;
c@225: 	
c@225:     case DF_PHASEDEV:
c@239: 	retVal = phaseDev( m_halfLength, m_thetaAngle);
c@225: 	break;
c@225: 	
c@225:     case DF_COMPLEXSD:
c@225: 	retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
c@225: 	break;
c@237: 
c@237:     case DF_BROADBAND:
c@239:         retVal = broadband( m_halfLength, m_magnitude);
c@239:         break;
c@225:     }
c@225: 	
c@225:     return retVal;
c@225: }
c@225: 
c@225: double DetectionFunction::HFC(unsigned int length, double *src)
c@225: {
c@225:     unsigned int i;
c@225:     double val = 0;
c@225: 
c@225:     for( i = 0; i < length; i++)
c@225:     {
c@225: 	val += src[ i ] * ( i + 1);
c@225:     }
c@225:     return val;
c@225: }
c@225: 
c@225: double DetectionFunction::specDiff(unsigned int length, double *src)
c@225: {
c@225:     unsigned int i;
c@225:     double val = 0.0;
c@225:     double temp = 0.0;
c@225:     double diff = 0.0;
c@225: 
c@225:     for( i = 0; i < length; i++)
c@225:     {
c@227: 	temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
c@225: 		
c@225: 	diff= sqrt(temp);
c@225: 
c@238:         // (See note in phaseDev below.)
c@238: 
c@238:         val += diff;
c@225: 
c@227: 	m_magHistory[ i ] = src[ i ];
c@225:     }
c@225: 
c@225:     return val;
c@225: }
c@225: 
c@225: 
c@239: double DetectionFunction::phaseDev(unsigned int length, double *srcPhase)
c@225: {
c@225:     unsigned int i;
c@225:     double tmpPhase = 0;
c@225:     double tmpVal = 0;
c@225:     double val = 0;
c@225: 
c@225:     double dev = 0;
c@225: 
c@225:     for( i = 0; i < length; i++)
c@225:     {
c@227: 	tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
c@225: 	dev = MathUtilities::princarg( tmpPhase );
c@238: 
c@238:         // A previous version of this code only counted the value here
c@238:         // if the magnitude exceeded 0.1.  My impression is that
c@238:         // doesn't greatly improve the results for "loud" music (so
c@238:         // long as the peak picker is reasonably sophisticated), but
c@238:         // does significantly damage its ability to work with quieter
c@238:         // music, so I'm removing it and counting the result always.
c@238:         // Same goes for the spectral difference measure above.
c@225: 		
c@238:         tmpVal  = fabs(dev);
c@238:         val += tmpVal ;
c@225: 
c@227: 	m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
c@227: 	m_phaseHistory[ i ] = srcPhase[ i ];
c@225:     }
c@225: 	
c@225: 	
c@225:     return val;
c@225: }
c@225: 
c@225: 
c@225: double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase)
c@225: {
c@225:     unsigned int i;
c@225:     double val = 0;
c@225:     double tmpPhase = 0;
c@225:     double tmpReal = 0;
c@225:     double tmpImag = 0;
c@225:    
c@225:     double dev = 0;
c@225:     ComplexData meas = ComplexData( 0, 0 );
c@227:     ComplexData j = ComplexData( 0, 1 );
c@225: 
c@225:     for( i = 0; i < length; i++)
c@225:     {
c@227: 	tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
c@225: 	dev= MathUtilities::princarg( tmpPhase );
c@225: 		
c@227: 	meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
c@225: 
c@225: 	tmpReal = real( meas );
c@225: 	tmpImag = imag( meas );
c@225: 
c@225: 	val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
c@225: 		
c@227: 	m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
c@227: 	m_phaseHistory[ i ] = srcPhase[ i ];
c@227: 	m_magHistory[ i ] = srcMagnitude[ i ];
c@225:     }
c@225: 
c@225:     return val;
c@225: }
c@225: 
c@239: double DetectionFunction::broadband(unsigned int length, double *src)
c@237: {
c@237:     double val = 0;
c@237:     for (unsigned int i = 0; i < length; ++i) {
c@239:         double sqrmag = src[i] * src[i];
c@237:         if (m_magHistory[i] > 0.0) {
c@237:             double diff = 10.0 * log10(sqrmag / m_magHistory[i]);
c@237:             if (diff > m_dbRise) val = val + 1;
c@237:         }
c@237:         m_magHistory[i] = sqrmag;
c@237:     }
c@237:     return val;
c@237: }        
c@237: 
c@225: double* DetectionFunction::getSpectrumMagnitude()
c@225: {
c@225:     return m_magnitude;
c@225: }
c@225: