DRC: math: Replace exponential function for performance

Replace exp_small_fixed() with sofm_exp_int32() for DRC performance. Included supporting modification to replace sofm_exp_int32() and moved exp_small_fixed for future usage. Signed-off-by: shastry <[email protected]>
ShriramShastry · Nov 3, 2023 · afb34d5 · afb34d5
1 parent 20705ad
commit afb34d5
Show file tree

Hide file tree

Showing 6 changed files with 57 additions and 42 deletions.
diff --git a/src/audio/Kconfig b/src/audio/Kconfig
@@ -280,6 +280,7 @@ config COMP_DRC
 	select NUMBERS_NORM
 	select MATH_DECIBELS
 	select COMP_BLOB
+	select MATH_EXP
 	default n
 	help
 	  Select for Dynamic Range Compressor (DRC) component. A DRC can be used

diff --git a/src/include/sof/math/decibels.h b/src/include/sof/math/decibels.h
@@ -17,6 +17,7 @@
 #define DB2LIN_FIXED_OUTPUT_QY 20
 
 int32_t exp_fixed(int32_t x); /* Input is Q5.27, output is Q12.20 */
+int32_t exp_small_fixed(int32_t x); /* Input is Q3.29, Output is Q9.23 */
 int32_t db2lin_fixed(int32_t x); /* Input is Q8.24, output is Q12.20 */
 
 #endif /* __SOF_MATH_DECIBELS_H__ */
diff --git a/src/math/Kconfig b/src/math/Kconfig
@@ -47,6 +47,15 @@ config MATH_EXP
 	  an input range of -5 to +5 gives positive numbers between 0.00673794699908547 and
 	  148.413159102577. The precision of this function is 1e-4.
 
+config MATH_EXP_SMALL_FXD
+	bool "Small Exponential functions"
+	default n
+	help
+	  By selecting this, the 32-bit exp_small_fixed() function can be used to calculate
+	  exponential values. With a mean thdn of -131.205(dBc), an exponential function with
+	  an input range of -2 to +2 gives positive numbers between 0.135335255 and
+	  7.38905609. The precision of this function is 1e-6.
+
 config NATURAL_LOGARITHM_FIXED
 	bool "Natural Logarithm function"
 	default n

diff --git a/src/math/decibels.c b/src/math/decibels.c
@@ -6,47 +6,14 @@
 
 #include <sof/audio/format.h>
 #include <sof/math/decibels.h>
+#include <sof/math/exp_fcn.h>
 #include <stdint.h>
 
 #define ONE_Q20         Q_CONVERT_FLOAT(1.0, 20)	  /* Use Q12.20 */
-#define ONE_Q23         Q_CONVERT_FLOAT(1.0, 23)	  /* Use Q9.23 */
 #define TWO_Q27         Q_CONVERT_FLOAT(2.0, 27)	  /* Use Q5.27 */
 #define MINUS_TWO_Q27   Q_CONVERT_FLOAT(-2.0, 27)	  /* Use Q5.27 */
 #define LOG10_DIV20_Q27 Q_CONVERT_FLOAT(0.1151292546, 27) /* Use Q5.27 */
 
-/* Exponent function for small values of x. This function calculates
- * fairly accurately exponent for x in range -2.0 .. +2.0. The iteration
- * uses first 11 terms of Taylor series approximation for exponent
- * function. With the current scaling the numerator just remains under
- * 64 bits with the 11 terms.
- *
- * See https://en.wikipedia.org/wiki/Exponential_function#Computation
- *
- * The input is Q3.29
- * The output is Q9.23
- */
-
-static int32_t exp_small_fixed(int32_t x)
-{
-	int64_t p;
-	int64_t num = Q_SHIFT_RND(x, 29, 23);
-	int32_t y0 = (int32_t)num;
-	int32_t den = 1;
-	int32_t inc;
-	int k;
-
-	/* Numerator is x^k, denominator is k! */
-	for (k = 2; k < 12; k++) {
-		p = num * x; /* Q9.23 x Q3.29 -> Q12.52 */
-		num = Q_SHIFT_RND(p, 52, 23);
-		den = den * k;
-		inc = (int32_t)(num / den);
-		y0 += inc;
-	}
-
-	return y0 + ONE_Q23;
-}
-
 /* Decibels to linear conversion: The function uses exp() to calculate
  * the linear value. The argument is multiplied by log(10)/20 to
  * calculate equivalent of 10^(db/20).
@@ -105,10 +72,10 @@ int32_t exp_fixed(int32_t x)
 		n++;
 	}
 
-	/* exp_small_fixed() input is Q3.29, while x1 is Q5.27
-	 * exp_small_fixed() output is Q9.23, while y0 is Q12.20
+	/* sofm_exp_int32() input is Q4.28, while x1 is Q5.27
+	 * sofm_exp_int32() output is Q9.23, while y0 is Q12.20
 	 */
-	y0 = Q_SHIFT_RND(exp_small_fixed(Q_SHIFT_LEFT(xs, 27, 29)), 23, 20);
+	y0 = Q_SHIFT_RND(sofm_exp_int32(Q_SHIFT_LEFT(xs, 27, 28)), 23, 20);
 	y = ONE_Q20;
 	for (i = 0; i < (1 << n); i++)
 		y = (int32_t)Q_MULTSR_32X32((int64_t)y, y0, 20, 20, 20);

