src/mathutils.h

   1 /*
   2   Copyright (C) 2003-2015 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
  23   Various math functions
  24
  25   \example test-mathutils.c
  26   \example test-mathutils-window.c
  27
  28  */
  29
  30 #ifndef AUBIO_MATHUTILS_H
  31 #define AUBIO_MATHUTILS_H
  32
  33 #include "fvec.h"
  34 #include "musicutils.h"
  35
  36 #ifdef __cplusplus
  37 extern "C" {
  38 #endif
  39
  40 /** compute the mean of a vector
  41
  42   \param s vector to compute mean from
  43   \return the mean of `v`
  44
  45 */
  46 smpl_t fvec_mean (fvec_t * s);
  47
  48 /** find the max of a vector
  49
  50   \param s vector to get the max from
  51
  52   \return the value of the minimum of v
  53
  54 */
  55 smpl_t fvec_max (fvec_t * s);
  56
  57 /** find the min of a vector
  58
  59   \param s vector to get the min from
  60
  61   \return the value of the maximum of v
  62
  63 */
  64 smpl_t fvec_min (fvec_t * s);
  65
  66 /** find the index of the min of a vector
  67
  68   \param s vector to get the index from
  69
  70   \return the index of the minimum element of v
  71
  72 */
  73 uint_t fvec_min_elem (fvec_t * s);
  74
  75 /** find the index of the max of a vector
  76
  77   \param s vector to get the index from
  78
  79   \return the index of the maximum element of v
  80
  81 */
  82 uint_t fvec_max_elem (fvec_t * s);
  83
  84 /** swap the left and right halves of a vector
  85
  86   This function swaps the left part of the signal with the right part of the
  87 signal. Therefore
  88
  89   \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
  90
  91   becomes
  92
  93   \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
  94
  95   This operation, known as 'fftshift' in the Matlab Signal Processing Toolbox,
  96 can be used before computing the FFT to simplify the phase relationship of the
  97 resulting spectrum. See Amalia de Götzen's paper referred to above.
  98
  99 */
 100 void fvec_shift (fvec_t * v);
 101
 102 /** swap the left and right halves of a vector
 103
 104   This function swaps the left part of the signal with the right part of the
 105 signal. Therefore
 106
 107   \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
 108
 109   becomes
 110
 111   \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
 112
 113   This operation, known as 'ifftshift' in the Matlab Signal Processing Toolbox,
 114 can be used after computing the inverse FFT to simplify the phase relationship
 115 of the resulting spectrum. See Amalia de Götzen's paper referred to above.
 116
 117 */
 118 void fvec_ishift (fvec_t * v);
 119
 120 /** push a new element to the end of a vector, erasing the first element and
 121  * sliding all others
 122
 123   \param in vector to push to
 124   \param new_elem new_element to add at the end of the vector
 125
 126   In numpy words, this is equivalent to: in = np.concatenate([in, [new_elem]])[1:]
 127
 128 */
 129 void fvec_push(fvec_t *in, smpl_t new_elem);
 130
 131 /** compute the sum of all elements of a vector
 132
 133   \param v vector to compute the sum of
 134
 135   \return the sum of v
 136
 137 */
 138 smpl_t fvec_sum (fvec_t * v);
 139
 140 /** compute the High Frequency Content of a vector
 141
 142   The High Frequency Content is defined as \f$ \sum_0^{N-1} (k+1) v[k] \f$.
 143
 144   \param v vector to get the energy from
 145
 146   \return the HFC of v
 147
 148 */
 149 smpl_t fvec_local_hfc (fvec_t * v);
 150
 151 /** computes the p-norm of a vector
 152
 153   Computes the p-norm of a vector for \f$ p = \alpha \f$
 154
 155   \f$ L^p = ||x||_p = (|x_1|^p + |x_2|^p + ... + |x_n|^p ) ^ \frac{1}{p} \f$
 156
 157   If p = 1, the result is the Manhattan distance.
 158
 159   If p = 2, the result is the Euclidean distance.
 160
 161   As p tends towards large values, \f$ L^p \f$ tends towards the maximum of the
 162 input vector.
 163
 164   References:
 165
 166     - <a href="http://en.wikipedia.org/wiki/Lp_space">\f$L^p\f$ space</a> on
 167   Wikipedia
 168
 169   \param v vector to compute norm from
 170   \param p order of the computed norm
 171
 172   \return the p-norm of v
 173
 174 */
 175 smpl_t fvec_alpha_norm (fvec_t * v, smpl_t p);
 176
 177 /**  alpha normalisation
 178
 179   This function divides all elements of a vector by the p-norm as computed by
 180 fvec_alpha_norm().
 181
 182   \param v vector to compute norm from
 183   \param p order of the computed norm
 184
 185 */
 186 void fvec_alpha_normalise (fvec_t * v, smpl_t p);
 187
 188 /** add a constant to each elements of a vector
 189
 190   \param v vector to add constant to
 191   \param c constant to add to v
 192
 193 */
 194 void fvec_add (fvec_t * v, smpl_t c);
 195
 196 /** remove the minimum value of the vector to each elements
 197
 198   \param v vector to remove minimum from
 199
 200 */
 201 void fvec_min_removal (fvec_t * v);
 202
 203 /** compute moving median threshold of a vector
 204
 205   This function computes the moving median threshold value of at the given
 206 position of a vector, taking the median among post elements before and up to
 207 pre elements after pos.
 208
 209   \param v input vector
 210   \param tmp temporary vector of length post+1+pre
 211   \param post length of causal part to take before pos
 212   \param pre length of anti-causal part to take after pos
 213   \param pos index to compute threshold for
 214
 215   \return moving median threshold value
 216
 217 */
 218 smpl_t fvec_moving_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre,
 219     uint_t pos);
 220
 221 /** apply adaptive threshold to a vector
 222
 223   For each points at position p of an input vector, this function remove the
 224 moving median threshold computed at p.
 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
 231 */
 232 void fvec_adapt_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre);
 233
 234 /** returns the median of a vector
 235
 236   The QuickSelect routine is based on the algorithm described in "Numerical
 237 recipes in C", Second Edition, Cambridge University Press, 1992, Section 8.5,
 238 ISBN 0-521-43108-5
 239
 240   This implementation of the QuickSelect routine is based on Nicolas
 241 Devillard's implementation, available at http://ndevilla.free.fr/median/median/
 242 and in the Public Domain.
 243
 244   \param v vector to get median from
 245
 246   \return the median of v
 247
 248 */
 249 smpl_t fvec_median (fvec_t * v);
 250
 251 /** finds exact peak index by quadratic interpolation
 252
 253   See [Quadratic Interpolation of Spectral
 254   Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
 255   by Julius O. Smith III
 256
 257   \f$ p_{frac} = \frac{1}{2} \frac {x[p-1] - x[p+1]} {x[p-1] - 2 x[p] + x[p+1]} \in [ -.5, .5] \f$
 258
 259   \param x vector to get the interpolated peak position from
 260   \param p index of the peak in vector `x`
 261   \return \f$ p + p_{frac} \f$ exact peak position of interpolated maximum or minimum
 262
 263 */
 264 smpl_t fvec_quadratic_peak_pos (const fvec_t * x, uint_t p);
 265
 266 /** finds magnitude of peak by quadratic interpolation
 267
 268   See [Quadratic Interpolation of Spectral
 269   Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
 270   by Julius O. Smith III
 271
 272   \param x vector to get the magnitude of the interpolated peak position from
 273   \param p index of the peak in vector `x`
 274   \return magnitude of interpolated peak
 275
 276 */
 277 smpl_t fvec_quadratic_peak_mag (fvec_t * x, smpl_t p);
 278
 279 /** Quadratic interpolation using Lagrange polynomial.
 280
 281   Inspired from ``Comparison of interpolation algorithms in real-time sound
 282 processing'', Vladimir Arnost,
 283
 284   \param s0,s1,s2 are 3 consecutive samples of a curve
 285   \param pf is the floating point index [0;2]
 286
 287   \return \f$ s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2); \f$
 288
 289 */
 290 smpl_t aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
 291
 292 /** return 1 if v[p] is a peak and positive, 0 otherwise
 293
 294   This function returns 1 if a peak is found at index p in the vector v. The
 295 peak is defined as follows:
 296
 297   - v[p] is positive
 298   - v[p-1] < v[p]
 299   - v[p] > v[p+1]
 300
 301   \param v input vector
 302   \param p position of supposed for peak
 303
 304   \return 1 if a peak is found, 0 otherwise
 305
 306 */
 307 uint_t fvec_peakpick (const fvec_t * v, uint_t p);
 308
 309 /** return 1 if a is a power of 2, 0 otherwise */
 310 uint_t aubio_is_power_of_two(uint_t a);
 311
 312 /** return the next power of power of 2 greater than a */
 313 uint_t aubio_next_power_of_two(uint_t a);
 314
 315 /** return the log2 factor of the given power of 2 value a */
 316 uint_t aubio_power_of_two_order(uint_t a);
 317
 318 /** compute normalised autocorrelation function
 319
 320   \param input vector to compute autocorrelation from
 321   \param output vector to store autocorrelation function to
 322
 323 */
 324 void aubio_autocorr (const fvec_t * input, fvec_t * output);
 325
 326 #ifdef __cplusplus
 327 }
 328 #endif
 329
 330 #endif /* AUBIO_MATHUTILS_H */