| OLD | NEW |
| 1 /* Copyright (c) 2002-2008 Jean-Marc Valin | 1 /* Copyright (c) 2002-2008 Jean-Marc Valin |
| 2 Copyright (c) 2007-2008 CSIRO | 2 Copyright (c) 2007-2008 CSIRO |
| 3 Copyright (c) 2007-2009 Xiph.Org Foundation | 3 Copyright (c) 2007-2009 Xiph.Org Foundation |
| 4 Written by Jean-Marc Valin */ | 4 Written by Jean-Marc Valin */ |
| 5 /** | 5 /** |
| 6 @file mathops.h | 6 @file mathops.h |
| 7 @brief Various math functions | 7 @brief Various math functions |
| 8 */ | 8 */ |
| 9 /* | 9 /* |
| 10 Redistribution and use in source and binary forms, with or without | 10 Redistribution and use in source and binary forms, with or without |
| (...skipping 26 matching lines...) Expand all Loading... |
| 37 #include "arch.h" | 37 #include "arch.h" |
| 38 #include "entcode.h" | 38 #include "entcode.h" |
| 39 #include "os_support.h" | 39 #include "os_support.h" |
| 40 | 40 |
| 41 /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is impor
tant */ | 41 /* Multiplies two 16-bit fractional values. Bit-exactness of this macro is impor
tant */ |
| 42 #define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>
15) | 42 #define FRAC_MUL16(a,b) ((16384+((opus_int32)(opus_int16)(a)*(opus_int16)(b)))>>
15) |
| 43 | 43 |
| 44 unsigned isqrt32(opus_uint32 _val); | 44 unsigned isqrt32(opus_uint32 _val); |
| 45 | 45 |
| 46 #ifndef OVERRIDE_CELT_MAXABS16 | 46 #ifndef OVERRIDE_CELT_MAXABS16 |
| 47 static inline opus_val32 celt_maxabs16(const opus_val16 *x, int len) | 47 static OPUS_INLINE opus_val32 celt_maxabs16(const opus_val16 *x, int len) |
| 48 { | 48 { |
| 49 int i; | 49 int i; |
| 50 opus_val16 maxval = 0; | 50 opus_val16 maxval = 0; |
| 51 opus_val16 minval = 0; | 51 opus_val16 minval = 0; |
| 52 for (i=0;i<len;i++) | 52 for (i=0;i<len;i++) |
| 53 { | 53 { |
| 54 maxval = MAX16(maxval, x[i]); | 54 maxval = MAX16(maxval, x[i]); |
| 55 minval = MIN16(minval, x[i]); | 55 minval = MIN16(minval, x[i]); |
| 56 } | 56 } |
| 57 return MAX32(EXTEND32(maxval),-EXTEND32(minval)); | 57 return MAX32(EXTEND32(maxval),-EXTEND32(minval)); |
| 58 } | 58 } |
| 59 #endif | 59 #endif |
| 60 | 60 |
| 61 #ifndef OVERRIDE_CELT_MAXABS32 | 61 #ifndef OVERRIDE_CELT_MAXABS32 |
| 62 #ifdef FIXED_POINT | 62 #ifdef FIXED_POINT |
| 63 static inline opus_val32 celt_maxabs32(const opus_val32 *x, int len) | 63 static OPUS_INLINE opus_val32 celt_maxabs32(const opus_val32 *x, int len) |
| 64 { | 64 { |
| 65 int i; | 65 int i; |
| 66 opus_val32 maxval = 0; | 66 opus_val32 maxval = 0; |
| 67 opus_val32 minval = 0; | 67 opus_val32 minval = 0; |
| 68 for (i=0;i<len;i++) | 68 for (i=0;i<len;i++) |
| 69 { | 69 { |
| 70 maxval = MAX32(maxval, x[i]); | 70 maxval = MAX32(maxval, x[i]); |
| 71 minval = MIN32(minval, x[i]); | 71 minval = MIN32(minval, x[i]); |
| 72 } | 72 } |
| 73 return MAX32(maxval, -minval); | 73 return MAX32(maxval, -minval); |
| (...skipping 14 matching lines...) Expand all Loading... |
| 88 #define celt_rcp(x) (1.f/(x)) | 88 #define celt_rcp(x) (1.f/(x)) |
| 89 #define celt_div(a,b) ((a)/(b)) | 89 #define celt_div(a,b) ((a)/(b)) |
| 90 #define frac_div32(a,b) ((float)(a)/(b)) | 90 #define frac_div32(a,b) ((float)(a)/(b)) |
| 91 | 91 |
| 92 #ifdef FLOAT_APPROX | 92 #ifdef FLOAT_APPROX |
| 93 | 93 |
| 94 /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 an
d an offset of 127 | 94 /* Note: This assumes radix-2 floating point with the exponent at bits 23..30 an
d an offset of 127 |
| 95 denorm, +/- inf and NaN are *not* handled */ | 95 denorm, +/- inf and NaN are *not* handled */ |
| 96 | 96 |
| 97 /** Base-2 log approximation (log2(x)). */ | 97 /** Base-2 log approximation (log2(x)). */ |
| 98 static inline float celt_log2(float x) | 98 static OPUS_INLINE float celt_log2(float x) |
| 99 { | 99 { |
| 100 int integer; | 100 int integer; |
| 101 float frac; | 101 float frac; |
| 102 union { | 102 union { |
| 103 float f; | 103 float f; |
| 104 opus_uint32 i; | 104 opus_uint32 i; |
| 105 } in; | 105 } in; |
| 106 in.f = x; | 106 in.f = x; |
| 107 integer = (in.i>>23)-127; | 107 integer = (in.i>>23)-127; |
| 108 in.i -= integer<<23; | 108 in.i -= integer<<23; |
| 109 frac = in.f - 1.5f; | 109 frac = in.f - 1.5f; |
| 110 frac = -0.41445418f + frac*(0.95909232f | 110 frac = -0.41445418f + frac*(0.95909232f |
| 111 + frac*(-0.33951290f + frac*0.16541097f)); | 111 + frac*(-0.33951290f + frac*0.16541097f)); |
| 112 return 1+integer+frac; | 112 return 1+integer+frac; |
| 113 } | 113 } |
| 114 | 114 |
| 115 /** Base-2 exponential approximation (2^x). */ | 115 /** Base-2 exponential approximation (2^x). */ |
| 116 static inline float celt_exp2(float x) | 116 static OPUS_INLINE float celt_exp2(float x) |
| 117 { | 117 { |
| 118 int integer; | 118 int integer; |
| 119 float frac; | 119 float frac; |
| 120 union { | 120 union { |
| 121 float f; | 121 float f; |
| 122 opus_uint32 i; | 122 opus_uint32 i; |
| 123 } res; | 123 } res; |
| 124 integer = floor(x); | 124 integer = floor(x); |
| 125 if (integer < -50) | 125 if (integer < -50) |
| 126 return 0; | 126 return 0; |
| (...skipping 11 matching lines...) Expand all Loading... |
| 138 #endif | 138 #endif |
| 139 | 139 |
| 140 #endif | 140 #endif |
| 141 | 141 |
| 142 #ifdef FIXED_POINT | 142 #ifdef FIXED_POINT |
| 143 | 143 |
| 144 #include "os_support.h" | 144 #include "os_support.h" |
| 145 | 145 |
| 146 #ifndef OVERRIDE_CELT_ILOG2 | 146 #ifndef OVERRIDE_CELT_ILOG2 |
| 147 /** Integer log in base2. Undefined for zero and negative numbers */ | 147 /** Integer log in base2. Undefined for zero and negative numbers */ |
| 148 static inline opus_int16 celt_ilog2(opus_int32 x) | 148 static OPUS_INLINE opus_int16 celt_ilog2(opus_int32 x) |
| 149 { | 149 { |
| 150 celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers"); | 150 celt_assert2(x>0, "celt_ilog2() only defined for strictly positive numbers"); |
| 151 return EC_ILOG(x)-1; | 151 return EC_ILOG(x)-1; |
| 152 } | 152 } |
| 153 #endif | 153 #endif |
| 154 | 154 |
| 155 | 155 |
| 156 /** Integer log in base2. Defined for zero, but not for negative numbers */ | 156 /** Integer log in base2. Defined for zero, but not for negative numbers */ |
| 157 static inline opus_int16 celt_zlog2(opus_val32 x) | 157 static OPUS_INLINE opus_int16 celt_zlog2(opus_val32 x) |
| 158 { | 158 { |
| 159 return x <= 0 ? 0 : celt_ilog2(x); | 159 return x <= 0 ? 0 : celt_ilog2(x); |
| 160 } | 160 } |
| 161 | 161 |
| 162 opus_val16 celt_rsqrt_norm(opus_val32 x); | 162 opus_val16 celt_rsqrt_norm(opus_val32 x); |
| 163 | 163 |
| 164 opus_val32 celt_sqrt(opus_val32 x); | 164 opus_val32 celt_sqrt(opus_val32 x); |
| 165 | 165 |
| 166 opus_val16 celt_cos_norm(opus_val32 x); | 166 opus_val16 celt_cos_norm(opus_val32 x); |
| 167 | 167 |
| 168 /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */ | 168 /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */ |
| 169 static inline opus_val16 celt_log2(opus_val32 x) | 169 static OPUS_INLINE opus_val16 celt_log2(opus_val32 x) |
| 170 { | 170 { |
| 171 int i; | 171 int i; |
| 172 opus_val16 n, frac; | 172 opus_val16 n, frac; |
| 173 /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605, | 173 /* -0.41509302963303146, 0.9609890551383969, -0.31836011537636605, |
| 174 0.15530808010959576, -0.08556153059057618 */ | 174 0.15530808010959576, -0.08556153059057618 */ |
| 175 static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545,
-1401}; | 175 static const opus_val16 C[5] = {-6801+(1<<(13-DB_SHIFT)), 15746, -5217, 2545,
-1401}; |
| 176 if (x==0) | 176 if (x==0) |
| 177 return -32767; | 177 return -32767; |
| 178 i = celt_ilog2(x); | 178 i = celt_ilog2(x); |
| 179 n = VSHR32(x,i-15)-32768-16384; | 179 n = VSHR32(x,i-15)-32768-16384; |
| 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])))))))); | 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])))))))); |
| 181 return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT); | 181 return SHL16(i-13,DB_SHIFT)+SHR16(frac,14-DB_SHIFT); |
| 182 } | 182 } |
| 183 | 183 |
| 184 /* | 184 /* |
| 185 K0 = 1 | 185 K0 = 1 |
| 186 K1 = log(2) | 186 K1 = log(2) |
| 187 K2 = 3-4*log(2) | 187 K2 = 3-4*log(2) |
| 188 K3 = 3*log(2) - 2 | 188 K3 = 3*log(2) - 2 |
| 189 */ | 189 */ |
| 190 #define D0 16383 | 190 #define D0 16383 |
| 191 #define D1 22804 | 191 #define D1 22804 |
| 192 #define D2 14819 | 192 #define D2 14819 |
| 193 #define D3 10204 | 193 #define D3 10204 |
| 194 | 194 |
| 195 static inline opus_val32 celt_exp2_frac(opus_val16 x) | 195 static OPUS_INLINE opus_val32 celt_exp2_frac(opus_val16 x) |
| 196 { | 196 { |
| 197 opus_val16 frac; | 197 opus_val16 frac; |
| 198 frac = SHL16(x, 4); | 198 frac = SHL16(x, 4); |
| 199 return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 ,
MULT16_16_Q15(D3,frac)))))); | 199 return ADD16(D0, MULT16_16_Q15(frac, ADD16(D1, MULT16_16_Q15(frac, ADD16(D2 ,
MULT16_16_Q15(D3,frac)))))); |
| 200 } | 200 } |
| 201 /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */ | 201 /** Base-2 exponential approximation (2^x). (Q10 input, Q16 output) */ |
| 202 static inline opus_val32 celt_exp2(opus_val16 x) | 202 static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x) |
| 203 { | 203 { |
| 204 int integer; | 204 int integer; |
| 205 opus_val16 frac; | 205 opus_val16 frac; |
| 206 integer = SHR16(x,10); | 206 integer = SHR16(x,10); |
| 207 if (integer>14) | 207 if (integer>14) |
| 208 return 0x7f000000; | 208 return 0x7f000000; |
| 209 else if (integer < -15) | 209 else if (integer < -15) |
| 210 return 0; | 210 return 0; |
| 211 frac = celt_exp2_frac(x-SHL16(integer,10)); | 211 frac = celt_exp2_frac(x-SHL16(integer,10)); |
| 212 return VSHR32(EXTEND32(frac), -integer-2); | 212 return VSHR32(EXTEND32(frac), -integer-2); |
| 213 } | 213 } |
| 214 | 214 |
| 215 opus_val32 celt_rcp(opus_val32 x); | 215 opus_val32 celt_rcp(opus_val32 x); |
| 216 | 216 |
| 217 #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b)) | 217 #define celt_div(a,b) MULT32_32_Q31((opus_val32)(a),celt_rcp(b)) |
| 218 | 218 |
| 219 opus_val32 frac_div32(opus_val32 a, opus_val32 b); | 219 opus_val32 frac_div32(opus_val32 a, opus_val32 b); |
| 220 | 220 |
| 221 #define M1 32767 | 221 #define M1 32767 |
| 222 #define M2 -21 | 222 #define M2 -21 |
| 223 #define M3 -11943 | 223 #define M3 -11943 |
| 224 #define M4 4936 | 224 #define M4 4936 |
| 225 | 225 |
| 226 /* Atan approximation using a 4th order polynomial. Input is in Q15 format | 226 /* Atan approximation using a 4th order polynomial. Input is in Q15 format |
| 227 and normalized by pi/4. Output is in Q15 format */ | 227 and normalized by pi/4. Output is in Q15 format */ |
| 228 static inline opus_val16 celt_atan01(opus_val16 x) | 228 static OPUS_INLINE opus_val16 celt_atan01(opus_val16 x) |
| 229 { | 229 { |
| 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))))))); | 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))))))); |
| 231 } | 231 } |
| 232 | 232 |
| 233 #undef M1 | 233 #undef M1 |
| 234 #undef M2 | 234 #undef M2 |
| 235 #undef M3 | 235 #undef M3 |
| 236 #undef M4 | 236 #undef M4 |
| 237 | 237 |
| 238 /* atan2() approximation valid for positive input values */ | 238 /* atan2() approximation valid for positive input values */ |
| 239 static inline opus_val16 celt_atan2p(opus_val16 y, opus_val16 x) | 239 static OPUS_INLINE opus_val16 celt_atan2p(opus_val16 y, opus_val16 x) |
| 240 { | 240 { |
| 241 if (y < x) | 241 if (y < x) |
| 242 { | 242 { |
| 243 opus_val32 arg; | 243 opus_val32 arg; |
| 244 arg = celt_div(SHL32(EXTEND32(y),15),x); | 244 arg = celt_div(SHL32(EXTEND32(y),15),x); |
| 245 if (arg >= 32767) | 245 if (arg >= 32767) |
| 246 arg = 32767; | 246 arg = 32767; |
| 247 return SHR16(celt_atan01(EXTRACT16(arg)),1); | 247 return SHR16(celt_atan01(EXTRACT16(arg)),1); |
| 248 } else { | 248 } else { |
| 249 opus_val32 arg; | 249 opus_val32 arg; |
| 250 arg = celt_div(SHL32(EXTEND32(x),15),y); | 250 arg = celt_div(SHL32(EXTEND32(x),15),y); |
| 251 if (arg >= 32767) | 251 if (arg >= 32767) |
| 252 arg = 32767; | 252 arg = 32767; |
| 253 return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1); | 253 return 25736-SHR16(celt_atan01(EXTRACT16(arg)),1); |
| 254 } | 254 } |
| 255 } | 255 } |
| 256 | 256 |
| 257 #endif /* FIXED_POINT */ | 257 #endif /* FIXED_POINT */ |
| 258 #endif /* MATHOPS_H */ | 258 #endif /* MATHOPS_H */ |
| OLD | NEW |