/*
- Copyright (C) 2003 Paul Brossier
+ Copyright (C) 2003-2009 Paul Brossier <piem@aubio.org>
- This program is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2 of the License, or
- (at your option) any later version.
+ This file is part of aubio.
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
+ aubio is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
- You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ aubio is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with aubio. If not, see <http://www.gnu.org/licenses/>.
*/
-/** @file
- * various math functions
- *
- * \todo multichannel (each function should return -or set- an array sized to
- * the number of channel in the input vector)
- *
- * \todo appropriate switches depending on types.h content
+/** \file
+
+ Various math functions
+
+ \example test-mathutils.c
+ \example test-mathutils-window.c
+
*/
#ifndef MATHUTILS_H
#define MATHUTILS_H
-#define PI (M_PI)
-#define TWO_PI (PI*2.)
-
-/* aliases to math.h functions */
-#define EXP expf
-#define COS cosf
-#define SIN sinf
-#define ABS fabsf
-#define POW powf
-#define SQRT sqrtf
-#define LOG10 log10f
-#define LOG logf
-#define FLOOR floorf
-#define TRUNC truncf
-
-/* aliases to complex.h functions */
-#if defined(WIN32)
-/* mingw32 does not know about c*f functions */
-#define EXPC cexp
-/** complex = CEXPC(complex) */
-#define CEXPC cexp
-/** sample = ARGC(complex) */
-#define ARGC carg
-/** sample = ABSC(complex) norm */
-#define ABSC cabs
-/** sample = REAL(complex) */
-#define REAL creal
-/** sample = IMAG(complex) */
-#define IMAG cimag
-#else
-/** sample = EXPC(complex) */
-#define EXPC cexpf
-/** complex = CEXPC(complex) */
-#define CEXPC cexp
-/** sample = ARGC(complex) */
-#define ARGC cargf
-/** sample = ABSC(complex) norm */
-#define ABSC cabsf
-/** sample = REAL(complex) */
-#define REAL crealf
-/** sample = IMAG(complex) */
-#define IMAG cimagf
+#include "fvec.h"
+#include "musicutils.h"
+
+#ifdef __cplusplus
+extern "C" {
#endif
-/* handy shortcuts */
-#define DB2LIN(g) (POW(10.0f,(g)*0.05f))
-#define LIN2DB(v) (20.0f*LOG10(v))
-#define SQR(_a) (_a*_a)
+/** compute the mean of a vector
-#define ELEM_SWAP(a,b) { register smpl_t t=(a);(a)=(b);(b)=t; }
+ \param s vector to compute mean from
+ \return the mean of `v`
-/** Window types
- *
- * inspired from
- *
- * - dafx : http://profs.sci.univr.it/%7Edafx/Final-Papers/ps/Bernardini.ps.gz
- * - freqtweak : http://freqtweak.sf.net/
- * - extace : http://extace.sf.net/
- */
+*/
+smpl_t fvec_mean (fvec_t * s);
-#ifdef __cplusplus
-extern "C" {
-#endif
+/** find the max of a vector
-typedef enum {
- aubio_win_rectangle,
- aubio_win_hamming,
- aubio_win_hanning,
- aubio_win_hanningz,
- aubio_win_blackman,
- aubio_win_blackman_harris,
- aubio_win_gaussian,
- aubio_win_welch,
- aubio_win_parzen
-} aubio_window_type;
-
-/** create window */
-void aubio_window(smpl_t *w, uint_t size, aubio_window_type wintype);
-
-/** principal argument
- *
- * mod(phase+PI,-TWO_PI)+PI
- */
-smpl_t aubio_unwrap2pi (smpl_t phase);
+ \param s vector to get the max from
-/** calculates the mean of a vector
- *
- * \bug mono
- */
-smpl_t vec_mean(fvec_t *s);
-/** returns the max of a vector
- *
- * \bug mono
- */
-smpl_t vec_max(fvec_t *s);
-/** returns the min of a vector
- *
- * \bug mono
- */
-smpl_t vec_min(fvec_t *s);
-/** returns the index of the min of a vector
- *
- * \bug mono
- */
-uint_t vec_min_elem(fvec_t *s);
-/** returns the index of the max of a vector
- *
- * \bug mono
- */
-uint_t vec_max_elem(fvec_t *s);
-/** implement 'fftshift' like function
- *
- * a[0]...,a[n/2],a[n/2+1],...a[n]
- *
- * becomes
- *
- * a[n/2+1],...a[n],a[0]...,a[n/2]
- */
-void vec_shift(fvec_t *s);
-/** returns sum */
-smpl_t vec_sum(fvec_t *s);
-/** returns energy
- *
- * \bug mono
- */
-smpl_t vec_local_energy(fvec_t * f);
-/** returns High Frequency Energy Content
- *
- * \bug mono */
-smpl_t vec_local_hfc(fvec_t * f);
-/** return alpha norm.
- *
- * alpha=2 means normalise variance.
- * alpha=1 means normalise abs value.
- * as alpha goes large, tends to normalisation
- * by max value.
- *
- * \bug should not use POW :(
- */
-smpl_t vec_alpha_norm(fvec_t * DF, smpl_t alpha);
-/* dc(min) removal */
-void vec_dc_removal(fvec_t * mag);
-/** alpha normalisation */
-void vec_alpha_normalise(fvec_t * mag, uint_t alpha);
-
-void vec_add(fvec_t * mag, smpl_t threshold);
-
-void vec_adapt_thres(fvec_t * vec, fvec_t * tmp,
- uint_t win_post, uint_t win_pre);
-/** adaptative thresholding
- *
- * y=fn_thresh(fn,x,post,pre)
- * compute adaptive threshold at each time
- * fn : a function name or pointer, eg 'median'
- * x: signal vector
- * post: window length, causal part
- * pre: window length, anti-causal part
- * Returns:
- * y: signal the same length as x
- *
- * Formerly median_thresh, used compute median over a
- * window of post+pre+1 samples, but now works with any
- * function that takes a vector or matrix and returns a
- * 'representative' value for each column, eg
- * medians=fn_thresh(median,x,8,8)
- * minima=fn_thresh(min,x,8,8)
- * see SPARMS for explanation of post and pre
- */
-smpl_t vec_moving_thres(fvec_t * vec, fvec_t * tmp,
- uint_t win_post, uint_t win_pre, uint_t win_pos);
-
-/** returns the median of the vector
- *
- * This Quickselect routine is based on the algorithm described in
- * "Numerical recipes in C", Second Edition,
- * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
- *
- * This code by Nicolas Devillard - 1998. Public domain,
- * available at http://ndevilla.free.fr/median/median/
- */
-smpl_t vec_median(fvec_t * input);
+ \return the value of the minimum of v
+
+*/
+smpl_t fvec_max (fvec_t * s);
+
+/** find the min of a vector
+
+ \param s vector to get the min from
+
+ \return the value of the maximum of v
+
+*/
+smpl_t fvec_min (fvec_t * s);
+
+/** find the index of the min of a vector
+
+ \param s vector to get the index from
+
+ \return the index of the minimum element of v
+
+*/
+uint_t fvec_min_elem (fvec_t * s);
+
+/** find the index of the max of a vector
+
+ \param s vector to get the index from
+
+ \return the index of the maximum element of v
+
+*/
+uint_t fvec_max_elem (fvec_t * s);
+
+/** swap the left and right halves of a vector
+
+ This function swaps the left part of the signal with the right part of the
+signal. Therefore
+
+ \f$ a[0], a[1], ..., a[\frac{N}{2}], a[\frac{N}{2}+1], ..., a[N-1], a[N] \f$
+
+ becomes
+
+ \f$ a[\frac{N}{2}+1], ..., a[N-1], a[N], a[0], a[1], ..., a[\frac{N}{2}] \f$
+
+ This operation, known as 'fftshift' in the Matlab Signal Processing Toolbox,
+can be used before computing the FFT to simplify the phase relationship of the
+resulting spectrum. See Amalia de Götzen's paper referred to above.
+
+*/
+void fvec_shift (fvec_t * v);
+
+/** compute the sum of all elements of a vector
+
+ \param v vector to compute the sum of
+
+ \return the sum of v
+
+*/
+smpl_t fvec_sum (fvec_t * v);
+
+/** compute the energy of a vector
+
+ This function compute the sum of the squared elements of a vector, normalised
+ by its length.
+
+ \param v vector to get the energy from
+
+ \return the energy of v
+
+*/
+smpl_t fvec_local_energy (fvec_t * v);
+
+/** compute the High Frequency Content of a vector
+
+ The High Frequency Content is defined as \f$ \sum_0^{N-1} (k+1) v[k] \f$.
+
+ \param v vector to get the energy from
+
+ \return the HFC of v
+
+*/
+smpl_t fvec_local_hfc (fvec_t * v);
+
+/** computes the p-norm of a vector
+
+ Computes the p-norm of a vector for \f$ p = \alpha \f$
-/** finds exact maximum position by quadratic interpolation*/
-smpl_t vec_quadint(fvec_t * x,uint_t pos);
+ \f$ L^p = ||x||_p = (|x_1|^p + |x_2|^p + ... + |x_n|^p ) ^ \frac{1}{p} \f$
+
+ If p = 1, the result is the Manhattan distance.
-/** finds exact minimum position by quadratic interpolation*/
-smpl_t vec_quadint_min(fvec_t * x,uint_t pos, uint_t span);
+ If p = 2, the result is the Euclidean distance.
+
+ As p tends towards large values, \f$ L^p \f$ tends towards the maximum of the
+input vector.
+
+ References:
+
+ - <a href="http://en.wikipedia.org/wiki/Lp_space">\f$L^p\f$ space</a> on
+ Wikipedia
+
+ \param v vector to compute norm from
+ \param p order of the computed norm
+
+ \return the p-norm of v
+
+*/
+smpl_t fvec_alpha_norm (fvec_t * v, smpl_t p);
+
+/** alpha normalisation
+
+ This function divides all elements of a vector by the p-norm as computed by
+fvec_alpha_norm().
+
+ \param v vector to compute norm from
+ \param p order of the computed norm
+
+*/
+void fvec_alpha_normalise (fvec_t * v, smpl_t p);
+
+/** add a constant to each elements of a vector
+
+ \param v vector to add constant to
+ \param c constant to add to v
+
+*/
+void fvec_add (fvec_t * v, smpl_t c);
+
+/** remove the minimum value of the vector to each elements
+
+ \param v vector to remove minimum from
+
+*/
+void fvec_min_removal (fvec_t * v);
+
+/** compute moving median threshold of a vector
+
+ This function computes the moving median threshold value of at the given
+position of a vector, taking the median among post elements before and up to
+pre elements after pos.
+
+ \param v input vector
+ \param tmp temporary vector of length post+1+pre
+ \param post length of causal part to take before pos
+ \param pre length of anti-causal part to take after pos
+ \param pos index to compute threshold for
+
+ \return moving median threshold value
+
+*/
+smpl_t fvec_moving_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre,
+ uint_t pos);
+
+/** apply adaptive threshold to a vector
+
+ For each points at position p of an input vector, this function remove the
+moving median threshold computed at p.
+
+ \param v input vector
+ \param tmp temporary vector of length post+1+pre
+ \param post length of causal part to take before pos
+ \param pre length of anti-causal part to take after pos
+
+*/
+void fvec_adapt_thres (fvec_t * v, fvec_t * tmp, uint_t post, uint_t pre);
+
+/** returns the median of a vector
+
+ The QuickSelect routine is based on the algorithm described in "Numerical
+recipes in C", Second Edition, Cambridge University Press, 1992, Section 8.5,
+ISBN 0-521-43108-5
+
+ This implementation of the QuickSelect routine is based on Nicolas
+Devillard's implementation, available at http://ndevilla.free.fr/median/median/
+and in the Public Domain.
+
+ \param v vector to get median from
+
+ \return the median of v
+
+*/
+smpl_t fvec_median (fvec_t * v);
+
+/** finds exact peak index by quadratic interpolation*/
+smpl_t fvec_quadint (fvec_t * x, uint_t pos);
+
+/** finds exact peak index by quadratic interpolation
+
+ See [Quadratic Interpolation of Spectral
+ Peaks](https://ccrma.stanford.edu/~jos/sasp/Quadratic_Peak_Interpolation.html),
+ by Julius O. Smith III
+
+ \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$
+
+ \param x vector to get the interpolated peak position from
+ \param p index of the peak in vector `x`
+ \return \f$ p + p_{frac} \f$ exact peak position of interpolated maximum or minimum
+
+*/
+smpl_t fvec_quadratic_peak_pos (fvec_t * x, uint_t p);
/** Quadratic interpolation using Lagrange polynomial.
- *
- * inspired from ``Comparison of interpolation algorithms in real-time sound
- * processing'', Vladimir Arnost,
- *
- * estimate = s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2);
- * where
- * \param s0,s1,s2 are 3 known points on the curve,
- * \param pf is the floating point index [0;2]
- */
-smpl_t aubio_quadfrac(smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
-
-/** returns 1 if X1 is a peak and positive */
-uint_t vec_peakpick(fvec_t * input, uint_t pos);
-
-smpl_t aubio_bintomidi(smpl_t bin, smpl_t samplerate, smpl_t fftsize);
-smpl_t aubio_miditobin(smpl_t midi, smpl_t samplerate, smpl_t fftsize);
-smpl_t aubio_bintofreq(smpl_t bin, smpl_t samplerate, smpl_t fftsize);
-smpl_t aubio_freqtobin(smpl_t freq, smpl_t samplerate, smpl_t fftsize);
-smpl_t aubio_freqtomidi(smpl_t freq);
-smpl_t aubio_miditofreq(smpl_t midi);
-
-uint_t aubio_silence_detection(fvec_t * ibuf, smpl_t threshold);
-smpl_t aubio_level_detection(fvec_t * ibuf, smpl_t threshold);
-/**
- * calculate normalised autocorrelation function
- */
-void aubio_autocorr(fvec_t * input, fvec_t * output);
-/**
- * clean up cached memory at the end of program
- *
- * use this function at the end of programs to purge all
- * cached memory. so far this function is only used to clean
- * fftw cache.
- */
-void aubio_cleanup(void);
+
+ Inspired from ``Comparison of interpolation algorithms in real-time sound
+processing'', Vladimir Arnost,
+
+ \param s0,s1,s2 are 3 consecutive samples of a curve
+ \param pf is the floating point index [0;2]
+
+ \return \f$ s0 + (pf/2.)*((pf-3.)*s0-2.*(pf-2.)*s1+(pf-1.)*s2); \f$
+
+*/
+smpl_t aubio_quadfrac (smpl_t s0, smpl_t s1, smpl_t s2, smpl_t pf);
+
+/** return 1 if v[p] is a peak and positive, 0 otherwise
+
+ This function returns 1 if a peak is found at index p in the vector v. The
+peak is defined as follows:
+
+ - v[p] is positive
+ - v[p-1] < v[p]
+ - v[p] > v[p+1]
+
+ \param v input vector
+ \param p position of supposed for peak
+
+ \return 1 if a peak is found, 0 otherwise
+
+*/
+uint_t fvec_peakpick (fvec_t * v, uint_t p);
+
+/** return 1 if a is a power of 2, 0 otherwise */
+uint_t aubio_is_power_of_two(uint_t a);
+
+/** return the next power of power of 2 greater than a */
+uint_t aubio_next_power_of_two(uint_t a);
+
+/** compute normalised autocorrelation function
+
+ \param input vector to compute autocorrelation from
+ \param output vector to store autocorrelation function to
+
+*/
+void aubio_autocorr (fvec_t * input, fvec_t * output);
#ifdef __cplusplus
}
projects / aubio.git / blobdiff