Mercurial > hg > pyin

annotate 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
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@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 172
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 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};
matthiasm@0 182
matthiasm@0 183 size_t nThreshold = 100;
matthiasm@0 184 int nThresholdInt = nThreshold;
matthiasm@0 185
matthiasm@0 186 for (int i = 0; i < nThresholdInt; ++i)
matthiasm@0 187 {
matthiasm@0 188 switch (prior) {
matthiasm@0 189 case 0:
matthiasm@0 190 distribution.push_back(uniformDist[i]);
matthiasm@0 191 break;
matthiasm@0 192 case 1:
matthiasm@0 193 distribution.push_back(betaDist1[i]);
matthiasm@0 194 break;
matthiasm@0 195 case 2:
matthiasm@0 196 distribution.push_back(betaDist2[i]);
matthiasm@0 197 break;
matthiasm@0 198 case 3:
matthiasm@0 199 distribution.push_back(betaDist3[i]);
matthiasm@0 200 break;
matthiasm@0 201 case 4:
matthiasm@0 202 distribution.push_back(betaDist4[i]);
matthiasm@0 203 break;
matthiasm@0 204 case 5:
matthiasm@0 205 distribution.push_back(single10[i]);
matthiasm@0 206 break;
matthiasm@0 207 case 6:
matthiasm@0 208 distribution.push_back(single15[i]);
matthiasm@0 209 break;
matthiasm@0 210 case 7:
matthiasm@0 211 distribution.push_back(single20[i]);
matthiasm@0 212 break;
matthiasm@0 213 default:
matthiasm@0 214 distribution.push_back(uniformDist[i]);
matthiasm@0 215 }
matthiasm@0 216 thresholds.push_back(0.01 + i*0.01);
matthiasm@0 217 }
matthiasm@0 218
matthiasm@0 219 // double minYin = 2936;
matthiasm@0 220 // for (size_t i = 2; i < yinBufferSize; ++i)
matthiasm@0 221 // {
matthiasm@0 222 // if (yinBuffer[i] < minYin)
matthiasm@0 223 // {
matthiasm@0 224 // minYin = yinBuffer[i];
matthiasm@0 225 // }
matthiasm@0 226 // }
matthiasm@0 227 // if (minYin < 0.01) std::cerr << "min Yin buffer element: " << minYin << std::endl;
matthiasm@0 228
matthiasm@0 229
matthiasm@0 230 int currThreshInd = nThreshold-1;
matthiasm@0 231 tau = 2;
matthiasm@0 232
matthiasm@0 233 // double factor = 1.0 / (0.25 * (nThresholdInt+1) * (nThresholdInt + 1)); // factor to scale down triangular weight
matthiasm@0 234 size_t minInd = 0;
matthiasm@0 235 float minVal = 42.f;
matthiasm@0 236 while (currThreshInd != -1 && tau < yinBufferSize)
matthiasm@0 237 {
matthiasm@0 238 if (yinBuffer[tau] < thresholds[currThreshInd])
matthiasm@0 239 {
matthiasm@0 240 while (tau + 1 < yinBufferSize && yinBuffer[tau+1] < yinBuffer[tau])
matthiasm@0 241 {
matthiasm@0 242 tau++;
matthiasm@0 243 }
matthiasm@0 244 // tau is now local minimum
matthiasm@0 245 // std::cerr << tau << " " << currThreshInd << " "<< thresholds[currThreshInd] << " " << distribution[currThreshInd] << std::endl;
matthiasm@0 246 if (yinBuffer[tau] < minVal && tau > 2){
matthiasm@0 247 minVal = yinBuffer[tau];
matthiasm@0 248 minInd = tau;
matthiasm@0 249 }
matthiasm@0 250 peakProb[tau] += distribution[currThreshInd];
matthiasm@0 251 currThreshInd--;
matthiasm@0 252 } else {
matthiasm@0 253 tau++;
matthiasm@0 254 }
matthiasm@0 255 }
matthiasm@0 256 double nonPeakProb = 1;
matthiasm@0 257 for (size_t i = 0; i < yinBufferSize; ++i)
matthiasm@0 258 {
matthiasm@0 259 nonPeakProb -= peakProb[i];
matthiasm@0 260 }
matthiasm@0 261 // std::cerr << nonPeakProb << std::endl;
matthiasm@0 262 if (minInd > 0)
matthiasm@0 263 {
matthiasm@0 264 // std::cerr << "min set " << minVal << " " << minInd << " " << nonPeakProb << std::endl;
matthiasm@0 265 peakProb[minInd] += nonPeakProb * minWeight;
matthiasm@0 266 }
matthiasm@0 267
matthiasm@0 268 return peakProb;
matthiasm@0 269 }
matthiasm@0 270
matthiasm@0 271 double
matthiasm@0 272 YinUtil::parabolicInterpolation(const double *yinBuffer, const size_t tau, const size_t yinBufferSize)
matthiasm@0 273 {
matthiasm@0 274 // this is taken almost literally from Joren Six's Java implementation
matthiasm@0 275 if (tau == yinBufferSize) // not valid anyway.
matthiasm@0 276 {
matthiasm@0 277 return static_cast<double>(tau);
matthiasm@0 278 }
matthiasm@0 279
matthiasm@0 280 double betterTau = 0.0;
matthiasm@0 281 size_t x0;
matthiasm@0 282 size_t x2;
matthiasm@0 283
matthiasm@0 284 if (tau < 1)
matthiasm@0 285 {
matthiasm@0 286 x0 = tau;
matthiasm@0 287 } else {
matthiasm@0 288 x0 = tau - 1;
matthiasm@0 289 }
matthiasm@0 290
matthiasm@0 291 if (tau + 1 < yinBufferSize)
matthiasm@0 292 {
matthiasm@0 293 x2 = tau + 1;
matthiasm@0 294 } else {
matthiasm@0 295 x2 = tau;
matthiasm@0 296 }
matthiasm@0 297
matthiasm@0 298 if (x0 == tau)
matthiasm@0 299 {
matthiasm@0 300 if (yinBuffer[tau] <= yinBuffer[x2])
matthiasm@0 301 {
matthiasm@0 302 betterTau = tau;
matthiasm@0 303 } else {
matthiasm@0 304 betterTau = x2;
matthiasm@0 305 }
matthiasm@0 306 }
matthiasm@0 307 else if (x2 == tau)
matthiasm@0 308 {
matthiasm@0 309 if (yinBuffer[tau] <= yinBuffer[x0])
matthiasm@0 310 {
matthiasm@0 311 betterTau = tau;
matthiasm@0 312 }
matthiasm@0 313 else
matthiasm@0 314 {
matthiasm@0 315 betterTau = x0;
matthiasm@0 316 }
matthiasm@0 317 }
matthiasm@0 318 else
matthiasm@0 319 {
matthiasm@0 320 float s0, s1, s2;
matthiasm@0 321 s0 = yinBuffer[x0];
matthiasm@0 322 s1 = yinBuffer[tau];
matthiasm@0 323 s2 = yinBuffer[x2];
matthiasm@0 324 // fixed AUBIO implementation, thanks to Karl Helgason:
matthiasm@0 325 // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1
matthiasm@0 326 betterTau = tau + (s2 - s0) / (2 * (2 * s1 - s2 - s0));
matthiasm@0 327
matthiasm@0 328 // std::cerr << tau << " --> " << betterTau << std::endl;
matthiasm@0 329
matthiasm@0 330 }
matthiasm@0 331 return betterTau;
matthiasm@0 332 }
matthiasm@0 333
matthiasm@0 334 double
matthiasm@0 335 YinUtil::sumSquare(const double *in, const size_t start, const size_t end)
matthiasm@0 336 {
matthiasm@0 337 double out = 0;
matthiasm@0 338 for (size_t i = start; i < end; ++i)
matthiasm@0 339 {
matthiasm@0 340 out += in[i] * in[i];
matthiasm@0 341 }
matthiasm@0 342 return out;
matthiasm@0 343 }