From 01beb9919b95479d8be040bec74abc5cc67a5e43 Mon Sep 17 00:00:00 2001 From: Mikko Tiusanen Date: Sun, 4 May 2014 01:18:52 +0300 Subject: Initial import. --- src/libspeex/math_approx.h | 332 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 332 insertions(+) create mode 100644 src/libspeex/math_approx.h (limited to 'src/libspeex/math_approx.h') diff --git a/src/libspeex/math_approx.h b/src/libspeex/math_approx.h new file mode 100644 index 0000000..9ca8307 --- /dev/null +++ b/src/libspeex/math_approx.h @@ -0,0 +1,332 @@ +/* Copyright (C) 2002 Jean-Marc Valin */ +/** + @file math_approx.h + @brief Various math approximation functions for Speex +*/ +/* + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + - Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + - Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + - Neither the name of the Xiph.org Foundation nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#ifndef MATH_APPROX_H +#define MATH_APPROX_H + +#include "arch.h" + +#ifndef FIXED_POINT + +#define spx_sqrt sqrt +#define spx_acos acos +#define spx_exp exp +#define spx_cos_norm(x) (cos((.5f*M_PI)*(x))) +#define spx_atan atan + +/** Generate a pseudo-random number */ +static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed) +{ + const unsigned int jflone = 0x3f800000; + const unsigned int jflmsk = 0x007fffff; + union {int i; float f;} ran; + *seed = 1664525 * *seed + 1013904223; + ran.i = jflone | (jflmsk & *seed); + ran.f -= 1.5; + return 3.4642*std*ran.f; +} + + +#endif + + +static inline spx_int16_t spx_ilog2(spx_uint32_t x) +{ + int r=0; + if (x>=(spx_int32_t)65536) + { + x >>= 16; + r += 16; + } + if (x>=256) + { + x >>= 8; + r += 8; + } + if (x>=16) + { + x >>= 4; + r += 4; + } + if (x>=4) + { + x >>= 2; + r += 2; + } + if (x>=2) + { + r += 1; + } + return r; +} + +static inline spx_int16_t spx_ilog4(spx_uint32_t x) +{ + int r=0; + if (x>=(spx_int32_t)65536) + { + x >>= 16; + r += 8; + } + if (x>=256) + { + x >>= 8; + r += 4; + } + if (x>=16) + { + x >>= 4; + r += 2; + } + if (x>=4) + { + r += 1; + } + return r; +} + +#ifdef FIXED_POINT + +/** Generate a pseudo-random number */ +static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed) +{ + spx_word32_t res; + *seed = 1664525 * *seed + 1013904223; + res = MULT16_16(EXTRACT16(SHR32(*seed,16)),std); + return EXTRACT16(PSHR32(SUB32(res, SHR32(res, 3)),14)); +} + +/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */ +/*#define C0 3634 +#define C1 21173 +#define C2 -12627 +#define C3 4215*/ + +/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */ +#define C0 3634 +#define C1 21173 +#define C2 -12627 +#define C3 4204 + +static inline spx_word16_t spx_sqrt(spx_word32_t x) +{ + int k; + spx_word32_t rt; + k = spx_ilog4(x)-6; + x = VSHR32(x, (k<<1)); + rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3))))))); + rt = VSHR32(rt,7-k); + return rt; +} + +/* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */ + + +#define A1 16469 +#define A2 2242 +#define A3 1486 + +static inline spx_word16_t spx_acos(spx_word16_t x) +{ + int s=0; + spx_word16_t ret; + spx_word16_t sq; + if (x<0) + { + s=1; + x = NEG16(x); + } + x = SUB16(16384,x); + + x = x >> 1; + sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3)))))); + ret = spx_sqrt(SHL32(EXTEND32(sq),13)); + + /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/ + if (s) + ret = SUB16(25736,ret); + return ret; +} + + +#define K1 8192 +#define K2 -4096 +#define K3 340 +#define K4 -10 + +static inline spx_word16_t spx_cos(spx_word16_t x) +{ + spx_word16_t x2; + + if (x<12868) + { + x2 = MULT16_16_P13(x,x); + return ADD32(K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); + } else { + x = SUB16(25736,x); + x2 = MULT16_16_P13(x,x); + return SUB32(-K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); + } +} + +#define L1 32767 +#define L2 -7651 +#define L3 8277 +#define L4 -626 + +static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x) +{ + spx_word16_t x2; + + x2 = MULT16_16_P15(x,x); + return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2)))))))); +} + +static inline spx_word16_t spx_cos_norm(spx_word32_t x) +{ + x = x&0x0001ffff; + if (x>SHL32(EXTEND32(1), 16)) + x = SUB32(SHL32(EXTEND32(1), 17),x); + if (x&0x00007fff) + { + if (x14) + return 0x7fffffff; + else if (integer < -15) + return 0; + frac = SHL16(x-SHL16(integer,11),3); + frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac)))))); + return VSHR32(EXTEND32(frac), -integer-2); +} + +/* Input in Q11 format, output in Q16 */ +static inline spx_word32_t spx_exp(spx_word16_t x) +{ + if (x>21290) + return 0x7fffffff; + else if (x<-21290) + return 0; + else + return spx_exp2(MULT16_16_P14(23637,x)); +} +#define M1 32767 +#define M2 -21 +#define M3 -11943 +#define M4 4936 + +static inline spx_word16_t spx_atan01(spx_word16_t x) +{ + return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x))))))); +} + +#undef M1 +#undef M2 +#undef M3 +#undef M4 + +/* Input in Q15, output in Q14 */ +static inline spx_word16_t spx_atan(spx_word32_t x) +{ + if (x <= 32767) + { + return SHR16(spx_atan01(x),1); + } else { + int e = spx_ilog2(x); + if (e>=29) + return 25736; + x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14))); + return SUB16(25736, SHR16(spx_atan01(x),1)); + } +} +#else + +#ifndef M_PI +#define M_PI 3.14159265358979323846 /* pi */ +#endif + +#define C1 0.9999932946f +#define C2 -0.4999124376f +#define C3 0.0414877472f +#define C4 -0.0012712095f + + +#define SPX_PI_2 1.5707963268 +static inline spx_word16_t spx_cos(spx_word16_t x) +{ + if (x