From 244d5a1cc1460d36dbd4c9b76d6ae49a5a4de0fb Mon Sep 17 00:00:00 2001 From: Pauli Virtanen Date: Sat, 26 Jul 2025 17:23:44 +0300 Subject: [PATCH] resample: use fixed point for resample phase and input rate If phase is float, calculations in impl_native_in_len/out_len can produce wrong results due to rounding error. It's probably better to not be in the business of predicting floating-point rounding, so replace this by fixed-point arithmetic. Also make sure `offset+1` cannot overflow data->filter array in do_resample_inter* due to float multiplication possibly rounding up. --- .../audioconvert/resample-native-impl.h | 45 ++++++++++++------- spa/plugins/audioconvert/resample-native.c | 40 +++++++++++------ spa/plugins/audioconvert/test-resample.c | 8 ++-- 3 files changed, 61 insertions(+), 32 deletions(-) diff --git a/spa/plugins/audioconvert/resample-native-impl.h b/spa/plugins/audioconvert/resample-native-impl.h index 023b350d7..534f36da6 100644 --- a/spa/plugins/audioconvert/resample-native-impl.h +++ b/spa/plugins/audioconvert/resample-native-impl.h @@ -12,6 +12,18 @@ typedef void (*resample_func_t)(struct resample *r, const void * SPA_RESTRICT src[], uint32_t ioffs, uint32_t *in_len, void * SPA_RESTRICT dst[], uint32_t ooffs, uint32_t *out_len); +#define FIXP_SHIFT 32 +#define FIXP_SCALE ((uint64_t)1 << FIXP_SHIFT) +#define FIXP_MASK (FIXP_SCALE - 1) +#define UINT32_TO_FIXP(v) ((struct fixp) { (uint64_t)((uint32_t)(v)) << FIXP_SHIFT }) +#define FLOAT_TO_FIXP(d) ((struct fixp) { (uint64_t)((d) * (float)FIXP_SCALE) }) +#define FIXP_TO_UINT32(f) ((f).value >> FIXP_SHIFT) +#define FIXP_TO_FLOAT(f) ((f).value / (float)FIXP_SCALE) + +struct fixp { + uint64_t value; +}; + struct resample_info { uint32_t format; resample_func_t process_copy; @@ -29,10 +41,10 @@ struct native_data { uint32_t n_phases; uint32_t in_rate; uint32_t out_rate; - float phase; + struct fixp phase; float pm; uint32_t inc; - uint32_t frac; + struct fixp frac; uint32_t filter_stride; uint32_t filter_stride_os; uint32_t gcd; @@ -86,25 +98,26 @@ DEFINE_RESAMPLER(full,arch) \ { \ struct native_data *data = r->data; \ uint32_t n_taps = data->n_taps, stride = data->filter_stride_os; \ - uint32_t index, phase, out_rate = data->out_rate; \ + uint32_t index; \ uint32_t c, o, olen = *out_len, ilen = *in_len; \ - uint32_t inc = data->inc, frac = data->frac, ch = r->channels; \ + uint32_t inc = data->inc, ch = r->channels; \ + uint64_t frac = data->frac.value, phase = data->phase.value; \ + uint64_t denom = UINT32_TO_FIXP(data->out_rate).value; \ \ index = ioffs; \ - phase = (uint32_t)data->phase; \ for (o = ooffs; o < olen && index + n_taps <= ilen; o++) { \ - float *filter = &data->filter[phase * stride]; \ + float *filter = &data->filter[(phase >> FIXP_SHIFT) * stride]; \ for (c = 0; c < ch; c++) { \ const float *s = src[c]; \ float *d = dst[c]; \ inner_product_##arch(&d[o], &s[index], \ filter, n_taps); \ } \ - INC(index, phase, out_rate); \ + INC(index, phase, denom); \ } \ *in_len = index; \ *out_len = o; \ - data->phase = phase; \ + data->phase.value = phase; \ } #define MAKE_RESAMPLER_INTER(arch) \ @@ -112,16 +125,18 @@ DEFINE_RESAMPLER(inter,arch) \ { \ struct native_data *data = r->data; \ uint32_t index, stride = data->filter_stride; \ - uint32_t n_taps = data->n_taps, out_rate = data->out_rate; \ + uint32_t n_taps = data->n_taps; \ uint32_t c, o, olen = *out_len, ilen = *in_len; \ - uint32_t inc = data->inc, frac = data->frac, ch = r->channels; \ - float phase, pm = data->pm; \ + uint32_t inc = data->inc, ch = r->channels; \ + uint32_t ph_max = data->n_phases - 1; \ + uint64_t frac = data->frac.value, phase = data->phase.value; \ + uint64_t denom = UINT32_TO_FIXP(data->out_rate).value; \ + float pm = data->pm; \ \ index = ioffs; \ - phase = data->phase; \ for (o = ooffs; o < olen && index + n_taps <= ilen; o++) { \ float ph = phase * pm; \ - uint32_t offset = (uint32_t)floorf(ph); \ + uint32_t offset = SPA_MIN((uint32_t)floorf(ph), ph_max); \ float *filter0 = &data->filter[(offset+0) * stride]; \ float *filter1 = &data->filter[(offset+1) * stride]; \ float pho = ph - offset; \ @@ -131,11 +146,11 @@ DEFINE_RESAMPLER(inter,arch) \ inner_product_ip_##arch(&d[o], &s[index], \ filter0, filter1, pho, n_taps); \ } \ - INC(index, phase, out_rate); \ + INC(index, phase, denom); \ } \ *in_len = index; \ *out_len = o; \ - data->phase = phase; \ + data->phase.value = phase; \ } diff --git a/spa/plugins/audioconvert/resample-native.c b/spa/plugins/audioconvert/resample-native.c index aab10efb9..d544b4cea 100644 --- a/spa/plugins/audioconvert/resample-native.c +++ b/spa/plugins/audioconvert/resample-native.c @@ -159,17 +159,31 @@ static void impl_native_update_rate(struct resample *r, double rate) data->func = data->info->process_full; } - data->in_rate = in_rate; if (data->out_rate != out_rate) { - data->phase = data->phase * out_rate / (float)data->out_rate; - data->out_rate = out_rate; + /* Cast to double to avoid overflows */ + data->phase.value = (uint64_t)(data->phase.value * (double)out_rate / data->out_rate); + if (data->phase.value >= UINT32_TO_FIXP(out_rate).value) + data->phase.value = UINT32_TO_FIXP(out_rate).value - 1; } - data->inc = data->in_rate / data->out_rate; - data->frac = data->in_rate % data->out_rate; - spa_log_trace_fp(r->log, "native %p: rate:%f in:%d out:%d phase:%f inc:%d frac:%d", r, - rate, r->i_rate, r->o_rate, data->phase, data->inc, data->frac); + data->in_rate = in_rate; + data->out_rate = out_rate; + data->inc = in_rate / out_rate; + data->frac = UINT32_TO_FIXP(in_rate % out_rate); + + spa_log_trace_fp(r->log, "native %p: rate:%f in:%d out:%d phase:%f inc:%d frac:%f", r, + rate, r->i_rate, r->o_rate, FIXP_TO_FLOAT(data->phase), + data->inc, FIXP_TO_FLOAT(data->frac)); +} + +static uint64_t fixp_floor_a_plus_bc(struct fixp a, uint32_t b, struct fixp c) +{ + /* (a + b*c) >> FIXP_SHIFT, with bigger overflow threshold */ + uint64_t hi, lo; + hi = (a.value >> FIXP_SHIFT) + b * (c.value >> FIXP_SHIFT); + lo = (a.value & FIXP_MASK) + b * (c.value & FIXP_MASK); + return hi + (lo >> FIXP_SHIFT); } static uint32_t impl_native_in_len(struct resample *r, uint32_t out_len) @@ -177,7 +191,7 @@ static uint32_t impl_native_in_len(struct resample *r, uint32_t out_len) struct native_data *data = r->data; uint32_t in_len; - in_len = (uint32_t)((data->phase + out_len * data->frac) / data->out_rate); + in_len = fixp_floor_a_plus_bc(data->phase, out_len, data->frac) / data->out_rate; in_len += out_len * data->inc + (data->n_taps - data->hist); spa_log_trace_fp(r->log, "native %p: hist:%d %d->%d", r, data->hist, out_len, in_len); @@ -191,7 +205,7 @@ static uint32_t impl_native_out_len(struct resample *r, uint32_t in_len) uint32_t out_len; in_len = in_len - SPA_MIN(in_len, data->n_taps - data->hist); - out_len = (uint32_t)(in_len * data->out_rate - data->phase); + out_len = in_len * data->out_rate - FIXP_TO_UINT32(data->phase); out_len = (out_len + data->in_rate - 1) / data->in_rate; spa_log_trace_fp(r->log, "native %p: hist:%d %d->%d", r, data->hist, in_len, out_len); @@ -300,7 +314,7 @@ static void impl_native_reset (struct resample *r) d->hist = d->n_taps - 1; else d->hist = d->n_taps / 2; - d->phase = 0; + d->phase.value = 0; } static uint32_t impl_native_delay (struct resample *r) @@ -315,13 +329,13 @@ static float impl_native_phase (struct resample *r) float pho = 0; if (d->func == d->info->process_full) { - pho = -(float)((int32_t)d->phase) / d->out_rate; + pho = -(float)FIXP_TO_UINT32(d->phase) / d->out_rate; /* XXX: this is how it seems to behave, but not clear why */ if (d->hist >= d->n_taps - 1) pho += 1.0f; } else if (d->func == d->info->process_inter) { - pho = -d->phase / d->out_rate; + pho = -FIXP_TO_FLOAT(d->phase) / d->out_rate; /* XXX: this is how it seems to behave, but not clear why */ if (d->hist >= d->n_taps - 1) @@ -388,7 +402,7 @@ int resample_native_init(struct resample *r) d->in_rate = in_rate; d->out_rate = out_rate; d->gcd = gcd; - d->pm = (float)n_phases / r->o_rate; + d->pm = (float)n_phases / r->o_rate / FIXP_SCALE; d->filter = SPA_PTROFF_ALIGN(d, sizeof(struct native_data), 64, float); d->hist_mem = SPA_PTROFF_ALIGN(d->filter, filter_size, 64, float); d->history = SPA_PTROFF(d->hist_mem, history_size, float*); diff --git a/spa/plugins/audioconvert/test-resample.c b/spa/plugins/audioconvert/test-resample.c index 0848d8f80..8c0a7de86 100644 --- a/spa/plugins/audioconvert/test-resample.c +++ b/spa/plugins/audioconvert/test-resample.c @@ -126,20 +126,20 @@ static void pull_blocks_out(struct resample *r, uint32_t first, uint32_t size, u } } -static void check_inout_len(struct resample *r, uint32_t first, uint32_t size, double rate, double phase) +static void check_inout_len(struct resample *r, uint32_t first, uint32_t size, double rate, float phase) { struct native_data *data = r->data; resample_reset(r); resample_update_rate(r, rate); if (phase != 0.0) - data->phase = (float)phase; + data->phase = FLOAT_TO_FIXP(phase); pull_blocks(r, first, size, 500); resample_reset(r); resample_update_rate(r, rate); if (phase != 0.0) - data->phase = (float)phase; + data->phase = FLOAT_TO_FIXP(phase); pull_blocks_out(r, first, size, 500); } @@ -237,7 +237,7 @@ static void test_inout_len(void) r.options = RESAMPLE_OPTION_PREFILL; resample_native_init(&r); - check_inout_len(&r, 64, 64, 1.0 + 1e-10, 7999.99); + check_inout_len(&r, 64, 64, 1.0 + 1e-10, 7999.99f); resample_free(&r); /* Test value of phase that overflows filter buffer due to floating point rounding