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 |