[650e39b] | 1 | /* |
---|
[a5f4b7d] | 2 | Copyright (C) 2003-2013 Paul Brossier <piem@aubio.org> |
---|
[650e39b] | 3 | |
---|
[e6a78ea] | 4 | This file is part of aubio. |
---|
[650e39b] | 5 | |
---|
[e6a78ea] | 6 | aubio is free software: you can redistribute it and/or modify |
---|
| 7 | it under the terms of the GNU General Public License as published by |
---|
| 8 | the Free Software Foundation, either version 3 of the License, or |
---|
| 9 | (at your option) any later version. |
---|
| 10 | |
---|
| 11 | aubio is distributed in the hope that it will be useful, |
---|
| 12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
---|
| 13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
---|
| 14 | GNU General Public License for more details. |
---|
| 15 | |
---|
| 16 | You should have received a copy of the GNU General Public License |
---|
| 17 | along with aubio. If not, see <http://www.gnu.org/licenses/>. |
---|
[650e39b] | 18 | |
---|
| 19 | */ |
---|
| 20 | |
---|
| 21 | #include "aubio_priv.h" |
---|
[6c7d49b] | 22 | #include "fvec.h" |
---|
| 23 | #include "cvec.h" |
---|
[650e39b] | 24 | #include "mathutils.h" |
---|
[32d6958] | 25 | #include "spectral/fft.h" |
---|
[2d8cffa] | 26 | #include "pitch/pitchyinfft.h" |
---|
[650e39b] | 27 | |
---|
[91879d9] | 28 | /** pitch yinfft structure */ |
---|
[fddfa64] | 29 | struct _aubio_pitchyinfft_t |
---|
| 30 | { |
---|
| 31 | fvec_t *win; /**< temporal weighting window */ |
---|
| 32 | fvec_t *winput; /**< windowed spectrum */ |
---|
| 33 | fvec_t *sqrmag; /**< square difference function */ |
---|
| 34 | fvec_t *weight; /**< spectral weighting window (psychoacoustic model) */ |
---|
[4a95f83] | 35 | fvec_t *fftout; /**< Fourier transform output */ |
---|
[fddfa64] | 36 | aubio_fft_t *fft; /**< fft object to compute square difference function */ |
---|
| 37 | fvec_t *yinfft; /**< Yin function */ |
---|
[22d33e2] | 38 | smpl_t tol; /**< Yin tolerance */ |
---|
[5284e0d] | 39 | smpl_t confidence; /**< confidence */ |
---|
[b4eadf0] | 40 | uint_t samplerate; /**< samplerate we got initialized with */ |
---|
[650e39b] | 41 | }; |
---|
| 42 | |
---|
[88fee8f] | 43 | static const smpl_t freqs[] = { |
---|
| 44 | 0., 20., 25., 31.5, 40., 50., 63., 80., 100., 125., |
---|
| 45 | 160., 200., 250., 315., 400., 500., 630., 800., 1000., 1250., |
---|
| 46 | 1600., 2000., 2500., 3150., 4000., 5000., 6300., 8000., 9000., 10000., |
---|
| 47 | 12500., 15000., 20000., 25100 |
---|
[fddfa64] | 48 | }; |
---|
[650e39b] | 49 | |
---|
[88fee8f] | 50 | static const smpl_t weight[] = { |
---|
| 51 | -75.8, -70.1, -60.8, -52.1, -44.2, -37.5, -31.3, -25.6, -20.9, -16.5, |
---|
| 52 | -12.6, -9.60, -7.00, -4.70, -3.00, -1.80, -0.80, -0.20, -0.00, 0.50, |
---|
| 53 | 1.60, 3.20, 5.40, 7.80, 8.10, 5.30, -2.40, -11.1, -12.8, -12.2, |
---|
| 54 | -7.40, -17.8, -17.8, -17.8 |
---|
[fddfa64] | 55 | }; |
---|
[650e39b] | 56 | |
---|
[fddfa64] | 57 | aubio_pitchyinfft_t * |
---|
[9c9202f] | 58 | new_aubio_pitchyinfft (uint_t samplerate, uint_t bufsize) |
---|
[650e39b] | 59 | { |
---|
[c21acb9] | 60 | uint_t i = 0, j = 1; |
---|
| 61 | smpl_t freq = 0, a0 = 0, a1 = 0, f0 = 0, f1 = 0; |
---|
[fddfa64] | 62 | aubio_pitchyinfft_t *p = AUBIO_NEW (aubio_pitchyinfft_t); |
---|
[168337e] | 63 | p->winput = new_fvec (bufsize); |
---|
| 64 | p->fft = new_aubio_fft (bufsize); |
---|
[0f2c1f4] | 65 | if (!p->fft) goto beach; |
---|
[4a95f83] | 66 | p->fftout = new_fvec (bufsize); |
---|
[168337e] | 67 | p->sqrmag = new_fvec (bufsize); |
---|
| 68 | p->yinfft = new_fvec (bufsize / 2 + 1); |
---|
[fddfa64] | 69 | p->tol = 0.85; |
---|
| 70 | p->win = new_aubio_window ("hanningz", bufsize); |
---|
[168337e] | 71 | p->weight = new_fvec (bufsize / 2 + 1); |
---|
[b4eadf0] | 72 | p->samplerate = samplerate; |
---|
[168337e] | 73 | for (i = 0; i < p->weight->length; i++) { |
---|
[9c9202f] | 74 | freq = (smpl_t) i / (smpl_t) bufsize *(smpl_t) samplerate; |
---|
[168337e] | 75 | while (freq > freqs[j]) { |
---|
| 76 | j += 1; |
---|
[10cf306] | 77 | } |
---|
[168337e] | 78 | a0 = weight[j - 1]; |
---|
| 79 | f0 = freqs[j - 1]; |
---|
| 80 | a1 = weight[j]; |
---|
| 81 | f1 = freqs[j]; |
---|
| 82 | if (f0 == f1) { // just in case |
---|
| 83 | p->weight->data[i] = a0; |
---|
| 84 | } else if (f0 == 0) { // y = ax+b |
---|
| 85 | p->weight->data[i] = (a1 - a0) / f1 * freq + a0; |
---|
| 86 | } else { |
---|
| 87 | p->weight->data[i] = (a1 - a0) / (f1 - f0) * freq + |
---|
| 88 | (a0 - (a1 - a0) / (f1 / f0 - 1.)); |
---|
| 89 | } |
---|
| 90 | while (freq > freqs[j]) { |
---|
| 91 | j += 1; |
---|
| 92 | } |
---|
| 93 | //AUBIO_DBG("%f\n",p->weight->data[i]); |
---|
| 94 | p->weight->data[i] = DB2LIN (p->weight->data[i]); |
---|
| 95 | //p->weight->data[i] = SQRT(DB2LIN(p->weight->data[i])); |
---|
[650e39b] | 96 | } |
---|
[b4eadf0] | 97 | |
---|
| 98 | // disable weighting |
---|
| 99 | fvec_set_all (p->weight, 1.0); |
---|
| 100 | |
---|
[650e39b] | 101 | return p; |
---|
[0f2c1f4] | 102 | |
---|
| 103 | beach: |
---|
| 104 | if (p->winput) del_fvec(p->winput); |
---|
| 105 | AUBIO_FREE(p); |
---|
| 106 | return NULL; |
---|
[650e39b] | 107 | } |
---|
| 108 | |
---|
[fddfa64] | 109 | void |
---|
[ce3ff2b] | 110 | aubio_pitchyinfft_do (aubio_pitchyinfft_t * p, const fvec_t * input, fvec_t * output) |
---|
[fddfa64] | 111 | { |
---|
[168337e] | 112 | uint_t tau, l; |
---|
[a5f4b7d] | 113 | uint_t length = p->fftout->length; |
---|
| 114 | fvec_t *fftout = p->fftout; |
---|
| 115 | fvec_t *yin = p->yinfft; |
---|
| 116 | smpl_t tmp = 0., sum = 0.; |
---|
[09558b8] | 117 | uint_t startbin = 0, endbin = 0, before = 0, after = 0; |
---|
[b4eadf0] | 118 | |
---|
[a5f4b7d] | 119 | // window the input |
---|
[c9e3a4e] | 120 | fvec_weighted_copy(input, p->win, p->winput); |
---|
[a5f4b7d] | 121 | // get the real / imag parts of its fft |
---|
| 122 | aubio_fft_do_complex (p->fft, p->winput, fftout); |
---|
| 123 | // get the squared magnitude spectrum, applying some weight |
---|
| 124 | p->sqrmag->data[0] = SQR(fftout->data[0]); |
---|
[8da0033] | 125 | p->sqrmag->data[0] *= p->weight->data[0]; |
---|
[4a95f83] | 126 | for (l = 1; l < length / 2; l++) { |
---|
[a5f4b7d] | 127 | p->sqrmag->data[l] = SQR(fftout->data[l]) + SQR(fftout->data[length - l]); |
---|
[168337e] | 128 | p->sqrmag->data[l] *= p->weight->data[l]; |
---|
[a5f4b7d] | 129 | p->sqrmag->data[length - l] = p->sqrmag->data[l]; |
---|
[168337e] | 130 | } |
---|
[a5f4b7d] | 131 | p->sqrmag->data[length / 2] = SQR(fftout->data[length / 2]); |
---|
[4a95f83] | 132 | p->sqrmag->data[length / 2] *= p->weight->data[length / 2]; |
---|
[a5f4b7d] | 133 | // get sum of weighted squared mags |
---|
[4a95f83] | 134 | for (l = 0; l < length / 2 + 1; l++) { |
---|
[168337e] | 135 | sum += p->sqrmag->data[l]; |
---|
| 136 | } |
---|
| 137 | sum *= 2.; |
---|
[a5f4b7d] | 138 | // get the real / imag parts of the fft of the squared magnitude |
---|
| 139 | aubio_fft_do_complex (p->fft, p->sqrmag, fftout); |
---|
[168337e] | 140 | yin->data[0] = 1.; |
---|
| 141 | for (tau = 1; tau < yin->length; tau++) { |
---|
[a5f4b7d] | 142 | // compute the square differences |
---|
| 143 | yin->data[tau] = sum - fftout->data[tau]; |
---|
| 144 | // and the cumulative mean normalized difference function |
---|
[168337e] | 145 | tmp += yin->data[tau]; |
---|
[e391790] | 146 | if (tmp != 0) { |
---|
| 147 | yin->data[tau] *= tau / tmp; |
---|
| 148 | } else { |
---|
| 149 | yin->data[tau] = 1.; |
---|
| 150 | } |
---|
[168337e] | 151 | } |
---|
[b4eadf0] | 152 | |
---|
| 153 | // calc min available confidence first |
---|
| 154 | tmp = fvec_min(yin); |
---|
| 155 | if (tmp > p->tol) { |
---|
[09558b8] | 156 | // give up - got no confident candidate at all |
---|
[168337e] | 157 | output->data[0] = 0.; |
---|
[b4eadf0] | 158 | return; |
---|
| 159 | } |
---|
| 160 | |
---|
| 161 | // choose lowest confident candidate first, to avoid choosing harmonics |
---|
| 162 | tau = 0; |
---|
| 163 | for (l = 1; l < yin->length; l++) { |
---|
| 164 | // is this candidate "roughly" as good as the lowest one? |
---|
| 165 | if (ABS (yin->data[l] - tmp) < 0.1) { |
---|
| 166 | tau = l; |
---|
| 167 | break; |
---|
| 168 | } |
---|
| 169 | } |
---|
| 170 | // find local min around current pick to sharpen the results |
---|
[09558b8] | 171 | before = 2, after = 5; |
---|
| 172 | // for low frequency, the peak may be large, we use a larger after |
---|
| 173 | if (tau > yin->length/2) { before = 1; after = 20; } |
---|
| 174 | startbin = tau > before ? tau - before : 0; |
---|
| 175 | endbin = tau < yin->length - 1 - after ? tau + after : yin->length - 1; |
---|
[b4eadf0] | 176 | tmp = yin->data[tau]; |
---|
| 177 | for (l = startbin; l < endbin; l++) { |
---|
| 178 | if (yin->data[l] < tmp ) { |
---|
| 179 | tmp = yin->data[l]; |
---|
| 180 | tau = l; |
---|
| 181 | } |
---|
[22d33e2] | 182 | } |
---|
[b4eadf0] | 183 | output->data[0] = fvec_quadratic_peak_pos(yin, tau); |
---|
[650e39b] | 184 | } |
---|
| 185 | |
---|
[fddfa64] | 186 | void |
---|
| 187 | del_aubio_pitchyinfft (aubio_pitchyinfft_t * p) |
---|
| 188 | { |
---|
| 189 | del_fvec (p->win); |
---|
| 190 | del_aubio_fft (p->fft); |
---|
| 191 | del_fvec (p->yinfft); |
---|
| 192 | del_fvec (p->sqrmag); |
---|
[4a95f83] | 193 | del_fvec (p->fftout); |
---|
[fddfa64] | 194 | del_fvec (p->winput); |
---|
| 195 | del_fvec (p->weight); |
---|
| 196 | AUBIO_FREE (p); |
---|
[650e39b] | 197 | } |
---|
[22d33e2] | 198 | |
---|
[5284e0d] | 199 | smpl_t |
---|
| 200 | aubio_pitchyinfft_get_confidence (aubio_pitchyinfft_t * o) { |
---|
| 201 | o->confidence = 1. - fvec_min (o->yinfft); |
---|
| 202 | return o->confidence; |
---|
| 203 | } |
---|
| 204 | |
---|
[fddfa64] | 205 | uint_t |
---|
| 206 | aubio_pitchyinfft_set_tolerance (aubio_pitchyinfft_t * p, smpl_t tol) |
---|
| 207 | { |
---|
[22d33e2] | 208 | p->tol = tol; |
---|
| 209 | return 0; |
---|
| 210 | } |
---|
| 211 | |
---|
[fddfa64] | 212 | smpl_t |
---|
| 213 | aubio_pitchyinfft_get_tolerance (aubio_pitchyinfft_t * p) |
---|
| 214 | { |
---|
[22d33e2] | 215 | return p->tol; |
---|
| 216 | } |
---|