From 9960a1d0e3050a9c88acdae7d709f021f12dd226 Mon Sep 17 00:00:00 2001 From: Uladzimir Karpenka Date: Thu, 9 Jul 2026 22:44:08 +0300 Subject: [PATCH] resample: vectorize the polyphase tap loop The tap loop wrapped the ring with a test on every tap: if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize; which made the address non-affine and stopped the vectorizer. Split the walk at the wrap point instead, so both halves are unit-stride, and split the complex ring into separate I/Q arrays so the taps load contiguously rather than through a de-interleaving ld2. The dot product also could not be vectorized as written: reassociating an fp reduction needs -ffast-math, which this library must not enable (it relies on IEEE semantics for 0/0 = NaN and x/0 = Inf, see linux_port.h). Carry four independent accumulator pairs instead, which both breaks the FMA dependency chain and lets the vectorizer in on any compiler. Measured on an Apple M1 Pro, 512-sample DSP buffers, best of 5: xresample, decimation to 48 kHz before after speedup 192k -> 48k (561 taps) 342.0 us 84.2 us 4.06x 384k -> 48k (1121 taps) 708.9 us 175.1 us 4.05x 576k -> 48k (1681 taps) 1074.4 us 266.9 us 4.03x 768k -> 48k (2241 taps) 1438.9 us 360.0 us 4.00x full xrxa() chain, 576k input 1148.2 us 337.6 us 3.40x 10.76% 3.16% of one core xresampleF (float, host audio) 2.6x - 3.4x Summation order changes, so the double path is not bit-identical: over 400 buffers of the full RX chain the worst deviation is 1.1e-12, an SNR of 251 dB. The float path is bit-identical, as the cast to float absorbs the difference. With the resampler bypassed (48k in, 48k out) the chain is unchanged bit-for-bit. Co-Authored-By: Claude Opus 4.8 --- comm.h | 8 +++ resample.c | 158 +++++++++++++++++++++++++++++++++++++++-------------- resample.h | 4 +- 3 files changed, 128 insertions(+), 42 deletions(-) diff --git a/comm.h b/comm.h index 430f2c7..ea1f729 100644 --- a/comm.h +++ b/comm.h @@ -145,6 +145,14 @@ warren@wpratt.com #define PI 3.1415926535897932 #define TWOPI 6.2831853071795864 +// Non-aliasing qualifier for DSP buffers. Spelled __restrict rather than +// restrict because the JNI translation unit is compiled as -std=gnu89. +#if defined(__GNUC__) || defined(__clang__) || defined(_MSC_VER) +#define WDSP_RESTRICT __restrict +#else +#define WDSP_RESTRICT +#endif + // miscellaneous typedef double complex[2]; #define PORT __declspec( dllexport ) diff --git a/resample.c b/resample.c index 76dd500..5d46926 100644 --- a/resample.c +++ b/resample.c @@ -32,6 +32,37 @@ warren@wpratt.com * * ************************************************************************************************/ +/* Accumulate n taps of a unit-stride complex dot product into *pI / *pQ. + + Four independent accumulator pairs are carried so the FMAs are not serialized + on a single dependency chain, and so the compiler is free to vectorize: a + plain 'I += h[j]*x[j]' reduction cannot be reassociated without -ffast-math, + which this library must not enable (it relies on IEEE semantics for 0/0 = NaN + and x/0 = Inf). Summation order therefore differs from a strict left-to-right + reduction, at the usual pairwise-summation accuracy gain. */ +static inline void resample_dot (const double* WDSP_RESTRICT hp, + const double* WDSP_RESTRICT xI, const double* WDSP_RESTRICT xQ, + int n, double* pI, double* pQ) +{ + double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0; + double q0 = 0.0, q1 = 0.0, q2 = 0.0, q3 = 0.0; + int j = 0; + for (; j <= n - 4; j += 4) + { + i0 += hp[j + 0] * xI[j + 0]; q0 += hp[j + 0] * xQ[j + 0]; + i1 += hp[j + 1] * xI[j + 1]; q1 += hp[j + 1] * xQ[j + 1]; + i2 += hp[j + 2] * xI[j + 2]; q2 += hp[j + 2] * xQ[j + 2]; + i3 += hp[j + 3] * xI[j + 3]; q3 += hp[j + 3] * xQ[j + 3]; + } + for (; j < n; j++) + { + i0 += hp[j] * xI[j]; + q0 += hp[j] * xQ[j]; + } + *pI += (i0 + i1) + (i2 + i3); + *pQ += (q0 + q1) + (q2 + q3); +} + void calc_resample (RESAMPLE a) { int x, y, z; @@ -71,7 +102,8 @@ void calc_resample (RESAMPLE a) for (k = 0; k < a->ncoef; k += a->L) a->h[i++] = impulse[j + k]; a->ringsize = a->cpp; - a->ring = (double *)malloc0(a->ringsize * sizeof(complex)); + a->ringI = (double *)malloc0(a->ringsize * sizeof(double)); + a->ringQ = (double *)malloc0(a->ringsize * sizeof(double)); a->idx_in = a->ringsize - 1; a->phnum = 0; _aligned_free(impulse); @@ -79,7 +111,8 @@ void calc_resample (RESAMPLE a) void decalc_resample (RESAMPLE a) { - _aligned_free(a->ring); + _aligned_free(a->ringQ); + _aligned_free(a->ringI); _aligned_free(a->h); } @@ -112,7 +145,8 @@ void destroy_resample (RESAMPLE a) PORT void flush_resample (RESAMPLE a) { - memset (a->ring, 0, a->ringsize * sizeof (complex)); + memset (a->ringI, 0, a->ringsize * sizeof (double)); + memset (a->ringQ, 0, a->ringsize * sizeof (double)); a->idx_in = a->ringsize - 1; a->phnum = 0; } @@ -123,40 +157,49 @@ int xresample (RESAMPLE a) int outsamps = 0; if (a->run) { - int i, j, n; - int idx_out; + int i, n1; double I, Q; - int cpp = a->cpp; + const int cpp = a->cpp; + const int ringsize = a->ringsize; + const int L = a->L; + const int M = a->M; + const int size = a->size; + const double* WDSP_RESTRICT h = a->h; + const double* WDSP_RESTRICT in = a->in; + double* WDSP_RESTRICT ringI = a->ringI; + double* WDSP_RESTRICT ringQ = a->ringQ; + double* WDSP_RESTRICT out = a->out; int idx_in = a->idx_in; - int ringsize = a->ringsize; - double* h = a->h; - double* ring = a->ring; + int phnum = a->phnum; - for (i = 0; i < a->size; i++) + for (i = 0; i < size; i++) { - ring[2 * idx_in + 0] = a->in[2 * i + 0]; - ring[2 * idx_in + 1] = a->in[2 * i + 1]; - while (a->phnum < a->L) + ringI[idx_in] = in[2 * i + 0]; + ringQ[idx_in] = in[2 * i + 1]; + while (phnum < L) { + const double* WDSP_RESTRICT hp = h + cpp * phnum; + /* The tap loop walks the ring forward from idx_in and wraps at + most once. Split it at the wrap point so both halves are + unit-stride: the wrap test that used to sit inside the loop + made the address non-affine and blocked vectorization. */ + if ((n1 = ringsize - idx_in) > cpp) n1 = cpp; I = 0.0; Q = 0.0; - n = cpp * a->phnum; - for (j = 0; j < cpp; j++) - { - if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize; - I += h[n + j] * ring[2 * idx_out + 0]; - Q += h[n + j] * ring[2 * idx_out + 1]; - } - a->out[2 * outsamps + 0] = I; - a->out[2 * outsamps + 1] = Q; + resample_dot (hp, ringI + idx_in, ringQ + idx_in, n1, &I, &Q); + if (n1 < cpp) + resample_dot (hp + n1, ringI, ringQ, cpp - n1, &I, &Q); + out[2 * outsamps + 0] = I; + out[2 * outsamps + 1] = Q; outsamps++; - a->phnum += a->M; + phnum += M; } - a->phnum -= a->L; - if (--idx_in < 0) idx_in = a->ringsize - 1; + phnum -= L; + if (--idx_in < 0) idx_in = ringsize - 1; } a->idx_in = idx_in; + a->phnum = phnum; } else if (a->in != a->out) memcpy (a->out, a->in, a->size * sizeof (complex)); @@ -240,6 +283,24 @@ void destroy_resampleV (void* ptr) * * ************************************************************************************************/ +/* Real-valued counterpart of resample_dot(). */ +static inline void resampleF_dot (const double* WDSP_RESTRICT hp, + const double* WDSP_RESTRICT x, int n, double* pI) +{ + double i0 = 0.0, i1 = 0.0, i2 = 0.0, i3 = 0.0; + int j = 0; + for (; j <= n - 4; j += 4) + { + i0 += hp[j + 0] * x[j + 0]; + i1 += hp[j + 1] * x[j + 1]; + i2 += hp[j + 2] * x[j + 2]; + i3 += hp[j + 3] * x[j + 3]; + } + for (; j < n; j++) + i0 += hp[j] * x[j]; + *pI += (i0 + i1) + (i2 + i3); +} + RESAMPLEF create_resampleF ( int run, int size, float* in, float* out, int in_rate, int out_rate) { RESAMPLEF a = (RESAMPLEF) malloc0 (sizeof (resampleF)); @@ -305,31 +366,46 @@ int xresampleF (RESAMPLEF a) int outsamps = 0; if (a->run) { - int i, j, n; - int idx_out; + int i; double I; - for (i = 0; i < a->size; i++) - { - a->ring[a->idx_in] = (double)a->in[i]; + const int cpp = a->cpp; + const int ringsize = a->ringsize; + const int L = a->L; + const int M = a->M; + const int size = a->size; + const double* WDSP_RESTRICT h = a->h; + const float* WDSP_RESTRICT in = a->in; + double* WDSP_RESTRICT ring = a->ring; + float* WDSP_RESTRICT out = a->out; + int idx_in = a->idx_in; + int phnum = a->phnum; + int n1; - while (a->phnum < a->L) + for (i = 0; i < size; i++) + { + ring[idx_in] = (double)in[i]; + + while (phnum < L) { + const double* WDSP_RESTRICT hp = h + cpp * phnum; + /* see resample_dot(): split at the ring wrap so both halves are + unit-stride, and carry independent accumulators */ + if ((n1 = ringsize - idx_in) > cpp) n1 = cpp; I = 0.0; - n = a->cpp * a->phnum; - for (j = 0; j < a->cpp; j++) - { - if ((idx_out = a->idx_in + j) >= a->ringsize) idx_out -= a->ringsize; - I += a->h[n + j] * a->ring[idx_out]; - } - a->out[outsamps] = (float)I; + resampleF_dot (hp, ring + idx_in, n1, &I); + if (n1 < cpp) + resampleF_dot (hp + n1, ring, cpp - n1, &I); + out[outsamps] = (float)I; outsamps++; - a->phnum += a->M; + phnum += M; } - a->phnum -= a->L; - if (--a->idx_in < 0) a->idx_in = a->ringsize - 1; + phnum -= L; + if (--idx_in < 0) idx_in = ringsize - 1; } + a->idx_in = idx_in; + a->phnum = phnum; } else if (a->in != a->out) memcpy (a->out, a->in, a->size * sizeof (float)); diff --git a/resample.h b/resample.h index 73c4b0b..aa25637 100644 --- a/resample.h +++ b/resample.h @@ -52,7 +52,9 @@ typedef struct _resample int M; // decimation factor double* h; // coefficients int ringsize; // number of complex pairs the ring buffer holds - double* ring; // ring buffer + double* ringI; // ring buffer, in-phase + double* ringQ; // ring buffer, quadrature (split from I so the tap loop + // reads unit-stride and vectorizes) int cpp; // coefficients of the phase int phnum; // phase number } resample, *RESAMPLE;