matthiasm@0: #include "YinUtil.h"
matthiasm@0: 
matthiasm@0: #include <vector>
matthiasm@0: 
matthiasm@0: #include <cstdio>
matthiasm@0: #include <cmath>
matthiasm@0: #include <algorithm>
matthiasm@0: 
matthiasm@0: #include <boost/math/distributions.hpp>
matthiasm@0: 
matthiasm@0: void 
matthiasm@0: YinUtil::fastDifference(const double *in, double *yinBuffer, const size_t yinBufferSize) 
matthiasm@0: {
matthiasm@0:     
matthiasm@0:     // DECLARE AND INITIALISE
matthiasm@0:     // initialisation of most of the arrays here was done in a separate function,
matthiasm@0:     // with all the arrays as members of the class... moved them back here.
matthiasm@0:     
matthiasm@0:     size_t frameSize = 2 * yinBufferSize;
matthiasm@0:     
matthiasm@0:     for (size_t j = 0; j < yinBufferSize; ++j)
matthiasm@0:     {
matthiasm@0:         yinBuffer[j] = 0.;
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     double *audioTransformedReal = new double[frameSize];
matthiasm@0:     double *audioTransformedImag = new double[frameSize];
matthiasm@0:     double *nullImag = new double[frameSize];
matthiasm@0:     double *kernel = new double[frameSize];
matthiasm@0:     double *kernelTransformedReal = new double[frameSize];
matthiasm@0:     double *kernelTransformedImag = new double[frameSize];
matthiasm@0:     double *yinStyleACFReal = new double[frameSize];
matthiasm@0:     double *yinStyleACFImag = new double[frameSize];
matthiasm@0:     double *powerTerms = new double[yinBufferSize];
matthiasm@0:     
matthiasm@0:     for (size_t j = 0; j < yinBufferSize; ++j)
matthiasm@0:     {
matthiasm@0:         powerTerms[j] = 0.;
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     for (size_t j = 0; j < frameSize; ++j)
matthiasm@0:     {
matthiasm@0:         nullImag[j] = 0.;
matthiasm@0:         audioTransformedReal[j] = 0.;
matthiasm@0:         audioTransformedImag[j] = 0.;
matthiasm@0:         kernel[j] = 0.;
matthiasm@0:         kernelTransformedReal[j] = 0.;
matthiasm@0:         kernelTransformedImag[j] = 0.;
matthiasm@0:         yinStyleACFReal[j] = 0.;
matthiasm@0:         yinStyleACFImag[j] = 0.;
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     // POWER TERM CALCULATION
matthiasm@0:     // ... for the power terms in equation (7) in the Yin paper
matthiasm@0:     powerTerms[0] = 0.0;
matthiasm@0:     for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0:         powerTerms[0] += in[j] * in[j];
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     // now iteratively calculate all others (saves a few multiplications)
matthiasm@0:     for (size_t tau = 1; tau < yinBufferSize; ++tau) {
matthiasm@0:         powerTerms[tau] = powerTerms[tau-1] - in[tau-1] * in[tau-1] + in[tau+yinBufferSize] * in[tau+yinBufferSize];  
matthiasm@0:     }
matthiasm@0: 
matthiasm@0:     // YIN-STYLE AUTOCORRELATION via FFT
matthiasm@0:     // 1. data
matthiasm@0:     Vamp::FFT::forward(frameSize, in, nullImag, audioTransformedReal, audioTransformedImag);
matthiasm@0:     
matthiasm@0:     // 2. half of the data, disguised as a convolution kernel
matthiasm@0:     for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0:         kernel[j] = in[yinBufferSize-1-j];
matthiasm@0:         kernel[j+yinBufferSize] = 0;
matthiasm@0:     }
matthiasm@0:     Vamp::FFT::forward(frameSize, kernel, nullImag, kernelTransformedReal, kernelTransformedImag);
matthiasm@0: 
matthiasm@0:     // 3. convolution via complex multiplication -- written into
matthiasm@0:     for (size_t j = 0; j < frameSize; ++j) {
matthiasm@0:         yinStyleACFReal[j] = audioTransformedReal[j]*kernelTransformedReal[j] - audioTransformedImag[j]*kernelTransformedImag[j]; // real
matthiasm@0:         yinStyleACFImag[j] = audioTransformedReal[j]*kernelTransformedImag[j] + audioTransformedImag[j]*kernelTransformedReal[j]; // imaginary
matthiasm@0:     }
matthiasm@0:     Vamp::FFT::inverse(frameSize, yinStyleACFReal, yinStyleACFImag, audioTransformedReal, audioTransformedImag);
matthiasm@0:     
matthiasm@0:     // CALCULATION OF difference function
matthiasm@0:     // ... according to (7) in the Yin paper.
matthiasm@0:     for (size_t j = 0; j < yinBufferSize; ++j) {
matthiasm@0:         // taking only the real part
matthiasm@0:         yinBuffer[j] = powerTerms[0] + powerTerms[j] - 2 * audioTransformedReal[j+yinBufferSize-1];
matthiasm@0:     }
matthiasm@0:     delete [] audioTransformedReal;
matthiasm@0:     delete [] audioTransformedImag;
matthiasm@0:     delete [] nullImag;
matthiasm@0:     delete [] kernel;
matthiasm@0:     delete [] kernelTransformedReal;
matthiasm@0:     delete [] kernelTransformedImag;
matthiasm@0:     delete [] yinStyleACFReal;
matthiasm@0:     delete [] yinStyleACFImag;
matthiasm@0:     delete [] powerTerms;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: void 
matthiasm@0: YinUtil::cumulativeDifference(double *yinBuffer, const size_t yinBufferSize)
matthiasm@0: {    
matthiasm@0:     size_t tau;
matthiasm@0:     
matthiasm@0:     yinBuffer[0] = 1;
matthiasm@0:     
matthiasm@0:     double runningSum = 0;
matthiasm@0:     
matthiasm@0:     for (tau = 1; tau < yinBufferSize; ++tau) {
matthiasm@0:         runningSum += yinBuffer[tau];
matthiasm@0:         if (runningSum == 0)
matthiasm@0:         {
matthiasm@0:             yinBuffer[tau] = 1;
matthiasm@0:         } else {
matthiasm@0:             yinBuffer[tau] *= tau / runningSum;
matthiasm@0:         }
matthiasm@0:     }    
matthiasm@0: }
matthiasm@0: 
matthiasm@0: int 
matthiasm@0: YinUtil::absoluteThreshold(const double *yinBuffer, const size_t yinBufferSize, const double thresh)
matthiasm@0: {
matthiasm@0:     size_t tau;
matthiasm@0:     size_t minTau = 0;
matthiasm@0:     double minVal = 1000.;
matthiasm@0:     
matthiasm@0:     // using Joren Six's "loop construct" from TarsosDSP
matthiasm@0:     tau = 2;
matthiasm@0:     while (tau < yinBufferSize)
matthiasm@0:     {
matthiasm@0:         if (yinBuffer[tau] < thresh)
matthiasm@0:         {
matthiasm@0:             while (tau+1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
matthiasm@0:             {
matthiasm@0:                 ++tau;
matthiasm@0:             }
matthiasm@0:             return tau;
matthiasm@0:         } else {
matthiasm@0:             if (yinBuffer[tau] < minVal)
matthiasm@0:             {
matthiasm@0:                 minVal = yinBuffer[tau];
matthiasm@0:                 minTau = tau;
matthiasm@0:             }
matthiasm@0:         }
matthiasm@0:         ++tau;
matthiasm@0:     }
matthiasm@0:     if (minTau > 0)
matthiasm@0:     {
matthiasm@0:         return -minTau;
matthiasm@0:     }
matthiasm@0:     return 0;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: 
matthiasm@0: std::vector<double>
matthiasm@0: YinUtil::yinProb(const double *yinBuffer, const size_t prior, const size_t yinBufferSize) 
matthiasm@0: {
matthiasm@0:     double minWeight = 0.01;
matthiasm@0:     size_t tau;
matthiasm@0:     std::vector<float> thresholds;
matthiasm@0:     std::vector<float> distribution;
matthiasm@0:     std::vector<double> peakProb = std::vector<double>(yinBufferSize);
matthiasm@0:     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:     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:     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:     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:     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:     // 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};
matthiasm@0:     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:     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:     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:     
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     // 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};
matthiasm@0:     
matthiasm@0:     size_t nThreshold = 100;
matthiasm@0:     int nThresholdInt = nThreshold;
matthiasm@0:     
matthiasm@0:     for (int i = 0; i < nThresholdInt; ++i)
matthiasm@0:     {
matthiasm@0:         switch (prior) {
matthiasm@0:             case 0:
matthiasm@0:                 distribution.push_back(uniformDist[i]);
matthiasm@0:                 break;
matthiasm@0:             case 1:
matthiasm@0:                 distribution.push_back(betaDist1[i]);
matthiasm@0:                 break;
matthiasm@0:             case 2:
matthiasm@0:                 distribution.push_back(betaDist2[i]);
matthiasm@0:                 break;
matthiasm@0:             case 3:
matthiasm@0:                 distribution.push_back(betaDist3[i]);
matthiasm@0:                 break;
matthiasm@0:             case 4:
matthiasm@0:                 distribution.push_back(betaDist4[i]);
matthiasm@0:                 break;
matthiasm@0:             case 5:
matthiasm@0:                 distribution.push_back(single10[i]);
matthiasm@0:                 break;
matthiasm@0:             case 6:
matthiasm@0:                 distribution.push_back(single15[i]);
matthiasm@0:                 break;
matthiasm@0:             case 7:
matthiasm@0:                 distribution.push_back(single20[i]);
matthiasm@0:                 break;
matthiasm@0:             default:
matthiasm@0:                 distribution.push_back(uniformDist[i]);
matthiasm@0:         }
matthiasm@0:         thresholds.push_back(0.01 + i*0.01);
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     // double minYin = 2936;
matthiasm@0:     // for (size_t i = 2; i < yinBufferSize; ++i)
matthiasm@0:     // {
matthiasm@0:     //     if (yinBuffer[i] < minYin)
matthiasm@0:     //     {
matthiasm@0:     //         minYin = yinBuffer[i];
matthiasm@0:     //     }
matthiasm@0:     // }
matthiasm@0:     // if (minYin < 0.01) std::cerr << "min Yin buffer element: " << minYin << std::endl;
matthiasm@0:     
matthiasm@0:     
matthiasm@0:     int currThreshInd = nThreshold-1;
matthiasm@0:     tau = 2;
matthiasm@0:     
matthiasm@0:     // double factor = 1.0 / (0.25 * (nThresholdInt+1) * (nThresholdInt + 1)); // factor to scale down triangular weight
matthiasm@0:     size_t minInd = 0;
matthiasm@0:     float minVal = 42.f;
matthiasm@0:     while (currThreshInd != -1 && tau < yinBufferSize)
matthiasm@0:     {
matthiasm@0:         if (yinBuffer[tau] < thresholds[currThreshInd])
matthiasm@0:         {
matthiasm@0:             while (tau + 1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
matthiasm@0:             {
matthiasm@0:                 tau++;
matthiasm@0:             }
matthiasm@0:             // tau is now local minimum
matthiasm@0:             // std::cerr << tau << " " << currThreshInd << " "<< thresholds[currThreshInd] << " " << distribution[currThreshInd] << std::endl;
matthiasm@0:             if (yinBuffer[tau] < minVal && tau > 2){
matthiasm@0:                 minVal = yinBuffer[tau];
matthiasm@0:                 minInd = tau;
matthiasm@0:             }
matthiasm@0:             peakProb[tau] += distribution[currThreshInd];
matthiasm@0:             currThreshInd--;
matthiasm@0:         } else {
matthiasm@0:             tau++;
matthiasm@0:         }
matthiasm@0:     }
matthiasm@0:     double nonPeakProb = 1;
matthiasm@0:     for (size_t i = 0; i < yinBufferSize; ++i)
matthiasm@0:     {
matthiasm@0:         nonPeakProb -= peakProb[i];
matthiasm@0:     }
matthiasm@0:     // std::cerr << nonPeakProb << std::endl;
matthiasm@0:     if (minInd > 0)
matthiasm@0:     {
matthiasm@0:         // std::cerr << "min set " << minVal << " " << minInd << " " << nonPeakProb << std::endl; 
matthiasm@0:         peakProb[minInd] += nonPeakProb * minWeight;
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     return peakProb;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: double
matthiasm@0: YinUtil::parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize) 
matthiasm@0: {
matthiasm@0:     // this is taken almost literally from Joren Six's Java implementation
matthiasm@0:     if (tau == yinBufferSize) // not valid anyway.
matthiasm@0:     {
matthiasm@0:         return static_cast<double>(tau);
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     double betterTau = 0.0;
matthiasm@0:     size_t x0;
matthiasm@0:     size_t x2;
matthiasm@0: 
matthiasm@0:     if (tau < 1) 
matthiasm@0:     {
matthiasm@0:         x0 = tau;
matthiasm@0:     } else {
matthiasm@0:         x0 = tau - 1;
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     if (tau + 1 < yinBufferSize) 
matthiasm@0:     {
matthiasm@0:         x2 = tau + 1;
matthiasm@0:     } else {
matthiasm@0:         x2 = tau;
matthiasm@0:     }
matthiasm@0:     
matthiasm@0:     if (x0 == tau) 
matthiasm@0:     {
matthiasm@0:         if (yinBuffer[tau] <= yinBuffer[x2]) 
matthiasm@0:         {
matthiasm@0:             betterTau = tau;
matthiasm@0:         } else {
matthiasm@0:             betterTau = x2;
matthiasm@0:         }
matthiasm@0:     } 
matthiasm@0:     else if (x2 == tau) 
matthiasm@0:     {
matthiasm@0:         if (yinBuffer[tau] <= yinBuffer[x0]) 
matthiasm@0:         {
matthiasm@0:             betterTau = tau;
matthiasm@0:         } 
matthiasm@0:         else 
matthiasm@0:         {
matthiasm@0:             betterTau = x0;
matthiasm@0:         }
matthiasm@0:     } 
matthiasm@0:     else 
matthiasm@0:     {
matthiasm@0:         float s0, s1, s2;
matthiasm@0:         s0 = yinBuffer[x0];
matthiasm@0:         s1 = yinBuffer[tau];
matthiasm@0:         s2 = yinBuffer[x2];
matthiasm@0:         // fixed AUBIO implementation, thanks to Karl Helgason:
matthiasm@0:         // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1
matthiasm@0:         betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
matthiasm@0:         
matthiasm@0:         // std::cerr << tau << " --> " << betterTau << std::endl;
matthiasm@0:         
matthiasm@0:     }
matthiasm@0:     return betterTau;
matthiasm@0: }
matthiasm@0: 
matthiasm@0: double 
matthiasm@0: YinUtil::sumSquare(const double *in, const size_t start, const size_t end)
matthiasm@0: {
matthiasm@0:     double out = 0;
matthiasm@0:     for (size_t i = start; i < end; ++i)
matthiasm@0:     {
matthiasm@0:         out += in[i] * in[i];
matthiasm@0:     }
matthiasm@0:     return out;
matthiasm@0: }