diff --git a/src/fastmath.h b/src/fastmath.h index f07f3c1..fe68f14 100644 --- a/src/fastmath.h +++ b/src/fastmath.h @@ -23,6 +23,51 @@ #include +/** + * IEEE 754 Floating point representation + */ + +#define LC3_IEEE754_SIGN_SHL (31) +#define LC3_IEEE754_SIGN_MASK (1 << LC3_IEEE754_SIGN_SHL) + +#define LC3_IEEE754_EXP_SHL (23) +#define LC3_IEEE754_EXP_MASK (0xff << LC3_IEEE754_EXP_SHL) +#define LC3_IEEE754_EXP_BIAS (127) + + +/** + * Fast multiply floating-point number by integral power of 2 + * x Operand, finite number + * exp Exponent such that 2^x is finite + * return 2^exp + */ +static inline float lc3_ldexpf(float _x, int exp) { + union { float f; uint32_t u; } x = { .f = _x }; + + if (x.u & LC3_IEEE754_EXP_MASK) + x.u += exp << LC3_IEEE754_EXP_SHL; + + return x.f; +} + +/** + * Fast convert floating-point number to fractional and integral components + * x Operand, finite number + * exp Return the exponent part + * return The normalized fraction in [0.5:1[ + */ +static inline float lc3_frexpf(float _x, int *exp) { + union { float f; uint32_t u; } x = { .f = _x }; + + int e = (x.u & LC3_IEEE754_EXP_MASK) >> LC3_IEEE754_EXP_SHL; + *exp = e - (LC3_IEEE754_EXP_BIAS - 1); + + x.u = (x.u & ~LC3_IEEE754_EXP_MASK) | + ((LC3_IEEE754_EXP_BIAS - 1) << LC3_IEEE754_EXP_SHL); + + return x.f; +} + /** * Fast 2^n approximation * x Operand, range -8 to 8 @@ -69,7 +114,7 @@ static inline float lc3_log2f(float x) static const float c[] = { -1.29479677, 5.11769018, -8.42295281, 8.10557963, -3.50567360 }; - x = frexpf(x, &e); + x = lc3_frexpf(x, &e); y = ( c[0]) * x; y = (y + c[1]) * x; diff --git a/src/lc3.c b/src/lc3.c index 9046abc..ac4427a 100644 --- a/src/lc3.c +++ b/src/lc3.c @@ -222,7 +222,7 @@ static void load_s24( for (int i = 0; i < ns; i++, pcm += stride) { xt[i] = *pcm >> 8; - xs[i] = ldexpf(*pcm, -8); + xs[i] = lc3_ldexpf(*pcm, -8); } } @@ -249,7 +249,7 @@ static void load_s24_3le( ((uint32_t)pcm[2] << 24) ; xt[i] = in >> 16; - xs[i] = ldexpf(in, -16); + xs[i] = lc3_ldexpf(in, -16); } } @@ -271,7 +271,7 @@ static void load_float( int ns = lc3_ns(dt, sr); for (int i = 0; i < ns; i++, pcm += stride) { - xs[i] = ldexpf(*pcm, 15); + xs[i] = lc3_ldexpf(*pcm, 15); xt[i] = LC3_SAT16((int32_t)xs[i]); } } @@ -500,8 +500,8 @@ static void store_s24( int ns = lc3_ns(dt, sr); for ( ; ns > 0; ns--, xs++, pcm += stride) { - int32_t s = *xs >= 0 ? (int32_t)(ldexpf(*xs, 8) + 0.5f) - : (int32_t)(ldexpf(*xs, 8) - 0.5f); + int32_t s = *xs >= 0 ? (int32_t)(lc3_ldexpf(*xs, 8) + 0.5f) + : (int32_t)(lc3_ldexpf(*xs, 8) - 0.5f); *pcm = LC3_SAT24(s); } } @@ -523,8 +523,8 @@ static void store_s24_3le( int ns = lc3_ns(dt, sr); for ( ; ns > 0; ns--, xs++, pcm += 3*stride) { - int32_t s = *xs >= 0 ? (int32_t)(ldexpf(*xs, 8) + 0.5f) - : (int32_t)(ldexpf(*xs, 8) - 0.5f); + int32_t s = *xs >= 0 ? (int32_t)(lc3_ldexpf(*xs, 8) + 0.5f) + : (int32_t)(lc3_ldexpf(*xs, 8) - 0.5f); s = LC3_SAT24(s); pcm[0] = (s >> 0) & 0xff; @@ -550,7 +550,7 @@ static void store_float( int ns = lc3_ns(dt, sr); for ( ; ns > 0; ns--, xs++, pcm += stride) { - float s = ldexpf(*xs, -15); + float s = lc3_ldexpf(*xs, -15); *pcm = fminf(fmaxf(s, -1.f), 1.f); } } diff --git a/src/spec.c b/src/spec.c index bc4d5d4..25a28d4 100644 --- a/src/spec.c +++ b/src/spec.c @@ -135,7 +135,7 @@ LC3_HOT static int estimate_gain( float x_max = sqrtf(x2_max); float nf = lc3_hr(sr) ? - ldexpf(x_max, -reg_bits) * lc3_exp2f(-low_bits) : 0; + lc3_ldexpf(x_max, -reg_bits) * lc3_exp2f(-low_bits) : 0; for (int i = 0; i < n4; i++) e[i].q16 = lc3_db_q16(fmaxf(e[i].f + nf, 1e-10f));