diff --git a/src/math/exp_fcn.c b/src/math/exp_fcn.c
@@ -6,6 +6,8 @@
  *
  */
 
+#include <sof/math/decibels.h>
+#include <sof/audio/format.h>
 #include <sof/math/exp_fcn.h>
 #include <sof/math/numbers.h>
 #include <sof/common.h>
@@ -20,7 +22,7 @@
 #define SOFM_BIT_MASK_Q62P2 0x4000000000000000LL
 #define SOFM_CONVERG_ERROR 28823037607936LL	// error smaller than 1e-4,1/2 ^ -44.7122876200884
 #define SOFM_BIT_MASK_LOW_Q27P5 0x8000000
-#define SOFM_QUOTIENT_SCALE BIT(30)
+#define SOFM_QUOTIENT_SCALE 0x400000000
 #define SOFM_TERMS_Q23P9 8388608
 #define SOFM_LSHIFT_BITS 8192
 
@@ -217,4 +219,38 @@ int32_t sofm_exp_int32(int32_t x)
 	}
 	return ts;
 }
+
+#define ONE_Q23         Q_CONVERT_FLOAT(1.0, 23)	  /* Use Q9.23 */
+
+/* Exponent function for small values of x. This function calculates
+ * fairly accurately exponent for x in range -2.0 .. +2.0. The iteration
+ * uses first 11 terms of Taylor series approximation for exponent
+ * function. With the current scaling the numerator just remains under
+ * 64 bits with the 11 terms.
+ *
+ * See https://en.wikipedia.org/wiki/Exponential_function#Computation
+ *
+ * The input is Q3.29
+ * The output is Q9.23
+ */
+int32_t exp_small_fixed(int32_t x)
+{
+	int64_t p;
+	int64_t num = Q_SHIFT_RND(x, 29, 23);
+	int32_t y0 = (int32_t)num;
+	int32_t den = 1;
+	int32_t inc;
+	int k;
+
+	/* Numerator is x^k, denominator is k! */
+	for (k = 2; k < 12; k++) {
+		p = num * x; /* Q9.23 x Q3.29 -> Q12.52 */
+		num = Q_SHIFT_RND(p, 52, 23);
+		den = den * k;
+		inc = (int32_t)(num / den);
+		y0 += inc;
+	}
+
+	return y0 + ONE_Q23;
+}
 #endif
diff --git a/src/math/exp_fcn_hifi.c b/src/math/exp_fcn_hifi.c
@@ -29,7 +29,7 @@
 #define SOFM_CONVERG_ERROR 28823037624320LL /* error smaller than 1e-4,1/2 ^ -44.7122876209085 */
 #define SOFM_BIT_MASK_LOW_Q27P5 0x0000000008000000
 #define SOFM_BIT_MASK_Q62P2 0x4000000000000000LL
-#define SOFM_QUOTIENT_SCALE BIT(30)
+#define SOFM_QUOTIENT_SCALE 0x400000000
 #define SOFM_TERMS_Q23P9 0x800000
 #define SOFM_LSHIFT_BITS 0x2000
 /*
@@ -89,7 +89,7 @@ static void mul_s64(ae_int64 in_0, ae_int64 in_1, ae_int64 *__restrict__ ptroutb
 		    ae_int64 *__restrict__ ptroutbitslo)
 {
 	ae_int64 producthihi, producthilo, productlolo;
-	ae_int64 producthi, productlo, product_hl_lh_h, product_hl_lh_l, carry;
+	ae_int64 producthi, product_hl_lh_h, product_hl_lh_l, carry;
 
 #if (SOFM_EXPONENTIAL_HIFI4 == 1 || SOFM_EXPONENTIAL_HIFI5 == 1)
 
@@ -130,6 +130,7 @@ static void mul_s64(ae_int64 in_0, ae_int64 in_1, ae_int64 *__restrict__ ptroutb
 	ae_int64 producthi_1c;
 	ae_int64 producthi_2c;
 	ae_int64 productlo_2c;
+	ae_int64 productlo;
 
 	ae_int64 s0 = AE_SRLI64(in_0, 63);
 	ae_int64 s1 = AE_SRLI64(in_1, 63);
@@ -194,7 +195,7 @@ static int64_t lomul_s64_sr_sat_near(int64_t a, int64_t b)
 	return AE_ADD64(temp, roundup);
 }
 
-static ae_int64 onebyfact_Q63[19] = {
+int64_t onebyfact_Q63[19] = {
 		4611686018427387904LL,
 		1537228672809129301LL,
 		384307168202282325LL,
@@ -239,7 +240,7 @@ int32_t sofm_exp_int32(int32_t x)
 	ae_int64 onebyfact;
 	ae_int64 temp;
 
-	ae_int64 *ponebyfact_Q63 = &onebyfact_Q63[0];
+	ae_int64 *ponebyfact_Q63 = (ae_int64 *)&onebyfact_Q63[0];
 	ae_int64 ts = SOFM_TERMS_Q23P9;
 	ae_int64 mp = (x + SOFM_LSHIFT_BITS) >> 14; /* x in Q50.14 */;
 	xtbool flag;