source: src/mathutils.c @ 28d8c4a

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

prefix mathutils function with aubio_

  • Property mode set to 100644
File size: 10.3 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/* see in mathutils.h for doc */
21
22#include "aubio_priv.h"
23#include "sample.h"
24#include "mathutils.h"
25
26void aubio_window(smpl_t *w, uint_t size, aubio_window_type wintype) {
27  uint_t i;
28  switch(wintype) {
29    case aubio_win_rectangle:
30      for (i=0;i<size;i++)
31        w[i] = 0.5; 
32      break;
33    case aubio_win_hamming:
34      for (i=0;i<size;i++)
35        w[i] = 0.54 - 0.46 * COS(TWO_PI * i / (size));
36      break;
37    case aubio_win_hanning:
38      for (i=0;i<size;i++)
39        w[i] = 0.5 - (0.5 * COS(TWO_PI * i / (size)));
40      break;
41    case aubio_win_hanningz:
42      for (i=0;i<size;i++)
43        w[i] = 0.5 * (1.0 - COS(TWO_PI * i / (size)));
44      break;
45    case aubio_win_blackman:
46      for (i=0;i<size;i++)
47        w[i] = 0.42
48          - 0.50 * COS(    TWO_PI*i/(size-1.0))
49          +     0.08 * COS(2.0*TWO_PI*i/(size-1.0));
50      break;
51    case aubio_win_blackman_harris:
52      for (i=0;i<size;i++)
53        w[i] = 0.35875 
54          - 0.48829 * COS(    TWO_PI*i/(size-1.0))
55          + 0.14128 * COS(2.0*TWO_PI*i/(size-1.0))
56          - 0.01168 * COS(3.0*TWO_PI*i/(size-1.0));
57      break;
58    case aubio_win_gaussian:
59      for (i=0;i<size;i++)
60        w[i] = EXP(- 1.0 / SQR(size) * SQR(2.0*i-size));
61      break;
62    case aubio_win_welch:
63      for (i=0;i<size;i++)
64        w[i] = 1.0 - SQR((2*i-size)/(size+1.0));
65      break;
66    case aubio_win_parzen:
67      for (i=0;i<size;i++)
68        w[i] = 1.0 - fabsf((2*i-size)/(size+1.0));
69      break;
70    default:
71      break;
72  }
73}
74
75
76smpl_t aubio_unwrap2pi(smpl_t phase) {
77  /* mod(phase+pi,-2pi)+pi */
78  return phase + TWO_PI * (1. + floorf(-(phase+PI)/TWO_PI));
79}
80
81
82smpl_t vec_mean(fvec_t *s) 
83{
84  uint_t i,j;
85  smpl_t tmp = 0.0f;
86  for (i=0; i < s->channels; i++)
87    for (j=0; j < s->length; j++)
88      tmp += s->data[i][j];
89  return tmp/(smpl_t)(s->length);
90}
91
92
93smpl_t vec_sum(fvec_t *s) 
94{
95  uint_t i,j;
96  smpl_t tmp = 0.0f;
97  for (i=0; i < s->channels; i++)
98    for (j=0; j < s->length; j++)
99      tmp += s->data[i][j];
100  return tmp;
101}
102
103
104smpl_t vec_max(fvec_t *s) 
105{
106  uint_t i,j;
107  smpl_t tmp = 0.0f;
108  for (i=0; i < s->channels; i++)
109    for (j=0; j < s->length; j++)
110      tmp = (tmp > s->data[i][j])? tmp : s->data[i][j];
111  return tmp;
112}
113
114smpl_t vec_min(fvec_t *s) 
115{
116  uint_t i,j;
117  smpl_t tmp = s->data[0][0];
118  for (i=0; i < s->channels; i++)
119    for (j=0; j < s->length; j++)
120      tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
121  return tmp;
122}
123
124
125uint_t vec_min_elem(fvec_t *s) 
126{
127  uint_t i,j=0, pos=0.;
128  smpl_t tmp = s->data[0][0];
129  for (i=0; i < s->channels; i++)
130    for (j=0; j < s->length; j++) {
131      pos = (tmp < s->data[i][j])? pos : j;
132      tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
133    }
134  return pos;
135}
136
137uint_t vec_max_elem(fvec_t *s) 
138{
139  uint_t i,j=0, pos=0.;
140  smpl_t tmp = 0.0f;
141  for (i=0; i < s->channels; i++)
142    for (j=0; j < s->length; j++) {
143      pos = (tmp > s->data[i][j])? pos : j;
144      tmp = (tmp > s->data[i][j])? tmp : s->data[i][j]  ;
145    }
146  return pos;
147}
148
149void vec_shift(fvec_t *s)
150{
151  uint_t i,j;
152  //smpl_t tmp = 0.0f;
153  for (i=0; i < s->channels; i++)
154    for (j=0; j < s->length / 2 ; j++) {
155      //tmp = s->data[i][j];
156      //s->data[i][j] = s->data[i][j+s->length/2];
157      //s->data[i][j+s->length/2] = tmp;
158      ELEM_SWAP(s->data[i][j],s->data[i][j+s->length/2]);
159    }
160}
161
162smpl_t vec_local_energy(fvec_t * f)
163{
164  smpl_t locE = 0.;
165  uint_t i,j;
166  for (i=0;i<f->channels;i++)
167    for (j=0;j<f->length;j++)
168      locE+=SQR(f->data[i][j]);
169  return locE;
170}
171
172smpl_t vec_local_hfc(fvec_t * f)
173{
174  smpl_t locE = 0.;
175  uint_t i,j;
176  for (i=0;i<f->channels;i++)
177    for (j=0;j<f->length;j++)
178      locE+=(i+1)*f->data[i][j];
179  return locE;
180}
181
182smpl_t vec_alpha_norm(fvec_t * DF, smpl_t alpha)
183{
184  smpl_t tmp = 0.;
185  uint_t i,j;
186  for (i=0;i<DF->channels;i++)
187    for (j=0;j<DF->length;j++)
188      tmp += POW(ABS(DF->data[i][j]),alpha);
189  return POW(tmp/DF->length,1./alpha);
190}
191
192
193void vec_dc_removal(fvec_t * mag)
194{
195    smpl_t mini = 0.;
196    uint_t length = mag->length, i=0, j;
197    mini = vec_min(mag);
198    for (j=0;j<length;j++) {
199      mag->data[i][j] -= mini;
200    }
201}
202
203
204void vec_alpha_normalise(fvec_t * mag, uint_t alpha)
205{
206  smpl_t alphan = 1.;
207  uint_t length = mag->length, i=0, j;
208  alphan = vec_alpha_norm(mag,alpha);
209  for (j=0;j<length;j++){
210    mag->data[i][j] /= alphan;
211  }
212}
213
214
215void vec_add(fvec_t * mag, smpl_t threshold) {
216  uint_t length = mag->length, i=0, j;
217  for (j=0;j<length;j++) {
218    mag->data[i][j] += threshold;
219  }
220}
221
222
223void vec_adapt_thres(fvec_t * vec, fvec_t * tmp, 
224    uint_t post, uint_t pre) 
225{
226  uint_t length = vec->length, i=0, j;
227  for (j=0;j<length;j++) {
228    vec->data[i][j] -= vec_moving_thres(vec, tmp, post, pre, j);
229  }
230}
231
232smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmpvec,
233    uint_t post, uint_t pre, uint_t pos) 
234{
235  smpl_t * medar = (smpl_t *)tmpvec->data[0];
236  uint_t k;
237  uint_t win_length =  post+pre+1;
238  uint_t length =  vec->length;
239  /* post part of the buffer does not exist */
240  if (pos<post+1) {
241    for (k=0;k<post+1-pos;k++) 
242      medar[k] = 0.; /* 0-padding at the beginning */
243    for (k=post+1-pos;k<win_length;k++)
244      medar[k] = vec->data[0][k+pos-post];
245  /* the buffer is fully defined */
246  } else if (pos+pre<length) {
247    for (k=0;k<win_length;k++)
248      medar[k] = vec->data[0][k+pos-post];
249  /* pre part of the buffer does not exist */
250  } else {
251    for (k=0;k<length-pos+post+1;k++)
252      medar[k] = vec->data[0][k+pos-post];
253    for (k=length-pos+post+1;k<win_length;k++) 
254      medar[k] = 0.; /* 0-padding at the end */
255  } 
256  return vec_median(tmpvec);
257}
258
259smpl_t vec_median(fvec_t * input) {
260  uint_t n = input->length;
261  smpl_t * arr = (smpl_t *) input->data[0];
262  uint_t low, high ;
263  uint_t median;
264  uint_t middle, ll, hh;
265
266  low = 0 ; high = n-1 ; median = (low + high) / 2;
267  for (;;) {
268    if (high <= low) /* One element only */
269      return arr[median] ;
270
271    if (high == low + 1) {  /* Two elements only */
272      if (arr[low] > arr[high])
273        ELEM_SWAP(arr[low], arr[high]) ;
274      return arr[median] ;
275    }
276
277    /* Find median of low, middle and high items; swap into position low */
278    middle = (low + high) / 2;
279    if (arr[middle] > arr[high])    ELEM_SWAP(arr[middle], arr[high]);
280    if (arr[low]    > arr[high])    ELEM_SWAP(arr[low],    arr[high]);
281    if (arr[middle] > arr[low])     ELEM_SWAP(arr[middle], arr[low]) ;
282
283    /* Swap low item (now in position middle) into position (low+1) */
284    ELEM_SWAP(arr[middle], arr[low+1]) ;
285
286    /* Nibble from each end towards middle, swapping items when stuck */
287    ll = low + 1;
288    hh = high;
289    for (;;) {
290      do ll++; while (arr[low] > arr[ll]) ;
291      do hh--; while (arr[hh]  > arr[low]) ;
292
293      if (hh < ll)
294        break;
295
296      ELEM_SWAP(arr[ll], arr[hh]) ;
297    }
298
299    /* Swap middle item (in position low) back into correct position */
300    ELEM_SWAP(arr[low], arr[hh]) ;
301
302    /* Re-set active partition */
303    if (hh <= median)
304      low = ll;
305    if (hh >= median)
306      high = hh - 1;
307  }
308}
309
310smpl_t vec_quadint(fvec_t * x,uint_t pos) {
311  uint_t span = 2;
312  smpl_t step = 1./200.;
313  /* hack : init resold to - something (in case x[pos+-span]<0)) */
314  smpl_t res, frac, s0, s1, s2, exactpos = (smpl_t)pos, resold = -1000.;
315  if ((pos > span) && (pos < x->length-span)) {
316    s0 = x->data[0][pos-span];
317    s1 = x->data[0][pos]     ;
318    s2 = x->data[0][pos+span];
319    /* increase frac */
320    for (frac = 0.; frac < 2.; frac = frac + step) {
321      res = aubio_quadfrac(s0, s1, s2, frac);
322      if (res > resold) 
323        resold = res;
324      else {                           
325        exactpos += (frac-step)*2. - 1.;
326        break;
327      }
328    }
329  }
330  return exactpos;
331}
332
333smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
334  smpl_t tmp = s0 + (pf/2.) * (pf * ( s0 - 2.*s1 + s2 ) - 3.*s0 + 4.*s1 - s2);
335  return tmp;
336}
337
338uint_t vec_peakpick(fvec_t * onset, uint_t pos) {
339        uint_t i=0, tmp=0;
340        /*for (i=0;i<onset->channels;i++)*/
341                tmp = (onset->data[i][pos] > onset->data[i][pos-1]
342                        &&  onset->data[i][pos] > onset->data[i][pos+1]
343                        &&      onset->data[i][pos] > 0.);
344        return tmp;
345}
346
347smpl_t aubio_freqtomidi(smpl_t freq) {
348  smpl_t midi = freq/6.875;
349  /* log(freq/A-2)/log(2) */
350  midi = LOG(midi)/0.69314718055995;
351  midi *= 12;
352  midi -= 3; 
353  return midi;
354}
355
356smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
357  smpl_t freq = samplerate/fftsize;
358  return freq*bin;
359}
360
361
362smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
363  smpl_t midi = aubio_bintofreq(bin,samplerate,fftsize);
364  return aubio_freqtomidi(midi);
365}
366
367
368
369/** returns 1 if wassilence is 0 and RMS(ibuf)<threshold
370 * \bug mono
371 */
372uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold) {
373  smpl_t loudness = 0;
374  uint_t i=0,j;
375  for (j=0;j<ibuf->length;j++) {
376    loudness += SQR(ibuf->data[i][j]);
377  }
378  loudness = SQRT(loudness);
379  loudness /= (smpl_t)ibuf->length;
380  loudness = LIN2DB(loudness);
381
382  return (loudness < threshold);
383}
384
385/** returns level log(RMS(ibuf)) if < threshold, 1 otherwise
386 * \bug mono
387 */
388smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold) {
389  smpl_t loudness = 0;
390  uint_t i=0,j;
391  for (j=0;j<ibuf->length;j++) {
392    loudness += SQR(ibuf->data[i][j]);
393  }
394  loudness = SQRT(loudness);
395  loudness /= (smpl_t)ibuf->length;
396  loudness = LIN2DB(loudness);
397
398  if (loudness < threshold)
399      return 1.;
400  else
401      return loudness;
402}
403
404void aubio_autocorr(fvec_t * input, fvec_t * output){
405        uint_t i = 0, j = 0, length = input->length;
406        smpl_t * data = input->data[0];
407        smpl_t * acf = output->data[0];
408        smpl_t tmp =0.;
409        for(i=0;i<length;i++){
410                for(j=i;j<length;j++){
411                        tmp += data[j-i]*data[j]; 
412                }
413                acf[i] = tmp /(smpl_t)(length-i);
414                tmp = 0.0;
415        }
416}
417
Note: See TracBrowser for help on using the repository browser.