source: src/mathutils.c @ 56ef7e1

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

Change peakpicker to match API specs, make quadint per channel

  • src/mathutils.c
    • add per channel mean and median
    • update moving thres and adapt_thres accordingly
    • change quadint unused span argument to a channel argument
  • src/onset/onset.c:
    • make wasonset a vector for multi channel, use new peakpicker
  • src/onset/peakpick.c:
    • update peakpicker do for multi channeling
  • src/pitch/: update use to fvec_quadint
  • src/tempo/beattracking.c: update calls to fvec_quadint
  • src/tempo/tempo.c: update peakpicker usage
  • tests/src/test-peakpick.c: update peakpicker usage
  • Property mode set to 100644
File size: 12.8 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/* see in mathutils.h for doc */
22
23#include "aubio_priv.h"
24#include "fvec.h"
25#include "mathutils.h"
26#include "config.h"
27
28
29/** Window types */
30typedef enum
31{
32  aubio_win_rectangle,
33  aubio_win_hamming,
34  aubio_win_hanning,
35  aubio_win_hanningz,
36  aubio_win_blackman,
37  aubio_win_blackman_harris,
38  aubio_win_gaussian,
39  aubio_win_welch,
40  aubio_win_parzen,
41  aubio_win_default = aubio_win_hanningz,
42} aubio_window_type;
43
44fvec_t *
45new_aubio_window (char_t * window_type, uint_t size)
46{
47  // create fvec of size x 1 channel
48  fvec_t * win = new_fvec( size, 1);
49  smpl_t * w = win->data[0];
50  uint_t i;
51  aubio_window_type wintype;
52  if (strcmp (window_type, "rectangle") == 0)
53      wintype = aubio_win_rectangle;
54  else if (strcmp (window_type, "hamming") == 0)
55      wintype = aubio_win_hamming;
56  else if (strcmp (window_type, "hanning") == 0)
57      wintype = aubio_win_hanning;
58  else if (strcmp (window_type, "hanningz") == 0)
59      wintype = aubio_win_hanningz;
60  else if (strcmp (window_type, "blackman") == 0)
61      wintype = aubio_win_blackman;
62  else if (strcmp (window_type, "blackman_harris") == 0)
63      wintype = aubio_win_blackman_harris;
64  else if (strcmp (window_type, "gaussian") == 0)
65      wintype = aubio_win_gaussian;
66  else if (strcmp (window_type, "welch") == 0)
67      wintype = aubio_win_welch;
68  else if (strcmp (window_type, "parzen") == 0)
69      wintype = aubio_win_parzen;
70  else if (strcmp (window_type, "default") == 0)
71      wintype = aubio_win_default;
72  else {
73      AUBIO_ERR ("unknown window type %s, using default.\n", window_type);
74      wintype = aubio_win_default;
75      return NULL;
76  }
77  switch(wintype) {
78    case aubio_win_rectangle:
79      for (i=0;i<size;i++)
80        w[i] = 0.5;
81      break;
82    case aubio_win_hamming:
83      for (i=0;i<size;i++)
84        w[i] = 0.54 - 0.46 * COS(TWO_PI * i / (size));
85      break;
86    case aubio_win_hanning:
87      for (i=0;i<size;i++)
88        w[i] = 0.5 - (0.5 * COS(TWO_PI * i / (size)));
89      break;
90    case aubio_win_hanningz:
91      for (i=0;i<size;i++)
92        w[i] = 0.5 * (1.0 - COS(TWO_PI * i / (size)));
93      break;
94    case aubio_win_blackman:
95      for (i=0;i<size;i++)
96        w[i] = 0.42
97          - 0.50 * COS(    TWO_PI*i/(size-1.0))
98          + 0.08 * COS(2.0*TWO_PI*i/(size-1.0));
99      break;
100    case aubio_win_blackman_harris:
101      for (i=0;i<size;i++)
102        w[i] = 0.35875
103          - 0.48829 * COS(    TWO_PI*i/(size-1.0))
104          + 0.14128 * COS(2.0*TWO_PI*i/(size-1.0))
105          - 0.01168 * COS(3.0*TWO_PI*i/(size-1.0));
106      break;
107    case aubio_win_gaussian:
108      for (i=0;i<size;i++)
109        w[i] = EXP(- 1.0 / SQR(size) * SQR(2.0*i-size));
110      break;
111    case aubio_win_welch:
112      for (i=0;i<size;i++)
113        w[i] = 1.0 - SQR((2*i-size)/(size+1.0));
114      break;
115    case aubio_win_parzen:
116      for (i=0;i<size;i++)
117        w[i] = 1.0 - ABS((2*i-size)/(size+1.0));
118      break;
119    default:
120      break;
121  }
122  return win;
123}
124
125smpl_t
126aubio_unwrap2pi (smpl_t phase)
127{
128  /* mod(phase+pi,-2pi)+pi */
129  return phase + TWO_PI * (1. + FLOOR (-(phase + PI) / TWO_PI));
130}
131
132smpl_t
133fvec_mean (fvec_t * s)
134{
135  uint_t i, j;
136  smpl_t tmp = 0.0;
137  for (i = 0; i < s->channels; i++)
138    for (j = 0; j < s->length; j++)
139      tmp += s->data[i][j];
140  return tmp / (smpl_t) (s->length);
141}
142
143smpl_t
144fvec_mean_channel (fvec_t * s, uint_t i)
145{
146  uint_t j;
147  smpl_t tmp = 0.0;
148  for (j = 0; j < s->length; j++)
149      tmp += s->data[i][j];
150  return tmp / (smpl_t) (s->length);
151}
152
153smpl_t
154fvec_sum (fvec_t * s)
155{
156  uint_t i, j;
157  smpl_t tmp = 0.0;
158  for (i = 0; i < s->channels; i++) {
159    for (j = 0; j < s->length; j++) {
160      tmp += s->data[i][j];
161    }
162  }
163  return tmp;
164}
165
166smpl_t
167fvec_max (fvec_t * s)
168{
169  uint_t i, j;
170  smpl_t tmp = 0.0;
171  for (i = 0; i < s->channels; i++) {
172    for (j = 0; j < s->length; j++) {
173      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
174    }
175  }
176  return tmp;
177}
178
179smpl_t
180fvec_min (fvec_t * s)
181{
182  uint_t i, j;
183  smpl_t tmp = s->data[0][0];
184  for (i = 0; i < s->channels; i++) {
185    for (j = 0; j < s->length; j++) {
186      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
187    }
188  }
189  return tmp;
190}
191
192uint_t
193fvec_min_elem (fvec_t * s)
194{
195  uint_t i, j, pos = 0.;
196  smpl_t tmp = s->data[0][0];
197  for (i = 0; i < s->channels; i++) {
198    for (j = 0; j < s->length; j++) {
199      pos = (tmp < s->data[i][j]) ? pos : j;
200      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
201    }
202  }
203  return pos;
204}
205
206uint_t
207fvec_max_elem (fvec_t * s)
208{
209  uint_t i, j, pos = 0;
210  smpl_t tmp = 0.0;
211  for (i = 0; i < s->channels; i++) {
212    for (j = 0; j < s->length; j++) {
213      pos = (tmp > s->data[i][j]) ? pos : j;
214      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
215    }
216  }
217  return pos;
218}
219
220void
221fvec_shift (fvec_t * s)
222{
223  uint_t i, j;
224  for (i = 0; i < s->channels; i++) {
225    for (j = 0; j < s->length / 2; j++) {
226      ELEM_SWAP (s->data[i][j], s->data[i][j + s->length / 2]);
227    }
228  }
229}
230
231smpl_t
232fvec_local_energy (fvec_t * f)
233{
234  smpl_t energy = 0.;
235  uint_t i, j;
236  for (i = 0; i < f->channels; i++) {
237    for (j = 0; j < f->length; j++) {
238      energy += SQR (f->data[i][j]);
239    }
240  }
241  return energy;
242}
243
244smpl_t
245fvec_local_hfc (fvec_t * v)
246{
247  smpl_t hfc = 0.;
248  uint_t i, j;
249  for (i = 0; i < v->channels; i++) {
250    for (j = 0; j < v->length; j++) {
251      hfc += (i + 1) * v->data[i][j];
252    }
253  }
254  return hfc;
255}
256
257void
258fvec_min_removal (fvec_t * v)
259{
260  smpl_t v_min = fvec_min (v);
261  fvec_add (v,  - v_min );
262}
263
264smpl_t
265fvec_alpha_norm (fvec_t * o, smpl_t alpha)
266{
267  uint_t i, j;
268  smpl_t tmp = 0.;
269  for (i = 0; i < o->channels; i++) {
270    for (j = 0; j < o->length; j++) {
271      tmp += POW (ABS (o->data[i][j]), alpha);
272    }
273  }
274  return POW (tmp / o->length, 1. / alpha);
275}
276
277void
278fvec_alpha_normalise (fvec_t * o, smpl_t alpha)
279{
280  uint_t i, j;
281  smpl_t norm = fvec_alpha_norm (o, alpha);
282  for (i = 0; i < o->channels; i++) {
283    for (j = 0; j < o->length; j++) {
284      o->data[i][j] /= norm;
285    }
286  }
287}
288
289void
290fvec_add (fvec_t * o, smpl_t val)
291{
292  uint_t i, j;
293  for (i = 0; i < o->channels; i++) {
294    for (j = 0; j < o->length; j++) {
295      o->data[i][j] += val;
296    }
297  }
298}
299
300void fvec_adapt_thres(fvec_t * vec, fvec_t * tmp,
301    uint_t post, uint_t pre, uint_t channel) {
302  uint_t length = vec->length, i=channel, j;
303  for (j=0;j<length;j++) {
304    vec->data[i][j] -= fvec_moving_thres(vec, tmp, post, pre, j, i);
305  }
306}
307
308smpl_t
309fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec,
310    uint_t post, uint_t pre, uint_t pos, uint_t channel)
311{
312  uint_t i = channel, k;
313  smpl_t *medar = (smpl_t *) tmpvec->data[i];
314  uint_t win_length = post + pre + 1;
315  uint_t length = vec->length;
316  /* post part of the buffer does not exist */
317  if (pos < post + 1) {
318    for (k = 0; k < post + 1 - pos; k++)
319      medar[k] = 0.;            /* 0-padding at the beginning */
320    for (k = post + 1 - pos; k < win_length; k++)
321      medar[k] = vec->data[0][k + pos - post];
322    /* the buffer is fully defined */
323  } else if (pos + pre < length) {
324    for (k = 0; k < win_length; k++)
325      medar[k] = vec->data[0][k + pos - post];
326    /* pre part of the buffer does not exist */
327  } else {
328    for (k = 0; k < length - pos + post; k++)
329      medar[k] = vec->data[0][k + pos - post];
330    for (k = length - pos + post; k < win_length; k++)
331      medar[k] = 0.;            /* 0-padding at the end */
332  }
333  return fvec_median_channel (tmpvec, i);
334}
335
336smpl_t fvec_median_channel (fvec_t * input, uint_t channel) {
337  uint_t n = input->length;
338  smpl_t * arr = (smpl_t *) input->data[channel];
339  uint_t low, high ;
340  uint_t median;
341  uint_t middle, ll, hh;
342
343  low = 0 ; high = n-1 ; median = (low + high) / 2;
344  for (;;) {
345    if (high <= low) /* One element only */
346      return arr[median] ;
347
348    if (high == low + 1) {  /* Two elements only */
349      if (arr[low] > arr[high])
350        ELEM_SWAP(arr[low], arr[high]) ;
351      return arr[median] ;
352    }
353
354    /* Find median of low, middle and high items; swap into position low */
355    middle = (low + high) / 2;
356    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]);
357    if (arr[low]    > arr[high])    ELEM_SWAP(arr[low],    arr[high]);
358    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
359
360    /* Swap low item (now in position middle) into position (low+1) */
361    ELEM_SWAP(arr[middle], arr[low+1]) ;
362
363    /* Nibble from each end towards middle, swapping items when stuck */
364    ll = low + 1;
365    hh = high;
366    for (;;) {
367      do ll++; while (arr[low] > arr[ll]) ;
368      do hh--; while (arr[hh]  > arr[low]) ;
369
370      if (hh < ll)
371        break;
372
373      ELEM_SWAP(arr[ll], arr[hh]) ;
374    }
375
376    /* Swap middle item (in position low) back into correct position */
377    ELEM_SWAP(arr[low], arr[hh]) ;
378
379    /* Re-set active partition */
380    if (hh <= median)
381      low = ll;
382    if (hh >= median)
383      high = hh - 1;
384  }
385}
386
387smpl_t fvec_quadint (fvec_t * x, uint_t pos, uint_t i) {
388  smpl_t s0, s1, s2;
389  uint_t x0 = (pos < 1) ? pos : pos - 1;
390  uint_t x2 = (pos + 1 < x->length) ? pos + 1 : pos;
391  if (x0 == pos) return (x->data[i][pos] <= x->data[i][x2]) ? pos : x2;
392  if (x2 == pos) return (x->data[i][pos] <= x->data[i][x0]) ? pos : x0;
393  s0 = x->data[i][x0];
394  s1 = x->data[i][pos];
395  s2 = x->data[i][x2];
396  return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
397}
398
399uint_t fvec_peakpick(fvec_t * onset, uint_t pos) {
400  uint_t i=0, tmp=0;
401  /*for (i=0;i<onset->channels;i++)*/
402  tmp = (onset->data[i][pos] > onset->data[i][pos-1]
403      &&  onset->data[i][pos] > onset->data[i][pos+1]
404      &&  onset->data[i][pos] > 0.);
405  return tmp;
406}
407
408smpl_t
409aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf)
410{
411  smpl_t tmp =
412      s0 + (pf / 2.) * (pf * (s0 - 2. * s1 + s2) - 3. * s0 + 4. * s1 - s2);
413  return tmp;
414}
415
416smpl_t
417aubio_freqtomidi (smpl_t freq)
418{
419  /* log(freq/A-2)/log(2) */
420  smpl_t midi = freq / 6.875;
421  midi = LOG (midi) / 0.69314718055995;
422  midi *= 12;
423  midi -= 3;
424  return midi;
425}
426
427smpl_t
428aubio_miditofreq (smpl_t midi)
429{
430  smpl_t freq = (midi + 3.) / 12.;
431  freq = EXP (freq * 0.69314718055995);
432  freq *= 6.875;
433  return freq;
434}
435
436smpl_t
437aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
438{
439  smpl_t freq = samplerate / fftsize;
440  return freq * bin;
441}
442
443smpl_t
444aubio_bintomidi (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
445{
446  smpl_t midi = aubio_bintofreq (bin, samplerate, fftsize);
447  return aubio_freqtomidi (midi);
448}
449
450smpl_t
451aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize)
452{
453  smpl_t bin = fftsize / samplerate;
454  return freq * bin;
455}
456
457smpl_t
458aubio_miditobin (smpl_t midi, smpl_t samplerate, smpl_t fftsize)
459{
460  smpl_t freq = aubio_miditofreq (midi);
461  return aubio_freqtobin (freq, samplerate, fftsize);
462}
463
464uint_t
465aubio_is_power_of_two (uint_t a)
466{
467  if ((a & (a - 1)) == 0) {
468    return 1;
469  } else {
470    return 0;
471  }
472}
473
474uint_t
475aubio_next_power_of_two (uint_t a)
476{
477  uint_t i;
478  a--;
479  for (i = 0; i < sizeof (uint_t) * CHAR_BIT; i++) {
480    a = a | a >> 1;
481  }
482  return a + 1;
483}
484
485smpl_t
486aubio_db_spl (fvec_t * o)
487{
488  smpl_t val = SQRT (fvec_local_energy (o));
489  val /= (smpl_t) o->length;
490  return LIN2DB (val);
491}
492
493uint_t
494aubio_silence_detection (fvec_t * o, smpl_t threshold)
495{
496  return (aubio_db_spl (o) < threshold);
497}
498
499smpl_t
500aubio_level_detection (fvec_t * o, smpl_t threshold)
501{
502  smpl_t db_spl = aubio_db_spl (o);
503  if (db_spl < threshold) {
504    return 1.;
505  } else {
506    return db_spl;
507  }
508}
509
510smpl_t
511aubio_zero_crossing_rate (fvec_t * input)
512{
513  uint_t i = 0, j;
514  uint_t zcr = 0;
515  for (j = 1; j < input->length; j++) {
516    // previous was strictly negative
517    if (input->data[i][j - 1] < 0.) {
518      // current is positive or null
519      if (input->data[i][j] >= 0.) {
520        zcr += 1;
521      }
522      // previous was positive or null
523    } else {
524      // current is strictly negative
525      if (input->data[i][j] < 0.) {
526        zcr += 1;
527      }
528    }
529  }
530  return zcr / (smpl_t) input->length;
531}
532
533void
534aubio_autocorr (fvec_t * input, fvec_t * output)
535{
536  uint_t i, j, k, length = input->length;
537  smpl_t *data, *acf;
538  smpl_t tmp = 0;
539  for (k = 0; k < input->channels; k++) {
540    data = input->data[k];
541    acf = output->data[k];
542    for (i = 0; i < length; i++) {
543      tmp = 0.;
544      for (j = i; j < length; j++) {
545        tmp += data[j - i] * data[j];
546      }
547      acf[i] = tmp / (smpl_t) (length - i);
548    }
549  }
550}
551
552void
553aubio_cleanup (void)
554{
555#if HAVE_FFTW3
556  fftw_cleanup ();
557#else
558#if HAVE_FFTW3F
559  fftwf_cleanup ();
560#endif
561#endif
562}
Note: See TracBrowser for help on using the repository browser.