varsamp: transpose the coefficients into phases, vectorize the tap loop

hshift() rebuilds the whole interpolated tap set on every output sample,
reading h[hidx + m*R] for m = 0..rsize-1. With R = 1024 (what rmatch asks
for) that strides 8 KB at a time through a 1.1 MB table, so every one of
the 2*rsize reads is its own cache line. It cost more than the filter it
was feeding: 258 ns per output sample against 146 ns for the tap loop.

Store the coefficients transposed instead, hp[p*rsize + m] = h[p + m*R],
so the two phases hshift() interpolates between are each contiguous. R+1
phases are needed since it reads hidx and hidx+1, and h_offset is kept in
[0,1) by the caller so hidx <= R-1. The untransposed h is freed; the
impulse cache hands back a copy, so varsamp owns it. Net memory is
unchanged.

The tap loop had the same wrap test per tap as resample.c did, so split it
at the wrap and carry four independent accumulator pairs; the ring is split
into I/Q so the taps load unit-stride.

Note a->hs is rewritten by hshift() inside the sample loop, so it must not
be hoisted behind a restrict pointer in xvarsamp().

Measured on an Apple M1 Pro, 512-sample buffers, best of 5:

    48k -> 48k   varmode=0    138580 ns -> 37370 ns   3.71x
    48k -> 48k   varmode=1    138390 ns -> 32333 ns   4.28x
    48k -> 44.1k varmode=1    141473 ns -> 38220 ns   3.70x
    44.1k -> 48k varmode=1    152307 ns -> 45263 ns   3.36x

hshift() is numerically identical -- same coefficients, different layout.
Only the reassociated tap sum rounds differently: worst deviation 7.2e-16
over four rate configurations, an SNR of 344 dB. Driven end to end through
rmatch's public API, output SNR is 306 dB and total energy matches.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

varsamp.c

@@ -59,10 +59,23 @@ void calc_varsamp (VARSAMP a)
 		fc_norm_low = a->fc_low / norm_rate;
 	a->rsize = (int)(140.0 * norm_rate / min_rate);
 	a->ncoef = a->rsize + 1;
-	a->ncoef += (a->R - 1) * (a->ncoef - 1);
+	a->ncoef += (a->R - 1) * (a->ncoef - 1);		// = R * rsize + 1
-	a->h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)a->R, 1, 0, (double)a->R * a->gain);
+	{
-	// print_impulse ("imp.txt", a->ncoef, a->h, 0, 0);
+		/* Store the coefficients transposed into phases.  hshift() walks
-	a->ring = (double *)malloc0(a->rsize * sizeof(complex));
+		   h[hidx + m*R] for m = 0..rsize-1, which strides by R doubles -- 8 KB
+		   at R = 1024 -- over a 1.1 MB table, so every tap is its own cache
+		   line.  Transposing makes each phase contiguous; hshift() interpolates
+		   between phases hidx and hidx+1, hence R+1 of them. */
+		int p, m;
+		double* h = fir_bandpass(a->ncoef, fc_norm_low, fc_norm_high, (double)a->R, 1, 0, (double)a->R * a->gain);
+		a->hp = (double *)malloc0 ((size_t)(a->R + 1) * a->rsize * sizeof (double));
+		for (p = 0; p <= a->R; p++)
+			for (m = 0; m < a->rsize; m++)
+				a->hp[(size_t)p * a->rsize + m] = h[p + (size_t)m * a->R];
+		_aligned_free (h);
+	}
+	a->ringI = (double *)malloc0(a->rsize * sizeof(double));
+	a->ringQ = (double *)malloc0(a->rsize * sizeof(double));
 	a->idx_in = a->rsize - 1;
 	a->h_offset = 0.0;
 	a->hs = (double *)malloc0 (a->rsize * sizeof (double));
