source: src/pitch/pitchmcomb.c @ e6a78ea

feature/autosinkfeature/cnnfeature/cnn_orgfeature/constantqfeature/crepefeature/crepe_orgfeature/pitchshiftfeature/pydocstringsfeature/timestretchfix/ffmpeg5pitchshiftsamplertimestretchyinfft+
Last change on this file since e6a78ea was e6a78ea, checked in by Paul Brossier <piem@piem.org>, 14 years ago

src: update all headers to GPLv3

  • Property mode set to 100644
File size: 14.4 KB
Line 
1/*
2  Copyright (C) 2003-2009 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#include "aubio_priv.h"
22#include "fvec.h"
23#include "cvec.h"
24#include "mathutils.h"
25#include "pitch/pitchmcomb.h"
26
27#define CAND_SWAP(a,b) { register aubio_spectralcandidate_t *t=(a);(a)=(b);(b)=t; }
28
29typedef struct _aubio_spectralpeak_t aubio_spectralpeak_t;
30typedef struct _aubio_spectralcandidate_t aubio_spectralcandidate_t;
31uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length);
32uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X);
33void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * oldmag);
34void aubio_pitchmcomb_combdet(aubio_pitchmcomb_t * p, fvec_t * newmag);
35/* not used but useful : sort by amplitudes (or anything else)
36 * sort_pitchpeak(peaks, length);
37 */
38/** spectral_peak comparison function (must return signed int) */
39static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y);
40/** sort spectral_peak against their mag */
41void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins);
42/** select the best candidates */
43uint_t aubio_pitch_cands(aubio_pitchmcomb_t * p, cvec_t * fftgrain, smpl_t * cands);
44
45/** sort spectral_candidate against their comb ene */
46void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins);
47/** sort spectral_candidate against their frequency */
48void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins);
49
50struct _aubio_pitchmcomb_t {
51  smpl_t threshold;                        /**< offset threshold [0.033 or 0.01]     */
52  smpl_t alpha;                            /**< normalisation exponent [9]           */
53  smpl_t cutoff;                           /**< low-pass filter cutoff [0.34, 1]     */
54  smpl_t tol;                              /**< tolerance [0.05]                     */
55  // smpl_t tau;                              /**< frequency precision [44100/4096]     */
56  uint_t win_post;                         /**< median filter window length          */
57  uint_t win_pre;                          /**< median filter window                 */
58  uint_t ncand;                            /**< maximum number of candidates (combs) */
59  uint_t npartials;                        /**< maximum number of partials per combs */
60  uint_t count;                            /**< picked picks                         */
61  uint_t goodcandidate;                    /**< best candidate                       */
62  uint_t spec_partition;                   /**< spectrum partition to consider       */
63  aubio_spectralpeak_t * peaks;            /**< up to length win/spec_partition      */
64  aubio_spectralcandidate_t ** candidates; /** up to five candidates                 */
65  /* some scratch pads */
66  /** \bug  (unnecessary copied from fftgrain?) */
67  fvec_t * newmag;                         /**< vec to store mag                     */
68  fvec_t * scratch;                        /**< vec to store modified mag            */
69  fvec_t * scratch2;                       /**< vec to compute moving median         */
70  fvec_t * theta;                         /**< vec to store phase                     */
71  smpl_t phasediff;
72  smpl_t phasefreq;
73  /** threshfn: name or handle of fn for computing adaptive threshold [median] */
74  /** aubio_thresholdfn_t thresholdfn; */
75  /** picker: name or handle of fn for picking event times [quadpick] */
76  /** aubio_pickerfn_t pickerfn; */
77};
78
79/** spectral peak object */
80struct _aubio_spectralpeak_t {
81  uint_t bin;     /**< bin [0-(length-1)] */
82  smpl_t ebin;    /**< estimated bin */
83  smpl_t mag;     /**< peak magnitude */
84};
85
86/** spectral candidates array object */
87struct _aubio_spectralcandidate_t {
88  smpl_t ebin;    /**< interpolated bin */
89  smpl_t * ecomb; /**< comb */
90  smpl_t ene;     /**< candidate energy */
91  smpl_t len;     /**< length */
92};
93
94
95void aubio_pitchmcomb_do (aubio_pitchmcomb_t * p, cvec_t * fftgrain, fvec_t * output) {
96  uint_t i,j;
97  smpl_t instfreq;
98  fvec_t * newmag = (fvec_t *)p->newmag;
99  //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
100  /* copy incoming grain to newmag */
101  for (i=0; i< fftgrain->channels; i++) {
102  for (j=0; j< newmag->length; j++)
103    newmag->data[0][j]=fftgrain->norm[i][j];
104  /* detect only if local energy > 10. */
105  //if (fvec_local_energy(newmag)>10.) {
106    //hfc = fvec_local_hfc(newmag); //not used
107    aubio_pitchmcomb_spectral_pp(p, newmag);
108    aubio_pitchmcomb_combdet(p,newmag);
109    //aubio_pitchmcomb_sort_cand_freq(p->candidates,p->ncand);
110    //return p->candidates[p->goodcandidate]->ebin;
111  j = (uint_t)FLOOR(p->candidates[p->goodcandidate]->ebin+.5);
112  instfreq  = aubio_unwrap2pi(fftgrain->phas[i][j]
113      - p->theta->data[i][j] - j*p->phasediff);
114  instfreq *= p->phasefreq;
115  /* store phase for next run */
116  for (j=0; j< p->theta->length; j++) {
117    p->theta->data[i][j]=fftgrain->phas[i][j];
118  }
119  //return p->candidates[p->goodcandidate]->ebin;
120  output->data[i][0] = FLOOR(p->candidates[p->goodcandidate]->ebin+.5) + instfreq;
121  /*} else {
122    return -1.;
123  }*/
124  }
125}
126
127uint_t aubio_pitch_cands(aubio_pitchmcomb_t * p, cvec_t * fftgrain,
128    smpl_t * cands) {
129  uint_t i=0,j;
130  uint_t k;
131  fvec_t * newmag = (fvec_t *)p->newmag;
132  aubio_spectralcandidate_t ** scands =
133    (aubio_spectralcandidate_t **)(p->candidates);
134  //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
135  /* copy incoming grain to newmag */
136  for (j=0; j< newmag->length; j++)
137    newmag->data[i][j]=fftgrain->norm[i][j];
138  /* detect only if local energy > 10. */
139  if (fvec_local_energy(newmag)>10.) {
140    /* hfc = fvec_local_hfc(newmag); do not use */
141    aubio_pitchmcomb_spectral_pp(p, newmag);
142    aubio_pitchmcomb_combdet(p,newmag);
143    aubio_pitchmcomb_sort_cand_freq(scands,p->ncand);
144    /* store ncand comb energies in cands[1:ncand] */
145    for (k = 0; k<p->ncand; k++)
146      cands[k] = p->candidates[k]->ene;
147    /* store ncand[end] freq in cands[end] */
148    cands[p->ncand] = p->candidates[p->ncand-1]->ebin;
149    return 1;
150  } else {
151    for (k = 0; k<p->ncand; k++)
152      cands[k] = 0;
153    return 0;
154  }
155}
156
157void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * newmag) {
158  fvec_t * mag = (fvec_t *)p->scratch;
159  fvec_t * tmp = (fvec_t *)p->scratch2;
160  uint_t i=0,j;
161  uint_t length = mag->length;
162  /* copy newmag to mag (scracth) */
163  for (j=0;j<length;j++) {
164    mag->data[i][j] = newmag->data[i][j];
165  }
166  fvec_min_removal(mag);              /* min removal          */
167  fvec_alpha_normalise(mag,p->alpha); /* alpha normalisation  */
168  /* skipped */                       /* low pass filtering   */
169  /** \bug fvec_moving_thres may write out of bounds */
170  fvec_adapt_thres(mag,tmp,p->win_post,p->win_pre,i); /* adaptative threshold */
171  fvec_add(mag,-p->threshold);        /* fixed threshold      */
172  {
173    aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
174    uint_t count;
175    /*  return bin and ebin */
176    count = aubio_pitchmcomb_quadpick(peaks,mag);
177    for (j=0;j<count;j++)
178      peaks[j].mag = newmag->data[i][peaks[j].bin];
179    /* reset non peaks */
180    for (j=count;j<length;j++)
181      peaks[j].mag = 0.;
182    p->peaks = peaks;
183    p->count = count;
184  }
185}
186
187void aubio_pitchmcomb_combdet(aubio_pitchmcomb_t * p, fvec_t * newmag) {
188  aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
189  aubio_spectralcandidate_t ** candidate =
190    (aubio_spectralcandidate_t **)p->candidates;
191
192  /* parms */
193  uint_t N = p->npartials; /* maximum number of partials to be considered 10 */
194  uint_t M = p->ncand;  /* maximum number of combs to be considered 5 */
195  uint_t length = newmag->length;
196  uint_t count = p->count;
197  uint_t k;
198  uint_t l;
199  uint_t d;
200  uint_t curlen = 0;
201
202  smpl_t delta2;
203  smpl_t xx;
204  uint_t position = 0;
205
206  uint_t root_peak = 0;
207  uint_t tmpl = 0;
208  smpl_t tmpene = 0.;
209
210  /* get the biggest peak in the spectrum */
211  root_peak = aubio_pitchmcomb_get_root_peak(peaks,count);
212  /* not enough partials in highest notes, could be forced */
213  //if (peaks[root_peak].ebin >= aubio_miditofreq(85.)/p->tau) N=2;
214  //if (peaks[root_peak].ebin >= aubio_miditofreq(90.)/p->tau) N=1;
215  /* now calculate the energy of each of the 5 combs */
216  for (l=0;l<M;l++) {
217    smpl_t scaler = (1./(l+1.));
218    candidate[l]->ene = 0.; /* reset ene and len sums */
219    candidate[l]->len = 0.;
220    candidate[l]->ebin=scaler*peaks[root_peak].ebin;
221    /* if less than N peaks available, curlen < N */
222    if (candidate[l]->ebin != 0.)
223      curlen = (uint_t)FLOOR(length/(candidate[l]->ebin));
224    curlen = (N < curlen )? N : curlen;
225    /* fill candidate[l]->ecomb[k] with (k+1)*candidate[l]->ebin */
226    for (k=0;k<curlen;k++)
227      candidate[l]->ecomb[k]=(candidate[l]->ebin)*(k+1.);
228    for (k=curlen;k<length;k++)
229      candidate[l]->ecomb[k]=0.;
230    /* for each in candidate[l]->ecomb[k] */
231    for (k=0;k<curlen;k++) {
232      xx = 100000.;
233      /** get the candidate->ecomb the closer to peaks.ebin
234       * (to cope with the inharmonicity)*/
235      for (d=0;d<count;d++) {
236        delta2 = ABS(candidate[l]->ecomb[k]-peaks[d].ebin);
237        if (delta2 <= xx) {
238          position = d;
239          xx = delta2;
240        }
241      }
242      /* for a Q factor of 17, maintaining "constant Q filtering",
243       * and sum energy and length over non null combs */
244      if ( 17. * xx < candidate[l]->ecomb[k] ) {
245        candidate[l]->ecomb[k]=peaks[position].ebin;
246        candidate[l]->ene += /* ecomb rounded to nearest int */
247          POW(newmag->data[0][(uint_t)FLOOR(candidate[l]->ecomb[k]+.5)],0.25);
248        candidate[l]->len += 1./curlen;
249      } else
250        candidate[l]->ecomb[k]=0.;
251    }
252    /* punishment */
253    /*if (candidate[l]->len<0.6)
254      candidate[l]->ene=0.; */
255    /* remember best candidate energy (in polyphonic, could check for
256     * tmpene*1.1 < candidate->ene to reduce jumps towards low frequencies) */
257    if (tmpene < candidate[l]->ene) {
258      tmpl = l;
259      tmpene = candidate[l]->ene;
260    }
261  }
262  //p->candidates=candidate;
263  //p->peaks=peaks;
264  p->goodcandidate = tmpl;
265}
266
267/** T=quadpick(X): return indices of elements of X which are peaks and positive
268 * exact peak positions are retrieved by quadratic interpolation
269 *
270 * \bug peak-picking too picky, sometimes counts too many peaks ?
271 */
272uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X){
273  uint_t i, j, ispeak, count = 0;
274  for (i=0;i<X->channels;i++)
275    for (j=1;j<X->length-1;j++) {
276      ispeak = fvec_peakpick(X,j);
277      if (ispeak) {
278        count += ispeak;
279        spectral_peaks[count-1].bin = j;
280        spectral_peaks[count-1].ebin = fvec_quadint(X, j, i) - 1.;
281      }
282    }
283  return count;
284}
285
286/* get predominant partial */
287uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length) {
288  uint_t i,pos=0;
289  smpl_t tmp = 0.;
290  for (i=0;i<length;i++)
291    if (tmp <= peaks[i].mag) {
292      pos = i;
293      tmp = peaks[i].mag;
294    }
295  return pos;
296}
297
298void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins) {
299  qsort(peaks, nbins, sizeof(aubio_spectralpeak_t),
300      aubio_pitchmcomb_sort_peak_comp);
301}
302static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y) {
303  return (((aubio_spectralpeak_t *)y)->mag - ((aubio_spectralpeak_t *)x)->mag);
304}
305
306
307void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
308  uint_t cur = 0;
309  uint_t run = 0;
310  for (cur=0;cur<nbins;cur++) {
311    run = cur + 1;
312    for (run=cur;run<nbins;run++) {
313      if(candidates[run]->ene > candidates[cur]->ene)
314        CAND_SWAP(candidates[run], candidates[cur]);
315    }
316  }
317}
318
319
320void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
321  uint_t cur = 0;
322  uint_t run = 0;
323  for (cur=0;cur<nbins;cur++) {
324    run = cur + 1;
325    for (run=cur;run<nbins;run++) {
326      if(candidates[run]->ebin < candidates[cur]->ebin)
327        CAND_SWAP(candidates[run], candidates[cur]);
328    }
329  }
330}
331
332aubio_pitchmcomb_t * new_aubio_pitchmcomb(uint_t bufsize, uint_t hopsize, uint_t channels) {
333  aubio_pitchmcomb_t * p = AUBIO_NEW(aubio_pitchmcomb_t);
334  /* bug: should check if size / 8 > post+pre+1 */
335  uint_t i, j;
336  uint_t spec_size;
337  p->spec_partition   = 4;
338  p->ncand            = 5;
339  p->npartials        = 5;
340  p->cutoff           = 1.;
341  p->threshold        = 0.01;
342  p->win_post         = 8;
343  p->win_pre          = 7;
344  // p->tau              = samplerate/bufsize;
345  p->alpha            = 9.;
346  p->goodcandidate    = 0;
347  p->phasefreq        = bufsize/hopsize/TWO_PI;
348  p->phasediff        = TWO_PI*hopsize/bufsize;
349  spec_size = bufsize/p->spec_partition;
350  //p->pickerfn = quadpick;
351  //p->biquad = new_biquad(0.1600,0.3200,0.1600, -0.5949, 0.2348);
352  /* allocate temp memory */
353  p->newmag     = new_fvec(spec_size,1);
354  /* array for median */
355  p->scratch    = new_fvec(spec_size,1);
356  /* array for phase */
357  p->theta      = new_fvec(spec_size,channels);
358  /* array for adaptative threshold */
359  p->scratch2   = new_fvec(p->win_post+p->win_pre+1,1);
360  /* array of spectral peaks */
361  p->peaks      = AUBIO_ARRAY(aubio_spectralpeak_t,spec_size);
362  for (i = 0; i < spec_size; i++) {
363    p->peaks[i].bin = 0.;
364    p->peaks[i].ebin = 0.;
365    p->peaks[i].mag = 0.;
366  }
367  /* array of pointers to spectral candidates */
368  p->candidates = AUBIO_ARRAY(aubio_spectralcandidate_t *,p->ncand);
369  for (i=0;i<p->ncand;i++) {
370    p->candidates[i] = AUBIO_NEW(aubio_spectralcandidate_t);
371    p->candidates[i]->ecomb = AUBIO_ARRAY(smpl_t, spec_size);
372    for (j=0; j < spec_size; j++) {
373      p->candidates[i]->ecomb[j] = 0.;
374    }
375    p->candidates[i]->ene = 0.;
376    p->candidates[i]->ebin = 0.;
377    p->candidates[i]->len = 0.;
378  }
379  return p;
380}
381
382
383void del_aubio_pitchmcomb (aubio_pitchmcomb_t *p) {
384  uint_t i;
385  del_fvec(p->newmag);
386  del_fvec(p->scratch);
387  del_fvec(p->theta);
388  del_fvec(p->scratch2);
389  AUBIO_FREE(p->peaks);
390  for (i=0;i<p->ncand;i++) {
391    AUBIO_FREE(p->candidates[i]->ecomb);
392    AUBIO_FREE(p->candidates[i]);
393  }
394  AUBIO_FREE(p->candidates);
395  AUBIO_FREE(p);
396}
Note: See TracBrowser for help on using the repository browser.