Update celt_log2_db fixed point to match new celt_log2 floating point.

Change-Id: Ibd0d16918272cb568923d384475e139dc312c61b
Signed-off-by: Jean-Marc Valin <jeanmarcv@google.com>

celt/mathops.h

@@ -325,8 +325,49 @@ static OPUS_INLINE opus_val32 celt_exp2(opus_val16 x)
 
 #ifdef ENABLE_QEXT
 
+/* Calculates the base-2 logarithm of a Q14 input value. The result is returned
+ * in Q(DB_SHIFT). If the input value is 0, the function will output -32.0f. */
 static OPUS_INLINE opus_val32 celt_log2_db(opus_val32 x) {
-   return (int)floor(.5 + (1<<DB_SHIFT) * 1.4426950409f*log(x/(float)(1<<14)));
+   /* Q30 */
+   static const opus_val32 log2_x_norm_coeff[8] = {
+      1073741824, 954437184, 858993472, 780903168,
+      715827904,  660764224, 613566784, 572662336};
+   /* Q24 */
+   static const opus_val32 log2_y_norm_coeff[8] = {
+      0,       2850868,  5401057,  7707983,
+      9814042, 11751428, 13545168, 15215099};
+   static const opus_val32 LOG2_COEFF_A0 = 1467383;     /* Q24 */
+   static const opus_val32 LOG2_COEFF_A1 = 182244800;   /* Q27 */
+   static const opus_val32 LOG2_COEFF_A2 = -21440512;   /* Q25 */
+   static const opus_val32 LOG2_COEFF_A3 = 107903336;   /* Q28 */
+   static const opus_val32 LOG2_COEFF_A4 = -610217024;  /* Q31 */
+
+   opus_int32 integer, norm_coeff_idx, tmp;
+   opus_val32 mantissa;
+   if (x==0) {
+      return -536870912; /* -32.0f */
+   }
+   integer =  SUB32(celt_ilog2(x), 14);  /* Q0 */
+   mantissa = VSHR32(x, integer + 14 - 29);  /* Q29 */
+   norm_coeff_idx = SHR32(mantissa, 29 - 3) & 0x7;
+   /* mantissa is in Q28 (29 + Q_NORM_CONST - 31 where Q_NORM_CONST is Q30)
+    * 285212672 (Q28) is 1.0625f. */
+   mantissa = SUB32(MULT32_32_Q31(mantissa, log2_x_norm_coeff[norm_coeff_idx]),
+                    285212672);
+
+   /* q_a3(Q28): q_mantissa + q_a4 - 31
+    * q_a2(Q25): q_mantissa + q_a3 - 31
+    * q_a1(Q27): q_mantissa + q_a2 - 31 + 5
+    * q_a0(Q24): q_mantissa + q_a1 - 31
+    * where  q_mantissa is Q28 */
+   /* Split evaluation in steps to avoid exploding macro expansion. */
+   tmp = MULT32_32_Q31(mantissa, LOG2_COEFF_A4);
+   tmp = MULT32_32_Q31(mantissa, ADD32(LOG2_COEFF_A3, tmp));
+   tmp = SHL32(MULT32_32_Q31(mantissa, ADD32(LOG2_COEFF_A2, tmp)), 5 /* SHL32 for LOG2_COEFF_A1 */);
+   tmp = MULT32_32_Q31(mantissa, ADD32(LOG2_COEFF_A1, tmp));
+   return ADD32(log2_y_norm_coeff[norm_coeff_idx],
+          ADD32(SHL32(integer, DB_SHIFT),
+          ADD32(LOG2_COEFF_A0, tmp)));
 }
 
 static OPUS_INLINE opus_val32 celt_exp2_db_frac(opus_val32 x)

celt/tests/test_unit_mathops.c

@@ -41,6 +41,8 @@
 
 #ifdef FIXED_POINT
 #define WORD "%d"
+#define FIX_INT_TO_DOUBLE(x,q) ((double)(x) / (double)(1L << q))
+#define DOUBLE_TO_FIX_INT(x,q) (((double)x * (double)(1L << q)))
 #else
 #define WORD "%f"
 #endif
@@ -204,6 +206,35 @@ void testexp2log2(void)
    fprintf (stdout, "celt_exp2, celt_log2 max_error: %15.25e\n", max_error);
 }
 #else
+
+void testlog2_db(void)
+{
+#if defined(ENABLE_QEXT)
+   /* celt_log2_db test */
+   float error = -1;
+   float max_error = -2;
+   float error_threshold = 2.e-07;
+   opus_int32 x = 0;
+   int q_input = 14;
+   for (x = 8; x < 1073741824; x += (x >> 3))
+   {
+      error = fabs((1.442695040888963387*log(FIX_INT_TO_DOUBLE(x, q_input))) -
+                   FIX_INT_TO_DOUBLE(celt_log2_db(x), DB_SHIFT));
+      if (error > max_error)
+      {
+         max_error = error;
+      }
+      if (error > error_threshold)
+      {
+         fprintf(stderr, "celt_log2_db failed: error: [%.5e > %.5e] (x = %f)\n",
+                 error, error_threshold, FIX_INT_TO_DOUBLE(x, DB_SHIFT));
+         ret = 1;
+      }
+   }
+   fprintf(stdout, "celt_log2_db max_error: %.7e\n", max_error);
+#endif  /* defined(ENABLE_QEXT) */
+}
+
 void testlog2(void)
 {
    opus_val32 x;
@@ -283,6 +314,7 @@ int main(void)
    testexp2log2();
 #ifdef FIXED_POINT
    testilog2();
+   testlog2_db();
 #endif
    return ret;
 }