diff --git a/src/audio/Kconfig b/src/audio/Kconfig index 1ea03b65e603..1e7e622dc891 100644 --- a/src/audio/Kconfig +++ b/src/audio/Kconfig @@ -278,8 +278,8 @@ config COMP_DRC bool "Dynamic Range Compressor component" select CORDIC_FIXED 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/audio/drc/drc_generic.c b/src/audio/drc/drc_generic.c index 472431a0e408..8acce480a238 100644 --- a/src/audio/drc/drc_generic.c +++ b/src/audio/drc/drc_generic.c @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -35,7 +36,7 @@ static int32_t knee_curveK(const struct sof_drc_params *p, int32_t x) * beta = -expf(k * linear_threshold) / k * gamma = -k * x */ - knee_exp_gamma = exp_fixed(Q_MULTSR_32X32((int64_t)x, -p->K, 31, 20, 27)); /* Q12.20 */ + knee_exp_gamma = sofm_exp_fixed(Q_MULTSR_32X32((int64_t)x, -p->K, 31, 20, 27)); /* Q12.20 */ return p->knee_alpha + Q_MULTSR_32X32((int64_t)p->knee_beta, knee_exp_gamma, 24, 20, 24); } @@ -65,8 +66,10 @@ static int32_t volume_gain(const struct sof_drc_params *p, int32_t x) * => y/x = ratio_base * x^(s - 1) * => y/x = ratio_base * e^(log(x) * (s - 1)) */ - exp_knee = exp_fixed(Q_MULTSR_32X32((int64_t)drc_log_fixed(Q_SHIFT_RND(x, 31, 26)), - (p->slope - ONE_Q30), 26, 30, 27)); /* Q12.20 */ + exp_knee = sofm_exp_fixed( + Q_MULTSR_32X32((int64_t) + drc_log_fixed(Q_SHIFT_RND(x, 31, 26)), + (p->slope - ONE_Q30), 26, 30, 27)); /* Q12.20 */ y = Q_MULTSR_32X32((int64_t)p->ratio_base, exp_knee, 30, 20, 30); } diff --git a/src/audio/drc/drc_hifi3.c b/src/audio/drc/drc_hifi3.c index 4613fff0bfac..82f403fc435f 100644 --- a/src/audio/drc/drc_hifi3.c +++ b/src/audio/drc/drc_hifi3.c @@ -10,6 +10,7 @@ #include #include #include +#include #include #include @@ -41,7 +42,7 @@ static int32_t knee_curveK(const struct sof_drc_params *p, int32_t x) * gamma = -k * x */ gamma = drc_mult_lshift(x, -p->K, drc_get_lshift(31, 20, 27)); - knee_exp_gamma = exp_fixed(gamma); + knee_exp_gamma = sofm_exp_fixed(gamma); knee_curve_k = drc_mult_lshift(p->knee_beta, knee_exp_gamma, drc_get_lshift(24, 20, 24)); knee_curve_k = AE_ADD32(knee_curve_k, p->knee_alpha); return knee_curve_k; @@ -77,7 +78,7 @@ static int32_t volume_gain(const struct sof_drc_params *p, int32_t x) tmp = AE_SRAI32R(x, 5); /* Q1.31 -> Q5.26 */ tmp = drc_log_fixed(tmp); /* Q6.26 */ tmp2 = AE_SUB32(p->slope, ONE_Q30); /* Q2.30 */ - exp_knee = exp_fixed(drc_mult_lshift(tmp, tmp2, drc_get_lshift(26, 30, 27))); + exp_knee = sofm_exp_fixed(drc_mult_lshift(tmp, tmp2, drc_get_lshift(26, 30, 27))); y = drc_mult_lshift(p->ratio_base, exp_knee, drc_get_lshift(30, 20, 30)); } diff --git a/src/audio/drc/drc_math_generic.c b/src/audio/drc/drc_math_generic.c index 0f927d86fe9e..4f78f0eb2bf8 100644 --- a/src/audio/drc/drc_math_generic.c +++ b/src/audio/drc/drc_math_generic.c @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -219,8 +220,6 @@ inline int32_t drc_inv_fixed(int32_t x, int32_t precision_x, int32_t precision_y #undef qc } -#endif /* DRC_GENERIC */ - /* * Input x is Q6.26; valid range: (0.0, 32.0); x <= 0 is not supported * y is Q2.30: (-2.0, 2.0) @@ -233,8 +232,8 @@ inline int32_t drc_pow_fixed(int32_t x, int32_t y) return 0; /* x^y = expf(y * log(x)) */ - return exp_fixed(q_mult(y, drc_log_fixed(x), 30, 26, 27)); + return sofm_exp_fixed(q_mult(y, drc_log_fixed(x), 30, 26, 27)); } - +#endif #undef q_multq #undef q_mult diff --git a/src/include/sof/math/decibels.h b/src/include/sof/math/decibels.h index a0336e9be9b9..c9465ebac362 100644 --- a/src/include/sof/math/decibels.h +++ b/src/include/sof/math/decibels.h @@ -11,12 +11,6 @@ #include -#define EXP_FIXED_INPUT_QY 27 -#define EXP_FIXED_OUTPUT_QY 20 -#define DB2LIN_FIXED_INPUT_QY 24 -#define DB2LIN_FIXED_OUTPUT_QY 20 - -int32_t exp_fixed(int32_t x); /* Input is Q5.27, output is Q12.20 */ int32_t db2lin_fixed(int32_t x); /* Input is Q8.24, output is Q12.20 */ #endif /* __SOF_MATH_DECIBELS_H__ */ diff --git a/src/include/sof/math/exp_fcn.h b/src/include/sof/math/exp_fcn.h index 74466f6bfb3a..b305c0b80ea2 100644 --- a/src/include/sof/math/exp_fcn.h +++ b/src/include/sof/math/exp_fcn.h @@ -26,6 +26,6 @@ #endif -int32_t sofm_exp_int32(int32_t x); - +int32_t sofm_exp_int32(int32_t x); /* Input is Q4.28, Output is Q9.23 */ +int32_t sofm_exp_fixed(int32_t x); /* Input is Q5.27, output is Q12.20 */ #endif 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/auditory/auditory.c b/src/math/auditory/auditory.c index 57c1fe0ee3f8..cbc52e39d0ba 100644 --- a/src/math/auditory/auditory.c +++ b/src/math/auditory/auditory.c @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -80,7 +81,7 @@ int16_t psy_mel_to_hz(int16_t mel) return 0; exp_arg = Q_MULTSR_32X32((int64_t)mel, ONE_OVER_MELDIV_Q31, 2, 31, 27); - exp = exp_fixed(exp_arg) - ONE_Q20; + exp = sofm_exp_fixed(exp_arg) - ONE_Q20; hz = Q_MULTSR_32X32((int64_t)exp, 700, 20, 0, 0); return hz; } diff --git a/src/math/decibels.c b/src/math/decibels.c index 560ce6e6964c..9f24d6f06036 100644 --- a/src/math/decibels.c +++ b/src/math/decibels.c @@ -6,47 +6,11 @@ #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). @@ -68,50 +32,5 @@ int32_t db2lin_fixed(int32_t db) /* Q8.24 x Q5.27, result needs to be Q5.27 */ arg = (int32_t)Q_MULTSR_32X32((int64_t)db, LOG10_DIV20_Q27, 24, 27, 27); - return exp_fixed(arg); -} - -/* Fixed point exponent function for approximate range -11.5 .. 7.6 - * that corresponds to decibels range -100 .. +66 dB. - * - * The functions uses rule exp(x) = exp(x/2) * exp(x/2) to reduce - * the input argument for private small value exp() function that is - * accurate with input range -2.0 .. +2.0. The number of possible - * divisions by 2 is computed into variable n. The returned value is - * exp()^(2^n). - * - * Input is Q5.27, -16.0 .. +16.0, but note the input range limitation - * Output is Q12.20, 0.0 .. +2048.0 - */ - -int32_t exp_fixed(int32_t x) -{ - int32_t xs; - int32_t y; - int32_t y0; - int i; - int n = 0; - - if (x < Q_CONVERT_FLOAT(-11.5, 27)) - return 0; - - if (x > Q_CONVERT_FLOAT(7.6245, 27)) - return INT32_MAX; - - /* x is Q5.27 */ - xs = x; - while (xs >= TWO_Q27 || xs <= MINUS_TWO_Q27) { - xs >>= 1; - 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 - */ - y0 = Q_SHIFT_RND(exp_small_fixed(Q_SHIFT_LEFT(xs, 27, 29)), 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); - - return y; + return sofm_exp_fixed(arg); } diff --git a/src/math/exp_fcn.c b/src/math/exp_fcn.c index ada4b573e4f1..eaeaf7758029 100644 --- a/src/math/exp_fcn.c +++ b/src/math/exp_fcn.c @@ -6,6 +6,9 @@ * */ +#include +#include +#include #include #include #include @@ -20,7 +23,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 @@ -176,6 +179,7 @@ static inline int64_t lomul_s64_sr_sat_near(int64_t a, int64_t b) *| 32 | 28 | 1 | 32 | 23 | 0 | 4.28 | 9.23 | *+------------------+-----------------+--------+--------+ */ +//static inline int32_t sofm_exp_int32(int32_t x) int32_t sofm_exp_int32(int32_t x) { uint64_t ou0Hi; @@ -217,4 +221,52 @@ int32_t sofm_exp_int32(int32_t x) } return ts; } + +#define ONE_Q20 Q_CONVERT_FLOAT(1.0, 20) /* Use Q12.20 */ +#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 */ + +/* Fixed point exponent function for approximate range -11.5 .. 7.6 + * that corresponds to decibels range -100 .. +66 dB. + * + * The functions uses rule exp(x) = exp(x/2) * exp(x/2) to reduce + * the input argument for private small value exp() function that is + * accurate with input range -2.0 .. +2.0. The number of possible + * divisions by 2 is computed into variable n. The returned value is + * exp()^(2^n). + * + * Input is Q5.27, -16.0 .. +16.0, but note the input range limitation + * Output is Q12.20, 0.0 .. +2048.0 + */ +int32_t sofm_exp_fixed(int32_t x) +{ + int32_t xs; + int32_t y; + int32_t y0; + int i; + int n = 0; + + if (x < Q_CONVERT_FLOAT(-11.5, 27)) + return 0; + + if (x > Q_CONVERT_FLOAT(7.6245, 27)) + return INT32_MAX; + + /* x is Q5.27 */ + xs = x; + while (xs >= TWO_Q27 || xs <= MINUS_TWO_Q27) { + xs >>= 1; + n++; + } + + /* 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(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); + + return y; +} #endif diff --git a/src/math/exp_fcn_hifi.c b/src/math/exp_fcn_hifi.c index 2a65d4932864..a4875fb0aa9b 100644 --- a/src/math/exp_fcn_hifi.c +++ b/src/math/exp_fcn_hifi.c @@ -6,13 +6,14 @@ * */ +#include +#include #include +#include #include #include #include #include -#include -#include #include #if defined(SOFM_EXPONENTIAL_HIFI3) || defined(SOFM_EXPONENTIAL_HIFI4) || \ @@ -26,10 +27,12 @@ #include #endif +#include + #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 +92,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 +133,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 +198,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, @@ -231,6 +235,7 @@ static ae_int64 onebyfact_Q63[19] = { *| 32 | 28 | 1 | 32 | 23 | 0 | 4.28 | 9.23 | *+------------------+-----------------+--------+--------+ */ + int32_t sofm_exp_int32(int32_t x) { ae_int64 outhi; @@ -239,7 +244,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; @@ -278,4 +283,124 @@ int32_t sofm_exp_int32(int32_t x) return AE_MOVAD32_L(AE_MOVINT32X2_FROMINT64(ts)); } + +/* Xtensa macros */ +/* Convert a float number to fractional Qnx.ny format. Note that there is no + * check for nx+ny number of bits to fit the word length of int. The parameter + * qy must be 31 or less. + */ +//#define XT_TRUNC_S TRUNC_S +//static int XT_Q_CONVERT_FLOAT(float a, int b) +//{ +// return XT_TRUNC_S(a, b);//note: b is const number +//} +/* Fractional multiplication with shift and round + * Note that the parameters px and py must be cast to (int64_t) if other type. + */ +static int XT_Q_MULTSR_32X32(int a, int b, int c, int d, int e) +{ + ae_int64 res; + int xt_o; + int shift; + + res = AE_MUL32_LL(a, b); + shift = XT_SUB(XT_ADD(c, d), XT_ADD(e, 1)); + res = AE_SRAA64(res, shift); + res = AE_ADD64(res, 1); + res = AE_SRAI64(res, 1); + xt_o = AE_MOVINT32_FROMINT64(res); + + return xt_o; + +} +/* A macro for Q-shifts */ +static int XT_Q_SHIFT_RND(int a, int b, int c) +{ + ae_int32 res; + int shift; + + shift = XT_SUB(b, XT_ADD(c, 1)); + res = AE_SRAA32(a, shift); + res = AE_ADD32(res, 1); + res = AE_SRAI32(res, 1); + + return res; +} +/* Alternative version since compiler does not allow (x >> -1) */ +static int XT_Q_SHIFT_LEFT(int a, int b, int c) +{ + ae_int32 xt_o; + int shift; + + shift = XT_SUB(b, c); + xt_o = AE_SLAA32(a, shift); + + return xt_o; + +} + +#define ONE_Q20 Q_CONVERT_FLOAT(1.0, 20) /* Use Q12.20 */ +#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 */ + + +/* Fixed point exponent function for approximate range -11.5 .. 7.6 + * that corresponds to decibels range -100 .. +66 dB. + * + * The functions uses rule exp(x) = exp(x/2) * exp(x/2) to reduce + * the input argument for private small value exp() function that is + * accurate with input range -2.0 .. +2.0. The number of possible + * divisions by 2 is computed into variable n. The returned value is + * exp()^(2^n). + * + * Input is Q5.27, -16.0 .. +16.0, but note the input range limitation + * Output is Q12.20, 0.0 .. +2048.0 + */ +int32_t sofm_exp_fixed(int32_t x) +{ + int32_t xs; + int32_t y; + int32_t y0; + int i; + int n = 0; + + if (x < Q_CONVERT_FLOAT(-11.5, 27)) + return 0; + + if (x > Q_CONVERT_FLOAT(7.6245, 27)) + return INT32_MAX; + + /* x is Q5.27 */ + xs = x; + while (xs >= TWO_Q27 || xs <= MINUS_TWO_Q27) { + xs >>= 1; + n++; + } + + /* 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 = XT_Q_SHIFT_RND(sofm_exp_int32(XT_Q_SHIFT_LEFT(xs, 27, 28)), 23, 20); + y = ONE_Q20; + for (i = 0; i < (1 << n); i++) + y = (int32_t)XT_Q_MULTSR_32X32((int64_t)y, y0, 20, 20, 20); + + return y; +} + +#define q_mult(a, b, qa, qb, qy) ((int32_t)XT_Q_MULTSR_32X32((int64_t)(a), b, qa, qb, qy)) +/* + * Input x is Q6.26; valid range: (0.0, 32.0); x <= 0 is not supported + * y is Q2.30: (-2.0, 2.0) + * Output is Q12.20: max 2048.0 + */ +int32_t drc_pow_fixed(int32_t x, int32_t y) +{ + /* Negative or zero input x is not supported, just return 0. */ + if (x <= 0) + return 0; + + /* x^y = expf(y * log(x)) */ + return sofm_exp_fixed(q_mult(y, drc_log_fixed(x), 30, 26, 27)); +} #endif diff --git a/src/math/window.c b/src/math/window.c index 4795112ffc0f..33bbcae38bab 100644 --- a/src/math/window.c +++ b/src/math/window.c @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -114,7 +115,7 @@ void win_povey_16b(int16_t win[], int length) /* Calculate x^0.85 as exp(0.85 * log(x)) */ x2 = (int32_t)(ln_int32((uint32_t)x1) >> 1) - WIN_LOG_2POW31_Q26; x3 = sat_int32(Q_MULTSR_32X32((int64_t)x2, WIN_085_Q31, 26, 31, 27)); /* Q5.27 */ - x4 = exp_fixed(x3); /* Q5.27 -> Q12.20 */ + x4 = sofm_exp_fixed(x3); /* Q5.27 -> Q12.20 */ /* Convert to Q1.15 */ win[n] = sat_int16(Q_SHIFT_RND(x4, 20, 15));