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 <noreply@anthropic.com>

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 )

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 i, n1;
-		int idx_out;
 		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;
+		int phnum = a->phnum;
-		double* h = a->h;
-		double* ring = a->ring;
 
-		for (i = 0; i < a->size; i++)
+		for (i = 0; i < size; i++)
 		{
-			ring[2 * idx_in + 0] = a->in[2 * i + 0];
+			ringI[idx_in] = in[2 * i + 0];
-			ring[2 * idx_in + 1] = a->in[2 * i + 1];
+			ringQ[idx_in] = in[2 * i + 1];
-			while (a->phnum < a->L)
+			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;
+				resample_dot (hp, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
-				for (j = 0; j < cpp; j++)
+				if (n1 < cpp)
-				{
+					resample_dot (hp + n1, ringI, ringQ, cpp - n1, &I, &Q);
-					if ((idx_out = idx_in + j) >= ringsize) idx_out -= ringsize;
+				out[2 * outsamps + 0] = I;
-					I += h[n + j] * ring[2 * idx_out + 0];
+				out[2 * outsamps + 1] = Q;
-					Q += h[n + j] * ring[2 * idx_out + 1];
-				}
-				a->out[2 * outsamps + 0] = I;
-				a->out[2 * outsamps + 1] = Q;
 				outsamps++;
-				a->phnum += a->M;
+				phnum += M;
 			}
-			a->phnum -= a->L;
+			phnum -= L;
-			if (--idx_in < 0) idx_in = a->ringsize - 1;
+			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 i;
-		int idx_out;
 		double I;
 
-		for (i = 0; i < a->size; i++)
+		const int cpp = a->cpp;
-		{
+		const int ringsize = a->ringsize;
-			a->ring[a->idx_in] = (double)a->in[i];
+		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;
+				resampleF_dot (hp, ring + idx_in, n1, &I);
-				for (j = 0; j < a->cpp; j++)
+				if (n1 < cpp)
-				{
+					resampleF_dot (hp + n1, ring, cpp - n1, &I);
-					if ((idx_out = a->idx_in + j) >= a->ringsize) idx_out -= a->ringsize;
+				out[outsamps] = (float)I;
-					I += a->h[n + j] * a->ring[idx_out];
-				}
-				a->out[outsamps] = (float)I;
 
 				outsamps++;
-				a->phnum += a->M;
+				phnum += M;
 			}
-			a->phnum -= a->L;
+			phnum -= L;
-			if (--a->idx_in < 0) a->idx_in = a->ringsize - 1;
+			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));

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;