Increase learning rate for some mis-adaptated conditions (lower bound on RER).
[speexdsp.git] / libspeex / pseudofloat.h
1 /* Copyright (C) 2005 Jean-Marc Valin */
2 /**
3    @file pseudofloat.h
4    @brief Pseudo-floating point
5  */
6 /*
7    Redistribution and use in source and binary forms, with or without
8    modification, are permitted provided that the following conditions
9    are met:
10    
11    - Redistributions of source code must retain the above copyright
12    notice, this list of conditions and the following disclaimer.
13    
14    - Redistributions in binary form must reproduce the above copyright
15    notice, this list of conditions and the following disclaimer in the
16    documentation and/or other materials provided with the distribution.
17    
18    - Neither the name of the Xiph.org Foundation nor the names of its
19    contributors may be used to endorse or promote products derived from
20    this software without specific prior written permission.
21    
22    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
25    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
26    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
27    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
28    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
29    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
30    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
31    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
32    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 */
34
35 #ifndef PSEUDOFLOAT_H
36 #define PSEUDOFLOAT_H
37
38 #include "misc.h"
39 #include "math_approx.h"
40 #include <math.h>
41
42 #ifdef FIXED_POINT
43
44 typedef struct {
45    spx_int16_t m;
46    spx_int16_t e;
47 } spx_float_t;
48
49 static const spx_float_t FLOAT_ZERO = {0,0};
50 static const spx_float_t FLOAT_ONE = {16384,-14};
51 static const spx_float_t FLOAT_HALF = {16384,-15};
52
53 #define MIN(a,b) ((a)<(b)?(a):(b))
54 static inline spx_float_t PSEUDOFLOAT(spx_int32_t x)
55 {
56    int e=0;
57    int sign=0;
58    if (x<0)
59    {
60       sign = 1;
61       x = -x;
62    }
63    if (x==0)
64    {
65       spx_float_t r = {0,0};
66       return r;
67    }
68    while (x>32767)
69    {
70       x >>= 1;
71       /*x *= .5;*/
72       e++;
73    }
74    while (x<16383)
75    {
76       x <<= 1;
77       /*x *= 2;*/
78       e--;
79    }
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
209 static inline spx_float_t FLOAT_SHL(spx_float_t a, int b)
210 {
211    spx_float_t r;
212    r.m = a.m;
213    r.e = a.e+b;
214    return r;
215 }
216
217 static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a)
218 {
219    if (a.e<0)
220       return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e);
221    else
222       return a.m<<a.e;
223 }
224
225 static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a)
226 {
227    if (a.e<0)
228       return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e;
229    else
230       return EXTEND32(a.m)<<a.e;
231 }
232
233 static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b)
234 {
235    if (a.e<-15)
236       return SHR32(MULT16_32_Q15(a.m, b),-a.e-15);
237    else
238       return SHL32(MULT16_32_Q15(a.m, b),15+a.e);
239 }
240
241 static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b)
242 {
243    int e=0;
244    spx_float_t r;
245    /* FIXME: Handle the sign */
246    if (a==0)
247    {
248       return FLOAT_ZERO;
249    }
250    while (a>32767)
251    {
252       a >>= 1;
253       e++;
254    }
255    while (a<16384)
256    {
257       a <<= 1;
258       e--;
259    }
260    while (b>32767)
261    {
262       b >>= 1;
263       e++;
264    }
265    while (b<16384)
266    {
267       b <<= 1;
268       e--;
269    }
270    r.m = MULT16_16_Q15(a,b);
271    r.e = e+15;
272    return r;
273 }
274
275 static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b)
276 {
277    int e=0;
278    spx_float_t r;
279    /* FIXME: Handle the sign */
280    if (a==0)
281    {
282       return FLOAT_ZERO;
283    }
284    while (a<SHL32(EXTEND32(b.m),14))
285    {
286       a <<= 1;
287       e--;
288    }
289    while (a>=SHL32(EXTEND32(b.m-1),15))
290    {
291       a >>= 1;
292       e++;
293    }
294    r.m = DIV32_16(a,b.m);
295    r.e = e-b.e;
296    return r;
297 }
298
299
300 static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b)
301 {
302    int e=0;
303    spx_float_t r;
304    /* FIXME: Handle the sign */
305    if (a==0)
306    {
307       return FLOAT_ZERO;
308    }
309    while (b>32767)
310    {
311       b >>= 1;
312       e--;
313    }
314    while (a<SHL32(b,14))
315    {
316       a <<= 1;
317       e--;
318    }
319    while (a>=SHL32(b-1,15))
320    {
321       a >>= 1;
322       e++;
323    }
324    r.m = DIV32_16(a,b);
325    r.e = e;
326    return r;
327 }
328
329 static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b)
330 {
331    int e=0;
332    spx_int32_t num;
333    spx_float_t r;
334    num = a.m;
335    while (a.m >= b.m)
336    {
337       e++;
338       a.m >>= 1;
339    }
340    num = num << (15-e);
341    r.m = DIV32_16(num,b.m);
342    r.e = a.e-b.e-15+e;
343    return r;
344 }
345
346 static inline spx_float_t FLOAT_SQRT(spx_float_t a)
347 {
348    spx_float_t r;
349    spx_int32_t m;
350    m = a.m << 14;
351    r.e = a.e - 14;
352    if (r.e & 1)
353    {
354       r.e -= 1;
355       m <<= 1;
356    }
357    r.e >>= 1;
358    r.m = spx_sqrt(m);
359    return r;
360 }
361
362 #else
363
364 #define spx_float_t float
365 #define FLOAT_ZERO 0.f
366 #define FLOAT_ONE 1.f
367 #define FLOAT_HALF 0.5f
368 #define PSEUDOFLOAT(x) (x)
369 #define FLOAT_MULT(a,b) ((a)*(b))
370 #define FLOAT_MUL32(a,b) ((a)*(b))
371 #define FLOAT_DIV32(a,b) ((a)/(b))
372 #define FLOAT_EXTRACT16(a) (a)
373 #define FLOAT_EXTRACT32(a) (a)
374 #define FLOAT_ADD(a,b) ((a)+(b))
375 #define FLOAT_SUB(a,b) ((a)-(b))
376 #define REALFLOAT(x) (x)
377 #define FLOAT_DIV32_FLOAT(a,b) ((a)/(b))
378 #define FLOAT_MUL32U(a,b) ((a)*(b))
379 #define FLOAT_SHL(a,b) (a)
380 #define FLOAT_LT(a,b) ((a)<(b))
381 #define FLOAT_GT(a,b) ((a)>(b))
382 #define FLOAT_DIVU(a,b) ((a)/(b))
383 #define FLOAT_SQRT(a) (spx_sqrt(a))
384
385 #endif
386
387 #endif