shithub: opus

Download patch

ref: 189acec5d36ae9a8031fb47b6aee8f3ef6b62f2f
parent: 385795ed7bbe22e9f0c4ddbc314fe16b88e6ae4e
author: Jean-Marc Valin <[email protected]>
date: Wed Mar 26 12:42:42 EDT 2008 

optimisation: defined a reciprocal square root (celt_rsqrt) for use in
find_spectral_pitch instead of using celt_rcp(celt_sqrt(x))


--- a/libcelt/mathops.h
+++ b/libcelt/mathops.h
@@ -74,7 +74,8 @@
 
 #ifndef FIXED_POINT
 
-#define celt_sqrt sqrt
+#define celt_sqrt(x) ((float)sqrt(x))
+#define celt_rsqrt(x) (1.f/celt_sqrt(x))
 #define celt_acos acos
 #define celt_exp exp
 #define celt_cos_norm(x) (cos((.5f*M_PI)*(x)))
@@ -115,6 +116,23 @@
 static inline celt_int16_t celt_zlog2(celt_word32_t x)
 {
    return EC_ILOG(x)-1;
+}
+
+/** Reciprocal sqrt approximation (Q30 input, Q0 output or equivalent) */
+static inline celt_word32_t celt_rsqrt(celt_word32_t x)
+{
+   int k;
+   celt_word16_t n;
+   celt_word32_t rt;
+   const celt_word16_t C[5] = {23126, -11496, 9812, -9097, 4100};
+   k = celt_ilog2(x)>>1;
+   x = VSHR32(x, (k-7)<<1);
+   /* Range of n is [-16384,32767] */
+   n = x-32768;
+   rt = 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])))))))));
+   rt = VSHR32(rt,k);
+   return rt;
 }
 
 /** Sqrt approximation (QX input, QX/2 output) */
--- a/libcelt/pitch.c
+++ b/libcelt/pitch.c
@@ -167,7 +167,8 @@
       celt_word32_t tmp;
       /*printf ("%d %d ", X[2*i]*X[2*i]+X[2*i+1]*X[2*i+1], Y[2*i]*Y[2*i]+Y[2*i+1]*Y[2*i+1]);*/
       /*n = DIV32_16(Q15ONE,celt_sqrt(EPSILON+curve[i]));*/
-      n = ROUND16(celt_rcp(celt_sqrt(EPSILON+curve[i])),16);
+      /*n = ROUND16(celt_rcp(celt_sqrt(EPSILON+curve[i])),16);*/
+      n = celt_rsqrt(EPSILON+curve[i]);
       /*printf ("%f ", n);*/
       tmp = X[2*i];
       X[2*i] = MULT16_32_Q15(n, ADD32(MULT16_16(X[2*i  ],Y[2*i  ]), MULT16_16(X[2*i+1],Y[2*i+1])));
--- a/tests/mathops-test.c
+++ b/tests/mathops-test.c
@@ -53,9 +53,28 @@
    }
 }
 
+void testrsqrt(void)
+{
+   celt_int32_t i;
+   for (i=1;i<=2000000;i++)
+   {
+      double ratio;
+      celt_word16_t val;
+      val = celt_rsqrt(i);
+      ratio = val*sqrt(i)/32768;
+      if (fabs(ratio - 1) > .05)
+      {
+         fprintf (stderr, "sqrt failed: sqrt(%d)="WORD" (ratio = %f)\n", i, val, ratio);
+         ret = 1;
+      }
+      i+= i>>10;
+   }
+}
+
 int main(void)
 {
    testdiv();
    testsqrt();
+   testrsqrt();
    return ret;
 }