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@225:     This file copyright 2006 Martin Gasser.
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 "ChangeDetectionFunction.h"
c@225: 
c@225: ChangeDetectionFunction::ChangeDetectionFunction(ChangeDFConfig config) :
cannam@482:     m_dFilterSigma(0.0), m_iFilterWidth(0)
c@225: {
cannam@482:     setFilterWidth(config.smoothingWidth);
c@225: }
c@225: 
c@225: ChangeDetectionFunction::~ChangeDetectionFunction()
c@225: {
c@225: }
c@225: 
c@225: void ChangeDetectionFunction::setFilterWidth(const int iWidth)
c@225: {
cannam@482:     m_iFilterWidth = iWidth*2+1;
cannam@482:         
cannam@482:     // it is assumed that the gaussian is 0 outside of +/- FWHM
cannam@482:     // => filter width = 2*FWHM = 2*2.3548*sigma
cannam@482:     m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548);
cannam@482:     m_vaGaussian.resize(m_iFilterWidth);
cannam@482:         
cannam@488:     double dScale = 1.0 / (m_dFilterSigma*sqrt(2*M_PI));
cannam@482:         
cannam@482:     for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++) {
cannam@482:         double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) );
cannam@482:         m_vaGaussian[x + (m_iFilterWidth-1)/2] = w;
cannam@482:     }
cannam@482:         
c@225: #ifdef DEBUG_CHANGE_DETECTION_FUNCTION
cannam@482:     std::cerr << "Filter sigma: " << m_dFilterSigma << std::endl;
cannam@482:     std::cerr << "Filter width: " << m_iFilterWidth << std::endl;
c@225: #endif
c@225: }
c@225: 
c@225: 
c@225: ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram)
c@225: {
cannam@482:     ChangeDistance retVal;
cannam@482:     retVal.resize(rTCSGram.getSize(), 0.0);
cannam@482:         
cannam@482:     TCSGram smoothedTCSGram;
c@225: 
cannam@482:     for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++) {
cannam@482:         
cannam@482:         int iSkipLower = 0;
cannam@482:         
cannam@482:         int iLowerPos = iPosition - (m_iFilterWidth-1)/2;
cannam@482:         int iUpperPos = iPosition + (m_iFilterWidth-1)/2;
cannam@482:         
cannam@482:         if (iLowerPos < 0) {
cannam@482:             iSkipLower = -iLowerPos;
cannam@482:             iLowerPos = 0;
cannam@482:         }
cannam@482:         
cannam@482:         if (iUpperPos >= rTCSGram.getSize()) {
cannam@482:             int iMaxIndex = rTCSGram.getSize() - 1;
cannam@482:             iUpperPos = iMaxIndex;
cannam@482:         }
cannam@482:         
cannam@482:         TCSVector smoothedVector;
c@225: 
cannam@482:         // for every bin of the vector, calculate the smoothed value
cannam@482:         for (int iPC = 0; iPC < 6; iPC++) {       
c@225: 
cannam@482:             size_t j = 0;
cannam@482:             double dSmoothedValue = 0.0;
cannam@482:             TCSVector rCV;
cannam@482:                 
cannam@482:             for (int i = iLowerPos; i <= iUpperPos; i++) {
cannam@482:                 rTCSGram.getTCSVector(i, rCV);
cannam@482:                 dSmoothedValue += m_vaGaussian[iSkipLower + j++] * rCV[iPC];
cannam@482:             }
c@225: 
cannam@482:             smoothedVector[iPC] = dSmoothedValue;
cannam@482:         }
cannam@482:                 
cannam@482:         smoothedTCSGram.addTCSVector(smoothedVector);
cannam@482:     }
c@225: 
cannam@482:     for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++) {
cannam@482:         
cannam@482:         /*
cannam@482:           TODO: calculate a confidence measure for the current estimation
cannam@482:           if the current estimate is not confident enough, look further into the future/the past
cannam@482:           e.g., High frequency content, zero crossing rate, spectral flatness
cannam@482:         */
cannam@482:                 
cannam@482:         TCSVector nextTCS;
cannam@482:         TCSVector previousTCS;
cannam@482:                 
cannam@482:         int iWindow = 1;
c@225: 
cannam@482:         // while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0)
cannam@482:         {
cannam@482:             smoothedTCSGram.getTCSVector(iPosition-iWindow, previousTCS);
cannam@482:             // std::cout << previousTCS.magnitude() << std::endl;
cannam@482:             iWindow++;
cannam@482:         }
cannam@482:                 
cannam@482:         iWindow = 1;
cannam@482:                 
cannam@482:         // while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) )
cannam@482:         {
cannam@482:             smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS);
cannam@482:             iWindow++;
cannam@482:         }
c@225: 
cannam@482:         double distance = 0.0;
cannam@482:         // Euclidean distance
cannam@482:         for (size_t j = 0; j < 6; j++) {
cannam@482:             distance += std::pow(nextTCS[j] - previousTCS[j], 2.0);
cannam@482:         }
cannam@482:         
cannam@482:         retVal[iPosition] = std::pow(distance, 0.5);
cannam@482:     }
cannam@482: 
cannam@482:     return retVal;
c@225: }