Commit 189acec5 authored by Jean-Marc Valin's avatar Jean-Marc Valin
Browse files

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

find_spectral_pitch instead of using celt_rcp(celt_sqrt(x))
parent 385795ed
......@@ -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)
{
......
......@@ -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])));
......
......@@ -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;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment