src/pitchfcomb.c

   1 /*
   2    Copyright (C) 2004, 2005  Mario Lang <mlang@delysid.org>
   3
   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.
   8
   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.
  13
  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.
  17
  18 */
  19
  20 /*
  21
  22    This file was taken from the tuneit project, in the file
  23    tuneit.c -- Detect fundamental frequency of a sound
  24    see http://delysid.org/tuneit.html
  25
  26    a fast harmonic comb filter algorithm for pitch tracking
  27
  28 */
  29
  30 #include "aubio_priv.h"
  31 #include "sample.h"
  32 #include "mathutils.h"
  33 #include "phasevoc.h"
  34 #include "pitchfcomb.h"
  35
  36 #define MAX_PEAKS 8
  37
  38 typedef struct {
  39   smpl_t freq;
  40   smpl_t db;
  41 } aubio_fpeak_t;
  42
  43 struct _aubio_pitchfcomb_t {
  44         uint_t fftSize;
  45         uint_t rate;
  46         cvec_t * fftOut;
  47         fvec_t * fftLastPhase;
  48         aubio_pvoc_t * pvoc;
  49 };
  50
  51 aubio_pitchfcomb_t * new_aubio_pitchfcomb (uint_t size, uint_t samplerate)
  52 {
  53   aubio_pitchfcomb_t * p = AUBIO_NEW(aubio_pitchfcomb_t);
  54   uint_t overlap_rate = 4;
  55   p->rate         = samplerate;
  56   p->fftSize      = size;
  57   p->fftOut       = new_cvec(size,1);
  58   p->fftLastPhase = new_fvec(size,1);
  59   p->pvoc = new_aubio_pvoc(size, size/overlap_rate, 1);
  60   return p;
  61 }
  62
  63 /* input must be stepsize long */
  64 smpl_t aubio_pitchfcomb_detect (aubio_pitchfcomb_t * p, fvec_t * input)
  65 {
  66   uint_t k, l, maxharm = 0, stepSize = input->length;
  67   smpl_t freqPerBin = p->rate/(smpl_t)p->fftSize,
  68     phaseDifference = TWO_PI*(smpl_t)stepSize/(smpl_t)p->fftSize;
  69   aubio_fpeak_t peaks[MAX_PEAKS];
  70
  71   for (k=0; k<MAX_PEAKS; k++) {
  72     peaks[k].db = -200.;
  73     peaks[k].freq = 0.;
  74   }
  75
  76   aubio_pvoc_do (p->pvoc, input, p->fftOut);
  77
  78   for (k=0; k<=p->fftSize; k++) {
  79     //long qpd;
  80     smpl_t
  81       magnitude = 20.*LOG10(2.*p->fftOut->norm[0][k]/(smpl_t)p->fftSize),
  82       phase     = p->fftOut->phas[0][k],
  83       tmp, freq;
  84
  85     /* compute phase difference */
  86     tmp = phase - p->fftLastPhase->data[0][k];
  87     p->fftLastPhase->data[0][k] = phase;
  88
  89     /* subtract expected phase difference */
  90     tmp -= (smpl_t)k*phaseDifference;
  91
  92     /* map delta phase into +/- Pi interval */
  93     tmp = aubio_unwrap2pi(tmp);
  94
  95     /* get deviation from bin frequency from the +/- Pi interval */
  96     tmp = p->fftSize/input->length*tmp/(TWO_PI);
  97
  98     /* compute the k-th partials' true frequency */
  99     freq = (smpl_t)k*freqPerBin + tmp*freqPerBin;
 100
 101     if (freq > 0.0 && magnitude > peaks[0].db && magnitude < 0) {
 102       memmove(peaks+1, peaks, sizeof(aubio_fpeak_t)*(MAX_PEAKS-1));
 103       peaks[0].freq = freq;
 104       peaks[0].db = magnitude;
 105     }
 106   }
 107
 108   k = 0;
 109   for (l=1; l<MAX_PEAKS && peaks[l].freq > 0.0; l++) {
 110     sint_t harmonic;
 111     for (harmonic=5; harmonic>1; harmonic--) {
 112       if (peaks[0].freq / peaks[l].freq < harmonic+.02 &&
 113         peaks[0].freq / peaks[l].freq > harmonic-.02) {
 114         if (harmonic > maxharm &&
 115           peaks[0].db < peaks[l].db/2) {
 116           maxharm = harmonic;
 117           k = l;
 118         }
 119       }
 120     }
 121   }
 122   /* quick hack to clean output a bit */
 123   if (peaks[k].freq > 10000) return 0.;
 124   return peaks[k].freq;
 125 }
 126
 127 void del_aubio_pitchfcomb (aubio_pitchfcomb_t * p)
 128 {
 129   del_cvec(p->fftOut);
 130   del_fvec(p->fftLastPhase);
 131   del_aubio_pvoc(p->pvoc);
 132   AUBIO_FREE(p);
 133 }
 134