From 34666474b8ed77ac88de81b23e6bcfc92209bc6d Mon Sep 17 00:00:00 2001 From: Michael Giacomelli Date: Mon, 5 Oct 2026 21:21:27 -0400 Subject: [PATCH] wma: fix dropouts in files with LSP exponents Where the LSP curve peaks, the value its -1/4 power is taken of is far below one step of 16.16 and was rounded to zero. pow_m1_4() of zero is enormous, the block's levels are taken relative to that maximum, and so the whole block came out about 46 dB down. In StutteringFile.wma (22 kHz mono, 20 kbps) this muted a block here and there, which is heard as a stutter. Pass the value to pow_m1_4() with all of its 38 fractional bits. Checked with perfsim (Sansa e200v1 build) against ffmpeg's decode of that file, in blocks of 46 ms: 11 of 841 were off by more than 1 dB, the worst by 10.6 dB; none are now, the worst by 0.4 dB. Four other files with LSP exponents are unchanged or closer, and files without them are byte-identical. Co-Authored-By: Claude Opus 5.5 --- lib/rbcodec/codecs/libwma/wmadeci.c | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/lib/rbcodec/codecs/libwma/wmadeci.c b/lib/rbcodec/codecs/libwma/wmadeci.c index be304f0b53..fd852aa3b7 100644 --- a/lib/rbcodec/codecs/libwma/wmadeci.c +++ b/lib/rbcodec/codecs/libwma/wmadeci.c @@ -587,7 +587,8 @@ int wma_decode_init(WMADecodeContext* s, asf_waveformatex_t *wfx) interpolation to reduce the mantissa table size at a small speed expense (linear interpolation approximately doubles the number of bits of precision). */ -static inline fixed32 pow_m1_4(WMADecodeContext *s, fixed64 x) +/* x^(-1/4) in 16.16, for x given with frac_bits fractional bits */ +static inline fixed32 pow_m1_4(WMADecodeContext *s, fixed64 x, int frac_bits) { union { float f; @@ -597,16 +598,16 @@ static inline fixed32 pow_m1_4(WMADecodeContext *s, fixed64 x) fixed32 a, b; int shift = 0; - /* x can be above what 16.16 holds; bring it into 32 bits to convert - it and put the shift back into the float's exponent */ + /* bring x into 32 bits to convert it, and put that shift and the + fractional bits into the float's exponent */ while (x >> 31) { x >>= 1; shift++; } - u.f = fixtof64((fixed32)x); + u.f = (float)(fixed32)x; if (x != 0) - u.v += shift << 23; + u.v += (shift - frac_bits) << 23; e = u.v >> 23; m = (u.v >> (23 - LSP_POW_BITS)) & ((1 << LSP_POW_BITS) - 1); /* build interpolation scale: 1 <= t < 2. */ @@ -689,10 +690,11 @@ static void wma_lsp_to_curve(WMADecodeContext *s, /* 2 in 5.27 format is 0x10000000. * These squares are far above what 16.16 holds where the curve - * is low, so the sum is taken in 64 bits. It is in 16.16. */ - v64 = ((fixed64)p * fixmul32b(p, (0x10000000 - w)) + - (fixed64)q * fixmul32b(q, (0x10000000 + w))) >> 22; - v = pow_m1_4(s, v64); + * is low, and far below one of its steps where the curve peaks, + * so the sum is taken in 64 bits, with 38 fractional bits. */ + v64 = (fixed64)p * fixmul32b(p, (0x10000000 - w)) + + (fixed64)q * fixmul32b(q, (0x10000000 + w)); + v = pow_m1_4(s, v64, 38); if (v > val_max) val_max = v; out[i] = v;