fixed-point: acos function approximated with fixed-point arithmetic
authorjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Fri, 7 Nov 2003 08:34:14 +0000 (08:34 +0000)
committerjm <jm@0101bb08-14d6-0310-b084-bc0e0c8e3800>
Fri, 7 Nov 2003 08:34:14 +0000 (08:34 +0000)
git-svn-id: http://svn.xiph.org/trunk/speex@5541 0101bb08-14d6-0310-b084-bc0e0c8e3800

libspeex/filters.c
libspeex/lsp.c
libspeex/math_approx.c
libspeex/math_approx.h
libspeex/nb_celp.c

index d4e8903..88e3919 100644 (file)
@@ -118,7 +118,7 @@ spx_word16_t compute_rms(spx_sig_t *x, int len)
    }
    
    /*FIXME: remove division*/
-   return SHR(SHL((spx_word32_t)sqroot(1+sum/len),(sig_shift+3)),SIG_SHIFT);
+   return SHR(SHL((spx_word32_t)spx_sqrt(1+sum/len),(sig_shift+3)),SIG_SHIFT);
 }
 
 
index f27f911..98854d6 100644 (file)
@@ -46,7 +46,7 @@ Modified by Jean-Marc Valin
 #include <math.h>
 #include "lsp.h"
 #include "stack_alloc.h"
-
+#include "math_approx.h"
 
 #ifndef M_PI
 #define M_PI           3.14159265358979323846  /* pi */
@@ -78,8 +78,8 @@ static spx_word16_t cos_32(spx_word16_t x)
 
 #define FREQ_SCALE 16384
 #define ANGLE2X(a) (SHL(cos_32(a),2))
-#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)
-
+//#define X2ANGLE(x) (acos(.00006103515625*(x))*LSP_SCALING)
+#define X2ANGLE(x) (spx_acos(x))
 /* maybe we could approximate acos as 
    sqrt(6.7349563814-5.213449731*sqrt(0.6688572663+x))
 */
index 9d42903..a40409f 100644 (file)
@@ -41,7 +41,7 @@
 #define C2 -12627
 #define C3 4215
 
-spx_word16_t sqroot(spx_word32_t x)
+spx_word16_t spx_sqrt(spx_word32_t x)
 {
    int k=0;
    spx_word32_t rt;
@@ -79,12 +79,47 @@ spx_word16_t sqroot(spx_word32_t x)
       k++;
       }
 #endif
+   while (x<4096)
+   {
+      x<<=2;
+      k--;
+   }
    rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3)))))));
-   rt <<= k;
+   if (k>0)
+      rt <<= k;
+   else
+      rt >>= -k;
    rt >>=7;
    return rt;
 }
 
 /* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */
 
+
+#define A1 16469
+#define A2 2242
+#define A3 1486
+
+spx_word16_t spx_acos(spx_word16_t x)
+{
+   int s=0;
+   spx_word16_t ret;
+   spx_word32_t sq;
+   if (x<0)
+   {
+      s=1;
+      x = -x;
+   }
+   x = 16384-x;
+   
+   x = x >> 1;
+   sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3))))));
+   ret = spx_sqrt(SHL(sq,13));
+   
+   /*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/
+   if (s)
+      ret = 25736-ret;
+   return ret;
+}
+
 #endif
index bf48900..1c84f5c 100644 (file)
 
 #include "misc.h"
 
-float speex_cos(float x);
-
 #ifdef FIXED_POINT
-spx_word16_t sqroot(spx_word32_t x);
+spx_word16_t spx_sqrt(spx_word32_t x);
+spx_word16_t spx_acos(spx_word16_t x);
 #else
-#define sqroot sqrt
+#define spx_sqrt sqrt
+#define spx_acos acos
 #endif
 
 #endif
index e2e848b..7af1aaf 100644 (file)
@@ -877,7 +877,7 @@ int nb_encode(void *state, short *in, SpeexBits *bits)
    /* The next frame will not be the first (Duh!) */
    st->first = 0;
 
-   {
+   if (0) {
       float ener=0, err=0;
       float snr;
       for (i=0;i<st->frameSize;i++)