libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     QM DSP Library
   5
   6     Centre for Digital Music, Queen Mary, University of London.
   7     This file copyright 2006 Martin Gasser.
   8
   9     This program is free software; you can redistribute it and/or
  10     modify it under the terms of the GNU General Public License as
  11     published by the Free Software Foundation; either version 2 of the
  12     License, or (at your option) any later version.  See the file
  13     COPYING included with this distribution for more information.
  14 */
  15
  16 #include "ChangeDetectionFunction.h"
  17
  18 #ifndef PI
  19 #define PI (3.14159265358979232846)
  20 #endif
  21
  22
  23
  24 ChangeDetectionFunction::ChangeDetectionFunction(ChangeDFConfig config) :
  25         m_dFilterSigma(0.0), m_iFilterWidth(0)
  26 {
  27         setFilterWidth(config.smoothingWidth);
  28 }
  29
  30 ChangeDetectionFunction::~ChangeDetectionFunction()
  31 {
  32 }
  33
  34 void ChangeDetectionFunction::setFilterWidth(const int iWidth)
  35 {
  36         m_iFilterWidth = iWidth*2+1;
  37
  38         // it is assumed that the gaussian is 0 outside of +/- FWHM
  39         // => filter width = 2*FWHM = 2*2.3548*sigma
  40         m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548);
  41         m_vaGaussian.resize(m_iFilterWidth);
  42
  43         double dScale = 1.0 / (m_dFilterSigma*sqrt(2*PI));
  44
  45         for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++)
  46         {
  47                 double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) );
  48                 m_vaGaussian[x + (m_iFilterWidth-1)/2] = w;
  49         }
  50
  51 #ifdef DEBUG_CHANGE_DETECTION_FUNCTION
  52         std::cerr << "Filter sigma: " << m_dFilterSigma << std::endl;
  53         std::cerr << "Filter width: " << m_iFilterWidth << std::endl;
  54 #endif
  55 }
  56
  57
  58 ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram)
  59 {
  60         ChangeDistance retVal;
  61         retVal.resize(rTCSGram.getSize(), 0.0);
  62
  63         TCSGram smoothedTCSGram;
  64
  65         for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
  66         {
  67                 int iSkipLower = 0;
  68
  69                 int iLowerPos = iPosition - (m_iFilterWidth-1)/2;
  70                 int iUpperPos = iPosition + (m_iFilterWidth-1)/2;
  71
  72                 if (iLowerPos < 0)
  73                 {
  74                         iSkipLower = -iLowerPos;
  75                         iLowerPos = 0;
  76                 }
  77
  78                 if (iUpperPos >= rTCSGram.getSize())
  79                 {
  80                         int iMaxIndex = rTCSGram.getSize() - 1;
  81                         iUpperPos = iMaxIndex;
  82                 }
  83
  84                 TCSVector smoothedVector;
  85
  86                 // for every bin of the vector, calculate the smoothed value
  87                 for (int iPC = 0; iPC < 6; iPC++)
  88                 {
  89                         size_t j = 0;
  90                         double dSmoothedValue = 0.0;
  91                         TCSVector rCV;
  92
  93                         for (int i = iLowerPos; i <= iUpperPos; i++)
  94                         {
  95                                 rTCSGram.getTCSVector(i, rCV);
  96                                 dSmoothedValue += m_vaGaussian[iSkipLower + j++] * rCV[iPC];
  97                         }
  98
  99                         smoothedVector[iPC] = dSmoothedValue;
 100                 }
 101
 102                 smoothedTCSGram.addTCSVector(smoothedVector);
 103         }
 104
 105         for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++)
 106         {
 107                 /*
 108                         TODO: calculate a confidence measure for the current estimation
 109                         if the current estimate is not confident enough, look further into the future/the past
 110                         e.g., High frequency content, zero crossing rate, spectral flatness
 111                 */
 112
 113                 TCSVector nextTCS;
 114                 TCSVector previousTCS;
 115
 116                 int iWindow = 1;
 117
 118                 // while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0)
 119                 {
 120                         smoothedTCSGram.getTCSVector(iPosition-iWindow, previousTCS);
 121                         // std::cout << previousTCS.magnitude() << std::endl;
 122                         iWindow++;
 123                 }
 124
 125                 iWindow = 1;
 126
 127                 // while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) )
 128                 {
 129                         smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS);
 130                         iWindow++;
 131                 }
 132
 133                 double distance = 0.0;
 134                 // Euclidean distance
 135                 for (size_t j = 0; j < 6; j++)
 136                 {
 137                         distance += std::pow(nextTCS[j] - previousTCS[j], 2.0);
 138                 }
 139
 140                 retVal[iPosition] = std::pow(distance, 0.5);
 141         }
 142
 143         return retVal;
 144 }