media/libopus/celt/mathops.h

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

michael@0 1 /* Copyright (c) 2002-2008 Jean-Marc Valin
michael@0 2 Copyright (c) 2007-2008 CSIRO
michael@0 3 Copyright (c) 2007-2009 Xiph.Org Foundation
michael@0 4 Written by Jean-Marc Valin */
michael@0 5 /**
michael@0 6 @file mathops.h
michael@0 7 @brief Various math functions
michael@0 8 */
michael@0 9 /*
michael@0 10 Redistribution and use in source and binary forms, with or without
michael@0 11 modification, are permitted provided that the following conditions
michael@0 12 are met:
michael@0 13
michael@0 14 - Redistributions of source code must retain the above copyright
michael@0 15 notice, this list of conditions and the following disclaimer.
michael@0 16
michael@0 17 - Redistributions in binary form must reproduce the above copyright
michael@0 18 notice, this list of conditions and the following disclaimer in the
michael@0 19 documentation and/or other materials provided with the distribution.
michael@0 20
michael@0 21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
michael@0 22 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
michael@0 23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
michael@0 24 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
michael@0 25 OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
michael@0 26 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
michael@0 27 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
michael@0 28 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
michael@0 29 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
michael@0 30 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
michael@0 31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0 32 */
michael@0 33
michael@0 34 #ifndef MATHOPS_H
michael@0 35 #define MATHOPS_H
michael@0 36
michael@0 37 #include "arch.h"
michael@0 38 #include "entcode.h"
michael@0 39 #include "os_support.h"
michael@0 40
michael@0 41 /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is important */
michael@0 42 #define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>15)
michael@0 43
michael@0 44 unsigned isqrt32(opus_uint32 _val);
michael@0 45
michael@0 46 #ifndef OVERRIDE_CELT_MAXABS16
michael@0 47 static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len)
michael@0 48 {
michael@0 49 int i;
michael@0 50 opus_val16 maxval = 0;
michael@0 51 opus_val16 minval = 0;
michael@0 52 for (i=0;i<len;i++)
michael@0 53 {
michael@0 54 maxval = MAX16(maxval, x[i]);
michael@0 55 minval = MIN16(minval, x[i]);
michael@0 56 }
michael@0 57 return MAX32(EXTEND32(maxval),-EXTEND32(minval));
michael@0 58 }
michael@0 59 #endif
michael@0 60
michael@0 61 #ifndef OVERRIDE_CELT_MAXABS32
michael@0 62 #ifdef FIXED_POINT
michael@0 63 static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len)
michael@0 64 {
michael@0 65 int i;
michael@0 66 opus_val32 maxval = 0;
michael@0 67 opus_val32 minval = 0;
michael@0 68 for (i=0;i<len;i++)
michael@0 69 {
michael@0 70 maxval = MAX32(maxval, x[i]);
michael@0 71 minval = MIN32(minval, x[i]);
michael@0 72 }
michael@0 73 return MAX32(maxval, -minval);
michael@0 74 }
michael@0 75 #else
michael@0 76 #define celt_maxabs32(x,len) celt_maxabs16(x,len)
michael@0 77 #endif
michael@0 78 #endif
michael@0 79
michael@0 80
michael@0 81 #ifndef FIXED_POINT
michael@0 82
michael@0 83 #define PI 3.141592653f
michael@0 84 #define celt_sqrt(x) ((float)sqrt(x))
michael@0 85 #define celt_rsqrt(x) (1.f/celt_sqrt(x))
michael@0 86 #define celt_rsqrt_norm(x) (celt_rsqrt(x))
michael@0 87 #define celt_cos_norm(x) ((float)cos((.5f*PI)*(x)))
michael@0 88 #define celt_rcp(x) (1.f/(x))
michael@0 89 #define celt_div(a,b) ((a)/(b))
michael@0 90 #define frac_div32(a,b) ((float)(a)/(b))
michael@0 91
michael@0 92 #ifdef FLOAT_APPROX
michael@0 93
michael@0 94 /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 and an offset of 127
michael@0 95 denorm, +/- inf and NaN are *not* handled */
michael@0 96
michael@0 97 /** Base-2 log approximation (log2(x)). */
michael@0 98 static OPUS_INLINE float celt_log2(float x)
michael@0 99 {
michael@0 100 int integer;
michael@0 101 float frac;
michael@0 102 union {
michael@0 103 float f;
michael@0 104 opus_uint32 i;
michael@0 105 } in;
michael@0 106 in.f = x;
michael@0 107 integer = (in.i>>23)-127;
michael@0 108 in.i -= integer<<23;
michael@0 109 frac = in.f - 1.5f;
michael@0 110 frac = -0.41445418f + frac*(0.95909232f
michael@0 111 + frac*(-0.33951290f + frac*0.16541097f));
michael@0 112 return 1+integer+frac;
michael@0 113 }
michael@0 114
michael@0 115 /** Base-2 exponential approximation (2^x). */
michael@0 116 static OPUS_INLINE float celt_exp2(float x)
michael@0 117 {
michael@0 118 int integer;
michael@0 119 float frac;
michael@0 120 union {
michael@0 121 float f;
michael@0 122 opus_uint32 i;
michael@0 123 } res;
michael@0 124 integer = floor(x);
michael@0 125 if (integer < -50)
michael@0 126 return 0;
michael@0 127 frac = x-integer;
michael@0 128 /* K0 = 1, K1 = log(2), K2 = 3-4*log(2), K3 = 3*log(2) - 2 */
michael@0 129 res.f = 0.99992522f + frac * (0.69583354f
michael@0 130 + frac * (0.22606716f + 0.078024523f*frac));
michael@0 131 res.i = (res.i + (integer<<23)) & 0x7fffffff;
michael@0 132 return res.f;
michael@0 133 }
michael@0 134
michael@0 135 #else
michael@0 136 #define celt_log2(x) ((float)(1.442695040888963387*log(x)))
michael@0 137 #define celt_exp2(x) ((float)exp(0.6931471805599453094*(x)))
michael@0 138 #endif
michael@0 139
michael@0 140 #endif
michael@0 141
michael@0 142 #ifdef FIXED_POINT
michael@0 143
michael@0 144 #include "os_support.h"
michael@0 145
michael@0 146 #ifndef OVERRIDE_CELT_ILOG2
michael@0 147 /** Integer log in base2. Undefined for zero and negative numbers */
michael@0 148 static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x)
michael@0 149 {
michael@0 150 celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers");
michael@0 151 return EC_ILOG(x)-1;
michael@0 152 }
michael@0 153 #endif
michael@0 154
michael@0 155
michael@0 156 /** Integer log in base2. Defined for zero, but not for negative numbers */
michael@0 157 static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x)
michael@0 158 {
michael@0 159 return x <= 0 ? 0 : celt_ilog2(x);
michael@0 160 }
michael@0 161
michael@0 162 opus_val16 celt_rsqrt_norm(opus_val32 x);
michael@0 163
michael@0 164 opus_val32 celt_sqrt(opus_val32 x);
michael@0 165
michael@0 166 opus_val16 celt_cos_norm(opus_val32 x);
michael@0 167
michael@0 168 /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
michael@0 169 static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
michael@0 170 {
michael@0 171 int i;
michael@0 172 opus_val16 n, frac;
michael@0 173 /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605,
michael@0 174 0.15530808010959576, -0.08556153059057618 */
michael@0 175 static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545, -1401};
michael@0 176 if (x==0)
michael@0 177 return -32767;
michael@0 178 i = celt_ilog2(x);
michael@0 179 n = VSHR32(x,i-15)-32768-16384;
michael@0 180 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]))))))));
michael@0 181 return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT);
michael@0 182 }
michael@0 183
michael@0 184 /*
michael@0 185 K0 = 1
michael@0 186 K1 = log(2)
michael@0 187 K2 = 3-4*log(2)
michael@0 188 K3 = 3*log(2) - 2
michael@0 189 */
michael@0 190 #define D0 16383
michael@0 191 #define D1 22804
michael@0 192 #define D2 14819
michael@0 193 #define D3 10204
michael@0 194
michael@0 195 static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x)
michael@0 196 {
michael@0 197 opus_val16 frac;
michael@0 198 frac = SHL16(x, 4);
michael@0 199 return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 , MULT16_16_Q15(D3,frac))))));
michael@0 200 }
michael@0 201 /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */
michael@0 202 static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
michael@0 203 {
michael@0 204 int integer;
michael@0 205 opus_val16 frac;
michael@0 206 integer = SHR16(x,10);
michael@0 207 if (integer>14)
michael@0 208 return 0x7f000000;
michael@0 209 else if (integer < -15)
michael@0 210 return 0;
michael@0 211 frac = celt_exp2_frac(x-SHL16(integer,10));
michael@0 212 return VSHR32(EXTEND32(frac), -integer-2);
michael@0 213 }
michael@0 214
michael@0 215 opus_val32 celt_rcp(opus_val32 x);
michael@0 216
michael@0 217 #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b))
michael@0 218
michael@0 219 opus_val32 frac_div32(opus_val32 a, opus_val32 b);
michael@0 220
michael@0 221 #define M1 32767
michael@0 222 #define M2 -21
michael@0 223 #define M3 -11943
michael@0 224 #define M4 4936
michael@0 225
michael@0 226 /* Atan approximation using a 4th order polynomial. Input is in Q15 format
michael@0 227 and normalized by pi/4. Output is in Q15 format */
michael@0 228 static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x)
michael@0 229 {
michael@0 230 return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x)))))));
michael@0 231 }
michael@0 232
michael@0 233 #undef M1
michael@0 234 #undef M2
michael@0 235 #undef M3
michael@0 236 #undef M4
michael@0 237
michael@0 238 /* atan2() approximation valid for positive input values */
michael@0 239 static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x)
michael@0 240 {
michael@0 241 if (y < x)
michael@0 242 {
michael@0 243 opus_val32 arg;
michael@0 244 arg = celt_div(SHL32(EXTEND32(y),15),x);
michael@0 245 if (arg >= 32767)
michael@0 246 arg = 32767;
michael@0 247 return SHR16(celt_atan01(EXTRACT16(arg)),1);
michael@0 248 } else {
michael@0 249 opus_val32 arg;
michael@0 250 arg = celt_div(SHL32(EXTEND32(x),15),y);
michael@0 251 if (arg >= 32767)
michael@0 252 arg = 32767;
michael@0 253 return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1);
michael@0 254 }
michael@0 255 }
michael@0 256
michael@0 257 #endif /* FIXED_POINT */
michael@0 258 #endif /* MATHOPS_H */

mercurial