2 Copyright (C) 2003-2015 Paul Brossier <piem@aubio.org>
4 This file is part of aubio.
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.
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.
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/>.
23 Various math functions
25 \example test-mathutils.c
26 \example test-mathutils-window.c
30 #ifndef AUBIO_MATHUTILS_H
31 #define AUBIO_MATHUTILS_H
34 #include "musicutils.h"
40 /** compute the mean of a vector
42 \param s vector to compute mean from
43 \return the mean of `v`
46 smpl_t fvec_mean (fvec_t * s);
48 /** find the max of a vector
50 \param s vector to get the max from
52 \return the value of the minimum of v
55 smpl_t fvec_max (fvec_t * s);
57 /** find the min of a vector
59 \param s vector to get the min from
61 \return the value of the maximum of v
64 smpl_t fvec_min (fvec_t * s);
66 /** find the index of the min of a vector
68 \param s vector to get the index from
70 \return the index of the minimum element of v
73 uint_t fvec_min_elem (fvec_t * s);
75 /** find the index of the max of a vector
77 \param s vector to get the index from
79 \return the index of the maximum element of v
82 uint_t fvec_max_elem (fvec_t * s);
84 /** swap the left and right halves of a vector
86 This function swaps the left part of the signal with the right part of the
89 \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
93 \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
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.
100 void fvec_shift (fvec_t * v);
102 /** swap the left and right halves of a vector
104 This function swaps the left part of the signal with the right part of the
107 \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
111 \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
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.
118 void fvec_ishift (fvec_t * v);
120 /** push a new element to the end of a vector, erasing the first element and
123 \param in vector to push to
124 \param new_elem new_element to add at the end of the vector
126 In numpy words, this is equivalent to: in = np.concatenate([in, [new_elem]])[1:]
129 void fvec_push(fvec_t *in, smpl_t new_elem);
131 /** compute the sum of all elements of a vector
133 \param v vector to compute the sum of
138 smpl_t fvec_sum (fvec_t * v);
140 /** compute the High Frequency Content of a vector
142 The High Frequency Content is defined as \f$ \sum_0^{N-1} (k+1) v[k] \f$.
144 \param v vector to get the energy from
149 smpl_t fvec_local_hfc (fvec_t * v);
151 /** computes the p-norm of a vector
153 Computes the p-norm of a vector for \f$ p = \alpha \f$
155 \f$ L^p = ||x||_p = (|x_1|^p + |x_2|^p + ... + |x_n|^p ) ^ \frac{1}{p} \f$
157 If p = 1, the result is the Manhattan distance.
159 If p = 2, the result is the Euclidean distance.
161 As p tends towards large values, \f$ L^p \f$ tends towards the maximum of the
166 - <a href="http://en.wikipedia.org/wiki/Lp_space">\f$L^p\f$ space</a> on
169 \param v vector to compute norm from
170 \param p order of the computed norm
172 \return the p-norm of v
175 smpl_t fvec_alpha_norm (fvec_t * v, smpl_t p);
177 /** alpha normalisation
179 This function divides all elements of a vector by the p-norm as computed by
182 \param v vector to compute norm from
183 \param p order of the computed norm
186 void fvec_alpha_normalise (fvec_t * v, smpl_t p);
188 /** add a constant to each elements of a vector
190 \param v vector to add constant to
191 \param c constant to add to v
194 void fvec_add (fvec_t * v, smpl_t c);
196 /** multiply each elements of a vector by a scalar
198 \param v vector to add constant to
199 \param s constant to scale v with
202 void fvec_mul (fvec_t * v, smpl_t s);
204 /** remove the minimum value of the vector to each elements
206 \param v vector to remove minimum from
209 void fvec_min_removal (fvec_t * v);
211 /** compute moving median threshold of a vector
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.
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
223 \return moving median threshold value
226 smpl_t fvec_moving_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre,
229 /** apply adaptive threshold to a vector
231 For each points at position p of an input vector, this function remove the
232 moving median threshold computed at p.
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
240 void fvec_adapt_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre);
242 /** returns the median of a vector
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,
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.
252 \param v vector to get median from
254 \return the median of v
257 smpl_t fvec_median (fvec_t * v);
259 /** finds exact peak index by quadratic interpolation
261 See [Quadratic Interpolation of Spectral
262 Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
263 by Julius O. Smith III
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$
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
272 smpl_t fvec_quadratic_peak_pos (const fvec_t * x, uint_t p);
274 /** finds magnitude of peak by quadratic interpolation
276 See [Quadratic Interpolation of Spectral
277 Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
278 by Julius O. Smith III
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
285 smpl_t fvec_quadratic_peak_mag (fvec_t * x, smpl_t p);
287 /** Quadratic interpolation using Lagrange polynomial.
289 Inspired from ``Comparison of interpolation algorithms in real-time sound
290 processing'', Vladimir Arnost,
292 \param s0,s1,s2 are 3 consecutive samples of a curve
293 \param pf is the floating point index [0;2]
295 \return \f$ s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2); \f$
298 smpl_t aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
300 /** return 1 if v[p] is a peak and positive, 0 otherwise
302 This function returns 1 if a peak is found at index p in the vector v. The
303 peak is defined as follows:
309 \param v input vector
310 \param p position of supposed for peak
312 \return 1 if a peak is found, 0 otherwise
315 uint_t fvec_peakpick (const fvec_t * v, uint_t p);
317 /** return 1 if a is a power of 2, 0 otherwise */
318 uint_t aubio_is_power_of_two(uint_t a);
320 /** return the next power of power of 2 greater than a */
321 uint_t aubio_next_power_of_two(uint_t a);
323 /** return the log2 factor of the given power of 2 value a */
324 uint_t aubio_power_of_two_order(uint_t a);
326 /** compute normalised autocorrelation function
328 \param input vector to compute autocorrelation from
329 \param output vector to store autocorrelation function to
332 void aubio_autocorr (const fvec_t * input, fvec_t * output);
338 #endif /* AUBIO_MATHUTILS_H */