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.
This commit is contained in:
Pauli Virtanen 2025-07-26 17:23:44 +03:00 committed by Wim Taymans
parent 3cade43cf3
commit 244d5a1cc1
3 changed files with 61 additions and 32 deletions

View file

@ -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; \
}

View file

@ -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*);

View file

@ -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