2 Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
4 This file is part of aubio.
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.
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.
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/>.
21 #include "aubio_priv.h"
24 #include "mathutils.h"
25 #include "pitch/pitchmcomb.h"
27 #define CAND_SWAP(a,b) { register aubio_spectralcandidate_t *t=(a);(a)=(b);(b)=t; }
29 typedef struct _aubio_spectralpeak_t aubio_spectralpeak_t;
30 typedef struct _aubio_spectralcandidate_t aubio_spectralcandidate_t;
31 uint_t aubio_pitchmcomb_get_root_peak (aubio_spectralpeak_t * peaks,
33 uint_t aubio_pitchmcomb_quadpick (aubio_spectralpeak_t * spectral_peaks,
35 void aubio_pitchmcomb_spectral_pp (aubio_pitchmcomb_t * p, fvec_t * oldmag);
36 void aubio_pitchmcomb_combdet (aubio_pitchmcomb_t * p, fvec_t * newmag);
37 /* not used but useful : sort by amplitudes (or anything else)
38 * sort_pitchpeak(peaks, length);
40 /** spectral_peak comparison function (must return signed int) */
41 static sint_t aubio_pitchmcomb_sort_peak_comp (const void *x, const void *y);
42 /** sort spectral_peak against their mag */
43 void aubio_pitchmcomb_sort_peak (aubio_spectralpeak_t * peaks, uint_t nbins);
44 /** select the best candidates */
45 uint_t aubio_pitch_cands (aubio_pitchmcomb_t * p, cvec_t * fftgrain,
48 /** sort spectral_candidate against their comb ene */
49 void aubio_pitchmcomb_sort_cand_ene (aubio_spectralcandidate_t ** candidates,
51 /** sort spectral_candidate against their frequency */
52 void aubio_pitchmcomb_sort_cand_freq (aubio_spectralcandidate_t ** candidates,
55 struct _aubio_pitchmcomb_t
57 smpl_t threshold; /**< offset threshold [0.033 or 0.01] */
58 smpl_t alpha; /**< normalisation exponent [9] */
59 smpl_t cutoff; /**< low-pass filter cutoff [0.34, 1] */
60 smpl_t tol; /**< tolerance [0.05] */
61 // smpl_t tau; /**< frequency precision [44100/4096] */
62 uint_t win_post; /**< median filter window length */
63 uint_t win_pre; /**< median filter window */
64 uint_t ncand; /**< maximum number of candidates (combs) */
65 uint_t npartials; /**< maximum number of partials per combs */
66 uint_t count; /**< picked picks */
67 uint_t goodcandidate; /**< best candidate */
68 uint_t spec_partition; /**< spectrum partition to consider */
69 aubio_spectralpeak_t *peaks; /**< up to length win/spec_partition */
70 aubio_spectralcandidate_t **candidates; /** up to five candidates */
71 /* some scratch pads */
72 /** \bug (unnecessary copied from fftgrain?) */
73 fvec_t *newmag; /**< vec to store mag */
74 fvec_t *scratch; /**< vec to store modified mag */
75 fvec_t *scratch2; /**< vec to compute moving median */
76 fvec_t *theta; /**< vec to store phase */
79 /** threshfn: name or handle of fn for computing adaptive threshold [median] */
80 /** aubio_thresholdfn_t thresholdfn; */
81 /** picker: name or handle of fn for picking event times [quadpick] */
82 /** aubio_pickerfn_t pickerfn; */
85 /** spectral peak object */
86 struct _aubio_spectralpeak_t
88 uint_t bin; /**< bin [0-(length-1)] */
89 smpl_t ebin; /**< estimated bin */
90 smpl_t mag; /**< peak magnitude */
93 /** spectral candidates array object */
94 struct _aubio_spectralcandidate_t
96 smpl_t ebin; /**< interpolated bin */
97 smpl_t *ecomb; /**< comb */
98 smpl_t ene; /**< candidate energy */
99 smpl_t len; /**< length */
104 aubio_pitchmcomb_do (aubio_pitchmcomb_t * p, cvec_t * fftgrain, fvec_t * output)
108 fvec_t *newmag = (fvec_t *) p->newmag;
109 //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
110 /* copy incoming grain to newmag */
111 for (i = 0; i < fftgrain->channels; i++) {
112 for (j = 0; j < newmag->length; j++)
113 newmag->data[0][j] = fftgrain->norm[i][j];
114 /* detect only if local energy > 10. */
115 //if (fvec_local_energy(newmag)>10.) {
116 //hfc = fvec_local_hfc(newmag); //not used
117 aubio_pitchmcomb_spectral_pp (p, newmag);
118 aubio_pitchmcomb_combdet (p, newmag);
119 //aubio_pitchmcomb_sort_cand_freq(p->candidates,p->ncand);
120 //return p->candidates[p->goodcandidate]->ebin;
121 j = (uint_t) FLOOR (p->candidates[p->goodcandidate]->ebin + .5);
122 instfreq = aubio_unwrap2pi (fftgrain->phas[i][j]
123 - p->theta->data[i][j] - j * p->phasediff);
124 instfreq *= p->phasefreq;
125 /* store phase for next run */
126 for (j = 0; j < p->theta->length; j++) {
127 p->theta->data[i][j] = fftgrain->phas[i][j];
129 //return p->candidates[p->goodcandidate]->ebin;
131 FLOOR (p->candidates[p->goodcandidate]->ebin + .5) + instfreq;
139 aubio_pitch_cands (aubio_pitchmcomb_t * p, cvec_t * fftgrain, smpl_t * cands)
143 fvec_t *newmag = (fvec_t *) p->newmag;
144 aubio_spectralcandidate_t **scands =
145 (aubio_spectralcandidate_t **) (p->candidates);
146 //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
147 /* copy incoming grain to newmag */
148 for (j = 0; j < newmag->length; j++)
149 newmag->data[i][j] = fftgrain->norm[i][j];
150 /* detect only if local energy > 10. */
151 if (fvec_local_energy (newmag) > 10.) {
152 /* hfc = fvec_local_hfc(newmag); do not use */
153 aubio_pitchmcomb_spectral_pp (p, newmag);
154 aubio_pitchmcomb_combdet (p, newmag);
155 aubio_pitchmcomb_sort_cand_freq (scands, p->ncand);
156 /* store ncand comb energies in cands[1:ncand] */
157 for (k = 0; k < p->ncand; k++)
158 cands[k] = p->candidates[k]->ene;
159 /* store ncand[end] freq in cands[end] */
160 cands[p->ncand] = p->candidates[p->ncand - 1]->ebin;
163 for (k = 0; k < p->ncand; k++)
170 aubio_pitchmcomb_spectral_pp (aubio_pitchmcomb_t * p, fvec_t * newmag)
172 fvec_t *mag = (fvec_t *) p->scratch;
173 fvec_t *tmp = (fvec_t *) p->scratch2;
175 uint_t length = mag->length;
176 /* copy newmag to mag (scracth) */
177 for (j = 0; j < length; j++) {
178 mag->data[i][j] = newmag->data[i][j];
180 fvec_min_removal (mag); /* min removal */
181 fvec_alpha_normalise (mag, p->alpha); /* alpha normalisation */
182 /* skipped *//* low pass filtering */
183 /** \bug fvec_moving_thres may write out of bounds */
184 fvec_adapt_thres (mag, tmp, p->win_post, p->win_pre, i); /* adaptative threshold */
185 fvec_add (mag, -p->threshold); /* fixed threshold */
187 aubio_spectralpeak_t *peaks = (aubio_spectralpeak_t *) p->peaks;
189 /* return bin and ebin */
190 count = aubio_pitchmcomb_quadpick (peaks, mag);
191 for (j = 0; j < count; j++)
192 peaks[j].mag = newmag->data[i][peaks[j].bin];
193 /* reset non peaks */
194 for (j = count; j < length; j++)
202 aubio_pitchmcomb_combdet (aubio_pitchmcomb_t * p, fvec_t * newmag)
204 aubio_spectralpeak_t *peaks = (aubio_spectralpeak_t *) p->peaks;
205 aubio_spectralcandidate_t **candidate =
206 (aubio_spectralcandidate_t **) p->candidates;
209 uint_t N = p->npartials; /* maximum number of partials to be considered 10 */
210 uint_t M = p->ncand; /* maximum number of combs to be considered 5 */
211 uint_t length = newmag->length;
212 uint_t count = p->count;
222 uint_t root_peak = 0;
226 /* get the biggest peak in the spectrum */
227 root_peak = aubio_pitchmcomb_get_root_peak (peaks, count);
228 /* not enough partials in highest notes, could be forced */
229 //if (peaks[root_peak].ebin >= aubio_miditofreq(85.)/p->tau) N=2;
230 //if (peaks[root_peak].ebin >= aubio_miditofreq(90.)/p->tau) N=1;
231 /* now calculate the energy of each of the 5 combs */
232 for (l = 0; l < M; l++) {
233 smpl_t scaler = (1. / (l + 1.));
234 candidate[l]->ene = 0.; /* reset ene and len sums */
235 candidate[l]->len = 0.;
236 candidate[l]->ebin = scaler * peaks[root_peak].ebin;
237 /* if less than N peaks available, curlen < N */
238 if (candidate[l]->ebin != 0.)
239 curlen = (uint_t) FLOOR (length / (candidate[l]->ebin));
240 curlen = (N < curlen) ? N : curlen;
241 /* fill candidate[l]->ecomb[k] with (k+1)*candidate[l]->ebin */
242 for (k = 0; k < curlen; k++)
243 candidate[l]->ecomb[k] = (candidate[l]->ebin) * (k + 1.);
244 for (k = curlen; k < length; k++)
245 candidate[l]->ecomb[k] = 0.;
246 /* for each in candidate[l]->ecomb[k] */
247 for (k = 0; k < curlen; k++) {
249 /** get the candidate->ecomb the closer to peaks.ebin
250 * (to cope with the inharmonicity)*/
251 for (d = 0; d < count; d++) {
252 delta2 = ABS (candidate[l]->ecomb[k] - peaks[d].ebin);
258 /* for a Q factor of 17, maintaining "constant Q filtering",
259 * and sum energy and length over non null combs */
260 if (17. * xx < candidate[l]->ecomb[k]) {
261 candidate[l]->ecomb[k] = peaks[position].ebin;
262 candidate[l]->ene += /* ecomb rounded to nearest int */
263 POW (newmag->data[0][(uint_t) FLOOR (candidate[l]->ecomb[k] + .5)],
265 candidate[l]->len += 1. / curlen;
267 candidate[l]->ecomb[k] = 0.;
270 /*if (candidate[l]->len<0.6)
271 candidate[l]->ene=0.; */
272 /* remember best candidate energy (in polyphonic, could check for
273 * tmpene*1.1 < candidate->ene to reduce jumps towards low frequencies) */
274 if (tmpene < candidate[l]->ene) {
276 tmpene = candidate[l]->ene;
279 //p->candidates=candidate;
281 p->goodcandidate = tmpl;
284 /** T=quadpick(X): return indices of elements of X which are peaks and positive
285 * exact peak positions are retrieved by quadratic interpolation
287 * \bug peak-picking too picky, sometimes counts too many peaks ?
290 aubio_pitchmcomb_quadpick (aubio_spectralpeak_t * spectral_peaks, fvec_t * X)
292 uint_t i, j, ispeak, count = 0;
293 for (i = 0; i < X->channels; i++)
294 for (j = 1; j < X->length - 1; j++) {
295 ispeak = fvec_peakpick (X, j);
298 spectral_peaks[count - 1].bin = j;
299 spectral_peaks[count - 1].ebin = fvec_quadint (X, j, i) - 1.;
305 /* get predominant partial */
307 aubio_pitchmcomb_get_root_peak (aubio_spectralpeak_t * peaks, uint_t length)
311 for (i = 0; i < length; i++)
312 if (tmp <= peaks[i].mag) {
320 aubio_pitchmcomb_sort_peak (aubio_spectralpeak_t * peaks, uint_t nbins)
322 qsort (peaks, nbins, sizeof (aubio_spectralpeak_t),
323 aubio_pitchmcomb_sort_peak_comp);
327 aubio_pitchmcomb_sort_peak_comp (const void *x, const void *y)
329 return (((aubio_spectralpeak_t *) y)->mag -
330 ((aubio_spectralpeak_t *) x)->mag);
335 aubio_pitchmcomb_sort_cand_ene (aubio_spectralcandidate_t ** candidates,
340 for (cur = 0; cur < nbins; cur++) {
342 for (run = cur; run < nbins; run++) {
343 if (candidates[run]->ene > candidates[cur]->ene)
344 CAND_SWAP (candidates[run], candidates[cur]);
351 aubio_pitchmcomb_sort_cand_freq (aubio_spectralcandidate_t ** candidates,
356 for (cur = 0; cur < nbins; cur++) {
358 for (run = cur; run < nbins; run++) {
359 if (candidates[run]->ebin < candidates[cur]->ebin)
360 CAND_SWAP (candidates[run], candidates[cur]);
366 new_aubio_pitchmcomb (uint_t bufsize, uint_t hopsize, uint_t channels)
368 aubio_pitchmcomb_t *p = AUBIO_NEW (aubio_pitchmcomb_t);
369 /* bug: should check if size / 8 > post+pre+1 */
372 p->spec_partition = 4;
379 // p->tau = samplerate/bufsize;
381 p->goodcandidate = 0;
382 p->phasefreq = bufsize / hopsize / TWO_PI;
383 p->phasediff = TWO_PI * hopsize / bufsize;
384 spec_size = bufsize / p->spec_partition;
385 //p->pickerfn = quadpick;
386 //p->biquad = new_biquad(0.1600,0.3200,0.1600, -0.5949, 0.2348);
387 /* allocate temp memory */
388 p->newmag = new_fvec (spec_size, 1);
389 /* array for median */
390 p->scratch = new_fvec (spec_size, 1);
391 /* array for phase */
392 p->theta = new_fvec (spec_size, channels);
393 /* array for adaptative threshold */
394 p->scratch2 = new_fvec (p->win_post + p->win_pre + 1, 1);
395 /* array of spectral peaks */
396 p->peaks = AUBIO_ARRAY (aubio_spectralpeak_t, spec_size);
397 for (i = 0; i < spec_size; i++) {
398 p->peaks[i].bin = 0.;
399 p->peaks[i].ebin = 0.;
400 p->peaks[i].mag = 0.;
402 /* array of pointers to spectral candidates */
403 p->candidates = AUBIO_ARRAY (aubio_spectralcandidate_t *, p->ncand);
404 for (i = 0; i < p->ncand; i++) {
405 p->candidates[i] = AUBIO_NEW (aubio_spectralcandidate_t);
406 p->candidates[i]->ecomb = AUBIO_ARRAY (smpl_t, spec_size);
407 for (j = 0; j < spec_size; j++) {
408 p->candidates[i]->ecomb[j] = 0.;
410 p->candidates[i]->ene = 0.;
411 p->candidates[i]->ebin = 0.;
412 p->candidates[i]->len = 0.;
419 del_aubio_pitchmcomb (aubio_pitchmcomb_t * p)
422 del_fvec (p->newmag);
423 del_fvec (p->scratch);
425 del_fvec (p->scratch2);
426 AUBIO_FREE (p->peaks);
427 for (i = 0; i < p->ncand; i++) {
428 AUBIO_FREE (p->candidates[i]->ecomb);
429 AUBIO_FREE (p->candidates[i]);
431 AUBIO_FREE (p->candidates);