Mercurial > hg > pyin

comparison YinUtil.cpp @ 0:99bac62ee2da

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added PYIN sources, should be compileable
author matthiasm
date Wed, 27 Nov 2013 11:59:49 +0000
parents
children 3dcef83df62a
comparison
equal deleted inserted replaced
-1:000000000000 0:99bac62ee2da
1 #include "YinUtil.h"
2
3 #include <vector>
4
5 #include <cstdio>
6 #include <cmath>
7 #include <algorithm>
8
9 #include <boost/math/distributions.hpp>
10
11 void
12 YinUtil::fastDifference(const double *in, double *yinBuffer, const size_t yinBufferSize)
13 {
14
15 // DECLARE AND INITIALISE
16 // initialisation of most of the arrays here was done in a separate function,
17 // with all the arrays as members of the class... moved them back here.
18
19 size_t frameSize = 2 * yinBufferSize;
20
21 for (size_t j = 0; j < yinBufferSize; ++j)
22 {
23 yinBuffer[j] = 0.;
24 }
25
26 double *audioTransformedReal = new double[frameSize];
27 double *audioTransformedImag = new double[frameSize];
28 double *nullImag = new double[frameSize];
29 double *kernel = new double[frameSize];
30 double *kernelTransformedReal = new double[frameSize];
31 double *kernelTransformedImag = new double[frameSize];
32 double *yinStyleACFReal = new double[frameSize];
33 double *yinStyleACFImag = new double[frameSize];
34 double *powerTerms = new double[yinBufferSize];
35
36 for (size_t j = 0; j < yinBufferSize; ++j)
37 {
38 powerTerms[j] = 0.;
39 }
40
41 for (size_t j = 0; j < frameSize; ++j)
42 {
43 nullImag[j] = 0.;
44 audioTransformedReal[j] = 0.;
45 audioTransformedImag[j] = 0.;
46 kernel[j] = 0.;
47 kernelTransformedReal[j] = 0.;
48 kernelTransformedImag[j] = 0.;
49 yinStyleACFReal[j] = 0.;
50 yinStyleACFImag[j] = 0.;
51 }
52
53 // POWER TERM CALCULATION
54 // ... for the power terms in equation (7) in the Yin paper
55 powerTerms[0] = 0.0;
56 for (size_t j = 0; j < yinBufferSize; ++j) {
57 powerTerms[0] += in[j] * in[j];
58 }
59
60 // now iteratively calculate all others (saves a few multiplications)
61 for (size_t tau = 1; tau < yinBufferSize; ++tau) {
62 powerTerms[tau] = powerTerms[tau-1] - in[tau-1] * in[tau-1] + in[tau+yinBufferSize] * in[tau+yinBufferSize];
63 }
64
65 // YIN-STYLE AUTOCORRELATION via FFT
66 // 1. data
67 Vamp::FFT::forward(frameSize, in, nullImag, audioTransformedReal, audioTransformedImag);
68
69 // 2. half of the data, disguised as a convolution kernel
70 for (size_t j = 0; j < yinBufferSize; ++j) {
71 kernel[j] = in[yinBufferSize-1-j];
72 kernel[j+yinBufferSize] = 0;
73 }
74 Vamp::FFT::forward(frameSize, kernel, nullImag, kernelTransformedReal, kernelTransformedImag);
75
76 // 3. convolution via complex multiplication -- written into
77 for (size_t j = 0; j < frameSize; ++j) {
78 yinStyleACFReal[j] = audioTransformedReal[j]*kernelTransformedReal[j] - audioTransformedImag[j]*kernelTransformedImag[j]; // real
79 yinStyleACFImag[j] = audioTransformedReal[j]*kernelTransformedImag[j] + audioTransformedImag[j]*kernelTransformedReal[j]; // imaginary
80 }
81 Vamp::FFT::inverse(frameSize, yinStyleACFReal, yinStyleACFImag, audioTransformedReal, audioTransformedImag);
82
83 // CALCULATION OF difference function
84 // ... according to (7) in the Yin paper.
85 for (size_t j = 0; j < yinBufferSize; ++j) {
86 // taking only the real part
87 yinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
88 }
89 delete [] audioTransformedReal;
90 delete [] audioTransformedImag;
91 delete [] nullImag;
92 delete [] kernel;
93 delete [] kernelTransformedReal;
94 delete [] kernelTransformedImag;
95 delete [] yinStyleACFReal;
96 delete [] yinStyleACFImag;
97 delete [] powerTerms;
98 }
99
100 void
101 YinUtil::cumulativeDifference(double *yinBuffer, const size_t yinBufferSize)
102 {
103 size_t tau;
104
105 yinBuffer[0] = 1;
106
107 double runningSum = 0;
108
109 for (tau = 1; tau < yinBufferSize; ++tau) {
110 runningSum += yinBuffer[tau];
111 if (runningSum == 0)
112 {
113 yinBuffer[tau] = 1;
114 } else {
115 yinBuffer[tau] *= tau / runningSum;
116 }
117 }
118 }
119
120 int
121 YinUtil::absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh)
122 {
123 size_t tau;
124 size_t minTau = 0;
125 double minVal = 1000.;
126
127 // using Joren Six's "loop construct" from TarsosDSP
128 tau = 2;
129 while (tau < yinBufferSize)
130 {
131 if (yinBuffer[tau] < thresh)
132 {
133 while (tau+1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
134 {
135 ++tau;
136 }
137 return tau;
138 } else {
139 if (yinBuffer[tau] < minVal)
140 {
141 minVal = yinBuffer[tau];
142 minTau = tau;
143 }
144 }
145 ++tau;
146 }
147 if (minTau > 0)
148 {
149 return -minTau;
150 }
151 return 0;
152 }
153
154
155 std::vector<double>
156 YinUtil::yinProb(const double *yinBuffer, const size_t prior, const size_t yinBufferSize)
157 {
158 double minWeight = 0.01;
159 size_t tau;
160 std::vector<float> thresholds;
161 std::vector<float> distribution;
162 std::vector<double> peakProb = std::vector<double>(yinBufferSize);
163 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};
164 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};
165 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};
166 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};
167 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};
168 // float betaDist4[100] = {0.0000119,0.0000470,0.0001048,0.0001843,0.0002850,0.0004061,0.0005469,0.0007066,0.0008846,0.0010801,0.0012924,0.0015208,0.0017645,0.0020229,0.0022952,0.0025807,0.0028787,0.0031885,0.0035093,0.0038404,0.0041811,0.0045307,0.0048884,0.0052536,0.0056256,0.0060035,0.0063867,0.0067744,0.0071660,0.0075608,0.0079579,0.0083567,0.0087564,0.0091564,0.0095560,0.0099543,0.0103507,0.0107444,0.0111348,0.0115212,0.0119027,0.0122787,0.0126484,0.0130112,0.0133663,0.0137131,0.0140506,0.0143784,0.0146956,0.0150015,0.0152954,0.0155766,0.0158443,0.0160979,0.0163366,0.0165597,0.0167665,0.0169563,0.0171282,0.0172817,0.0174160,0.0175304,0.0176241,0.0176965,0.0177468,0.0177743,0.0177782,0.0177579,0.0177127,0.0176418,0.0175444,0.0174200,0.0172677,0.0170868,0.0168767,0.0166365,0.0163657,0.0160634,0.0157289,0.0153615,0.0149606,0.0145253,0.0140550,0.0135489,0.0130063,0.0124265,0.0118088,0.0111525,0.0104568,0.0097210,0.0089444,0.0081263,0.0072659,0.0063626,0.0054155,0.0044241,0.0033876,0.0023052,0.0011762,0.0000000};
169 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};
170 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};
171 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};
172
173 // float betaDist4[100] = {0.00001,0.00005,0.00010,0.00018,0.00029,0.00041,0.00055,0.00071,0.00088,0.00108,0.00129,0.00152,0.00176,0.00202,0.00230,0.00258,0.00288,0.00319,0.00351,0.00384,0.00418,0.00453,0.00489,0.00525,0.00563,0.00600,0.00639,0.00677,0.00717,0.00756,0.00796,0.00836,0.00876,0.00916,0.00956,0.00995,0.01035,0.01074,0.01113,0.01152,0.01190,0.01228,0.01265,0.01301,0.01337,0.01371,0.01405,0.01438,0.01470,0.01500,0.01530,0.01558,0.01584,0.01610,0.01634,0.01656,0.01677,0.01696,0.01713,0.01728,0.01742,0.01753,0.01762,0.01770,0.01775,0.01777,0.01778,0.01776,0.01771,0.01764,0.01754,0.01742,0.01727,0.01709,0.01688,0.01664,0.01637,0.01606,0.01573,0.01536,0.01496,0.01453,0.01405,0.01355,0.01301,0.01243,0.01181,0.01115,0.01046,0.00972,0.00894,0.00813,0.00727,0.00636,0.00542,0.00442,0.00339,0.00231,0.00118,0.00000};
174 // float uniformDist[40] = {0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500,0.02500};
175 // float betaDist1[40] = {0.00632,0.02052,0.03731,0.05327,0.06644,0.07587,0.08133,0.08305,0.08153,0.07743,0.07144,0.06421,0.05633,0.04831,0.04052,0.03326,0.02671,0.02098,0.01611,0.01209,0.00884,0.00629,0.00436,0.00292,0.00189,0.00118,0.00070,0.00040,0.00021,0.00011,0.00005,0.00002,0.00001,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000,0.00000};
176 // float betaDist2[40] = {0.00148,0.00535,0.01081,0.01722,0.02404,0.03083,0.03724,0.04301,0.04794,0.05191,0.05485,0.05672,0.05756,0.05740,0.05633,0.05443,0.05183,0.04864,0.04499,0.04102,0.03683,0.03256,0.02832,0.02419,0.02028,0.01664,0.01334,0.01042,0.00789,0.00577,0.00404,0.00269,0.00168,0.00096,0.00049,0.00021,0.00007,0.00001,0.00000,0.00000};
177 // float betaDist3[40] = {0.00045,0.00169,0.00361,0.00608,0.00897,0.01219,0.01563,0.01920,0.02280,0.02637,0.02981,0.03308,0.03609,0.03882,0.04120,0.04320,0.04479,0.04594,0.04664,0.04688,0.04664,0.04594,0.04479,0.04320,0.04120,0.03882,0.03609,0.03308,0.02981,0.02637,0.02280,0.01920,0.01563,0.01219,0.00897,0.00608,0.00361,0.00169,0.00045,0.00000};
178 // float betaDist4[40] = {0.00018,0.00071,0.00156,0.00270,0.00410,0.00574,0.00758,0.00961,0.01178,0.01407,0.01646,0.01891,0.02140,0.02390,0.02638,0.02882,0.03118,0.03343,0.03556,0.03752,0.03930,0.04086,0.04218,0.04323,0.04397,0.04439,0.04445,0.04413,0.04339,0.04221,0.04057,0.03842,0.03576,0.03253,0.02873,0.02432,0.01926,0.01355,0.00713,0.00000};
179 // float single10[40] = {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};
180 // float single15[40] = {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};
181 // float single20[40] = {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};
182
183 size_t nThreshold = 100;
184 int nThresholdInt = nThreshold;
185
186 for (int i = 0; i < nThresholdInt; ++i)
187 {
188 switch (prior) {
189 case 0:
190 distribution.push_back(uniformDist[i]);
191 break;
192 case 1:
193 distribution.push_back(betaDist1[i]);
194 break;
195 case 2:
196 distribution.push_back(betaDist2[i]);
197 break;
198 case 3:
199 distribution.push_back(betaDist3[i]);
200 break;
201 case 4:
202 distribution.push_back(betaDist4[i]);
203 break;
204 case 5:
205 distribution.push_back(single10[i]);
206 break;
207 case 6:
208 distribution.push_back(single15[i]);
209 break;
210 case 7:
211 distribution.push_back(single20[i]);
212 break;
213 default:
214 distribution.push_back(uniformDist[i]);
215 }
216 thresholds.push_back(0.01 + i*0.01);
217 }
218
219 // double minYin = 2936;
220 // for (size_t i = 2; i < yinBufferSize; ++i)
221 // {
222 // if (yinBuffer[i] < minYin)
223 // {
224 // minYin = yinBuffer[i];
225 // }
226 // }
227 // if (minYin < 0.01) std::cerr << "min Yin buffer element: " << minYin << std::endl;
228
229
230 int currThreshInd = nThreshold-1;
231 tau = 2;
232
233 // double factor = 1.0 / (0.25 * (nThresholdInt+1) * (nThresholdInt + 1)); // factor to scale down triangular weight
234 size_t minInd = 0;
235 float minVal = 42.f;
236 while (currThreshInd != -1 && tau < yinBufferSize)
237 {
238 if (yinBuffer[tau] < thresholds[currThreshInd])
239 {
240 while (tau + 1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
241 {
242 tau++;
243 }
244 // tau is now local minimum
245 // std::cerr << tau << " " << currThreshInd << " "<< thresholds[currThreshInd] << " " << distribution[currThreshInd] << std::endl;
246 if (yinBuffer[tau] < minVal && tau > 2){
247 minVal = yinBuffer[tau];
248 minInd = tau;
249 }
250 peakProb[tau] += distribution[currThreshInd];
251 currThreshInd--;
252 } else {
253 tau++;
254 }
255 }
256 double nonPeakProb = 1;
257 for (size_t i = 0; i < yinBufferSize; ++i)
258 {
259 nonPeakProb -= peakProb[i];
260 }
261 // std::cerr << nonPeakProb << std::endl;
262 if (minInd > 0)
263 {
264 // std::cerr << "min set " << minVal << " " << minInd << " " << nonPeakProb << std::endl;
265 peakProb[minInd] += nonPeakProb * minWeight;
266 }
267
268 return peakProb;
269 }
270
271 double
272 YinUtil::parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize)
273 {
274 // this is taken almost literally from Joren Six's Java implementation
275 if (tau == yinBufferSize) // not valid anyway.
276 {
277 return static_cast<double>(tau);
278 }
279
280 double betterTau = 0.0;
281 size_t x0;
282 size_t x2;
283
284 if (tau < 1)
285 {
286 x0 = tau;
287 } else {
288 x0 = tau - 1;
289 }
290
291 if (tau + 1 < yinBufferSize)
292 {
293 x2 = tau + 1;
294 } else {
295 x2 = tau;
296 }
297
298 if (x0 == tau)
299 {
300 if (yinBuffer[tau] <= yinBuffer[x2])
301 {
302 betterTau = tau;
303 } else {
304 betterTau = x2;
305 }
306 }
307 else if (x2 == tau)
308 {
309 if (yinBuffer[tau] <= yinBuffer[x0])
310 {
311 betterTau = tau;
312 }
313 else
314 {
315 betterTau = x0;
316 }
317 }
318 else
319 {
320 float s0, s1, s2;
321 s0 = yinBuffer[x0];
322 s1 = yinBuffer[tau];
323 s2 = yinBuffer[x2];
324 // fixed AUBIO implementation, thanks to Karl Helgason:
325 // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1
326 betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
327
328 // std::cerr << tau << " --> " << betterTau << std::endl;
329
330 }
331 return betterTau;
332 }
333
334 double
335 YinUtil::sumSquare(const double *in, const size_t start, const size_t end)
336 {
337 double out = 0;
338 for (size_t i = start; i < end; ++i)
339 {
340 out += in[i] * in[i];
341 }
342 return out;
343 }