-/* Copyright (C) 2002-2008 Jean-Marc Valin */
+/* Copyright (c) 2002-2008 Jean-Marc Valin
+ Copyright (c) 2007-2008 CSIRO
+ Copyright (c) 2007-2009 Xiph.Org Foundation
+ Written by Jean-Marc Valin */
/**
@file mathops.h
@brief Various math functions
#include "entcode.h"
#include "os_support.h"
-#ifndef OVERRIDE_CELT_ILOG2
-/** Integer log in base2. Undefined for zero and negative numbers */
-static inline celt_int16_t celt_ilog2(celt_word32_t x)
-{
- celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
- return EC_ILOG(x)-1;
-}
-#endif
-
-#ifndef OVERRIDE_FIND_MAX16
-static inline int find_max16(celt_word16_t *x, int len)
-{
- celt_word16_t max_corr=-VERY_LARGE16;
- int i, id = 0;
- for (i=0;i<len;i++)
- {
- if (x[i] > max_corr)
- {
- id = i;
- max_corr = x[i];
- }
- }
- return id;
-}
-#endif
-
-#ifndef OVERRIDE_FIND_MAX32
-static inline int find_max32(celt_word32_t *x, int len)
-{
- celt_word32_t max_corr=-VERY_LARGE32;
- int i, id = 0;
- for (i=0;i<len;i++)
- {
- if (x[i] > max_corr)
- {
- id = i;
- max_corr = x[i];
- }
- }
- return id;
-}
-#endif
-
-#define FRAC_MUL16(a,b) ((16384+((celt_int32_t)(celt_int16_t)(a)*(celt_int16_t)(b)))>>15)
-static inline celt_int16_t bitexact_cos(celt_int16_t x)
-{
- celt_int32_t tmp;
- celt_int16_t x2;
- tmp = (4096+((celt_int32_t)(x)*(x)))>>13;
- if (tmp > 32767)
- tmp = 32767;
- x2 = tmp;
- x2 = (32767-x2) + FRAC_MUL16(x2, (-7651 + FRAC_MUL16(x2, (8277 + FRAC_MUL16(-626, x2)))));
- if (x2 > 32766)
- x2 = 32766;
- return 1+x2;
-}
+/* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
+#define FRAC_MUL16(a,b) ((16384+((celt_int32)(celt_int16)(a)*(celt_int16)(b)))>>15)
+unsigned isqrt32(celt_uint32 _val);
#ifndef FIXED_POINT
#define celt_sqrt(x) ((float)sqrt(x))
-#define celt_psqrt(x) ((float)sqrt(x))
#define celt_rsqrt(x) (1.f/celt_sqrt(x))
+#define celt_rsqrt_norm(x) (celt_rsqrt(x))
#define celt_acos acos
#define celt_exp exp
-#define celt_cos_norm(x) (cos((.5f*M_PI)*(x)))
+#define celt_cos_norm(x) ((float)cos((.5f*M_PI)*(x)))
#define celt_atan atan
#define celt_rcp(x) (1.f/(x))
#define celt_div(a,b) ((a)/(b))
+#define frac_div32(a,b) ((float)(a)/(b))
#ifdef FLOAT_APPROX
float frac;
union {
float f;
- celt_uint32_t i;
+ celt_uint32 i;
} in;
in.f = x;
integer = (in.i>>23)-127;
in.i -= integer<<23;
- frac = in.f - 1.5;
- /* -0.41446 0.96093 -0.33981 0.15600 */
- frac = -0.41446 + frac*(0.96093 + frac*(-0.33981 + frac*0.15600));
- /*printf ("%f %d %f %f %f %f\n", x, integer, in.f, frac, log2(x), 1+integer+frac);*/
+ frac = in.f - 1.5f;
+ frac = -0.41445418f + frac*(0.95909232f
+ + frac*(-0.33951290f + frac*0.16541097f));
return 1+integer+frac;
}
float frac;
union {
float f;
- celt_uint32_t i;
+ celt_uint32 i;
} res;
integer = floor(x);
if (integer < -50)
return 0;
frac = x-integer;
/* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
- res.f = 1.f + frac * (0.696147f + frac * (0.224411f + 0.079442f*frac));
+ res.f = 0.99992522f + frac * (0.69583354f
+ + frac * (0.22606716f + 0.078024523f*frac));
res.i = (res.i + (integer<<23)) & 0x7fffffff;
- /*printf ("%f %f %f %f\n", x, frac, exp2(x), res.f);*/
return res.f;
}
#else
-#define celt_log2(x) (1.442695040888963387*log(x))
-#define celt_exp2(x) (exp(0.6931471805599453094*(x)))
+#define celt_log2(x) ((float)(1.442695040888963387*log(x)))
+#define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
#endif
#endif
#include "os_support.h"
+#ifndef OVERRIDE_CELT_ILOG2
+/** Integer log in base2. Undefined for zero and negative numbers */
+static inline celt_int16 celt_ilog2(celt_int32 x)
+{
+ celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
+ return EC_ILOG(x)-1;
+}
+#endif
+
+
#ifndef OVERRIDE_CELT_MAXABS16
-static inline celt_word16_t celt_maxabs16(celt_word16_t *x, int len)
+static inline celt_word16 celt_maxabs16(celt_word16 *x, int len)
{
int i;
- celt_word16_t maxval = 0;
+ celt_word16 maxval = 0;
for (i=0;i<len;i++)
maxval = MAX16(maxval, ABS16(x[i]));
return maxval;
#endif
/** Integer log in base2. Defined for zero, but not for negative numbers */
-static inline celt_int16_t celt_zlog2(celt_word32_t x)
+static inline celt_int16 celt_zlog2(celt_word32 x)
{
return x <= 0 ? 0 : celt_ilog2(x);
}
-/** Reciprocal sqrt approximation (Q30 input, Q0 output or equivalent) */
-static inline celt_word32_t celt_rsqrt(celt_word32_t x)
-{
- int k;
- celt_word16_t n;
- celt_word32_t rt;
- const celt_word16_t C[5] = {23126, -11496, 9812, -9097, 4100};
- k = celt_ilog2(x)>>1;
- x = VSHR32(x, (k-7)<<1);
- /* Range of n is [-16384,32767] */
- n = x-32768;
- rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
- MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
- rt = VSHR32(rt,k);
- return rt;
-}
-
-/** Sqrt approximation (QX input, QX/2 output) */
-static inline celt_word32_t celt_sqrt(celt_word32_t x)
-{
- int k;
- celt_word16_t n;
- celt_word32_t rt;
- const celt_word16_t C[5] = {23174, 11584, -3011, 1570, -557};
- if (x==0)
- return 0;
- k = (celt_ilog2(x)>>1)-7;
- x = VSHR32(x, (k<<1));
- n = x-32768;
- rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
- MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
- rt = VSHR32(rt,7-k);
- return rt;
-}
-
-/** Sqrt approximation (QX input, QX/2 output) that assumes that the input is
- strictly positive */
-static inline celt_word32_t celt_psqrt(celt_word32_t x)
-{
- int k;
- celt_word16_t n;
- celt_word32_t rt;
- const celt_word16_t C[5] = {23174, 11584, -3011, 1570, -557};
- k = (celt_ilog2(x)>>1)-7;
- x = VSHR32(x, (k<<1));
- n = x-32768;
- rt = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
- MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
- rt = VSHR32(rt,7-k);
- return rt;
-}
-
-#define L1 32767
-#define L2 -7651
-#define L3 8277
-#define L4 -626
+celt_word16 celt_rsqrt_norm(celt_word32 x);
-static inline celt_word16_t _celt_cos_pi_2(celt_word16_t x)
-{
- celt_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
- ))))))));
-}
+celt_word32 celt_sqrt(celt_word32 x);
-#undef L1
-#undef L2
-#undef L3
-#undef L4
+celt_word16 celt_cos_norm(celt_word32 x);
-static inline celt_word16_t celt_cos_norm(celt_word32_t x)
-{
- x = x&0x0001ffff;
- if (x>SHL32(EXTEND32(1), 16))
- x = SUB32(SHL32(EXTEND32(1), 17),x);
- if (x&0x00007fff)
- {
- if (x<SHL32(EXTEND32(1), 15))
- {
- return _celt_cos_pi_2(EXTRACT16(x));
- } else {
- return NEG32(_celt_cos_pi_2(EXTRACT16(65536-x)));
- }
- } else {
- if (x&0x0000ffff)
- return 0;
- else if (x&0x0001ffff)
- return -32767;
- else
- return 32767;
- }
-}
-static inline celt_word16_t celt_log2(celt_word32_t x)
+static inline celt_word16 celt_log2(celt_word32 x)
{
int i;
- celt_word16_t n, frac;
- /*-0.41446 0.96093 -0.33981 0.15600 */
- const celt_word16_t C[4] = {-6791, 7872, -1392, 319};
+ celt_word16 n, frac;
+ /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
+ 0.15530808010959576, -0.08556153059057618 */
+ static const celt_word16 C[5] = {-6801+(1<<13-DB_SHIFT), 15746, -5217, 2545, -1401};
if (x==0)
return -32767;
i = celt_ilog2(x);
n = VSHR32(x,i-15)-32768-16384;
- frac = ADD16(C[0], MULT16_16_Q14(n, ADD16(C[1], MULT16_16_Q14(n, ADD16(C[2], MULT16_16_Q14(n, (C[3])))))));
- /*printf ("%d %d %d %d\n", x, n, ret, SHL16(i-13,8)+SHR16(ret,14-8));*/
- return SHL16(i-13,8)+SHR16(frac,14-8);
+ frac = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2], MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, C[4]))))))));
+ return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
}
/*
K2 = 3-4*log(2)
K3 = 3*log(2) - 2
*/
-#define D0 16384
-#define D1 11356
-#define D2 3726
-#define D3 1301
+#define D0 16383
+#define D1 22804
+#define D2 14819
+#define D3 10204
/** Base-2 exponential approximation (2^x). (Q11 input, Q16 output) */
-static inline celt_word32_t celt_exp2(celt_word16_t x)
+static inline celt_word32 celt_exp2(celt_word16 x)
{
int integer;
- celt_word16_t frac;
+ celt_word16 frac;
integer = SHR16(x,11);
if (integer>14)
return 0x7f000000;
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))))));
+ frac = ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
return VSHR32(EXTEND32(frac), -integer-2);
}
-/** Reciprocal approximation (Q15 input, Q16 output) */
-static inline celt_word32_t celt_rcp(celt_word32_t x)
-{
- int i;
- celt_word16_t n, frac;
- const celt_word16_t C[5] = {21848, -7251, 2403, -934, 327};
- celt_assert2(x>0, "celt_rcp() only defined for positive values");
- i = celt_ilog2(x);
- n = VSHR32(x,i-16)-SHL32(EXTEND32(3),15);
- frac = ADD16(C[0], MULT16_16_Q15(n, ADD16(C[1], MULT16_16_Q15(n, ADD16(C[2],
- MULT16_16_Q15(n, ADD16(C[3], MULT16_16_Q15(n, (C[4])))))))));
- return VSHR32(EXTEND32(frac),i-16);
-}
+celt_word32 celt_rcp(celt_word32 x);
-#define celt_div(a,b) MULT32_32_Q31((celt_word32_t)(a),celt_rcp(b))
+#define celt_div(a,b) MULT32_32_Q31((celt_word32)(a),celt_rcp(b))
+celt_word32 frac_div32(celt_word32 a, celt_word32 b);
#define M1 32767
#define M2 -21
#define M3 -11943
#define M4 4936
-static inline celt_word16_t celt_atan01(celt_word16_t x)
+/* Atan approximation using a 4th order polynomial. Input is in Q15 format
+ and normalized by pi/4. Output is in Q15 format */
+static inline celt_word16 celt_atan01(celt_word16 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 M3
#undef M4
-static inline celt_word16_t celt_atan2p(celt_word16_t y, celt_word16_t x)
+/* atan2() approximation valid for positive input values */
+static inline celt_word16 celt_atan2p(celt_word16 y, celt_word16 x)
{
if (y < x)
{
- celt_word32_t arg;
+ celt_word32 arg;
arg = celt_div(SHL32(EXTEND32(y),15),x);
if (arg >= 32767)
arg = 32767;
return SHR16(celt_atan01(EXTRACT16(arg)),1);
} else {
- celt_word32_t arg;
+ celt_word32 arg;
arg = celt_div(SHL32(EXTEND32(x),15),y);
if (arg >= 32767)
arg = 32767;
xiph / opus.git / blobdiff