Changeset eb7f743 for src/mathutils.c


Ignore:
Timestamp:
Oct 2, 2009, 6:11:07 AM (10 years ago)
Author:
Paul Brossier <piem@piem.org>
Branches:
feature/autosink, feature/constantq, feature/pitchshift, feature/pydocstrings, feature/timestretch, master, pitchshift, sampler, timestretch, yinfft+
Children:
352fd5f
Parents:
1498ced
Message:

src/mathutils.{c,h}: loop over channels when possible, improve documentation, indent

File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/mathutils.c

    r1498ced reb7f743  
    2626#include "config.h"
    2727
    28 fvec_t * new_aubio_window(uint_t size, aubio_window_type wintype) {
     28fvec_t *
     29new_aubio_window (uint_t size, aubio_window_type wintype)
     30{
    2931  // create fvec of size x 1 channel
    3032  fvec_t * win = new_fvec( size, 1);
     
    7981}
    8082
    81 smpl_t aubio_unwrap2pi(smpl_t phase) {
     83smpl_t
     84aubio_unwrap2pi (smpl_t phase)
     85{
    8286  /* mod(phase+pi,-2pi)+pi */
    83   return phase + TWO_PI * (1. + FLOOR(-(phase+PI)/TWO_PI));
    84 }
    85 
    86 smpl_t fvec_mean(fvec_t *s) {
    87   uint_t i,j;
     87  return phase + TWO_PI * (1. + FLOOR (-(phase + PI) / TWO_PI));
     88}
     89
     90smpl_t
     91fvec_mean (fvec_t * s)
     92{
     93  uint_t i, j;
    8894  smpl_t tmp = 0.0;
    89   for (i=0; i < s->channels; i++)
    90     for (j=0; j < s->length; j++)
     95  for (i = 0; i < s->channels; i++)
     96    for (j = 0; j < s->length; j++)
    9197      tmp += s->data[i][j];
    92   return tmp/(smpl_t)(s->length);
    93 }
    94 
    95 smpl_t fvec_sum(fvec_t *s) {
    96   uint_t i,j;
     98  return tmp / (smpl_t) (s->length);
     99}
     100
     101smpl_t
     102fvec_sum (fvec_t * s)
     103{
     104  uint_t i, j;
    97105  smpl_t tmp = 0.0;
    98   for (i=0; i < s->channels; i++)
    99     for (j=0; j < s->length; j++)
     106  for (i = 0; i < s->channels; i++) {
     107    for (j = 0; j < s->length; j++) {
    100108      tmp += s->data[i][j];
     109    }
     110  }
    101111  return tmp;
    102112}
    103113
    104 smpl_t fvec_max(fvec_t *s) {
    105   uint_t i,j;
     114smpl_t
     115fvec_max (fvec_t * s)
     116{
     117  uint_t i, j;
    106118  smpl_t tmp = 0.0;
    107   for (i=0; i < s->channels; i++)
    108     for (j=0; j < s->length; j++)
    109       tmp = (tmp > s->data[i][j])? tmp : s->data[i][j];
     119  for (i = 0; i < s->channels; i++) {
     120    for (j = 0; j < s->length; j++) {
     121      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
     122    }
     123  }
    110124  return tmp;
    111125}
    112126
    113 smpl_t fvec_min(fvec_t *s) {
    114   uint_t i,j;
     127smpl_t
     128fvec_min (fvec_t * s)
     129{
     130  uint_t i, j;
    115131  smpl_t tmp = s->data[0][0];
    116   for (i=0; i < s->channels; i++)
    117     for (j=0; j < s->length; j++)
    118       tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
     132  for (i = 0; i < s->channels; i++) {
     133    for (j = 0; j < s->length; j++) {
     134      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
     135    }
     136  }
    119137  return tmp;
    120138}
    121139
    122 uint_t fvec_min_elem(fvec_t *s) {
    123   uint_t i,j=0, pos=0.;
     140uint_t
     141fvec_min_elem (fvec_t * s)
     142{
     143  uint_t i, j = 0, pos = 0.;
    124144  smpl_t tmp = s->data[0][0];
    125   for (i=0; i < s->channels; i++)
    126     for (j=0; j < s->length; j++) {
    127       pos = (tmp < s->data[i][j])? pos : j;
    128       tmp = (tmp < s->data[i][j])? tmp : s->data[i][j]  ;
    129     }
     145  for (i = 0; i < s->channels; i++) {
     146    for (j = 0; j < s->length; j++) {
     147      pos = (tmp < s->data[i][j]) ? pos : j;
     148      tmp = (tmp < s->data[i][j]) ? tmp : s->data[i][j];
     149    }
     150  }
    130151  return pos;
    131152}
    132153
    133 uint_t fvec_max_elem(fvec_t *s) {
    134   uint_t i,j=0, pos=0.;
     154uint_t
     155fvec_max_elem (fvec_t * s)
     156{
     157  uint_t i, j, pos;
    135158  smpl_t tmp = 0.0;
    136   for (i=0; i < s->channels; i++)
    137     for (j=0; j < s->length; j++) {
    138       pos = (tmp > s->data[i][j])? pos : j;
    139       tmp = (tmp > s->data[i][j])? tmp : s->data[i][j]  ;
    140     }
     159  for (i = 0; i < s->channels; i++) {
     160    for (j = 0; j < s->length; j++) {
     161      pos = (tmp > s->data[i][j]) ? pos : j;
     162      tmp = (tmp > s->data[i][j]) ? tmp : s->data[i][j];
     163    }
     164  }
    141165  return pos;
    142166}
    143167
    144 void fvec_shift(fvec_t *s) {
    145   uint_t i,j;
    146   //smpl_t tmp = 0.0;
    147   for (i=0; i < s->channels; i++)
    148     for (j=0; j < s->length / 2 ; j++) {
    149       //tmp = s->data[i][j];
    150       //s->data[i][j] = s->data[i][j+s->length/2];
    151       //s->data[i][j+s->length/2] = tmp;
    152       ELEM_SWAP(s->data[i][j],s->data[i][j+s->length/2]);
    153     }
    154 }
    155 
    156 smpl_t fvec_local_energy(fvec_t * f) {
    157   smpl_t locE = 0.;
    158   uint_t i,j;
    159   for (i=0;i<f->channels;i++)
    160     for (j=0;j<f->length;j++)
    161       locE+=SQR(f->data[i][j]);
    162   return locE;
    163 }
    164 
    165 smpl_t fvec_local_hfc(fvec_t * f) {
    166   smpl_t locE = 0.;
    167   uint_t i,j;
    168   for (i=0;i<f->channels;i++)
    169     for (j=0;j<f->length;j++)
    170       locE+=(i+1)*f->data[i][j];
    171   return locE;
    172 }
    173 
    174 smpl_t fvec_alpha_norm(fvec_t * DF, smpl_t alpha) {
     168void
     169fvec_shift (fvec_t * s)
     170{
     171  uint_t i, j;
     172  for (i = 0; i < s->channels; i++) {
     173    for (j = 0; j < s->length / 2; j++) {
     174      ELEM_SWAP (s->data[i][j], s->data[i][j + s->length / 2]);
     175    }
     176  }
     177}
     178
     179smpl_t
     180fvec_local_energy (fvec_t * f)
     181{
     182  smpl_t energy = 0.;
     183  uint_t i, j;
     184  for (i = 0; i < f->channels; i++) {
     185    for (j = 0; j < f->length; j++) {
     186      energy += SQR (f->data[i][j]);
     187    }
     188  }
     189  return energy;
     190}
     191
     192smpl_t
     193fvec_local_hfc (fvec_t * v)
     194{
     195  smpl_t hfc = 0.;
     196  uint_t i, j;
     197  for (i = 0; i < v->channels; i++) {
     198    for (j = 0; j < v->length; j++) {
     199      hfc += (i + 1) * v->data[i][j];
     200    }
     201  }
     202  return hfc;
     203}
     204
     205void
     206fvec_min_removal (fvec_t * v)
     207{
     208  smpl_t v_min = fvec_min (v);
     209  fvec_add (v,  - v_min );
     210}
     211
     212smpl_t
     213fvec_alpha_norm (fvec_t * o, smpl_t alpha)
     214{
     215  uint_t i, j;
    175216  smpl_t tmp = 0.;
    176   uint_t i,j;
    177   for (i=0;i<DF->channels;i++)
    178     for (j=0;j<DF->length;j++)
    179       tmp += POW(ABS(DF->data[i][j]),alpha);
    180   return POW(tmp/DF->length,1./alpha);
    181 }
    182 
    183 void
    184 fvec_min_removal (fvec_t * o)
    185 {
    186   uint_t i, j;
    187   smpl_t mini = fvec_min (o);
    188217  for (i = 0; i < o->channels; i++) {
    189218    for (j = 0; j < o->length; j++) {
    190       o->data[i][j] -= mini;
    191     }
    192   }
    193 }
    194 
    195 void fvec_alpha_normalise(fvec_t * mag, uint_t alpha) {
    196   smpl_t alphan = 1.;
    197   uint_t length = mag->length, i=0, j;
    198   alphan = fvec_alpha_norm(mag,alpha);
    199   for (j=0;j<length;j++){
    200     mag->data[i][j] /= alphan;
    201   }
    202 }
    203 
    204 void fvec_add(fvec_t * mag, smpl_t threshold) {
    205   uint_t length = mag->length, i=0, j;
    206   for (j=0;j<length;j++) {
    207     mag->data[i][j] += threshold;
     219      tmp += POW (ABS (o->data[i][j]), alpha);
     220    }
     221  }
     222  return POW (tmp / o->length, 1. / alpha);
     223}
     224
     225void
     226fvec_alpha_normalise (fvec_t * o, smpl_t alpha)
     227{
     228  uint_t i, j;
     229  smpl_t norm = fvec_alpha_norm (o, alpha);
     230  for (i = 0; i < o->channels; i++) {
     231    for (j = 0; j < o->length; j++) {
     232      o->data[i][j] /= norm;
     233    }
     234  }
     235}
     236
     237void
     238fvec_add (fvec_t * o, smpl_t val)
     239{
     240  uint_t i, j;
     241  for (i = 0; i < o->channels; i++) {
     242    for (j = 0; j < o->length; j++) {
     243      o->data[i][j] += val;
     244    }
    208245  }
    209246}
     
    217254}
    218255
    219 smpl_t fvec_moving_thres(fvec_t * vec, fvec_t * tmpvec,
    220     uint_t post, uint_t pre, uint_t pos) {
    221   smpl_t * medar = (smpl_t *)tmpvec->data[0];
     256smpl_t
     257fvec_moving_thres (fvec_t * vec, fvec_t * tmpvec,
     258    uint_t post, uint_t pre, uint_t pos)
     259{
     260  smpl_t *medar = (smpl_t *) tmpvec->data[0];
    222261  uint_t k;
    223   uint_t win_length =  post+pre+1;
    224   uint_t length =  vec->length;
     262  uint_t win_length = post + pre + 1;
     263  uint_t length = vec->length;
    225264  /* post part of the buffer does not exist */
    226   if (pos<post+1) {
    227     for (k=0;k<post+1-pos;k++)
    228       medar[k] = 0.; /* 0-padding at the beginning */
    229     for (k=post+1-pos;k<win_length;k++)
    230       medar[k] = vec->data[0][k+pos-post];
    231   /* the buffer is fully defined */
    232   } else if (pos+pre<length) {
    233     for (k=0;k<win_length;k++)
    234       medar[k] = vec->data[0][k+pos-post];
    235   /* pre part of the buffer does not exist */
     265  if (pos < post + 1) {
     266    for (k = 0; k < post + 1 - pos; k++)
     267      medar[k] = 0.;            /* 0-padding at the beginning */
     268    for (k = post + 1 - pos; k < win_length; k++)
     269      medar[k] = vec->data[0][k + pos - post];
     270    /* the buffer is fully defined */
     271  } else if (pos + pre < length) {
     272    for (k = 0; k < win_length; k++)
     273      medar[k] = vec->data[0][k + pos - post];
     274    /* pre part of the buffer does not exist */
    236275  } else {
    237     for (k=0;k<length-pos+post;k++)
    238       medar[k] = vec->data[0][k+pos-post];
    239     for (k=length-pos+post;k<win_length;k++)
    240       medar[k] = 0.; /* 0-padding at the end */
    241   }
    242   return fvec_median(tmpvec);
     276    for (k = 0; k < length - pos + post; k++)
     277      medar[k] = vec->data[0][k + pos - post];
     278    for (k = length - pos + post; k < win_length; k++)
     279      medar[k] = 0.;            /* 0-padding at the end */
     280  }
     281  return fvec_median (tmpvec);
    243282}
    244283
     
    301340  if (x2 == pos) return (x->data[0][pos] <= x->data[0][x0]) ? pos : x0;
    302341  s0 = x->data[0][x0];
    303   s1 = x->data[0][pos]     ;
     342  s1 = x->data[0][pos];
    304343  s2 = x->data[0][x2];
    305344  return pos + 0.5 * (s2 - s0 ) / (s2 - 2.* s1 + s0);
    306 }
    307 
    308 smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf) {
    309   smpl_t tmp = s0 + (pf/2.) * (pf * ( s0 - 2.*s1 + s2 ) - 3.*s0 + 4.*s1 - s2);
    310   return tmp;
    311345}
    312346
     
    320354}
    321355
    322 smpl_t aubio_freqtomidi(smpl_t freq) {
     356smpl_t
     357aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf)
     358{
     359  smpl_t tmp =
     360      s0 + (pf / 2.) * (pf * (s0 - 2. * s1 + s2) - 3. * s0 + 4. * s1 - s2);
     361  return tmp;
     362}
     363
     364smpl_t
     365aubio_freqtomidi (smpl_t freq)
     366{
    323367  /* log(freq/A-2)/log(2) */
    324   smpl_t midi = freq/6.875;
    325   midi = LOG(midi)/0.69314718055995;
     368  smpl_t midi = freq / 6.875;
     369  midi = LOG (midi) / 0.69314718055995;
    326370  midi *= 12;
    327371  midi -= 3;
     
    329373}
    330374
    331 smpl_t aubio_miditofreq(smpl_t midi) {
    332   smpl_t freq = (midi+3.)/12.;
    333   freq = EXP(freq*0.69314718055995);
     375smpl_t
     376aubio_miditofreq (smpl_t midi)
     377{
     378  smpl_t freq = (midi + 3.) / 12.;
     379  freq = EXP (freq * 0.69314718055995);
    334380  freq *= 6.875;
    335381  return freq;
    336382}
    337383
    338 smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
    339   smpl_t freq = samplerate/fftsize;
    340   return freq*bin;
    341 }
    342 
    343 smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize) {
    344   smpl_t midi = aubio_bintofreq(bin,samplerate,fftsize);
    345   return aubio_freqtomidi(midi);
    346 }
    347 
    348 smpl_t aubio_freqtobin(smpl_t freq, smpl_t samplerate, smpl_t fftsize) {
    349   smpl_t bin = fftsize/samplerate;
    350   return freq*bin;
    351 }
    352 
    353 smpl_t aubio_miditobin(smpl_t midi, smpl_t samplerate, smpl_t fftsize) {
    354   smpl_t freq = aubio_miditofreq(midi);
    355   return aubio_freqtobin(freq,samplerate,fftsize);
    356 }
    357 
    358 /** returns 1 if wassilence is 0 and RMS(ibuf)<threshold
    359  * \bug mono
    360  */
    361 uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold) {
    362   smpl_t loudness = 0;
    363   uint_t i=0,j;
    364   for (j=0;j<ibuf->length;j++) {
    365     loudness += SQR(ibuf->data[i][j]);
    366   }
    367   loudness = SQRT(loudness);
    368   loudness /= (smpl_t)ibuf->length;
    369   loudness = LIN2DB(loudness);
    370 
    371   return (loudness < threshold);
    372 }
    373 
    374 /** returns level log(RMS(ibuf)) if < threshold, 1 otherwise
    375  * \bug mono
    376  */
    377 smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold) {
    378   smpl_t loudness = 0;
    379   uint_t i=0,j;
    380   for (j=0;j<ibuf->length;j++) {
    381     loudness += SQR(ibuf->data[i][j]);
    382   }
    383   loudness = SQRT(loudness);
    384   loudness /= (smpl_t)ibuf->length;
    385   loudness = LIN2DB(loudness);
    386 
    387   if (loudness < threshold)
     384smpl_t
     385aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
     386{
     387  smpl_t freq = samplerate / fftsize;
     388  return freq * bin;
     389}
     390
     391smpl_t
     392aubio_bintomidi (smpl_t bin, smpl_t samplerate, smpl_t fftsize)
     393{
     394  smpl_t midi = aubio_bintofreq (bin, samplerate, fftsize);
     395  return aubio_freqtomidi (midi);
     396}
     397
     398smpl_t
     399aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize)
     400{
     401  smpl_t bin = fftsize / samplerate;
     402  return freq * bin;
     403}
     404
     405smpl_t
     406aubio_miditobin (smpl_t midi, smpl_t samplerate, smpl_t fftsize)
     407{
     408  smpl_t freq = aubio_miditofreq (midi);
     409  return aubio_freqtobin (freq, samplerate, fftsize);
     410}
     411
     412smpl_t
     413aubio_db_spl (fvec_t * o)
     414{
     415  smpl_t val = SQRT (fvec_local_energy (o));
     416  val /= (smpl_t) o->length;
     417  return LIN2DB (val);
     418}
     419
     420uint_t
     421aubio_silence_detection (fvec_t * o, smpl_t threshold)
     422{
     423  return (aubio_db_spl (o) < threshold);
     424}
     425
     426smpl_t
     427aubio_level_detection (fvec_t * o, smpl_t threshold)
     428{
     429  smpl_t db_spl = aubio_db_spl (o);
     430  if (db_spl < threshold) {
    388431    return 1.;
    389   else
    390     return loudness;
    391 }
    392 
    393 smpl_t aubio_zero_crossing_rate(fvec_t * input) {
    394   uint_t i=0,j;
     432  } else {
     433    return db_spl;
     434  }
     435}
     436
     437smpl_t
     438aubio_zero_crossing_rate (fvec_t * input)
     439{
     440  uint_t i = 0, j;
    395441  uint_t zcr = 0;
    396   for ( j = 1; j < input->length; j++ ) {
     442  for (j = 1; j < input->length; j++) {
    397443    // previous was strictly negative
    398     if( input->data[i][j-1] < 0. ) {
     444    if (input->data[i][j - 1] < 0.) {
    399445      // current is positive or null
    400       if ( input->data[i][j] >= 0. ) {
     446      if (input->data[i][j] >= 0.) {
    401447        zcr += 1;
    402448      }
    403     // previous was positive or null
     449      // previous was positive or null
    404450    } else {
    405451      // current is strictly negative
    406       if ( input->data[i][j] < 0. ) {
     452      if (input->data[i][j] < 0.) {
    407453        zcr += 1;
    408454      }
    409455    }
    410456  }
    411   return zcr/(smpl_t)input->length;
    412 }
    413 
    414 void aubio_autocorr(fvec_t * input, fvec_t * output) {
    415   uint_t i = 0, j = 0, length = input->length;
    416   smpl_t * data = input->data[0];
    417   smpl_t * acf = output->data[0];
    418   smpl_t tmp =0.;
    419   for(i=0;i<length;i++){
    420     for(j=i;j<length;j++){
    421       tmp += data[j-i]*data[j];
    422     }
    423     acf[i] = tmp /(smpl_t)(length-i);
    424     tmp = 0.0;
    425   }
    426 }
    427 
    428 void aubio_cleanup(void) {
     457  return zcr / (smpl_t) input->length;
     458}
     459
     460void
     461aubio_autocorr (fvec_t * input, fvec_t * output)
     462{
     463  uint_t i, j, k, length = input->length;
     464  smpl_t *data, *acf;
     465  smpl_t tmp = 0;
     466  for (k = 0; k < input->channels; k++) {
     467    data = input->data[k];
     468    acf = output->data[k];
     469    for (i = 0; i < length; i++) {
     470      tmp = 0.;
     471      for (j = i; j < length; j++) {
     472        tmp += data[j - i] * data[j];
     473      }
     474      acf[i] = tmp / (smpl_t) (length - i);
     475    }
     476  }
     477}
     478
     479void
     480aubio_cleanup (void)
     481{
    429482#if HAVE_FFTW3
    430   fftw_cleanup();
     483  fftw_cleanup ();
    431484#else
    432485#if HAVE_FFTW3F
    433   fftwf_cleanup();
     486  fftwf_cleanup ();
    434487#endif
    435488#endif
Note: See TracChangeset for help on using the changeset viewer.