optimisation: defined a reciprocal square root (celt_rsqrt) for use in

find_spectral_pitch instead of using celt_rcp(celt_sqrt(x))

3 files changed, 40 insertions, 2 deletions

diff --git a/libcelt/mathops.h b/libcelt/mathops.h
index bb58207..2218429 100644
--- a/libcelt/mathops.h
+++ b/libcelt/mathops.h
@@ -74,7 +74,8 @@ static inline int find_max32(celt_word32_t *x, int len)
 
 #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)))
@@ -117,6 +118,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) */
 static inline celt_word32_t celt_sqrt(celt_word32_t x)
 {
diff --git a/libcelt/pitch.c b/libcelt/pitch.c
index 01eec28..2dab60f 100644
--- a/libcelt/pitch.c
+++ b/libcelt/pitch.c
@@ -167,7 +167,8 @@ void find_spectral_pitch(kiss_fftr_cfg fft, const struct PsyDecay *decay, const
       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])));
diff --git a/tests/mathops-test.c b/tests/mathops-test.c
index dc208da..b7ffdac 100644
--- a/tests/mathops-test.c
+++ b/tests/mathops-test.c
@@ -53,9 +53,28 @@ void testsqrt(void)
    }
 }
 
+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;
 }