@@ -72,8 +85,9 @@ void calc_varsamp (VARSAMP a)
 void decalc_varsamp (VARSAMP a)
 {
 	_aligned_free (a->hs);
-	_aligned_free (a->ring);
+	_aligned_free (a->ringQ);
-	_aligned_free (a->h);
+	_aligned_free (a->ringI);
+	_aligned_free (a->hp);
 }
 
 VARSAMP create_varsamp ( int run, int size, double* in, double* out,
@@ -105,22 +119,60 @@ void destroy_varsamp (VARSAMP a)
 
 void flush_varsamp (VARSAMP a)
 {
-	memset (a->ring, 0, a->rsize * sizeof (complex));
+	memset (a->ringI, 0, a->rsize * sizeof (double));
+	memset (a->ringQ, 0, a->rsize * sizeof (double));
 	a->idx_in = a->rsize - 1;
 	a->h_offset = 0.0;
 	a->isamps = 0.0;
 }
 
+/* Accumulate n taps of a unit-stride complex dot product into *pI / *pQ.
+
+   Four independent accumulator pairs keep the FMAs off a single dependency
+   chain and let the vectorizer in: an '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). */
+static inline void varsamp_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 hshift (VARSAMP a)
 {
-	int i, j, k;
+	int m;
 	int hidx;
 	double frac, pos;
+	const int rsize = a->rsize;
+	const double* WDSP_RESTRICT h0;
+	const double* WDSP_RESTRICT h1;
+	double* WDSP_RESTRICT hs = a->hs;
+	/* h_offset is normalized to [0,1) by the caller, so hidx is in [0, R-1]
+	   and phase hidx+1 <= R exists. */
 	pos = (double)a->R * a->h_offset;
 	hidx = (int)(pos);
 	frac = pos - (double)hidx;
-	for (i = a->rsize - 1, j = hidx, k = hidx + 1; i >= 0; i--, j += a->R, k += a->R)
+	h0 = a->hp + (size_t)hidx * rsize;
-		a->hs[i] = a->h[j] + frac * (a->h[k] - a->h[j]);
+	h1 = h0 + rsize;
+	for (m = 0; m < rsize; m++)
+		hs[rsize - 1 - m] = h0[m] + frac * (h1[m] - h0[m]);
 }
 
 int xvarsamp (VARSAMP a, double var)
@@ -140,13 +192,21 @@ int xvarsamp (VARSAMP a, double var)
 	else			a->dicvar = 0.0;
 	if (a->run)
 	{
-		int i, j;
+		int i, n1;
-		int idx_out;
 		double I, Q;
+		const int rsize = a->rsize;
+		/* a->hs is rewritten by hshift() on every output sample, so it must not
+		   be hoisted behind a restrict pointer here; varsamp_dot() re-reads it. */
+		const double* in = a->in;
+		double* out = a->out;
+		double* WDSP_RESTRICT ringI = a->ringI;
+		double* WDSP_RESTRICT ringQ = a->ringQ;
+		int idx_in = a->idx_in;
+
 		for (i = 0; i < a->size; i++)
 		{
-			a->ring[2 * a->idx_in + 0] = a->in[2 * i + 0];
+			ringI[idx_in] = in[2 * i + 0];
-			a->ring[2 * a->idx_in + 1] = a->in[2 * i + 1];
+			ringQ[idx_in] = in[2 * i + 1];
 			a->inv_cvar += a->dicvar;
 			picvar = (uint64_t*)(&a->inv_cvar);
 			N = *picvar & 0xffffffffffff0000;
@@ -160,20 +220,21 @@ int xvarsamp (VARSAMP a, double var)
 				a->h_offset += a->delta;
 				while (a->h_offset >= 1.0) a->h_offset -= 1.0;
 				while (a->h_offset <  0.0) a->h_offset += 1.0;
-				for (j = 0; j < a->rsize; j++)
+				/* the ring wraps at most once over rsize taps; split it so both
-				{
+				   halves are unit-stride */
-					if ((idx_out = a->idx_in + j) >= a->rsize) idx_out -= a->rsize;
+				n1 = rsize - idx_in;
-					I += a->hs[j] * a->ring[2 * idx_out + 0];
+				varsamp_dot (a->hs, ringI + idx_in, ringQ + idx_in, n1, &I, &Q);
-					Q += a->hs[j] * a->ring[2 * idx_out + 1];
+				if (n1 < rsize)
-				}
+					varsamp_dot (a->hs + n1, ringI, ringQ, rsize - n1, &I, &Q);
-				a->out[2 * outsamps + 0] = I;
+				out[2 * outsamps + 0] = I;
-				a->out[2 * outsamps + 1] = Q;
+				out[2 * outsamps + 1] = Q;
 				outsamps++;
 				a->isamps += a->inv_cvar;
 			}
 			a->isamps -= 1.0;
-			if (--a->idx_in < 0) a->idx_in = a->rsize - 1;
+			if (--idx_in < 0) idx_in = rsize - 1;
 		}
+		a->idx_in = idx_in;
 	}
 	else if (a->in != a->out)
 		memcpy (a->out, a->in, a->size * sizeof (complex));

varsamp.h

@@ -41,9 +41,13 @@ typedef struct _varsamp
 	double gain;
 	int idx_in;
 	int ncoef;
-	double* h;
+	double* hp;			// coefficients, polyphase: hp[p * rsize + m] = h[p + m * R],
+						//   p = 0..R.  hshift() then reads two adjacent phases
+						//   contiguously instead of striding by R.
 	int rsize;
-	double* ring;
+	double* ringI;		// ring buffer, in-phase
+	double* ringQ;		// ring buffer, quadrature (split from I so the tap loop
+						//   reads unit-stride and vectorizes)
 	double var;
 	int varmode;
 	double cvar;