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 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 }