qm-dsp  1.8
ChangeDetectionFunction.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 copyright 2006 Martin Gasser.
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 "ChangeDetectionFunction.h"
00017 
00018 #ifndef PI
00019 #define PI (3.14159265358979232846)
00020 #endif
00021 
00022 
00023 
00024 ChangeDetectionFunction::ChangeDetectionFunction(ChangeDFConfig config) :
00025         m_dFilterSigma(0.0), m_iFilterWidth(0)
00026 {
00027         setFilterWidth(config.smoothingWidth);
00028 }
00029 
00030 ChangeDetectionFunction::~ChangeDetectionFunction()
00031 {
00032 }
00033 
00034 void ChangeDetectionFunction::setFilterWidth(const int iWidth)
00035 {
00036         m_iFilterWidth = iWidth*2+1;
00037         
00038         // it is assumed that the gaussian is 0 outside of +/- FWHM
00039         // => filter width = 2*FWHM = 2*2.3548*sigma
00040         m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548);
00041         m_vaGaussian.resize(m_iFilterWidth);
00042         
00043         double dScale = 1.0 / (m_dFilterSigma*sqrt(2*PI));
00044         
00045         for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++)
00046         {
00047                 double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) );
00048                 m_vaGaussian[x + (m_iFilterWidth-1)/2] = w;
00049         }
00050         
00051 #ifdef DEBUG_CHANGE_DETECTION_FUNCTION
00052         std::cerr << "Filter sigma: " << m_dFilterSigma << std::endl;
00053         std::cerr << "Filter width: " << m_iFilterWidth << std::endl;
00054 #endif
00055 }
00056 
00057 
00058 ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram)
00059 {
00060         ChangeDistance retVal;
00061         retVal.resize(rTCSGram.getSize(), 0.0);
00062         
00063         TCSGram smoothedTCSGram;
00064 
00065         for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
00066         {
00067                 int iSkipLower = 0;
00068         
00069                 int iLowerPos = iPosition - (m_iFilterWidth-1)/2;
00070                 int iUpperPos = iPosition + (m_iFilterWidth-1)/2;
00071         
00072                 if (iLowerPos < 0)
00073                 {
00074                         iSkipLower = -iLowerPos;
00075                         iLowerPos = 0;
00076                 }
00077         
00078                 if (iUpperPos >= rTCSGram.getSize())
00079                 {
00080                         int iMaxIndex = rTCSGram.getSize() - 1;
00081                         iUpperPos = iMaxIndex;
00082                 }
00083         
00084                 TCSVector smoothedVector;
00085 
00086                 // for every bin of the vector, calculate the smoothed value
00087                 for (int iPC = 0; iPC < 6; iPC++)
00088                 {       
00089                         size_t j = 0;
00090                         double dSmoothedValue = 0.0;
00091                         TCSVector rCV;
00092                 
00093                         for (int i = iLowerPos; i <= iUpperPos; i++)
00094                         {
00095                                 rTCSGram.getTCSVector(i, rCV);
00096                                 dSmoothedValue += m_vaGaussian[iSkipLower + j++] * rCV[iPC];
00097                         }
00098 
00099                         smoothedVector[iPC] = dSmoothedValue;
00100                 }
00101                 
00102                 smoothedTCSGram.addTCSVector(smoothedVector);
00103         }
00104 
00105         for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
00106         {
00107                 /*
00108                         TODO: calculate a confidence measure for the current estimation
00109                         if the current estimate is not confident enough, look further into the future/the past
00110                         e.g., High frequency content, zero crossing rate, spectral flatness
00111                 */
00112                 
00113                 TCSVector nextTCS;
00114                 TCSVector previousTCS;
00115                 
00116                 int iWindow = 1;
00117 
00118                 // while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0)
00119                 {
00120                         smoothedTCSGram.getTCSVector(iPosition-iWindow, previousTCS);
00121                         // std::cout << previousTCS.magnitude() << std::endl;
00122                         iWindow++;
00123                 }
00124                 
00125                 iWindow = 1;
00126                 
00127                 // while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) )
00128                 {
00129                         smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS);
00130                         iWindow++;
00131                 }
00132 
00133                 double distance = 0.0;
00134                 // Euclidean distance
00135                 for (size_t j = 0; j < 6; j++)
00136                 {
00137                         distance += std::pow(nextTCS[j] - previousTCS[j], 2.0);
00138                 }
00139         
00140                 retVal[iPosition] = std::pow(distance, 0.5);
00141         }
00142 
00143         return retVal;
00144 }