2 Copyright (C) 2003 Paul Brossier
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20 #include "aubio_priv.h"
22 #include "mathutils.h"
23 #include "pitchmcomb.h"
25 #define CAND_SWAP(a,b) { register aubio_spectralcandidate_t *t=(a);(a)=(b);(b)=t; }
27 typedef struct _aubio_spectralpeak_t aubio_spectralpeak_t;
28 typedef struct _aubio_spectralcandidate_t aubio_spectralcandidate_t;
29 uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length);
30 uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X);
31 void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * oldmag);
32 void aubio_pitchmcomb_combdet(aubio_pitchmcomb_t * p, fvec_t * newmag);
33 /* not used but useful : sort by amplitudes (or anything else)
34 * sort_pitchpeak(peaks, length);
36 /** spectral_peak comparison function (must return signed int) */
37 static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y);
38 /** sort spectral_peak against their mag */
39 void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins);
41 /** sort spectral_candidate against their comb ene */
42 void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins);
43 /** sort spectral_candidate against their frequency */
44 void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins);
46 struct _aubio_pitchmcomb_t {
47 smpl_t threshold; /**< offset threshold [0.033 or 0.01] */
48 smpl_t alpha; /**< normalisation exponent [9] */
49 smpl_t cutoff; /**< low-pass filter cutoff [0.34, 1] */
50 smpl_t tol; /**< tolerance [0.05] */
51 smpl_t tau; /**< frequency precision [44100/4096] */
52 uint_t win_post; /**< median filter window length */
53 uint_t win_pre; /**< median filter window */
54 uint_t ncand; /**< maximum number of candidates (combs) */
55 uint_t npartials; /**< maximum number of partials per combs */
56 uint_t count; /**< picked picks */
57 uint_t goodcandidate; /**< best candidate */
58 uint_t spec_partition; /**< spectrum partition to consider */
59 aubio_spectralpeak_t * peaks; /**< up to length win/spec_partition */
60 aubio_spectralcandidate_t ** candidates; /** up to five candidates */
61 /* some scratch pads */
62 /** \bug (unnecessary copied from fftgrain?) */
63 fvec_t * newmag; /**< vec to store mag */
64 fvec_t * scratch; /**< vec to store modified mag */
65 fvec_t * scratch2; /**< vec to compute moving median */
66 fvec_t * theta; /**< vec to store phase */
69 /** threshfn: name or handle of fn for computing adaptive threshold [median] */
70 /** aubio_thresholdfn_t thresholdfn; */
71 /** picker: name or handle of fn for picking event times [quadpick] */
72 /** aubio_pickerfn_t pickerfn; */
75 /** spectral peak object */
76 struct _aubio_spectralpeak_t {
77 uint_t bin; /**< bin [0-(length-1)] */
78 smpl_t ebin; /**< estimated bin */
79 smpl_t mag; /**< peak magnitude */
82 /** spectral candidates array object */
83 struct _aubio_spectralcandidate_t {
84 smpl_t ebin; /**< interpolated bin */
85 smpl_t * ecomb; /**< comb */
86 smpl_t ene; /**< candidate energy */
87 smpl_t len; /**< length */
91 smpl_t aubio_pitchmcomb_detect(aubio_pitchmcomb_t * p, cvec_t * fftgrain) {
94 fvec_t * newmag = (fvec_t *)p->newmag;
95 //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
96 /* copy incoming grain to newmag */
97 for (j=0; j< newmag->length; j++)
98 newmag->data[i][j]=fftgrain->norm[i][j];
99 /* detect only if local energy > 10. */
100 //if (vec_local_energy(newmag)>10.) {
101 //hfc = vec_local_hfc(newmag); //not used
102 aubio_pitchmcomb_spectral_pp(p, newmag);
103 aubio_pitchmcomb_combdet(p,newmag);
104 //aubio_pitchmcomb_sort_cand_freq(p->candidates,p->ncand);
105 //return p->candidates[p->goodcandidate]->ebin;
106 j = (uint_t)FLOOR(p->candidates[p->goodcandidate]->ebin+.5);
107 instfreq = aubio_unwrap2pi(fftgrain->phas[0][j]
108 - p->theta->data[0][j] - j*p->phasediff);
109 instfreq *= p->phasefreq;
110 /* store phase for next run */
111 for (j=0; j< p->theta->length; j++) {
112 p->theta->data[i][j]=fftgrain->phas[i][j];
114 //return p->candidates[p->goodcandidate]->ebin;
115 return FLOOR(p->candidates[p->goodcandidate]->ebin+.5) + instfreq;
121 uint_t aubio_pitch_cands(aubio_pitchmcomb_t * p, cvec_t * fftgrain,
125 fvec_t * newmag = (fvec_t *)p->newmag;
126 aubio_spectralcandidate_t ** scands =
127 (aubio_spectralcandidate_t **)(p->candidates);
128 //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
129 /* copy incoming grain to newmag */
130 for (j=0; j< newmag->length; j++)
131 newmag->data[i][j]=fftgrain->norm[i][j];
132 /* detect only if local energy > 10. */
133 if (vec_local_energy(newmag)>10.) {
134 /* hfc = vec_local_hfc(newmag); do not use */
135 aubio_pitchmcomb_spectral_pp(p, newmag);
136 aubio_pitchmcomb_combdet(p,newmag);
137 aubio_pitchmcomb_sort_cand_freq(scands,p->ncand);
138 /* store ncand comb energies in cands[1:ncand] */
139 for (k = 0; k<p->ncand; k++)
140 cands[k] = p->candidates[k]->ene;
141 /* store ncand[end] freq in cands[end] */
142 cands[p->ncand] = p->candidates[p->ncand-1]->ebin;
145 for (k = 0; k<p->ncand; k++)
151 void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * newmag) {
152 fvec_t * mag = (fvec_t *)p->scratch;
153 fvec_t * tmp = (fvec_t *)p->scratch2;
155 uint_t length = mag->length;
156 /* copy newmag to mag (scracth) */
157 for (j=0;j<length;j++) {
158 mag->data[i][j] = newmag->data[i][j];
160 vec_dc_removal(mag); /* dc removal */
161 vec_alpha_normalise(mag,p->alpha); /* alpha normalisation */
162 /* skipped */ /* low pass filtering */
163 /** \bug vec_moving_thres may write out of bounds */
164 vec_adapt_thres(mag,tmp,p->win_post,p->win_pre); /* adaptative threshold */
165 vec_add(mag,-p->threshold); /* fixed threshold */
167 aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
169 /* return bin and ebin */
170 count = aubio_pitchmcomb_quadpick(peaks,mag);
171 for (j=0;j<count;j++)
172 peaks[j].mag = newmag->data[i][peaks[j].bin];
173 /* reset non peaks */
174 for (j=count;j<length;j++)
181 void aubio_pitchmcomb_combdet(aubio_pitchmcomb_t * p, fvec_t * newmag) {
182 aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
183 aubio_spectralcandidate_t ** candidate =
184 (aubio_spectralcandidate_t **)p->candidates;
187 uint_t N = p->npartials; /* maximum number of partials to be considered 10 */
188 uint_t M = p->ncand; /* maximum number of combs to be considered 5 */
189 uint_t length = newmag->length;
190 uint_t count = p->count;
200 uint_t root_peak = 0;
204 /* get the biggest peak in the spectrum */
205 root_peak = aubio_pitchmcomb_get_root_peak(peaks,count);
206 /* not enough partials in highest notes, could be forced */
207 //if (peaks[root_peak].ebin >= aubio_miditofreq(85.)/p->tau) N=2;
208 //if (peaks[root_peak].ebin >= aubio_miditofreq(90.)/p->tau) N=1;
209 /* now calculate the energy of each of the 5 combs */
211 smpl_t scaler = (1./(l+1.));
212 candidate[l]->ene = 0.; /* reset ene and len sums */
213 candidate[l]->len = 0.;
214 candidate[l]->ebin=scaler*peaks[root_peak].ebin;
215 /* if less than N peaks available, curlen < N */
216 curlen = (uint_t)FLOOR(length/(candidate[l]->ebin));
217 curlen = (N < curlen )? N : curlen;
218 /* fill candidate[l]->ecomb[k] with (k+1)*candidate[l]->ebin */
219 for (k=0;k<curlen;k++)
220 candidate[l]->ecomb[k]=(candidate[l]->ebin)*(k+1.);
221 for (k=curlen;k<length;k++)
222 candidate[l]->ecomb[k]=0.;
223 /* for each in candidate[l]->ecomb[k] */
224 for (k=0;k<curlen;k++) {
226 /** get the candidate->ecomb the closer to peaks.ebin
227 * (to cope with the inharmonicity)*/
228 for (d=0;d<count;d++) {
229 delta2 = ABS(candidate[l]->ecomb[k]-peaks[d].ebin);
235 /* for a Q factor of 17, maintaining "constant Q filtering",
236 * and sum energy and length over non null combs */
237 if ( 17. * xx < candidate[l]->ecomb[k] ) {
238 candidate[l]->ecomb[k]=peaks[position].ebin;
239 candidate[l]->ene += /* ecomb rounded to nearest int */
240 POW(newmag->data[0][(uint_t)FLOOR(candidate[l]->ecomb[k]+.5)],0.25);
241 candidate[l]->len += 1./curlen;
243 candidate[l]->ecomb[k]=0.;
246 /*if (candidate[l]->len<0.6)
247 candidate[l]->ene=0.; */
248 /* remember best candidate energy (in polyphonic, could check for
249 * tmpene*1.1 < candidate->ene to reduce jumps towards low frequencies) */
250 if (tmpene < candidate[l]->ene) {
252 tmpene = candidate[l]->ene;
255 //p->candidates=candidate;
257 p->goodcandidate = tmpl;
260 /** T=quadpick(X): return indices of elements of X which are peaks and positive
261 * exact peak positions are retrieved by quadratic interpolation
263 * \bug peak-picking too picky, sometimes counts too many peaks ?
265 uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X){
266 uint_t i, j, ispeak, count = 0;
267 for (i=0;i<X->channels;i++)
268 for (j=1;j<X->length-1;j++) {
269 ispeak = vec_peakpick(X,j);
272 spectral_peaks[count-1].bin = j;
273 spectral_peaks[count-1].ebin = vec_quadint(X,j) - 1.;
279 /* get predominant partial */
280 uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length) {
283 for (i=0;i<length;i++)
284 if (tmp <= peaks[i].mag) {
291 void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins) {
292 qsort(peaks, nbins, sizeof(aubio_spectralpeak_t),
293 aubio_pitchmcomb_sort_peak_comp);
295 static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y) {
296 return (((aubio_spectralpeak_t *)y)->mag - ((aubio_spectralpeak_t *)x)->mag);
300 void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
303 for (cur=0;cur<nbins;cur++) {
305 for (run=cur;run<nbins;run++) {
306 if(candidates[run]->ene > candidates[cur]->ene)
307 CAND_SWAP(candidates[run], candidates[cur]);
313 void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
316 for (cur=0;cur<nbins;cur++) {
318 for (run=cur;run<nbins;run++) {
319 if(candidates[run]->ebin < candidates[cur]->ebin)
320 CAND_SWAP(candidates[run], candidates[cur]);
325 aubio_pitchmcomb_t * new_aubio_pitchmcomb(uint_t bufsize, uint_t hopsize, uint_t channels, uint_t samplerate) {
326 aubio_pitchmcomb_t * p = AUBIO_NEW(aubio_pitchmcomb_t);
327 /* bug: should check if size / 8 > post+pre+1 */
330 p->spec_partition = 4;
337 p->tau = samplerate/bufsize;
339 p->goodcandidate = 0;
340 p->phasefreq = bufsize/hopsize/TWO_PI;
341 p->phasediff = TWO_PI*hopsize/bufsize;
342 spec_size = bufsize/p->spec_partition;
343 //p->pickerfn = quadpick;
344 //p->biquad = new_biquad(0.1600,0.3200,0.1600, -0.5949, 0.2348);
345 /* allocate temp memory */
346 p->newmag = new_fvec(spec_size,channels);
347 /* array for median */
348 p->scratch = new_fvec(spec_size,channels);
349 /* array for phase */
350 p->theta = new_fvec(spec_size,channels);
351 /* array for adaptative threshold */
352 p->scratch2 = new_fvec(p->win_post+p->win_pre+1,channels);
353 /* array of spectral peaks */
354 p->peaks = AUBIO_ARRAY(aubio_spectralpeak_t,spec_size);
355 /* array of pointers to spectral candidates */
356 p->candidates = AUBIO_ARRAY(aubio_spectralcandidate_t *,p->ncand);
357 for (i=0;i<p->ncand;i++) {
358 p->candidates[i] = AUBIO_NEW(aubio_spectralcandidate_t);
359 p->candidates[i]->ecomb = AUBIO_ARRAY(smpl_t, spec_size);
365 void del_aubio_pitchmcomb (aubio_pitchmcomb_t *p) {
368 del_fvec(p->scratch);
369 del_fvec(p->scratch2);
370 AUBIO_FREE(p->peaks);
371 for (i=0;i<p->ncand;i++) {
372 AUBIO_FREE(p->candidates[i]);
374 AUBIO_FREE(p->candidates);