Mercurial > hg > pyin

annotate YinUtil.cpp @ 1:3dcef83df62a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
notes work again, now based on the hard PYIN estimate
author matthiasm
date Wed, 27 Nov 2013 15:31:47 +0000
parents 99bac62ee2da
children 5945b8905d1f
rev   line source
matthiasm@0 1 #include "YinUtil.h"
matthiasm@0 2
matthiasm@0 3 #include <vector>
matthiasm@0 4
matthiasm@0 5 #include <cstdio>
matthiasm@0 6 #include <cmath>
matthiasm@0 7 #include <algorithm>
matthiasm@0 8
matthiasm@0 9 #include <boost/math/distributions.hpp>
matthiasm@0 10
matthiasm@0 11 void
matthiasm@0 12 YinUtil::fastDifference(const double *in, double *yinBuffer, const size_t yinBufferSize)
matthiasm@0 13 {
matthiasm@0 14
matthiasm@0 15 // DECLARE AND INITIALISE
matthiasm@0 16 // initialisation of most of the arrays here was done in a separate function,
matthiasm@0 17 // with all the arrays as members of the class... moved them back here.
matthiasm@0 18
matthiasm@0 19 size_t frameSize = 2 * yinBufferSize;
matthiasm@0 20
matthiasm@0 21 for (size_t j = 0; j < yinBufferSize; ++j)
matthiasm@0 22 {
matthiasm@0 23 yinBuffer[j] = 0.;
matthiasm@0 24 }
matthiasm@0 25
matthiasm@0 26 double *audioTransformedReal = new double[frameSize];
matthiasm@0 27 double *audioTransformedImag = new double[frameSize];
matthiasm@0 28 double *nullImag = new double[frameSize];
matthiasm@0 29 double *kernel = new double[frameSize];
matthiasm@0 30 double *kernelTransformedReal = new double[frameSize];
matthiasm@0 31 double *kernelTransformedImag = new double[frameSize];
matthiasm@0 32 double *yinStyleACFReal = new double[frameSize];
matthiasm@0 33 double *yinStyleACFImag = new double[frameSize];
matthiasm@0 34 double *powerTerms = new double[yinBufferSize];
matthiasm@0 35
matthiasm@0 36 for (size_t j = 0; j < yinBufferSize; ++j)
matthiasm@0 37 {
matthiasm@0 38 powerTerms[j] = 0.;
matthiasm@0 39 }
matthiasm@0 40
matthiasm@0 41 for (size_t j = 0; j < frameSize; ++j)
matthiasm@0 42 {
matthiasm@0 43 nullImag[j] = 0.;
matthiasm@0 44 audioTransformedReal[j] = 0.;
matthiasm@0 45 audioTransformedImag[j] = 0.;
matthiasm@0 46 kernel[j] = 0.;
matthiasm@0 47 kernelTransformedReal[j] = 0.;
matthiasm@0 48 kernelTransformedImag[j] = 0.;
matthiasm@0 49 yinStyleACFReal[j] = 0.;
matthiasm@0 50 yinStyleACFImag[j] = 0.;
matthiasm@0 51 }
matthiasm@0 52
matthiasm@0 53 // POWER TERM CALCULATION
matthiasm@0 54 // ... for the power terms in equation (7) in the Yin paper
matthiasm@0 55 powerTerms[0] = 0.0;
matthiasm@0 56 for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0 57 powerTerms[0] += in[j] * in[j];
matthiasm@0 58 }
matthiasm@0 59
matthiasm@0 60 // now iteratively calculate all others (saves a few multiplications)
matthiasm@0 61 for (size_t tau = 1; tau < yinBufferSize; ++tau) {
matthiasm@0 62 powerTerms[tau] = powerTerms[tau-1] - in[tau-1] * in[tau-1] + in[tau+yinBufferSize] * in[tau+yinBufferSize];
matthiasm@0 63 }
matthiasm@0 64
matthiasm@0 65 // YIN-STYLE AUTOCORRELATION via FFT
matthiasm@0 66 // 1. data
matthiasm@0 67 Vamp::FFT::forward(frameSize, in, nullImag, audioTransformedReal, audioTransformedImag);
matthiasm@0 68
matthiasm@0 69 // 2. half of the data, disguised as a convolution kernel
matthiasm@0 70 for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0 71 kernel[j] = in[yinBufferSize-1-j];
matthiasm@0 72 kernel[j+yinBufferSize] = 0;
matthiasm@0 73 }
matthiasm@0 74 Vamp::FFT::forward(frameSize, kernel, nullImag, kernelTransformedReal, kernelTransformedImag);
matthiasm@0 75
matthiasm@0 76 // 3. convolution via complex multiplication -- written into
matthiasm@0 77 for (size_t j = 0; j < frameSize; ++j) {
matthiasm@0 78 yinStyleACFReal[j] = audioTransformedReal[j]*kernelTransformedReal[j] - audioTransformedImag[j]*kernelTransformedImag[j]; // real
matthiasm@0 79 yinStyleACFImag[j] = audioTransformedReal[j]*kernelTransformedImag[j] + audioTransformedImag[j]*kernelTransformedReal[j]; // imaginary
matthiasm@0 80 }
matthiasm@0 81 Vamp::FFT::inverse(frameSize, yinStyleACFReal, yinStyleACFImag, audioTransformedReal, audioTransformedImag);
matthiasm@0 82
matthiasm@0 83 // CALCULATION OF difference function
matthiasm@0 84 // ... according to (7) in the Yin paper.
matthiasm@0 85 for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0 86 // taking only the real part
matthiasm@0 87 yinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
matthiasm@0 88 }
matthiasm@0 89 delete [] audioTransformedReal;
matthiasm@0 90 delete [] audioTransformedImag;
matthiasm@0 91 delete [] nullImag;
matthiasm@0 92 delete [] kernel;
matthiasm@0 93 delete [] kernelTransformedReal;
matthiasm@0 94 delete [] kernelTransformedImag;
matthiasm@0 95 delete [] yinStyleACFReal;
matthiasm@0 96 delete [] yinStyleACFImag;
matthiasm@0 97 delete [] powerTerms;
matthiasm@0 98 }
matthiasm@0 99
matthiasm@0 100 void
matthiasm@0 101 YinUtil::cumulativeDifference(double *yinBuffer, const size_t yinBufferSize)
matthiasm@0 102 {
matthiasm@0 103 size_t tau;
matthiasm@0 104
matthiasm@0 105 yinBuffer[0] = 1;
matthiasm@0 106
matthiasm@0 107 double runningSum = 0;
matthiasm@0 108
matthiasm@0 109 for (tau = 1; tau < yinBufferSize; ++tau) {
matthiasm@0 110 runningSum += yinBuffer[tau];
matthiasm@0 111 if (runningSum == 0)
matthiasm@0 112 {
matthiasm@0 113 yinBuffer[tau] = 1;
matthiasm@0 114 } else {
matthiasm@0 115 yinBuffer[tau] *= tau / runningSum;
matthiasm@0 116 }
matthiasm@0 117 }
matthiasm@0 118 }
matthiasm@0 119
matthiasm@0 120 int
matthiasm@0 121 YinUtil::absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh)
matthiasm@0 122 {
matthiasm@0 123 size_t tau;
matthiasm@0 124 size_t minTau = 0;
matthiasm@0 125 double minVal = 1000.;
matthiasm@0 126
matthiasm@0 127 // using Joren Six's "loop construct" from TarsosDSP
matthiasm@0 128 tau = 2;
matthiasm@0 129 while (tau < yinBufferSize)
matthiasm@0 130 {
matthiasm@0 131 if (yinBuffer[tau] < thresh)
matthiasm@0 132 {
matthiasm@0 133 while (tau+1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
matthiasm@0 134 {
matthiasm@0 135 ++tau;
matthiasm@0 136 }
matthiasm@0 137 return tau;
matthiasm@0 138 } else {
matthiasm@0 139 if (yinBuffer[tau] < minVal)
matthiasm@0 140 {
matthiasm@0 141 minVal = yinBuffer[tau];
matthiasm@0 142 minTau = tau;
matthiasm@0 143 }
matthiasm@0 144 }
matthiasm@0 145 ++tau;
matthiasm@0 146 }
matthiasm@0 147 if (minTau > 0)
matthiasm@0 148 {
matthiasm@0 149 return -minTau;
matthiasm@0 150 }
matthiasm@0 151 return 0;
matthiasm@0 152 }
matthiasm@0 153
matthiasm@0 154
matthiasm@0 155 std::vector<double>
matthiasm@0 156 YinUtil::yinProb(const double *yinBuffer, const size_t prior, const size_t yinBufferSize)
matthiasm@0 157 {
matthiasm@0 158 double minWeight = 0.01;
matthiasm@0 159 size_t tau;
matthiasm@0 160 std::vector<float> thresholds;
matthiasm@0 161 std::vector<float> distribution;
matthiasm@0 162 std::vector<double> peakProb = std::vector<double>(yinBufferSize);
matthiasm@1 163 // TODO: make the distributions below part of a class, so they don't have to
matthiasm@1 164 // be allocated every time.
matthiasm@0 165 float uniformDist[100] = {0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000,0.0100000};
matthiasm@0 166 float betaDist1[100] = {0.028911,0.048656,0.061306,0.068539,0.071703,0.071877,0.069915,0.066489,0.062117,0.057199,0.052034,0.046844,0.041786,0.036971,0.032470,0.028323,0.024549,0.021153,0.018124,0.015446,0.013096,0.011048,0.009275,0.007750,0.006445,0.005336,0.004397,0.003606,0.002945,0.002394,0.001937,0.001560,0.001250,0.000998,0.000792,0.000626,0.000492,0.000385,0.000300,0.000232,0.000179,0.000137,0.000104,0.000079,0.000060,0.000045,0.000033,0.000024,0.000018,0.000013,0.000009,0.000007,0.000005,0.000003,0.000002,0.000002,0.000001,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
matthiasm@0 167 float betaDist2[100] = {0.012614,0.022715,0.030646,0.036712,0.041184,0.044301,0.046277,0.047298,0.047528,0.047110,0.046171,0.044817,0.043144,0.041231,0.039147,0.036950,0.034690,0.032406,0.030133,0.027898,0.025722,0.023624,0.021614,0.019704,0.017900,0.016205,0.014621,0.013148,0.011785,0.010530,0.009377,0.008324,0.007366,0.006497,0.005712,0.005005,0.004372,0.003806,0.003302,0.002855,0.002460,0.002112,0.001806,0.001539,0.001307,0.001105,0.000931,0.000781,0.000652,0.000542,0.000449,0.000370,0.000303,0.000247,0.000201,0.000162,0.000130,0.000104,0.000082,0.000065,0.000051,0.000039,0.000030,0.000023,0.000018,0.000013,0.000010,0.000007,0.000005,0.000004,0.000003,0.000002,0.000001,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
matthiasm@0 168 float betaDist3[100] = {0.006715,0.012509,0.017463,0.021655,0.025155,0.028031,0.030344,0.032151,0.033506,0.034458,0.035052,0.035331,0.035332,0.035092,0.034643,0.034015,0.033234,0.032327,0.031314,0.030217,0.029054,0.027841,0.026592,0.025322,0.024042,0.022761,0.021489,0.020234,0.019002,0.017799,0.016630,0.015499,0.014409,0.013362,0.012361,0.011407,0.010500,0.009641,0.008830,0.008067,0.007351,0.006681,0.006056,0.005475,0.004936,0.004437,0.003978,0.003555,0.003168,0.002814,0.002492,0.002199,0.001934,0.001695,0.001481,0.001288,0.001116,0.000963,0.000828,0.000708,0.000603,0.000511,0.000431,0.000361,0.000301,0.000250,0.000206,0.000168,0.000137,0.000110,0.000088,0.000070,0.000055,0.000043,0.000033,0.000025,0.000019,0.000014,0.000010,0.000007,0.000005,0.000004,0.000002,0.000002,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
matthiasm@0 169 float betaDist4[100] = {0.003996,0.007596,0.010824,0.013703,0.016255,0.018501,0.020460,0.022153,0.023597,0.024809,0.025807,0.026607,0.027223,0.027671,0.027963,0.028114,0.028135,0.028038,0.027834,0.027535,0.027149,0.026687,0.026157,0.025567,0.024926,0.024240,0.023517,0.022763,0.021983,0.021184,0.020371,0.019548,0.018719,0.017890,0.017062,0.016241,0.015428,0.014627,0.013839,0.013068,0.012315,0.011582,0.010870,0.010181,0.009515,0.008874,0.008258,0.007668,0.007103,0.006565,0.006053,0.005567,0.005107,0.004673,0.004264,0.003880,0.003521,0.003185,0.002872,0.002581,0.002312,0.002064,0.001835,0.001626,0.001434,0.001260,0.001102,0.000959,0.000830,0.000715,0.000612,0.000521,0.000440,0.000369,0.000308,0.000254,0.000208,0.000169,0.000136,0.000108,0.000084,0.000065,0.000050,0.000037,0.000027,0.000019,0.000014,0.000009,0.000006,0.000004,0.000002,0.000001,0.000001,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000,0.000000};
matthiasm@0 170 float single10[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
matthiasm@0 171 float single15[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
matthiasm@0 172 float single20[100] = {0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,1.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
matthiasm@0 173
matthiasm@0 174 size_t nThreshold = 100;
matthiasm@0 175 int nThresholdInt = nThreshold;
matthiasm@0 176
matthiasm@0 177 for (int i = 0; i < nThresholdInt; ++i)
matthiasm@0 178 {
matthiasm@0 179 switch (prior) {
matthiasm@0 180 case 0:
matthiasm@0 181 distribution.push_back(uniformDist[i]);
matthiasm@0 182 break;
matthiasm@0 183 case 1:
matthiasm@0 184 distribution.push_back(betaDist1[i]);
matthiasm@0 185 break;
matthiasm@0 186 case 2:
matthiasm@0 187 distribution.push_back(betaDist2[i]);
matthiasm@0 188 break;
matthiasm@0 189 case 3:
matthiasm@0 190 distribution.push_back(betaDist3[i]);
matthiasm@0 191 break;
matthiasm@0 192 case 4:
matthiasm@0 193 distribution.push_back(betaDist4[i]);
matthiasm@0 194 break;
matthiasm@0 195 case 5:
matthiasm@0 196 distribution.push_back(single10[i]);
matthiasm@0 197 break;
matthiasm@0 198 case 6:
matthiasm@0 199 distribution.push_back(single15[i]);
matthiasm@0 200 break;
matthiasm@0 201 case 7:
matthiasm@0 202 distribution.push_back(single20[i]);
matthiasm@0 203 break;
matthiasm@0 204 default:
matthiasm@0 205 distribution.push_back(uniformDist[i]);
matthiasm@0 206 }
matthiasm@0 207 thresholds.push_back(0.01 + i*0.01);
matthiasm@0 208 }
matthiasm@0 209
matthiasm@0 210 // double minYin = 2936;
matthiasm@0 211 // for (size_t i = 2; i < yinBufferSize; ++i)
matthiasm@0 212 // {
matthiasm@0 213 // if (yinBuffer[i] < minYin)
matthiasm@0 214 // {
matthiasm@0 215 // minYin = yinBuffer[i];
matthiasm@0 216 // }
matthiasm@0 217 // }
matthiasm@0 218 // if (minYin < 0.01) std::cerr << "min Yin buffer element: " << minYin << std::endl;
matthiasm@0 219
matthiasm@0 220
matthiasm@0 221 int currThreshInd = nThreshold-1;
matthiasm@0 222 tau = 2;
matthiasm@0 223
matthiasm@0 224 // double factor = 1.0 / (0.25 * (nThresholdInt+1) * (nThresholdInt + 1)); // factor to scale down triangular weight
matthiasm@0 225 size_t minInd = 0;
matthiasm@0 226 float minVal = 42.f;
matthiasm@0 227 while (currThreshInd != -1 && tau < yinBufferSize)
matthiasm@0 228 {
matthiasm@0 229 if (yinBuffer[tau] < thresholds[currThreshInd])
matthiasm@0 230 {
matthiasm@0 231 while (tau + 1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
matthiasm@0 232 {
matthiasm@0 233 tau++;
matthiasm@0 234 }
matthiasm@0 235 // tau is now local minimum
matthiasm@0 236 // std::cerr << tau << " " << currThreshInd << " "<< thresholds[currThreshInd] << " " << distribution[currThreshInd] << std::endl;
matthiasm@0 237 if (yinBuffer[tau] < minVal && tau > 2){
matthiasm@0 238 minVal = yinBuffer[tau];
matthiasm@0 239 minInd = tau;
matthiasm@0 240 }
matthiasm@0 241 peakProb[tau] += distribution[currThreshInd];
matthiasm@0 242 currThreshInd--;
matthiasm@0 243 } else {
matthiasm@0 244 tau++;
matthiasm@0 245 }
matthiasm@0 246 }
matthiasm@0 247 double nonPeakProb = 1;
matthiasm@0 248 for (size_t i = 0; i < yinBufferSize; ++i)
matthiasm@0 249 {
matthiasm@0 250 nonPeakProb -= peakProb[i];
matthiasm@0 251 }
matthiasm@0 252 // std::cerr << nonPeakProb << std::endl;
matthiasm@0 253 if (minInd > 0)
matthiasm@0 254 {
matthiasm@0 255 // std::cerr << "min set " << minVal << " " << minInd << " " << nonPeakProb << std::endl;
matthiasm@0 256 peakProb[minInd] += nonPeakProb * minWeight;
matthiasm@0 257 }
matthiasm@0 258
matthiasm@0 259 return peakProb;
matthiasm@0 260 }
matthiasm@0 261
matthiasm@0 262 double
matthiasm@0 263 YinUtil::parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize)
matthiasm@0 264 {
matthiasm@0 265 // this is taken almost literally from Joren Six's Java implementation
matthiasm@0 266 if (tau == yinBufferSize) // not valid anyway.
matthiasm@0 267 {
matthiasm@0 268 return static_cast<double>(tau);
matthiasm@0 269 }
matthiasm@0 270
matthiasm@0 271 double betterTau = 0.0;
matthiasm@0 272 size_t x0;
matthiasm@0 273 size_t x2;
matthiasm@0 274
matthiasm@0 275 if (tau < 1)
matthiasm@0 276 {
matthiasm@0 277 x0 = tau;
matthiasm@0 278 } else {
matthiasm@0 279 x0 = tau - 1;
matthiasm@0 280 }
matthiasm@0 281
matthiasm@0 282 if (tau + 1 < yinBufferSize)
matthiasm@0 283 {
matthiasm@0 284 x2 = tau + 1;
matthiasm@0 285 } else {
matthiasm@0 286 x2 = tau;
matthiasm@0 287 }
matthiasm@0 288
matthiasm@0 289 if (x0 == tau)
matthiasm@0 290 {
matthiasm@0 291 if (yinBuffer[tau] <= yinBuffer[x2])
matthiasm@0 292 {
matthiasm@0 293 betterTau = tau;
matthiasm@0 294 } else {
matthiasm@0 295 betterTau = x2;
matthiasm@0 296 }
matthiasm@0 297 }
matthiasm@0 298 else if (x2 == tau)
matthiasm@0 299 {
matthiasm@0 300 if (yinBuffer[tau] <= yinBuffer[x0])
matthiasm@0 301 {
matthiasm@0 302 betterTau = tau;
matthiasm@0 303 }
matthiasm@0 304 else
matthiasm@0 305 {
matthiasm@0 306 betterTau = x0;
matthiasm@0 307 }
matthiasm@0 308 }
matthiasm@0 309 else
matthiasm@0 310 {
matthiasm@0 311 float s0, s1, s2;
matthiasm@0 312 s0 = yinBuffer[x0];
matthiasm@0 313 s1 = yinBuffer[tau];
matthiasm@0 314 s2 = yinBuffer[x2];
matthiasm@0 315 // fixed AUBIO implementation, thanks to Karl Helgason:
matthiasm@0 316 // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1
matthiasm@0 317 betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
matthiasm@0 318
matthiasm@0 319 // std::cerr << tau << " --> " << betterTau << std::endl;
matthiasm@0 320
matthiasm@0 321 }
matthiasm@0 322 return betterTau;
matthiasm@0 323 }
matthiasm@0 324
matthiasm@0 325 double
matthiasm@0 326 YinUtil::sumSquare(const double *in, const size_t start, const size_t end)
matthiasm@0 327 {
matthiasm@0 328 double out = 0;
matthiasm@0 329 for (size_t i = start; i < end; ++i)
matthiasm@0 330 {
matthiasm@0 331 out += in[i] * in[i];
matthiasm@0 332 }
matthiasm@0 333 return out;
matthiasm@0 334 }