Add the fixed-point cos(PI*0.5*x) function (cos_norm2).

Change-Id: I08be579a4224b7b09057e2f6b2771d0a4c3567ca Signed-off-by: Jean-Marc Valin <[email protected]>
xiph · Dec 19, 2024 · ec1e887 · ec1e887
1 parent 57a17aa
commit ec1e887
Show file tree

Hide file tree

Showing 4 changed files with 63 additions and 7 deletions.
diff --git a/celt/bands.c b/celt/bands.c
@@ -1051,8 +1051,8 @@ static unsigned quant_partition(struct band_ctx *ctx, celt_norm *X,
 # ifdef ENABLE_QEXT
       (void)imid;
       (void)iside;
-      mid = (int)MIN32(2147483647, floor(.5+2147483648*celt_cos_norm2(sctx.itheta_q30*(1.f/(1<<30)))));
-      side = (int)MIN32(2147483647, floor(.5+2147483648*celt_cos_norm2(1.f-sctx.itheta_q30*(1.f/(1<<30)))));
+      mid = celt_cos_norm32(sctx.itheta_q30);
+      side = celt_cos_norm32((1<<30)-sctx.itheta_q30);
 # else
       mid = SHL32(EXTEND32(imid), 16);
       side = SHL32(EXTEND32(iside), 16);
@@ -1261,8 +1261,8 @@ static unsigned cubic_quant_partition(struct band_ctx *ctx, celt_norm *X, int N,
       delta = (N0-1) * 23 * ((itheta_q30>>16)-8192) >> (17-BITRES);
 
 #ifdef FIXED_POINT
-      g1 = (int)MIN32(2147483647, floor(.5+2147483648*celt_cos_norm2(itheta_q30*(1.f/(1<<30)))));
-      g2 = (int)MIN32(2147483647, floor(.5+2147483648*celt_cos_norm2(1.f-itheta_q30*(1.f/(1<<30)))));
+      g1 = celt_cos_norm32(itheta_q30);
+      g2 = celt_cos_norm32((1<<30)-itheta_q30);
 #else
       g1 = celt_cos_norm2(itheta_q30*(1.f/(1<<30)));
       g2 = celt_cos_norm2(1.f-itheta_q30*(1.f/(1<<30)));
@@ -1459,8 +1459,8 @@ static unsigned quant_band_stereo(struct band_ctx *ctx, celt_norm *X, celt_norm
 # ifdef ENABLE_QEXT
    (void)imid;
    (void)iside;
-   mid = (int)MIN32(2147483647, floor(.5+2147483648*cos(.5*M_PI*sctx.itheta_q30*(1.f/(1<<30)))));
-   side = (int)MIN32(2147483647, floor(.5+2147483648*sin(.5*M_PI*sctx.itheta_q30*(1.f/(1<<30)))));
+   mid = celt_cos_norm32(sctx.itheta_q30);
+   side = celt_cos_norm32((1<<30)-sctx.itheta_q30);
 # else
    mid = SHL32(EXTEND32(imid), 16);
    side = SHL32(EXTEND32(iside), 16);

diff --git a/celt/mathops.c b/celt/mathops.c
@@ -199,6 +199,26 @@ opus_val16 celt_cos_norm(opus_val32 x)
    }
 }
 
+/* Calculates the cosine of (PI*0.5*x) where the input x ranges from -1 to 1 and
+ * is in Q30 format. The output will also be in Q31 format. */
+opus_val32 celt_cos_norm32(opus_val32 x)
+{
+   static const opus_val32 COS_NORM_COEFF_A0 = 134217720;   /* Q27 */
+   static const opus_val32 COS_NORM_COEFF_A1 = -662336704;  /* Q29 */
+   static const opus_val32 COS_NORM_COEFF_A2 = 544710848;   /* Q31 */
+   static const opus_val32 COS_NORM_COEFF_A3 = -178761936;  /* Q33 */
+   static const opus_val32 COS_NORM_COEFF_A4 = 29487206;    /* Q35 */
+   opus_int32 x_sq_q29, tmp;
+   /* The expected x is in the range of [-1.0f, 1.0f] */
+   celt_sig_assert((x >= -1073741824) && (x <= 1073741824));
+   x_sq_q29 = MULT32_32_Q31(x, x);
+   /* Split evaluation in steps to avoid exploding macro expansion. */
+   tmp = ADD32(COS_NORM_COEFF_A3, MULT32_32_Q31(x_sq_q29, COS_NORM_COEFF_A4));
+   tmp = ADD32(COS_NORM_COEFF_A2, MULT32_32_Q31(x_sq_q29, tmp));
+   tmp = ADD32(COS_NORM_COEFF_A1, MULT32_32_Q31(x_sq_q29, tmp));
+   return SHL32(ADD32(COS_NORM_COEFF_A0, MULT32_32_Q31(x_sq_q29, tmp)), 4);
+}
+
 /** Reciprocal approximation (Q15 input, Q16 output) */
 opus_val32 celt_rcp(opus_val32 x)
 {

diff --git a/celt/mathops.h b/celt/mathops.h
@@ -371,6 +371,8 @@ opus_val32 celt_sqrt(opus_val32 x);
 
 opus_val16 celt_cos_norm(opus_val32 x);
 
+opus_val32 celt_cos_norm32(opus_val32 x);
+
 /** Base-2 logarithm approximation (log2(x)). (Q14 input, Q10 output) */
 static OPUS_INLINE opus_val16 celt_log2(opus_val32 x)
 {

diff --git a/celt/tests/test_unit_mathops.c b/celt/tests/test_unit_mathops.c
@@ -519,6 +519,40 @@ void testatan(void)
    testatan_norm();
    testatan2p_norm();
 }
+
+void test_cos_norm32(void)
+{
+   float error = -1;
+   float max_error = -2;
+   float error_threshold = 1e-07;
+   float fx = 0;
+   opus_int32 x = 0;
+   int q_input = 30;
+   int q_output = 31;
+   for (fx = -1.0f; fx <= 1.0f; fx += 0.007f)
+   {
+      x = DOUBLE_TO_FIX_INT(fx, q_input);
+      error = fabs(cos(1.5707963267948966 * FIX_INT_TO_DOUBLE(x, q_input)) -
+                   FIX_INT_TO_DOUBLE(celt_cos_norm32(x), q_output));
+      if (error > max_error)
+      {
+         max_error = error;
+      }
+      if (error > error_threshold)
+      {
+         fprintf(stderr,
+                 "celt_cos_norm32 failed: error: [%.5e > %.5e] (x = %f)\n",
+                 error, error_threshold, FIX_INT_TO_DOUBLE(x, DB_SHIFT));
+         ret = 1;
+      }
+   }
+   fprintf(stdout, "celt_cos_norm32 max_error: %.7e\n", max_error);
+}
+
+void test_cos(void)
+{
+   test_cos_norm32();
+}
 #endif
 
 int main(void)
@@ -530,14 +564,14 @@ int main(void)
    testlog2();
    testexp2();
    testexp2log2();
+   test_cos();
 #ifdef FIXED_POINT
    testilog2();
    testlog2_db();
    testexp2_db();
    testrsqrt();
    testatan();
 #else
-   test_cos();
    test_atan2();
 #endif
    return ret;