source: src/pitch/pitchmcomb.c @ 9ee1ee9

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

src/pitch/pitchmcomb.c: fix candidates sorting function, really comparing current to next

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