source: src/pitchmcomb.c @ d3efa4e

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

fixed memory freeing of pitchmcomb

  • Property mode set to 100644
File size: 12.8 KB
Line 
1/*
2   Copyright (C) 2003 Paul Brossier
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#include "aubio_priv.h"
21#include "sample.h"
22#include "mathutils.h"
23#include "pitchmcomb.h"
24
25#define CAND_SWAP(a,b) { register aubio_spectralcandidate_t *t=(a);(a)=(b);(b)=t; }
26
27typedef struct _aubio_spectralpeak_t aubio_spectralpeak_t;
28typedef struct _aubio_spectralcandidate_t aubio_spectralcandidate_t;
29uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length);
30uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X);
31void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * oldmag);
32void 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);
35 */
36/** spectral_peak comparison function (must return signed int) */
37static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y);
38/** sort spectral_peak against their mag */
39void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins);
40
41/** sort spectral_candidate against their comb ene */
42void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins);
43/** sort spectral_candidate against their frequency */
44void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins);
45
46struct _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  /** threshfn: name or handle of fn for computing adaptive threshold [median] */
67  /** aubio_thresholdfn_t thresholdfn; */
68  /** picker: name or handle of fn for picking event times [quadpick] */
69  /** aubio_pickerfn_t pickerfn; */
70};
71
72/** spectral peak object */
73struct _aubio_spectralpeak_t {
74  uint_t bin;     /**< bin [0-(length-1)] */
75  smpl_t ebin;    /**< estimated bin */
76  smpl_t mag;     /**< peak magnitude */
77};
78
79/** spectral candidates array object */
80struct _aubio_spectralcandidate_t {
81  smpl_t ebin;    /**< interpolated bin */
82  smpl_t * ecomb; /**< comb */
83  smpl_t ene;     /**< candidate energy */
84  smpl_t len;     /**< length */
85};
86
87
88smpl_t aubio_pitchmcomb_detect(aubio_pitchmcomb_t * p, cvec_t * fftgrain) {
89  uint_t i=0,j;
90  fvec_t * newmag = (fvec_t *)p->newmag;
91  //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
92  /* copy incoming grain to newmag */
93  for (j=0; j< newmag->length; j++)
94    newmag->data[i][j]=fftgrain->norm[i][j];
95  /* detect only if local energy > 10. */ 
96  if (vec_local_energy(newmag)>10.) {
97    //hfc = vec_local_hfc(newmag); //not used
98    aubio_pitchmcomb_spectral_pp(p, newmag);
99    aubio_pitchmcomb_combdet(p,newmag);
100    // print picked candidate
101    //AUBIO_DBG("%f\n",p->candidates[p->goodcandidate]->ebin);
102    //AUBIO_DBG("%f\n",p->candidates[p->goodcandidate]->ebin);
103    return p->candidates[p->goodcandidate]->ebin;
104  } else {
105    return -1.;
106  }
107}
108
109uint_t aubio_pitch_cands(aubio_pitchmcomb_t * p, cvec_t * fftgrain, 
110    smpl_t * cands) {
111  uint_t i=0,j;
112  uint_t k;
113  fvec_t * newmag = (fvec_t *)p->newmag;
114  aubio_spectralcandidate_t ** scands = 
115    (aubio_spectralcandidate_t **)(p->candidates);
116  //smpl_t hfc; //fe=instfreq(theta1,theta,ops); //theta1=theta;
117  /* copy incoming grain to newmag */
118  for (j=0; j< newmag->length; j++)
119    newmag->data[i][j]=fftgrain->norm[i][j];
120  /* detect only if local energy > 10. */ 
121  if (vec_local_energy(newmag)>10.)     {
122    /* hfc = vec_local_hfc(newmag); do not use */
123    aubio_pitchmcomb_spectral_pp(p, newmag);
124    aubio_pitchmcomb_combdet(p,newmag);
125    aubio_pitchmcomb_sort_cand_freq(scands,p->ncand);
126    /* store ncand comb energies in cands[1:ncand] */ 
127    for (k = 0; k<p->ncand; k++) 
128      cands[k] = p->candidates[k]->ene;
129    /* store ncand[end] freq in cands[end] */ 
130    cands[p->ncand] = p->candidates[p->ncand-1]->ebin;
131    return 1;
132  } else {
133    for (k = 0; k<p->ncand; k++)
134      cands[k] = 0;
135    return 0;
136  }
137}
138
139void aubio_pitchmcomb_spectral_pp(aubio_pitchmcomb_t * p, fvec_t * newmag) {
140  fvec_t * mag = (fvec_t *)p->scratch;
141  //fvec_t * tmp = (fvec_t *)p->scratch2;
142  uint_t i=0,j;
143  uint_t length = mag->length;
144  /* copy newmag to mag (scracth) */
145  for (j=0;j<length;j++) {
146    mag->data[i][j] = newmag->data[i][j]; 
147  }
148  vec_dc_removal(mag);               /* dc removal           */
149  vec_alpha_normalise(mag,p->alpha); /* alpha normalisation  */
150  /* skipped */                      /* low pass filtering   */
151  /** \bug: vec_movind_thres writes out of bounds */
152  //vec_adapt_thres(mag,tmp,p->win_post,p->win_pre); /* adaptative threshold */
153  vec_add(mag,-p->threshold);        /* fixed threshold      */
154  {
155    aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
156    uint_t count;
157    /*  return bin and ebin */
158    count = aubio_pitchmcomb_quadpick(peaks,mag);
159    for (j=0;j<count;j++) 
160      peaks[j].mag = newmag->data[i][peaks[j].bin];
161    /* reset non peaks */
162    for (j=count;j<length;j++)
163      peaks[j].mag = 0.;
164    p->peaks = peaks;
165    p->count = count;
166  }
167}
168
169void aubio_pitchmcomb_combdet(aubio_pitchmcomb_t * p, fvec_t * newmag) {
170  aubio_spectralpeak_t * peaks = (aubio_spectralpeak_t *)p->peaks;
171  aubio_spectralcandidate_t ** candidate = 
172    (aubio_spectralcandidate_t **)p->candidates;
173
174  /* parms */
175  uint_t N = p->npartials; /* maximum number of partials to be considered 10 */
176  uint_t M = p->ncand;  /* maximum number of combs to be considered 5 */
177  uint_t length = newmag->length;
178  uint_t count = p->count;
179  uint_t k;
180  uint_t l;
181  uint_t d;
182  uint_t curlen;
183
184  smpl_t delta2;
185  smpl_t xx;
186  uint_t position = 0;
187
188  uint_t root_peak = 0;
189  uint_t tmpl = 0;
190  smpl_t tmpene = 0.;
191
192  /* get the biggest peak in the spectrum */
193  root_peak = aubio_pitchmcomb_get_root_peak(peaks,count);
194  /* now calculate the energy of each of the 5 combs */
195  for (l=0;l<M;l++) {
196    smpl_t scaler = (1./(l+1.));
197    candidate[l]->ene = 0.; /* reset ene and len sums */
198    candidate[l]->len = 0.;
199    candidate[l]->ebin=scaler*peaks[root_peak].ebin;
200    /* if less than N peaks available, curlen < N */
201    curlen = (uint_t)FLOOR(length/(candidate[l]->ebin));
202    curlen = (N < curlen )? N : curlen;
203    /* fill candidate[l]->ecomb[k] with (k+1)*candidate[l]->ebin */
204    for (k=0;k<curlen;k++)
205      candidate[l]->ecomb[k]=(candidate[l]->ebin)*(k+1.);
206    for (k=curlen;k<length;k++)
207      candidate[l]->ecomb[k]=0.;
208    /* for each in candidate[l]->ecomb[k] */
209    for (k=0;k<curlen;k++) {
210      xx = 10000.;
211      /** get the candidate->ecomb the closer to peaks.ebin
212       * (to cope with the inharmonicity)*/
213      for (d=0;d<count;d++) { 
214        delta2 = ABS(candidate[l]->ecomb[k]-peaks[d].ebin);
215        if (delta2 <= xx) {
216          position = d;
217          xx = delta2;
218        }
219      }
220      /* for a Q factor of 17, maintaining "constant Q filtering",
221       * and sum energy and length over non null combs */
222      if ( 17. * xx < candidate[l]->ecomb[k] ) {
223        candidate[l]->ecomb[k]=peaks[position].ebin;
224        candidate[l]->ene += /* ecomb rounded to nearest int */
225          SQR(newmag->data[0][(uint_t)FLOOR(candidate[l]->ecomb[k]-.5)]);
226        candidate[l]->len += 1./curlen;
227      } else
228        candidate[l]->ecomb[k]=0.;
229    }
230    /* punishment */
231    if (candidate[l]->len<0.6)
232      candidate[l]->ene=0.; 
233    /* remember best candidate energy */
234    if (tmpene <= candidate[l]->ene) {
235      tmpl = l;
236      tmpene = candidate[l]->ene;
237    }
238  }
239  //p->candidates=candidate;
240  //p->peaks=peaks;
241  p->goodcandidate = tmpl;
242}
243
244/** T=quadpick(X): return indices of elements of X which are peaks and positive
245 * exact peak positions are retrieved by quadratic interpolation
246 *
247 * \bug peak-picking too picky, sometimes counts too many peaks ?
248 */
249uint_t aubio_pitchmcomb_quadpick(aubio_spectralpeak_t * spectral_peaks, fvec_t * X){
250  uint_t i, j, ispeak, count = 0;
251  for (i=0;i<X->channels;i++)
252    for (j=1;j<X->length-1;j++) {
253      ispeak = vec_peakpick(X,j);
254      if (ispeak) {
255        count += ispeak;
256        spectral_peaks[count-1].bin = j;
257        spectral_peaks[count-1].ebin = vec_quadint(X,j);
258      }
259    }
260  return count;
261}
262
263/* get predominant partial */
264uint_t aubio_pitchmcomb_get_root_peak(aubio_spectralpeak_t * peaks, uint_t length) {
265  uint_t i,pos=0;
266  smpl_t tmp = 0.;
267  for (i=0;i<length;i++)
268    if (tmp <= peaks[i].mag) {
269      pos = i;
270      tmp = peaks[i].mag;
271    }
272  return pos;
273}
274
275void aubio_pitchmcomb_sort_peak(aubio_spectralpeak_t * peaks, uint_t nbins) {
276  qsort(peaks, nbins, sizeof(aubio_spectralpeak_t), 
277      aubio_pitchmcomb_sort_peak_comp);
278}
279static sint_t aubio_pitchmcomb_sort_peak_comp(const void *x, const void *y) {
280  return (((aubio_spectralpeak_t *)y)->mag - ((aubio_spectralpeak_t *)x)->mag);
281}
282
283
284void aubio_pitchmcomb_sort_cand_ene(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
285  uint_t cur = 0;
286  uint_t run = 0;
287  for (cur=0;cur<nbins;cur++) {
288    run = cur + 1;
289    for (run=cur;run<nbins;run++) {
290      if(candidates[run]->ene > candidates[cur]->ene)
291        CAND_SWAP(candidates[run], candidates[cur]);
292    }
293  }
294}
295
296
297void aubio_pitchmcomb_sort_cand_freq(aubio_spectralcandidate_t ** candidates, uint_t nbins) {
298  uint_t cur = 0;
299  uint_t run = 0;
300  for (cur=0;cur<nbins;cur++) {
301    run = cur + 1;
302    for (run=cur;run<nbins;run++) {
303      if(candidates[run]->ebin < candidates[cur]->ebin)
304        CAND_SWAP(candidates[run], candidates[cur]);
305    }
306  }
307}
308
309aubio_pitchmcomb_t * new_aubio_pitchmcomb(uint_t size, uint_t channels) {
310  aubio_pitchmcomb_t * p = AUBIO_NEW(aubio_pitchmcomb_t);
311  /** \bug should check if size / 8 > post+pre+1 */
312  uint_t i;
313  uint_t spec_size;
314  p->spec_partition   = 2;
315  p->ncand            = 5;
316  p->npartials        = 10;
317  p->cutoff           = 1.;
318  p->threshold        = 0.01;
319  p->win_post         = 8;
320  p->win_pre          = 7;
321  p->tau              = 44100./size;
322  p->alpha            = 9.;
323  p->goodcandidate    = 0;
324  spec_size = size/p->spec_partition;
325  //p->pickerfn = quadpick;
326  //p->biquad = new_biquad(0.1600,0.3200,0.1600, -0.5949, 0.2348);
327  /* allocate temp memory */
328  p->newmag     = new_fvec(spec_size,channels);
329  /* array for median */
330  p->scratch    = new_fvec(spec_size,channels);
331  /* array for adaptative threshold */
332  p->scratch2   = new_fvec(p->win_post+p->win_pre+1,channels);
333  /* array of spectral peaks */
334  p->peaks      = AUBIO_ARRAY(aubio_spectralpeak_t,spec_size);
335  /* array of pointers to spectral candidates */
336  p->candidates = AUBIO_ARRAY(aubio_spectralcandidate_t *,p->ncand);
337  for (i=0;i<p->ncand;i++) {
338    p->candidates[i] = AUBIO_NEW(aubio_spectralcandidate_t);
339    p->candidates[i]->ecomb = AUBIO_ARRAY(smpl_t, spec_size);
340  }
341  return p;
342}
343
344
345void del_aubio_pitchmcomb (aubio_pitchmcomb_t *p) {
346  uint_t i;
347  del_fvec(p->newmag);
348  del_fvec(p->scratch);
349  del_fvec(p->scratch2);
350  AUBIO_FREE(p->peaks);
351  for (i=0;i<p->ncand;i++) {
352    AUBIO_FREE(p->candidates[i]);
353  }
354  AUBIO_FREE(p->candidates);
355  AUBIO_FREE(p);
356}
Note: See TracBrowser for help on using the repository browser.