Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 1 | /* Copyright (C) 2002 Jean-Marc Valin |
| 2 | File: math_approx.c |
| 3 | Various math approximation functions for Speex |
| 4 | |
| 5 | Redistribution and use in source and binary forms, with or without |
| 6 | modification, are permitted provided that the following conditions |
| 7 | are met: |
| 8 | |
| 9 | - Redistributions of source code must retain the above copyright |
| 10 | notice, this list of conditions and the following disclaimer. |
| 11 | |
| 12 | - Redistributions in binary form must reproduce the above copyright |
| 13 | notice, this list of conditions and the following disclaimer in the |
| 14 | documentation and/or other materials provided with the distribution. |
| 15 | |
| 16 | - Neither the name of the Xiph.org Foundation nor the names of its |
| 17 | contributors may be used to endorse or promote products derived from |
| 18 | this software without specific prior written permission. |
| 19 | |
| 20 | THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| 21 | ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| 22 | LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| 23 | A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
| 24 | CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| 25 | EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| 26 | PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
| 27 | PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
| 28 | LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
| 29 | NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| 30 | SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 31 | */ |
| 32 | |
| 33 | #ifdef HAVE_CONFIG_H |
| 34 | #include "config.h" |
| 35 | #endif |
| 36 | |
| 37 | #include "math_approx.h" |
| 38 | #include "misc.h" |
| 39 | |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 40 | spx_int16_t spx_ilog2(spx_uint32_t x) |
| 41 | { |
| 42 | int r=0; |
| 43 | if (x>=(spx_int32_t)65536) |
| 44 | { |
| 45 | x >>= 16; |
| 46 | r += 16; |
| 47 | } |
| 48 | if (x>=256) |
| 49 | { |
| 50 | x >>= 8; |
| 51 | r += 8; |
| 52 | } |
| 53 | if (x>=16) |
| 54 | { |
| 55 | x >>= 4; |
| 56 | r += 4; |
| 57 | } |
| 58 | if (x>=4) |
| 59 | { |
| 60 | x >>= 2; |
| 61 | r += 2; |
| 62 | } |
| 63 | if (x>=2) |
| 64 | { |
| 65 | r += 1; |
| 66 | } |
| 67 | return r; |
| 68 | } |
| 69 | |
| 70 | spx_int16_t spx_ilog4(spx_uint32_t x) |
| 71 | { |
| 72 | int r=0; |
| 73 | if (x>=(spx_int32_t)65536) |
| 74 | { |
| 75 | x >>= 16; |
| 76 | r += 8; |
| 77 | } |
| 78 | if (x>=256) |
| 79 | { |
| 80 | x >>= 8; |
| 81 | r += 4; |
| 82 | } |
| 83 | if (x>=16) |
| 84 | { |
| 85 | x >>= 4; |
| 86 | r += 2; |
| 87 | } |
| 88 | if (x>=4) |
| 89 | { |
| 90 | r += 1; |
| 91 | } |
| 92 | return r; |
| 93 | } |
| 94 | |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 95 | #ifdef FIXED_POINT |
| 96 | |
| 97 | /* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */ |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 98 | /*#define C0 3634 |
| 99 | #define C1 21173 |
| 100 | #define C2 -12627 |
| 101 | #define C3 4215*/ |
| 102 | |
| 103 | /* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */ |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 104 | #define C0 3634 |
| 105 | #define C1 21173 |
| 106 | #define C2 -12627 |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 107 | #define C3 4204 |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 108 | |
| 109 | spx_word16_t spx_sqrt(spx_word32_t x) |
| 110 | { |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 111 | int k; |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 112 | spx_word32_t rt; |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 113 | k = spx_ilog4(x)-6; |
| 114 | x = VSHR32(x, (k<<1)); |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 115 | rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3))))))); |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 116 | rt = VSHR32(rt,7-k); |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 117 | return rt; |
| 118 | } |
| 119 | |
| 120 | /* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */ |
| 121 | |
| 122 | |
| 123 | #define A1 16469 |
| 124 | #define A2 2242 |
| 125 | #define A3 1486 |
| 126 | |
| 127 | spx_word16_t spx_acos(spx_word16_t x) |
| 128 | { |
| 129 | int s=0; |
| 130 | spx_word16_t ret; |
| 131 | spx_word16_t sq; |
| 132 | if (x<0) |
| 133 | { |
| 134 | s=1; |
| 135 | x = NEG16(x); |
| 136 | } |
| 137 | x = SUB16(16384,x); |
| 138 | |
| 139 | x = x >> 1; |
| 140 | sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3)))))); |
| 141 | ret = spx_sqrt(SHL32(EXTEND32(sq),13)); |
| 142 | |
| 143 | /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/ |
| 144 | if (s) |
| 145 | ret = SUB16(25736,ret); |
| 146 | return ret; |
| 147 | } |
| 148 | |
| 149 | |
| 150 | #define K1 8192 |
| 151 | #define K2 -4096 |
| 152 | #define K3 340 |
| 153 | #define K4 -10 |
| 154 | |
| 155 | spx_word16_t spx_cos(spx_word16_t x) |
| 156 | { |
| 157 | spx_word16_t x2; |
| 158 | |
| 159 | if (x<12868) |
| 160 | { |
| 161 | x2 = MULT16_16_P13(x,x); |
| 162 | return ADD32(K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); |
| 163 | } else { |
| 164 | x = SUB16(25736,x); |
| 165 | x2 = MULT16_16_P13(x,x); |
| 166 | return SUB32(-K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); |
| 167 | } |
| 168 | } |
| 169 | |
Benny Prijono | d2990b9 | 2006-11-23 10:19:46 +0000 | [diff] [blame] | 170 | #define L1 32767 |
| 171 | #define L2 -7651 |
| 172 | #define L3 8277 |
| 173 | #define L4 -626 |
| 174 | |
| 175 | static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x) |
| 176 | { |
| 177 | spx_word16_t x2; |
| 178 | |
| 179 | x2 = MULT16_16_P15(x,x); |
| 180 | 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)))))))); |
| 181 | } |
| 182 | |
| 183 | spx_word16_t spx_cos_norm(spx_word32_t x) |
| 184 | { |
| 185 | x = x&0x0001ffff; |
| 186 | if (x>SHL32(EXTEND32(1), 16)) |
| 187 | x = SUB32(SHL32(EXTEND32(1), 17),x); |
| 188 | if (x&0x00007fff) |
| 189 | { |
| 190 | if (x<SHL32(EXTEND32(1), 15)) |
| 191 | { |
| 192 | return _spx_cos_pi_2(EXTRACT16(x)); |
| 193 | } else { |
| 194 | return NEG32(_spx_cos_pi_2(EXTRACT16(65536-x))); |
| 195 | } |
| 196 | } else { |
| 197 | if (x&0x0000ffff) |
| 198 | return 0; |
| 199 | else if (x&0x0001ffff) |
| 200 | return -32767; |
| 201 | else |
| 202 | return 32767; |
| 203 | } |
| 204 | } |
| 205 | |
| 206 | /* |
| 207 | K0 = 1 |
| 208 | K1 = log(2) |
| 209 | K2 = 3-4*log(2) |
| 210 | K3 = 3*log(2) - 2 |
| 211 | */ |
| 212 | #define D0 16384 |
| 213 | #define D1 11356 |
| 214 | #define D2 3726 |
| 215 | #define D3 1301 |
| 216 | /* Input in Q11 format, output in Q16 */ |
| 217 | static spx_word32_t spx_exp2(spx_word16_t x) |
| 218 | { |
| 219 | int integer; |
| 220 | spx_word16_t frac; |
| 221 | integer = SHR16(x,11); |
| 222 | if (integer>14) |
| 223 | return 0x7fffffff; |
| 224 | else if (integer < -15) |
| 225 | return 0; |
| 226 | frac = SHL16(x-SHL16(integer,11),3); |
| 227 | frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac)))))); |
| 228 | return VSHR32(EXTEND32(frac), -integer-2); |
| 229 | } |
| 230 | |
| 231 | /* Input in Q11 format, output in Q16 */ |
| 232 | spx_word32_t spx_exp(spx_word16_t x) |
| 233 | { |
| 234 | if (x>21290) |
| 235 | return 0x7fffffff; |
| 236 | else if (x<-21290) |
| 237 | return 0; |
| 238 | else |
| 239 | return spx_exp2(MULT16_16_P14(23637,x)); |
| 240 | } |
| 241 | #define M1 32767 |
| 242 | #define M2 -21 |
| 243 | #define M3 -11943 |
| 244 | #define M4 4936 |
| 245 | |
| 246 | static inline spx_word16_t spx_atan01(spx_word16_t x) |
| 247 | { |
| 248 | return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x))))))); |
| 249 | } |
| 250 | |
| 251 | /* Input in Q15, output in Q14 */ |
| 252 | spx_word16_t spx_atan(spx_word32_t x) |
| 253 | { |
| 254 | if (x <= 32767) |
| 255 | { |
| 256 | return SHR16(spx_atan01(x),1); |
| 257 | } else { |
| 258 | int e = spx_ilog2(x); |
| 259 | if (e>=29) |
| 260 | return 25736; |
| 261 | x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14))); |
| 262 | return SUB16(25736, SHR16(spx_atan01(x),1)); |
| 263 | } |
| 264 | } |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 265 | #else |
| 266 | |
| 267 | #ifndef M_PI |
| 268 | #define M_PI 3.14159265358979323846 /* pi */ |
| 269 | #endif |
| 270 | |
| 271 | #define C1 0.9999932946f |
| 272 | #define C2 -0.4999124376f |
| 273 | #define C3 0.0414877472f |
| 274 | #define C4 -0.0012712095f |
| 275 | |
| 276 | |
| 277 | #define SPX_PI_2 1.5707963268 |
| 278 | spx_word16_t spx_cos(spx_word16_t x) |
| 279 | { |
| 280 | if (x<SPX_PI_2) |
| 281 | { |
| 282 | x *= x; |
| 283 | return C1 + x*(C2+x*(C3+C4*x)); |
| 284 | } else { |
| 285 | x = M_PI-x; |
| 286 | x *= x; |
| 287 | return NEG16(C1 + x*(C2+x*(C3+C4*x))); |
| 288 | } |
| 289 | } |
| 290 | |
Benny Prijono | eb30bf5 | 2006-03-04 20:43:52 +0000 | [diff] [blame] | 291 | #endif |