source: src/pitch/pitchyinfast.c @ 5e499fcd

feature/crepe_org
Last change on this file since 5e499fcd was f73f3fb, checked in by Paul Brossier <piem@piem.org>, 6 years ago

[pitch] prevent null pointer dereference in yinfast

  • Property mode set to 100644
File size: 5.7 KB
Line 
1/*
2  Copyright (C) 2003-2017 Paul Brossier <piem@aubio.org>
3
4  This file is part of aubio.
5
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/>.
18
19*/
20
21/* This algorithm was developed by A. de Cheveigné and H. Kawahara and
22 * published in:
23 *
24 * de Cheveigné, A., Kawahara, H. (2002) "YIN, a fundamental frequency
25 * estimator for speech and music", J. Acoust. Soc. Am. 111, 1917-1930.
26 *
27 * see http://recherche.ircam.fr/equipes/pcm/pub/people/cheveign.html
28 */
29
30#include "aubio_priv.h"
31#include "fvec.h"
32#include "mathutils.h"
33#include "cvec.h"
34#include "spectral/fft.h"
35#include "pitch/pitchyinfast.h"
36
37struct _aubio_pitchyinfast_t
38{
39  fvec_t *yin;
40  smpl_t tol;
41  uint_t peak_pos;
42  fvec_t *tmpdata;
43  fvec_t *sqdiff;
44  fvec_t *kernel;
45  fvec_t *samples_fft;
46  fvec_t *kernel_fft;
47  aubio_fft_t *fft;
48};
49
50aubio_pitchyinfast_t *
51new_aubio_pitchyinfast (uint_t bufsize)
52{
53  aubio_pitchyinfast_t *o = AUBIO_NEW (aubio_pitchyinfast_t);
54  o->yin = new_fvec (bufsize / 2);
55  o->tmpdata = new_fvec (bufsize);
56  o->sqdiff = new_fvec (bufsize / 2);
57  o->kernel = new_fvec (bufsize);
58  o->samples_fft = new_fvec (bufsize);
59  o->kernel_fft = new_fvec (bufsize);
60  o->fft = new_aubio_fft (bufsize);
61  if (!o->yin || !o->tmpdata || !o->tmpdata || !o->sqdiff
62      || !o->kernel || !o->samples_fft || !o->kernel || !o->fft)
63  {
64    del_aubio_pitchyinfast(o);
65    return NULL;
66  }
67  o->tol = 0.15;
68  o->peak_pos = 0;
69  return o;
70}
71
72void
73del_aubio_pitchyinfast (aubio_pitchyinfast_t * o)
74{
75  if (o->yin)
76    del_fvec (o->yin);
77  if (o->tmpdata)
78    del_fvec (o->tmpdata);
79  if (o->sqdiff)
80    del_fvec (o->sqdiff);
81  if (o->kernel)
82    del_fvec (o->kernel);
83  if (o->samples_fft)
84    del_fvec (o->samples_fft);
85  if (o->kernel_fft)
86    del_fvec (o->kernel_fft);
87  if (o->fft)
88    del_aubio_fft (o->fft);
89  AUBIO_FREE (o);
90}
91
92/* all the above in one */
93void
94aubio_pitchyinfast_do (aubio_pitchyinfast_t * o, const fvec_t * input, fvec_t * out)
95{
96  const smpl_t tol = o->tol;
97  fvec_t* yin = o->yin;
98  const uint_t length = yin->length;
99  uint_t B = o->tmpdata->length;
100  uint_t W = o->yin->length; // B / 2
101  fvec_t tmp_slice, kernel_ptr;
102  uint_t tau;
103  sint_t period;
104  smpl_t tmp2 = 0.;
105
106  // compute r_t(0) + r_t+tau(0)
107  {
108    fvec_t *squares = o->tmpdata;
109    fvec_weighted_copy(input, input, squares);
110#if 0
111    for (tau = 0; tau < W; tau++) {
112      tmp_slice.data = squares->data + tau;
113      tmp_slice.length = W;
114      o->sqdiff->data[tau] = fvec_sum(&tmp_slice);
115    }
116#else
117    tmp_slice.data = squares->data;
118    tmp_slice.length = W;
119    o->sqdiff->data[0] = fvec_sum(&tmp_slice);
120    for (tau = 1; tau < W; tau++) {
121      o->sqdiff->data[tau] = o->sqdiff->data[tau-1];
122      o->sqdiff->data[tau] -= squares->data[tau-1];
123      o->sqdiff->data[tau] += squares->data[W+tau-1];
124    }
125#endif
126    fvec_add(o->sqdiff, o->sqdiff->data[0]);
127  }
128  // compute r_t(tau) = -2.*ifft(fft(samples)*fft(samples[W-1::-1]))
129  {
130    fvec_t *compmul = o->tmpdata;
131    fvec_t *rt_of_tau = o->samples_fft;
132    aubio_fft_do_complex(o->fft, input, o->samples_fft);
133    // build kernel, take a copy of first half of samples
134    tmp_slice.data = input->data;
135    tmp_slice.length = W;
136    kernel_ptr.data = o->kernel->data + 1;
137    kernel_ptr.length = W;
138    fvec_copy(&tmp_slice, &kernel_ptr);
139    // reverse them
140    fvec_rev(&kernel_ptr);
141    // compute fft(kernel)
142    aubio_fft_do_complex(o->fft, o->kernel, o->kernel_fft);
143    // compute complex product
144    compmul->data[0]  = o->kernel_fft->data[0] * o->samples_fft->data[0];
145    for (tau = 1; tau < W; tau++) {
146      compmul->data[tau]    = o->kernel_fft->data[tau] * o->samples_fft->data[tau];
147      compmul->data[tau]   -= o->kernel_fft->data[B-tau] * o->samples_fft->data[B-tau];
148    }
149    compmul->data[W]    = o->kernel_fft->data[W] * o->samples_fft->data[W];
150    for (tau = 1; tau < W; tau++) {
151      compmul->data[B-tau]  = o->kernel_fft->data[B-tau] * o->samples_fft->data[tau];
152      compmul->data[B-tau] += o->kernel_fft->data[tau] * o->samples_fft->data[B-tau];
153    }
154    // compute inverse fft
155    aubio_fft_rdo_complex(o->fft, compmul, rt_of_tau);
156    // compute square difference r_t(tau) = sqdiff - 2 * r_t_tau[W-1:-1]
157    for (tau = 0; tau < W; tau++) {
158      yin->data[tau] = o->sqdiff->data[tau] - 2. * rt_of_tau->data[tau+W];
159    }
160  }
161
162  // now build yin and look for first minimum
163  fvec_zeros(out);
164  yin->data[0] = 1.;
165  for (tau = 1; tau < length; tau++) {
166    tmp2 += yin->data[tau];
167    if (tmp2 != 0) {
168      yin->data[tau] *= tau / tmp2;
169    } else {
170      yin->data[tau] = 1.;
171    }
172    period = tau - 3;
173    if (tau > 4 && (yin->data[period] < tol) &&
174        (yin->data[period] < yin->data[period + 1])) {
175      o->peak_pos = (uint_t)period;
176      out->data[0] = fvec_quadratic_peak_pos (yin, o->peak_pos);
177      return;
178    }
179  }
180  // use global minimum
181  o->peak_pos = (uint_t)fvec_min_elem (yin);
182  out->data[0] = fvec_quadratic_peak_pos (yin, o->peak_pos);
183}
184
185smpl_t
186aubio_pitchyinfast_get_confidence (aubio_pitchyinfast_t * o) {
187  return 1. - o->yin->data[o->peak_pos];
188}
189
190uint_t
191aubio_pitchyinfast_set_tolerance (aubio_pitchyinfast_t * o, smpl_t tol)
192{
193  o->tol = tol;
194  return 0;
195}
196
197smpl_t
198aubio_pitchyinfast_get_tolerance (aubio_pitchyinfast_t * o)
199{
200  return o->tol;
201}
Note: See TracBrowser for help on using the repository browser.