OLD | NEW |
(Empty) | |
| 1 /* Copyright (C) 2005 Jean-Marc Valin */ |
| 2 /** |
| 3 @file pseudofloat.h |
| 4 @brief Pseudo-floating point |
| 5 * This header file provides a lightweight floating point type for |
| 6 * use on fixed-point platforms when a large dynamic range is |
| 7 * required. The new type is not compatible with the 32-bit IEEE format, |
| 8 * it is not even remotely as accurate as 32-bit floats, and is not |
| 9 * even guaranteed to produce even remotely correct results for code |
| 10 * other than Speex. It makes all kinds of shortcuts that are acceptable |
| 11 * for Speex, but may not be acceptable for your application. You're |
| 12 * quite welcome to reuse this code and improve it, but don't assume |
| 13 * it works out of the box. Most likely, it doesn't. |
| 14 */ |
| 15 /* |
| 16 Redistribution and use in source and binary forms, with or without |
| 17 modification, are permitted provided that the following conditions |
| 18 are met: |
| 19 |
| 20 - Redistributions of source code must retain the above copyright |
| 21 notice, this list of conditions and the following disclaimer. |
| 22 |
| 23 - Redistributions in binary form must reproduce the above copyright |
| 24 notice, this list of conditions and the following disclaimer in the |
| 25 documentation and/or other materials provided with the distribution. |
| 26 |
| 27 - Neither the name of the Xiph.org Foundation nor the names of its |
| 28 contributors may be used to endorse or promote products derived from |
| 29 this software without specific prior written permission. |
| 30 |
| 31 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| 32 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| 33 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| 34 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
| 35 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| 36 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| 37 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
| 38 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
| 39 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
| 40 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| 41 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| 42 */ |
| 43 |
| 44 #ifndef PSEUDOFLOAT_H |
| 45 #define PSEUDOFLOAT_H |
| 46 |
| 47 #include "arch.h" |
| 48 #include "os_support.h" |
| 49 #include "math_approx.h" |
| 50 #include <math.h> |
| 51 |
| 52 #ifdef FIXED_POINT |
| 53 |
| 54 typedef struct { |
| 55 spx_int16_t m; |
| 56 spx_int16_t e; |
| 57 } spx_float_t; |
| 58 |
| 59 static const spx_float_t FLOAT_ZERO = {0,0}; |
| 60 static const spx_float_t FLOAT_ONE = {16384,-14}; |
| 61 static const spx_float_t FLOAT_HALF = {16384,-15}; |
| 62 |
| 63 #define MIN(a,b) ((a)<(b)?(a):(b)) |
| 64 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x) |
| 65 { |
| 66 int e=0; |
| 67 int sign=0; |
| 68 if (x<0) |
| 69 { |
| 70 sign = 1; |
| 71 x = -x; |
| 72 } |
| 73 if (x==0) |
| 74 { |
| 75 spx_float_t r = {0,0}; |
| 76 return r; |
| 77 } |
| 78 e = spx_ilog2(ABS32(x))-14; |
| 79 x = VSHR32(x, e); |
| 80 if (sign) |
| 81 { |
| 82 spx_float_t r; |
| 83 r.m = -x; |
| 84 r.e = e; |
| 85 return r; |
| 86 } |
| 87 else |
| 88 { |
| 89 spx_float_t r; |
| 90 r.m = x; |
| 91 r.e = e; |
| 92 return r; |
| 93 } |
| 94 } |
| 95 |
| 96 |
| 97 static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b) |
| 98 { |
| 99 spx_float_t r; |
| 100 if (a.m==0) |
| 101 return b; |
| 102 else if (b.m==0) |
| 103 return a; |
| 104 if ((a).e > (b).e) |
| 105 { |
| 106 r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1)); |
| 107 r.e = (a).e+1; |
| 108 } |
| 109 else |
| 110 { |
| 111 r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1)); |
| 112 r.e = (b).e+1; |
| 113 } |
| 114 if (r.m>0) |
| 115 { |
| 116 if (r.m<16384) |
| 117 { |
| 118 r.m<<=1; |
| 119 r.e-=1; |
| 120 } |
| 121 } else { |
| 122 if (r.m>-16384) |
| 123 { |
| 124 r.m<<=1; |
| 125 r.e-=1; |
| 126 } |
| 127 } |
| 128 /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
| 129 return r; |
| 130 } |
| 131 |
| 132 static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b) |
| 133 { |
| 134 spx_float_t r; |
| 135 if (a.m==0) |
| 136 return b; |
| 137 else if (b.m==0) |
| 138 return a; |
| 139 if ((a).e > (b).e) |
| 140 { |
| 141 r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1)); |
| 142 r.e = (a).e+1; |
| 143 } |
| 144 else |
| 145 { |
| 146 r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1); |
| 147 r.e = (b).e+1; |
| 148 } |
| 149 if (r.m>0) |
| 150 { |
| 151 if (r.m<16384) |
| 152 { |
| 153 r.m<<=1; |
| 154 r.e-=1; |
| 155 } |
| 156 } else { |
| 157 if (r.m>-16384) |
| 158 { |
| 159 r.m<<=1; |
| 160 r.e-=1; |
| 161 } |
| 162 } |
| 163 /*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
| 164 return r; |
| 165 } |
| 166 |
| 167 static inline int FLOAT_LT(spx_float_t a, spx_float_t b) |
| 168 { |
| 169 if (a.m==0) |
| 170 return b.m>0; |
| 171 else if (b.m==0) |
| 172 return a.m<0; |
| 173 if ((a).e > (b).e) |
| 174 return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1)); |
| 175 else |
| 176 return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1)); |
| 177 |
| 178 } |
| 179 |
| 180 static inline int FLOAT_GT(spx_float_t a, spx_float_t b) |
| 181 { |
| 182 return FLOAT_LT(b,a); |
| 183 } |
| 184 |
| 185 static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b) |
| 186 { |
| 187 spx_float_t r; |
| 188 r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15); |
| 189 r.e = (a).e+(b).e+15; |
| 190 if (r.m>0) |
| 191 { |
| 192 if (r.m<16384) |
| 193 { |
| 194 r.m<<=1; |
| 195 r.e-=1; |
| 196 } |
| 197 } else { |
| 198 if (r.m>-16384) |
| 199 { |
| 200 r.m<<=1; |
| 201 r.e-=1; |
| 202 } |
| 203 } |
| 204 /*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
| 205 return r; |
| 206 } |
| 207 |
| 208 static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b) |
| 209 { |
| 210 spx_float_t r; |
| 211 r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15); |
| 212 r.e = (a).e+(b).e+15; |
| 213 return r; |
| 214 } |
| 215 |
| 216 |
| 217 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b) |
| 218 { |
| 219 spx_float_t r; |
| 220 r.m = a.m; |
| 221 r.e = a.e+b; |
| 222 return r; |
| 223 } |
| 224 |
| 225 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a) |
| 226 { |
| 227 if (a.e<0) |
| 228 return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e); |
| 229 else |
| 230 return a.m<<a.e; |
| 231 } |
| 232 |
| 233 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a) |
| 234 { |
| 235 if (a.e<0) |
| 236 return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e; |
| 237 else |
| 238 return EXTEND32(a.m)<<a.e; |
| 239 } |
| 240 |
| 241 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b) |
| 242 { |
| 243 return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15); |
| 244 } |
| 245 |
| 246 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b) |
| 247 { |
| 248 int e1, e2; |
| 249 spx_float_t r; |
| 250 if (a==0 || b==0) |
| 251 { |
| 252 return FLOAT_ZERO; |
| 253 } |
| 254 e1 = spx_ilog2(ABS32(a)); |
| 255 a = VSHR32(a, e1-14); |
| 256 e2 = spx_ilog2(ABS32(b)); |
| 257 b = VSHR32(b, e2-14); |
| 258 r.m = MULT16_16_Q15(a,b); |
| 259 r.e = e1+e2-13; |
| 260 return r; |
| 261 } |
| 262 |
| 263 /* Do NOT attempt to divide by a negative number */ |
| 264 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b) |
| 265 { |
| 266 int e=0; |
| 267 spx_float_t r; |
| 268 if (a==0) |
| 269 { |
| 270 return FLOAT_ZERO; |
| 271 } |
| 272 e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15; |
| 273 a = VSHR32(a, e); |
| 274 if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15)) |
| 275 { |
| 276 a >>= 1; |
| 277 e++; |
| 278 } |
| 279 r.m = DIV32_16(a,b.m); |
| 280 r.e = e-b.e; |
| 281 return r; |
| 282 } |
| 283 |
| 284 |
| 285 /* Do NOT attempt to divide by a negative number */ |
| 286 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b) |
| 287 { |
| 288 int e0=0,e=0; |
| 289 spx_float_t r; |
| 290 if (a==0) |
| 291 { |
| 292 return FLOAT_ZERO; |
| 293 } |
| 294 if (b>32767) |
| 295 { |
| 296 e0 = spx_ilog2(b)-14; |
| 297 b = VSHR32(b, e0); |
| 298 e0 = -e0; |
| 299 } |
| 300 e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15; |
| 301 a = VSHR32(a, e); |
| 302 if (ABS32(a)>=SHL32(EXTEND32(b-1),15)) |
| 303 { |
| 304 a >>= 1; |
| 305 e++; |
| 306 } |
| 307 e += e0; |
| 308 r.m = DIV32_16(a,b); |
| 309 r.e = e; |
| 310 return r; |
| 311 } |
| 312 |
| 313 /* Do NOT attempt to divide by a negative number */ |
| 314 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) |
| 315 { |
| 316 int e=0; |
| 317 spx_int32_t num; |
| 318 spx_float_t r; |
| 319 if (b.m<=0) |
| 320 { |
| 321 speex_warning_int("Attempted to divide by", b.m); |
| 322 return FLOAT_ONE; |
| 323 } |
| 324 num = a.m; |
| 325 a.m = ABS16(a.m); |
| 326 while (a.m >= b.m) |
| 327 { |
| 328 e++; |
| 329 a.m >>= 1; |
| 330 } |
| 331 num = num << (15-e); |
| 332 r.m = DIV32_16(num,b.m); |
| 333 r.e = a.e-b.e-15+e; |
| 334 return r; |
| 335 } |
| 336 |
| 337 static inline spx_float_t FLOAT_SQRT(spx_float_t a) |
| 338 { |
| 339 spx_float_t r; |
| 340 spx_int32_t m; |
| 341 m = SHL32(EXTEND32(a.m), 14); |
| 342 r.e = a.e - 14; |
| 343 if (r.e & 1) |
| 344 { |
| 345 r.e -= 1; |
| 346 m <<= 1; |
| 347 } |
| 348 r.e >>= 1; |
| 349 r.m = spx_sqrt(m); |
| 350 return r; |
| 351 } |
| 352 |
| 353 #else |
| 354 |
| 355 #define spx_float_t float |
| 356 #define FLOAT_ZERO 0.f |
| 357 #define FLOAT_ONE 1.f |
| 358 #define FLOAT_HALF 0.5f |
| 359 #define PSEUDOFLOAT(x) (x) |
| 360 #define FLOAT_MULT(a,b) ((a)*(b)) |
| 361 #define FLOAT_AMULT(a,b) ((a)*(b)) |
| 362 #define FLOAT_MUL32(a,b) ((a)*(b)) |
| 363 #define FLOAT_DIV32(a,b) ((a)/(b)) |
| 364 #define FLOAT_EXTRACT16(a) (a) |
| 365 #define FLOAT_EXTRACT32(a) (a) |
| 366 #define FLOAT_ADD(a,b) ((a)+(b)) |
| 367 #define FLOAT_SUB(a,b) ((a)-(b)) |
| 368 #define REALFLOAT(x) (x) |
| 369 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b)) |
| 370 #define FLOAT_MUL32U(a,b) ((a)*(b)) |
| 371 #define FLOAT_SHL(a,b) (a) |
| 372 #define FLOAT_LT(a,b) ((a)<(b)) |
| 373 #define FLOAT_GT(a,b) ((a)>(b)) |
| 374 #define FLOAT_DIVU(a,b) ((a)/(b)) |
| 375 #define FLOAT_SQRT(a) (spx_sqrt(a)) |
| 376 |
| 377 #endif |
| 378 |
| 379 #endif |
OLD | NEW |