Showing the changes to TestPlugin/Source/Yin.cpp (1 of 15 files). Show all
Unified · Split
| @@ -0,0 +1,233 @@ | ||
| 1 | +/* | |
| 2 | + ============================================================================== | |
| 3 | + | |
| 4 | + Yin.c | |
| 5 | + Created: 17 Jan 2017 1:20:23pm | |
| 6 | + Author: Michael R Mulshine | |
| 7 | + | |
| 8 | + ============================================================================== | |
| 9 | +*/ | |
| 10 | + | |
| 11 | +#include "Yin.h" | |
| 12 | + | |
| 13 | +#include <stdint.h> /* For standard interger types (int16_t) */ | |
| 14 | +#include <stdlib.h> /* For call to malloc */ | |
| 15 | + | |
| 16 | +/* ------------------------------------------------------------------------------------------ | |
| 17 | +--------------------------------------------------------------------------- PRIVATE FUNCTIONS | |
| 18 | +-------------------------------------------------------------------------------------------*/ | |
| 19 | + | |
| 20 | +/** | |
| 21 | + * Step 1: Calculates the squared difference of the signal with a shifted version of itself. | |
| 22 | + * @param buffer Buffer of samples to process. | |
| 23 | + * | |
| 24 | + * This is the Yin algorithms tweak on autocorellation. Read http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf | |
| 25 | + * for more details on what is in here and why it's done this way. | |
| 26 | + */ | |
| 27 | +void Yin_difference(Yin *yin, int16_t* buffer){ | |
| 28 | + int16_t i; | |
| 29 | + int16_t tau; | |
| 30 | + float delta; | |
| 31 | + | |
| 32 | + /* Calculate the difference for difference shift values (tau) for the half of the samples */ | |
| 33 | + for(tau = 0 ; tau < yin->halfBufferSize; tau++){ | |
| 34 | + | |
| 35 | + /* Take the difference of the signal with a shifted version of itself, then square it. | |
| 36 | + * (This is the Yin algorithm's tweak on autocorellation) */ | |
| 37 | + for(i = 0; i < yin->halfBufferSize; i++){ | |
| 38 | + delta = buffer[i] - buffer[i + tau]; | |
| 39 | + yin->yinBuffer[tau] += delta * delta; | |
| 40 | + } | |
| 41 | + } | |
| 42 | +} | |
| 43 | + | |
| 44 | + | |
| 45 | +/** | |
| 46 | + * Step 2: Calculate the cumulative mean on the normalised difference calculated in step 1 | |
| 47 | + * @param yin #Yin structure with information about the signal | |
| 48 | + * | |
| 49 | + * This goes through the Yin autocorellation values and finds out roughly where shift is which | |
| 50 | + * produced the smallest difference | |
| 51 | + */ | |
| 52 | +void Yin_cumulativeMeanNormalizedDifference(Yin *yin){ | |
| 53 | + int16_t tau; | |
| 54 | + float runningSum = 0; | |
| 55 | + yin->yinBuffer[0] = 1; | |
| 56 | + | |
| 57 | + /* Sum all the values in the autocorellation buffer and nomalise the result, replacing | |
| 58 | + * the value in the autocorellation buffer with a cumulative mean of the normalised difference */ | |
| 59 | + for (tau = 1; tau < yin->halfBufferSize; tau++) { | |
| 60 | + runningSum += yin->yinBuffer[tau]; | |
| 61 | + yin->yinBuffer[tau] *= tau / runningSum; | |
| 62 | + } | |
| 63 | +} | |
| 64 | + | |
| 65 | +/** | |
| 66 | + * Step 3: Search through the normalised cumulative mean array and find values that are over the threshold | |
| 67 | + * @return Shift (tau) which caused the best approximate autocorellation. -1 if no suitable value is found over the threshold. | |
| 68 | + */ | |
| 69 | +int16_t Yin_absoluteThreshold(Yin *yin){ | |
| 70 | + int16_t tau; | |
| 71 | + | |
| 72 | + /* Search through the array of cumulative mean values, and look for ones that are over the threshold | |
| 73 | + * The first two positions in yinBuffer are always so start at the third (index 2) */ | |
| 74 | + for (tau = 2; tau < yin->halfBufferSize ; tau++) { | |
| 75 | + if (yin->yinBuffer[tau] < yin->threshold) { | |
| 76 | + while (tau + 1 < yin->halfBufferSize && yin->yinBuffer[tau + 1] < yin->yinBuffer[tau]) { | |
| 77 | + tau++; | |
| 78 | + } | |
| 79 | + /* found tau, exit loop and return | |
| 80 | + * store the probability | |
| 81 | + * From the YIN paper: The yin->threshold determines the list of | |
| 82 | + * candidates admitted to the set, and can be interpreted as the | |
| 83 | + * proportion of aperiodic power tolerated | |
| 84 | + * within a periodic signal. | |
| 85 | + * | |
| 86 | + * Since we want the periodicity and and not aperiodicity: | |
| 87 | + * periodicity = 1 - aperiodicity */ | |
| 88 | + yin->probability = 1 - yin->yinBuffer[tau]; | |
| 89 | + break; | |
| 90 | + } | |
| 91 | + } | |
| 92 | + | |
| 93 | + /* if no pitch found, tau => -1 */ | |
| 94 | + if (tau == yin->halfBufferSize || yin->yinBuffer[tau] >= yin->threshold) { | |
| 95 | + tau = -1; | |
| 96 | + yin->probability = 0; | |
| 97 | + } | |
| 98 | + | |
| 99 | + return tau; | |
| 100 | +} | |
| 101 | + | |
| 102 | +/** | |
| 103 | + * Step 5: Interpolate the shift value (tau) to improve the pitch estimate. | |
| 104 | + * @param yin [description] | |
| 105 | + * @param tauEstimate [description] | |
| 106 | + * @return [description] | |
| 107 | + * | |
| 108 | + * The 'best' shift value for autocorellation is most likely not an interger shift of the signal. | |
| 109 | + * As we only autocorellated using integer shifts we should check that there isn't a better fractional | |
| 110 | + * shift value. | |
| 111 | + */ | |
| 112 | +float Yin_parabolicInterpolation(Yin *yin, int16_t tauEstimate) { | |
| 113 | + float betterTau; | |
| 114 | + int16_t x0; | |
| 115 | + int16_t x2; | |
| 116 | + | |
| 117 | + /* Calculate the first polynomial coeffcient based on the current estimate of tau */ | |
| 118 | + if (tauEstimate < 1) { | |
| 119 | + x0 = tauEstimate; | |
| 120 | + } | |
| 121 | + else { | |
| 122 | + x0 = tauEstimate - 1; | |
| 123 | + } | |
| 124 | + | |
| 125 | + /* Calculate the second polynomial coeffcient based on the current estimate of tau */ | |
| 126 | + if (tauEstimate + 1 < yin->halfBufferSize) { | |
| 127 | + x2 = tauEstimate + 1; | |
| 128 | + } | |
| 129 | + else { | |
| 130 | + x2 = tauEstimate; | |
| 131 | + } | |
| 132 | + | |
| 133 | + /* Algorithm to parabolically interpolate the shift value tau to find a better estimate */ | |
| 134 | + if (x0 == tauEstimate) { | |
| 135 | + if (yin->yinBuffer[tauEstimate] <= yin->yinBuffer[x2]) { | |
| 136 | + betterTau = tauEstimate; | |
| 137 | + } | |
| 138 | + else { | |
| 139 | + betterTau = x2; | |
| 140 | + } | |
| 141 | + } | |
| 142 | + else if (x2 == tauEstimate) { | |
| 143 | + if (yin->yinBuffer[tauEstimate] <= yin->yinBuffer[x0]) { | |
| 144 | + betterTau = tauEstimate; | |
| 145 | + } | |
| 146 | + else { | |
| 147 | + betterTau = x0; | |
| 148 | + } | |
| 149 | + } | |
| 150 | + else { | |
| 151 | + float s0, s1, s2; | |
| 152 | + s0 = yin->yinBuffer[x0]; | |
| 153 | + s1 = yin->yinBuffer[tauEstimate]; | |
| 154 | + s2 = yin->yinBuffer[x2]; | |
| 155 | + // fixed AUBIO implementation, thanks to Karl Helgason: | |
| 156 | + // (2.0f * s1 - s2 - s0) was incorrectly multiplied with -1 | |
| 157 | + betterTau = tauEstimate + (s2 - s0) / (2 * (2 * s1 - s2 - s0)); | |
| 158 | + } | |
| 159 | + | |
| 160 | + | |
| 161 | + return betterTau; | |
| 162 | +} | |
| 163 | + | |
| 164 | + | |
| 165 | + | |
| 166 | + | |
| 167 | + | |
| 168 | +/* ------------------------------------------------------------------------------------------ | |
| 169 | +---------------------------------------------------------------------------- PUBLIC FUNCTIONS | |
| 170 | +-------------------------------------------------------------------------------------------*/ | |
| 171 | + | |
| 172 | + | |
| 173 | + | |
| 174 | +/** | |
| 175 | + * Initialise the Yin pitch detection object | |
| 176 | + * @param yin Yin pitch detection object to initialise | |
| 177 | + * @param bufferSize Length of the audio buffer to analyse | |
| 178 | + * @param threshold Allowed uncertainty (e.g 0.05 will return a pitch with ~95% probability) | |
| 179 | + */ | |
| 180 | +void Yin_init(Yin *yin, int16_t bufferSize, float threshold){ | |
| 181 | + /* Initialise the fields of the Yin structure passed in */ | |
| 182 | + yin->bufferSize = bufferSize; | |
| 183 | + yin->halfBufferSize = bufferSize / 2; | |
| 184 | + yin->probability = 0.0; | |
| 185 | + yin->threshold = threshold; | |
| 186 | + | |
| 187 | + /* Allocate the autocorellation buffer and initialise it to zero */ | |
| 188 | + yin->yinBuffer = (float *) malloc(sizeof(float)* yin->halfBufferSize); | |
| 189 | + | |
| 190 | + int16_t i; | |
| 191 | + for(i = 0; i < yin->halfBufferSize; i++){ | |
| 192 | + yin->yinBuffer[i] = 0; | |
| 193 | + } | |
| 194 | +} | |
| 195 | + | |
| 196 | +/** | |
| 197 | + * Runs the Yin pitch detection algortihm | |
| 198 | + * @param yin Initialised Yin object | |
| 199 | + * @param buffer Buffer of samples to analyse | |
| 200 | + * @return Fundamental frequency of the signal in Hz. Returns -1 if pitch can't be found | |
| 201 | + */ | |
| 202 | +float Yin_getPitch(Yin *yin, int16_t* buffer){ | |
| 203 | + int16_t tauEstimate = -1; | |
| 204 | + float pitchInHertz = -1; | |
| 205 | + | |
| 206 | + /* Step 1: Calculates the squared difference of the signal with a shifted version of itself. */ | |
| 207 | + Yin_difference(yin, buffer); | |
| 208 | + | |
| 209 | + /* Step 2: Calculate the cumulative mean on the normalised difference calculated in step 1 */ | |
| 210 | + Yin_cumulativeMeanNormalizedDifference(yin); | |
| 211 | + | |
| 212 | + /* Step 3: Search through the normalised cumulative mean array and find values that are over the threshold */ | |
| 213 | + tauEstimate = Yin_absoluteThreshold(yin); | |
| 214 | + | |
| 215 | + /* Step 5: Interpolate the shift value (tau) to improve the pitch estimate. */ | |
| 216 | + if(tauEstimate != -1){ | |
| 217 | + pitchInHertz = YIN_SAMPLING_RATE / Yin_parabolicInterpolation(yin, tauEstimate); | |
| 218 | + } | |
| 219 | + | |
| 220 | + return pitchInHertz; | |
| 221 | +} | |
| 222 | + | |
| 223 | +/** | |
| 224 | + * Certainty of the pitch found | |
| 225 | + * @param yin Yin object that has been run over a buffer | |
| 226 | + * @return Returns the certainty of the note found as a decimal (i.e 0.3 is 30%) | |
| 227 | + */ | |
| 228 | +float Yin_getProbability(Yin *yin){ | |
| 229 | + return yin->probability; | |
| 230 | +} | |
| 231 | + | |
| 232 | + | |
| 233 | + | |