src/mathutils.h

   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 /** @file
  21  *  various math functions
  22  *
  23  *  \todo multichannel (each function should return -or set- an array sized to
  24  *  the number of channel in the input vector)
  25  *
  26  *  \todo appropriate switches depending on types.h content
  27  */
  28
  29 #ifndef MATHUTILS_H
  30 #define MATHUTILS_H
  31
  32 #define PI                              (M_PI)
  33 #define TWO_PI          (PI*2.)
  34
  35 /* aliases to math.h functions */
  36 #define EXP                             expf
  37 #define COS                             cosf
  38 #define SIN                             sinf
  39 #define ABS                             fabsf
  40 #define POW                             powf
  41 #define SQRT                    sqrtf
  42 #define LOG10                   log10f
  43 #define LOG                       logf
  44 #define FLOOR                   floorf
  45 #define TRUNC                   truncf
  46
  47 /* aliases to complex.h functions */
  48 #if defined(WIN32)
  49 /* mingw32 does not know about c*f functions */
  50 #define EXPC                    cexp
  51 /** complex = CEXPC(complex) */
  52 #define CEXPC                   cexp
  53 /** sample = ARGC(complex) */
  54 #define ARGC                    carg
  55 /** sample = ABSC(complex) norm */
  56 #define ABSC                    cabs
  57 /** sample = REAL(complex) */
  58 #define REAL                    creal
  59 /** sample = IMAG(complex) */
  60 #define IMAG                    cimag
  61 #else
  62 /** sample = EXPC(complex) */
  63 #define EXPC                    cexpf
  64 /** complex = CEXPC(complex) */
  65 #define CEXPC                   cexp
  66 /** sample = ARGC(complex) */
  67 #define ARGC                    cargf
  68 /** sample = ABSC(complex) norm */
  69 #define ABSC                    cabsf
  70 /** sample = REAL(complex) */
  71 #define REAL                    crealf
  72 /** sample = IMAG(complex) */
  73 #define IMAG                    cimagf
  74 #endif
  75
  76 /* handy shortcuts */
  77 #define DB2LIN(g) (POW(10.0f,(g)*0.05f))
  78 #define LIN2DB(v) (20.0f*LOG10(v))
  79 #define SQR(_a)         (_a*_a)
  80
  81 #define ELEM_SWAP(a,b) { register smpl_t t=(a);(a)=(b);(b)=t; }
  82
  83 /** Window types
  84  *
  85  * inspired from
  86  *
  87  *  - dafx : http://profs.sci.univr.it/%7Edafx/Final-Papers/ps/Bernardini.ps.gz
  88  *  - freqtweak : http://freqtweak.sf.net/
  89  *  - extace : http://extace.sf.net/
  90  */
  91
  92 #ifdef __cplusplus
  93 extern "C" {
  94 #endif
  95
  96 typedef enum {
  97         rectangle,
  98         hamming,
  99         hanning,
 100         hanningz,
 101         blackman,
 102         blackman_harris,
 103         gaussian,
 104         welch,
 105         parzen
 106 } window_type_t;
 107
 108 /** create window */
 109 void window(smpl_t *w, uint_t size, window_type_t wintype);
 110
 111 /** principal argument
 112  *
 113  * mod(phase+PI,-TWO_PI)+PI
 114  */
 115 smpl_t unwrap2pi (smpl_t phase);
 116
 117 /** calculates the mean of a vector
 118  *
 119  * \bug mono
 120  */
 121 smpl_t vec_mean(fvec_t *s);
 122 /** returns the max of a vector
 123  *
 124  * \bug mono
 125  */
 126 smpl_t vec_max(fvec_t *s);
 127 /** returns the min of a vector
 128  *
 129  * \bug mono
 130  */
 131 smpl_t vec_min(fvec_t *s);
 132 /** returns the index of the min of a vector
 133  *
 134  * \bug mono
 135  */
 136 uint_t vec_min_elem(fvec_t *s);
 137 /** returns the index of the max of a vector
 138  *
 139  * \bug mono
 140  */
 141 uint_t vec_max_elem(fvec_t *s);
 142 /** implement 'fftshift' like function
 143  *
 144  * a[0]...,a[n/2],a[n/2+1],...a[n]
 145  *
 146  * becomes
 147  *
 148  * a[n/2+1],...a[n],a[0]...,a[n/2]
 149  */
 150 void vec_shift(fvec_t *s);
 151 /** returns sum */
 152 smpl_t vec_sum(fvec_t *s);
 153 /** returns energy
 154  *
 155  * \bug mono
 156  */
 157 smpl_t vec_local_energy(fvec_t * f);
 158 /** returns High Frequency Energy Content
 159  *
 160  * \bug mono */
 161 smpl_t vec_local_hfc(fvec_t * f);
 162 /** return alpha norm.
 163  *
 164  * alpha=2 means normalise variance.
 165  * alpha=1 means normalise abs value.
 166  * as alpha goes large, tends to normalisation
 167  * by max value.
 168  *
 169  * \bug should not use POW :(
 170  */
 171 smpl_t vec_alpha_norm(fvec_t * DF, smpl_t alpha);
 172 /*  dc(min) removal */
 173 void vec_dc_removal(fvec_t * mag);
 174 /**  alpha normalisation */
 175 void vec_alpha_normalise(fvec_t * mag, uint_t alpha);
 176
 177 void vec_add(fvec_t * mag, smpl_t threshold);
 178
 179 void vec_adapt_thres(fvec_t * vec, fvec_t * tmp,
 180     uint_t post, uint_t pre);
 181 /**  adaptative thresholding
 182  *
 183  * y=fn_thresh(fn,x,post,pre)
 184  * compute adaptive threshold at each time
 185  *    fn : a function name or pointer, eg 'median'
 186  *    x:   signal vector
 187  *    post: window length, causal part
 188  *    pre: window length, anti-causal part
 189  * Returns:
 190  *    y:   signal the same length as x
 191  *
 192  * Formerly median_thresh, used compute median over a
 193  * window of post+pre+1 samples, but now works with any
 194  * function that takes a vector or matrix and returns a
 195  * 'representative' value for each column, eg
 196  *    medians=fn_thresh(median,x,8,8)
 197  *    minima=fn_thresh(min,x,8,8)
 198  * see SPARMS for explanation of post and pre
 199  */
 200 smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmp,
 201     uint_t post, uint_t pre, uint_t pos);
 202
 203 /** returns the median of the vector
 204  *
 205  *  This Quickselect routine is based on the algorithm described in
 206  *  "Numerical recipes in C", Second Edition,
 207  *  Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
 208  *
 209  *  This code by Nicolas Devillard - 1998. Public domain,
 210  *  available at http://ndevilla.free.fr/median/median/
 211  */
 212 smpl_t vec_median(fvec_t * input);
 213
 214 /** finds exact peak position by quadratic interpolation*/
 215 smpl_t vec_quadint(fvec_t * x,uint_t pos);
 216
 217 /** Quadratic interpolation using Lagrange polynomial.
 218  *
 219  * inspired from ``Comparison of interpolation algorithms in real-time sound
 220  * processing'', Vladimir Arnost,
 221  *
 222  * estimate = s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2);
 223  *    where
 224  *    \param s0,s1,s2 are 3 known points on the curve,
 225  *    \param pf is the floating point index [0;2]
 226  */
 227 smpl_t quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
 228
 229 /** returns 1 if X1 is a peak and positive */
 230 uint_t vec_peakpick(fvec_t * input, uint_t pos);
 231
 232 smpl_t bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize);
 233 smpl_t bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize);
 234 smpl_t freqtomidi(smpl_t freq);
 235
 236 uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold);
 237 smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold);
 238
 239 #ifdef __cplusplus
 240 }
 241 #endif
 242
 243 #endif
 244