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 /** multiply each elements of a vector by a scalar
 197
 198   \param v vector to add constant to
 199   \param s constant to scale v with
 200
 201 */
 202 void fvec_mul (fvec_t * v, smpl_t s);
 203
 204 /** remove the minimum value of the vector to each elements
 205
 206   \param v vector to remove minimum from
 207
 208 */
 209 void fvec_min_removal (fvec_t * v);
 210
 211 /** compute moving median threshold of a vector
 212
 213   This function computes the moving median threshold value of at the given
 214 position of a vector, taking the median among post elements before and up to
 215 pre elements after pos.
 216
 217   \param v input vector
 218   \param tmp temporary vector of length post+1+pre
 219   \param post length of causal part to take before pos
 220   \param pre length of anti-causal part to take after pos
 221   \param pos index to compute threshold for
 222
 223   \return moving median threshold value
 224
 225 */
 226 smpl_t fvec_moving_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre,
 227     uint_t pos);
 228
 229 /** apply adaptive threshold to a vector
 230
 231   For each points at position p of an input vector, this function remove the
 232 moving median threshold computed at p.
 233
 234   \param v input vector
 235   \param tmp temporary vector of length post+1+pre
 236   \param post length of causal part to take before pos
 237   \param pre length of anti-causal part to take after pos
 238
 239 */
 240 void fvec_adapt_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre);
 241
 242 /** returns the median of a vector
 243
 244   The QuickSelect routine is based on the algorithm described in "Numerical
 245 recipes in C", Second Edition, Cambridge University Press, 1992, Section 8.5,
 246 ISBN 0-521-43108-5
 247
 248   This implementation of the QuickSelect routine is based on Nicolas
 249 Devillard's implementation, available at http://ndevilla.free.fr/median/median/
 250 and in the Public Domain.
 251
 252   \param v vector to get median from
 253
 254   \return the median of v
 255
 256 */
 257 smpl_t fvec_median (fvec_t * v);
 258
 259 /** finds exact peak index by quadratic interpolation
 260
 261   See [Quadratic Interpolation of Spectral
 262   Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
 263   by Julius O. Smith III
 264
 265   \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$
 266
 267   \param x vector to get the interpolated peak position from
 268   \param p index of the peak in vector `x`
 269   \return \f$ p + p_{frac} \f$ exact peak position of interpolated maximum or minimum
 270
 271 */
 272 smpl_t fvec_quadratic_peak_pos (const fvec_t * x, uint_t p);
 273
 274 /** finds magnitude of peak by quadratic interpolation
 275
 276   See [Quadratic Interpolation of Spectral
 277   Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
 278   by Julius O. Smith III
 279
 280   \param x vector to get the magnitude of the interpolated peak position from
 281   \param p index of the peak in vector `x`
 282   \return magnitude of interpolated peak
 283
 284 */
 285 smpl_t fvec_quadratic_peak_mag (fvec_t * x, smpl_t p);
 286
 287 /** Quadratic interpolation using Lagrange polynomial.
 288
 289   Inspired from ``Comparison of interpolation algorithms in real-time sound
 290 processing'', Vladimir Arnost,
 291
 292   \param s0,s1,s2 are 3 consecutive samples of a curve
 293   \param pf is the floating point index [0;2]
 294
 295   \return \f$ s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2); \f$
 296
 297 */
 298 smpl_t aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
 299
 300 /** return 1 if v[p] is a peak and positive, 0 otherwise
 301
 302   This function returns 1 if a peak is found at index p in the vector v. The
 303 peak is defined as follows:
 304
 305   - v[p] is positive
 306   - v[p-1] < v[p]
 307   - v[p] > v[p+1]
 308
 309   \param v input vector
 310   \param p position of supposed for peak
 311
 312   \return 1 if a peak is found, 0 otherwise
 313
 314 */
 315 uint_t fvec_peakpick (const fvec_t * v, uint_t p);
 316
 317 /** return 1 if a is a power of 2, 0 otherwise */
 318 uint_t aubio_is_power_of_two(uint_t a);
 319
 320 /** return the next power of power of 2 greater than a */
 321 uint_t aubio_next_power_of_two(uint_t a);
 322
 323 /** return the log2 factor of the given power of 2 value a */
 324 uint_t aubio_power_of_two_order(uint_t a);
 325
 326 /** compute normalised autocorrelation function
 327
 328   \param input vector to compute autocorrelation from
 329   \param output vector to store autocorrelation function to
 330
 331 */
 332 void aubio_autocorr (const fvec_t * input, fvec_t * output);
 333
 334 #ifdef __cplusplus
 335 }
 336 #endif
 337
 338 #endif /* AUBIO_MATHUTILS_H */