src/mathutils.h

   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 /** @file
  22  *  various math functions
  23  */
  24
  25 #ifndef MATHUTILS_H
  26 #define MATHUTILS_H
  27
  28 #ifdef __cplusplus
  29 extern "C" {
  30 #endif
  31
  32 /** Window types
  33
  34   References:
  35
  36     - <a href="http://en.wikipedia.org/wiki/Window_function">Window
  37 function</a> on Wikipedia
  38     - Amalia de Götzen, Nicolas Bernardini, and Daniel Arfib. Traditional (?)
  39 implementations of a phase vocoder: the tricks of the trade. In Proceedings of
  40 the International Conference on Digital Audio Effects (DAFx-00), pages 37–44,
  41 Uni- versity of Verona, Italy, 2000.
  42   (<a href="http://profs.sci.univr.it/%7Edafx/Final-Papers/ps/Bernardini.ps.gz">
  43   ps.gz</a>)
  44
  45 */
  46 typedef enum
  47 {
  48   aubio_win_rectangle,
  49   aubio_win_hamming,
  50   aubio_win_hanning,
  51   aubio_win_hanningz,
  52   aubio_win_blackman,
  53   aubio_win_blackman_harris,
  54   aubio_win_gaussian,
  55   aubio_win_welch,
  56   aubio_win_parzen
  57 } aubio_window_type;
  58
  59 /** create window */
  60 fvec_t *new_aubio_window (uint_t size, aubio_window_type wintype);
  61
  62 /** compute the principal argument
  63
  64   This function maps the input phase to its corresponding value wrapped in the
  65 range \f$ [-\pi, \pi] \f$.
  66
  67   \param phase unwrapped phase to map to the unit circle
  68
  69   \return equivalent phase wrapped to the unit circle
  70
  71 */
  72 smpl_t aubio_unwrap2pi (smpl_t phase);
  73
  74 /** compute the mean of a vector
  75
  76   \param s vector to compute norm from
  77
  78   \return the mean of v
  79
  80 */
  81 smpl_t fvec_mean (fvec_t * s);
  82
  83 /** find the max of a vector
  84
  85   \param s vector to get the max from
  86
  87   \return the value of the minimum of v
  88
  89 */
  90 smpl_t fvec_max (fvec_t * s);
  91
  92 /** find the min of a vector
  93
  94   \param s vector to get the min from
  95
  96   \return the value of the maximum of v
  97
  98 */
  99 smpl_t fvec_min (fvec_t * s);
 100
 101 /** find the index of the min of a vector
 102
 103   \param s vector to get the index from
 104
 105   \return the index of the minimum element of v
 106
 107 */
 108 uint_t fvec_min_elem (fvec_t * s);
 109
 110 /** find the index of the max of a vector
 111
 112   \param s vector to get the index from
 113
 114   \return the index of the maximum element of v
 115
 116 */
 117 uint_t fvec_max_elem (fvec_t * s);
 118
 119 /** swap the left and right halves of a vector
 120
 121   This function swaps the left part of the signal with the right part of the
 122 signal. Therefore
 123
 124   \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
 125
 126   becomes
 127
 128   \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
 129
 130   This operation, known as 'fftshift' in the Matlab Signal Processing Toolbox,
 131 can be used before computing the FFT to simplify the phase relationship of the
 132 resulting spectrum. See Amalia de Götzen's paper referred to above.
 133
 134 */
 135 void fvec_shift (fvec_t * v);
 136
 137 /** compute the sum of all elements of a vector
 138
 139   \param v vector to compute the sum of
 140
 141   \return the sum of v
 142
 143 */
 144 smpl_t fvec_sum (fvec_t * v);
 145
 146 /** compute the energy of a vector
 147
 148   This function compute the sum of the squared elements of a vector.
 149
 150   \param v vector to get the energy from
 151
 152   \return the energy of v
 153
 154 */
 155 smpl_t fvec_local_energy (fvec_t * v);
 156
 157 /** compute the High Frequency Content of a vector
 158
 159   The High Frequency Content is defined as \f$ \sum_0^{N-1} (k+1) v[k] \f$.
 160
 161   \param v vector to get the energy from
 162
 163   \return the HFC of v
 164
 165 */
 166 smpl_t fvec_local_hfc (fvec_t * v);
 167
 168 /** computes the p-norm of a vector
 169
 170   Computes the p-norm of a vector for \f$ p = \alpha \f$
 171
 172   \f$ L^p = ||x||_p = (|x_1|^p + |x_2|^p + ... + |x_n|^p ) ^ \frac{1}{p} \f$
 173
 174   If p = 1, the result is the Manhattan distance.
 175
 176   If p = 2, the result is the Euclidean distance.
 177
 178   As p tends towards large values, \f$ L^p \f$ tends towards the maximum of the
 179 input vector.
 180
 181   References:
 182
 183     - <a href="http://en.wikipedia.org/wiki/Lp_space">\f$L^p\f$ space</a> on
 184   Wikipedia
 185
 186   \param v vector to compute norm from
 187   \param p order of the computed norm
 188
 189   \return the p-norm of v
 190
 191 */
 192 smpl_t fvec_alpha_norm (fvec_t * v, smpl_t p);
 193
 194 /**  alpha normalisation
 195
 196   This function divides all elements of a vector by the p-norm as computed by
 197 fvec_alpha_norm().
 198
 199   \param v vector to compute norm from
 200   \param p order of the computed norm
 201
 202 */
 203 void fvec_alpha_normalise (fvec_t * v, smpl_t p);
 204
 205 /** add a constant to each elements of a vector
 206
 207   \param v vector to add constant to
 208   \param c constant to add to v
 209
 210 */
 211 void fvec_add (fvec_t * v, smpl_t c);
 212
 213 /** remove the minimum value of the vector to each elements
 214
 215   \param v vector to remove minimum from
 216
 217 */
 218 void fvec_min_removal (fvec_t * v);
 219
 220 /** compute moving median theshold of a vector
 221
 222   This function computes the moving median threshold value of at the given
 223 position of a vector, taking the median amongs post elements before and up to
 224 pre elements after pos.
 225
 226   \param v input vector
 227   \param tmp temporary vector of length post+1+pre
 228   \param post length of causal part to take before pos
 229   \param pre length of anti-causal part to take after pos
 230   \param pos index to compute threshold for
 231
 232   \return moving median threshold value
 233
 234 */
 235 smpl_t fvec_moving_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre,
 236     uint_t pos);
 237
 238 /** apply adaptive threshold to a vector
 239
 240   For each points at position p of an input vector, this function remove the
 241 moving median threshold computed at p.
 242
 243   \param v input vector
 244   \param tmp temporary vector of length post+1+pre
 245   \param post length of causal part to take before pos
 246   \param pre length of anti-causal part to take after pos
 247
 248 */
 249 void fvec_adapt_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre);
 250
 251 /** returns the median of a vector
 252
 253   The QuickSelect routine is based on the algorithm described in "Numerical
 254 recipes in C", Second Edition, Cambridge University Press, 1992, Section 8.5,
 255 ISBN 0-521-43108-5
 256
 257   This implementation of the QuickSelect routine is based on Nicolas
 258 Devillard's implementation, available at http://ndevilla.free.fr/median/median/
 259 and in the Public Domain.
 260
 261   \param v vector to get median from
 262
 263   \return the median of v
 264
 265 */
 266 smpl_t fvec_median (fvec_t * v);
 267
 268 /** finds exact peak index by quadratic interpolation*/
 269 smpl_t fvec_quadint (fvec_t * x, uint_t pos, uint_t span);
 270
 271 /** Quadratic interpolation using Lagrange polynomial.
 272
 273   Inspired from ``Comparison of interpolation algorithms in real-time sound
 274 processing'', Vladimir Arnost,
 275
 276   \param s0,s1,s2 are 3 consecutive samples of a curve
 277   \param pf is the floating point index [0;2]
 278
 279   \return s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2);
 280
 281 */
 282 smpl_t aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
 283
 284 /** return 1 if v[p] is a peak and positive, 0 otherwise
 285
 286   This function returns 1 if a peak is found at index p in the vector v. The
 287 peak is defined as follows:
 288
 289   - v[p] is positive
 290   - v[p-1] < v[p]
 291   - v[p] > v[p+1]
 292
 293   \param v input vector
 294   \param p position of supposed for peak
 295
 296   \return 1 if a peak is found, 0 otherwise
 297
 298 */
 299 uint_t fvec_peakpick (fvec_t * v, uint_t p);
 300
 301 /** convert frequency bin to midi value */
 302 smpl_t aubio_bintomidi (smpl_t bin, smpl_t samplerate, smpl_t fftsize);
 303
 304 /** convert midi value to frequency bin */
 305 smpl_t aubio_miditobin (smpl_t midi, smpl_t samplerate, smpl_t fftsize);
 306
 307 /** convert frequency bin to frequency (Hz) */
 308 smpl_t aubio_bintofreq (smpl_t bin, smpl_t samplerate, smpl_t fftsize);
 309
 310 /** convert frequency (Hz) to frequency bin */
 311 smpl_t aubio_freqtobin (smpl_t freq, smpl_t samplerate, smpl_t fftsize);
 312
 313 /** convert frequency (Hz) to midi value (0-128) */
 314 smpl_t aubio_freqtomidi (smpl_t freq);
 315
 316 /** convert midi value (0-128) to frequency (Hz) */
 317 smpl_t aubio_miditofreq (smpl_t midi);
 318
 319 /** compute sound pressure level (SPL) in dB
 320
 321   This quantity is often wrongly called 'loudness'.
 322
 323   \param v vector to compute dB SPL from
 324
 325   \return level of v in dB SPL
 326
 327 */
 328 smpl_t aubio_db_spl (fvec_t * v);
 329
 330 /** check if buffer level in dB SPL is under a given threshold
 331
 332   \param v vector to get level from
 333   \param threshold threshold in dB SPL
 334
 335   \return 0 if level is under the given threshold, 1 otherwise
 336
 337 */
 338 uint_t aubio_silence_detection (fvec_t * v, smpl_t threshold);
 339
 340 /** get buffer level if level >= threshold, 1. otherwise
 341
 342   \param v vector to get level from
 343   \param threshold threshold in dB SPL
 344
 345   \return level in dB SPL if level >= threshold, 1. otherwise
 346
 347 */
 348 smpl_t aubio_level_detection (fvec_t * v, smpl_t threshold);
 349
 350 /** compute normalised autocorrelation function
 351
 352   \param input vector to compute autocorrelation from
 353   \param output vector to store autocorrelation function to
 354
 355 */
 356 void aubio_autocorr (fvec_t * input, fvec_t * output);
 357
 358 /** zero-crossing rate (ZCR)
 359
 360   The zero-crossing rate is the number of times a signal changes sign,
 361   divided by the length of this signal.
 362
 363   \param v vector to compute ZCR from
 364
 365   \return zero-crossing rate of v
 366
 367 */
 368 smpl_t aubio_zero_crossing_rate (fvec_t * v);
 369
 370 /** clean up cached memory at the end of program
 371
 372   This function should be used at the end of programs to purge all cached
 373   memory. So far it is only useful to clean FFTW's cache.
 374
 375 */
 376 void aubio_cleanup (void);
 377
 378 #ifdef __cplusplus
 379 }
 380 #endif
 381
 382 #endif
 383