From afb34d59e73dc68721680e19d5e41484561d7a45 Mon Sep 17 00:00:00 2001 From: shastry Date: Fri, 3 Nov 2023 11:12:27 +0530 Subject: [PATCH] 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 --- src/audio/Kconfig | 1 + src/include/sof/math/decibels.h | 1 + src/math/Kconfig | 9 ++++++++ src/math/decibels.c | 41 ++++----------------------------- src/math/exp_fcn.c | 38 +++++++++++++++++++++++++++++- src/math/exp_fcn_hifi.c | 9 ++++---- 6 files changed, 57 insertions(+), 42 deletions(-) diff --git a/src/audio/Kconfig b/src/audio/Kconfig index b05957a835d3..d13ec89b6404 100644 --- 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 index a0336e9be9b9..b882973a3b02 100644 --- 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 index 78a67ac868d1..766dde075360 100644 --- 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 index 560ce6e6964c..38696a72b81f 100644 --- a/src/math/decibels.c +++ b/src/math/decibels.c @@ -6,47 +6,14 @@ #include #include +#include #include #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 index ada4b573e4f1..39dcadfbadd2 100644 --- a/src/math/exp_fcn.c +++ b/src/math/exp_fcn.c @@ -6,6 +6,8 @@ * */ +#include +#include #include #include #include @@ -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 index 2a65d4932864..6784621d77c1 100644 --- 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;