src/pitch/pitchmcomb.c

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