From f39eea6b01d43bc8b04399cc3f54b36f3c8483a6 Mon Sep 17 00:00:00 2001 From: Uladzimir Karpenka Date: Fri, 10 Jul 2026 00:31:32 +0300 Subject: [PATCH] phrot, cfcomp: cascade in a register; drop idiv from cfcomp's ring walks xphrot() kept x0[] and y0[] in the struct, but neither carried state between samples: x0[n] was always the previous stage's output and y0[n] this stage's. Cascade that single value in a register and leave x1/y1 as the real filter state. Stores through out[] may alias the struct's doubles, so hoist the coefficients as well. Eight first-order sections then run without touching memory except for the state. xcfcomp() is structurally the same overlap-add loop as xemnr(), and had the same defect: four ring indices advanced with '% size' per step, ~5100 integer divisions per call at fsize = 2048, which a profile showed dominating the block (1820 samples in xcfcomp against 362 in calc_mask and ~520 in the FFTs). The indices step by one, and iasize >= fsize and oasize >= incr always hold, so they wrap at most once per loop: walk contiguous runs and wrap between them. Measured in situ on an Apple M1 Pro, 512-sample buffers, cost of turning the block on, best of 5: phrot 17219 ns -> 7047 ns 2.44x cfcomp 28918 ns -> 13637 ns 2.12x Both are bit-identical: phrot preserves the operation order, and cfcomp only changes integer index arithmetic. The RX chain is unchanged bit-for-bit. Co-Authored-By: Claude Opus 4.8 --- cfcomp.c | 77 ++++++++++++++++++++++++++++++++++++++------------------ iir.c | 32 +++++++++++++++-------- 2 files changed, 75 insertions(+), 34 deletions(-) diff --git a/cfcomp.c b/cfcomp.c index 37b5e38..95772f8 100644 --- a/cfcomp.c +++ b/cfcomp.c @@ -336,18 +336,36 @@ void xcfcomp (CFCOMP a, int pos) if (a->run && pos == a->position) { int i, j, k, sbuff, sbegin; - for (i = 0; i < 2 * a->bsize; i += 2) + /* Each ring index below steps by one and, since iasize >= fsize and + oasize >= incr always hold, wraps at most once per loop. The '% size' + per step was therefore an integer division for nothing -- about 5100 + of them per call at fsize = 2048. Walk contiguous runs instead. */ + const int iasize = a->iasize; + const int oasize = a->oasize; + const int fsize = a->fsize; + const int incr = a->incr; + const int bsize = a->bsize; + const int ovrlp = a->ovrlp; + const double pregain = a->pregain; + const double postgain = a->postgain; + + for (i = 0, j = a->iainidx; i < 2 * bsize; i += 2) { - a->inaccum[a->iainidx] = a->in[i]; - a->iainidx = (a->iainidx + 1) % a->iasize; + a->inaccum[j] = a->in[i]; + if (++j == iasize) j = 0; } - a->nsamps += a->bsize; - while (a->nsamps >= a->fsize) + a->iainidx = j; + a->nsamps += bsize; + while (a->nsamps >= fsize) { - for (i = 0, j = a->iaoutidx; i < a->fsize; i++, j = (j + 1) % a->iasize) - a->forfftin[i] = a->pregain * a->window[i] * a->inaccum[j]; - a->iaoutidx = (a->iaoutidx + a->incr) % a->iasize; - a->nsamps -= a->incr; + int n1 = iasize - a->iaoutidx; + if (n1 > fsize) n1 = fsize; + for (i = 0; i < n1; i++) + a->forfftin[i] = pregain * a->window[i] * a->inaccum[a->iaoutidx + i]; + for (; i < fsize; i++) + a->forfftin[i] = pregain * a->window[i] * a->inaccum[i - n1]; + if ((a->iaoutidx += incr) >= iasize) a->iaoutidx -= iasize; + a->nsamps -= incr; fftw_execute (a->Rfor); calc_mask(a); for (i = 0; i < a->msize; i++) @@ -356,29 +374,40 @@ void xcfcomp (CFCOMP a, int pos) a->revfftin[2 * i + 1] = a->mask[i] * a->forfftout[2 * i + 1]; } fftw_execute (a->Rrev); - for (i = 0; i < a->fsize; i++) - a->save[a->saveidx][i] = a->postgain * a->window[i] * a->revfftout[i]; - for (i = a->ovrlp; i > 0; i--) + for (i = 0; i < fsize; i++) + a->save[a->saveidx][i] = postgain * a->window[i] * a->revfftout[i]; + for (i = ovrlp; i > 0; i--) { - sbuff = (a->saveidx + i) % a->ovrlp; - sbegin = a->incr * (a->ovrlp - i); - for (j = sbegin, k = a->oainidx; j < a->incr + sbegin; j++, k = (k + 1) % a->oasize) + const double* WDSP_RESTRICT sv; + double* WDSP_RESTRICT oa = a->outaccum; + int m1; + sbuff = (a->saveidx + i) % ovrlp; + sbegin = incr * (ovrlp - i); + sv = a->save[sbuff] + sbegin; + m1 = oasize - a->oainidx; + if (m1 > incr) m1 = incr; + k = a->oainidx; + if (i == ovrlp) { - if ( i == a->ovrlp) - a->outaccum[k] = a->save[sbuff][j]; - else - a->outaccum[k] += a->save[sbuff][j]; + for (j = 0; j < m1; j++) oa[k + j] = sv[j]; + for (; j < incr; j++) oa[j - m1] = sv[j]; + } + else + { + for (j = 0; j < m1; j++) oa[k + j] += sv[j]; + for (; j < incr; j++) oa[j - m1] += sv[j]; } } - a->saveidx = (a->saveidx + 1) % a->ovrlp; - a->oainidx = (a->oainidx + a->incr) % a->oasize; + if (++a->saveidx == ovrlp) a->saveidx = 0; + if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize; } - for (i = 0; i < a->bsize; i++) + for (i = 0, k = a->oaoutidx; i < bsize; i++) { - a->out[2 * i + 0] = a->outaccum[a->oaoutidx]; + a->out[2 * i + 0] = a->outaccum[k]; a->out[2 * i + 1] = 0.0; - a->oaoutidx = (a->oaoutidx + 1) % a->oasize; + if (++k == oasize) k = 0; } + a->oaoutidx = k; } else if (a->out != a->in) memcpy (a->out, a->in, a->bsize * sizeof (complex)); diff --git a/iir.c b/iir.c index 1df7bf3..d0081d9 100644 --- a/iir.c +++ b/iir.c @@ -616,19 +616,31 @@ void xphrot (PHROT a) if (a->run) { int i, n; - for (i = 0; i < a->size; i++) + const int size = a->size; + const int nstages = a->nstages; + const double b0 = a->b0, b1 = a->b1, a1 = a->a1; + /* in and out are the same buffer in TXA, so neither may be restrict */ + const double* in = a->in; + double* out = a->out; + /* x0[]/y0[] never carried state between samples: x0[n] was only ever the + previous stage's output and y0[n] this stage's. Keep that single value + in a register and cascade it, leaving x1/y1 as the actual filter state. + Stores through out[] could alias the struct's doubles, so hoist the + coefficients too. */ + double* WDSP_RESTRICT x1 = a->x1; + double* WDSP_RESTRICT y1 = a->y1; + + for (i = 0; i < size; i++) { - a->x0[0] = a->in[2 * i + 0]; - for (n = 0; n < a->nstages; n++) + double v = in[2 * i + 0]; + for (n = 0; n < nstages; n++) { - if (n > 0) a->x0[n] = a->y0[n - 1]; - a->y0[n] = a->b0 * a->x0[n] - + a->b1 * a->x1[n] - - a->a1 * a->y1[n]; - a->y1[n] = a->y0[n]; - a->x1[n] = a->x0[n]; + double y = b0 * v + b1 * x1[n] - a1 * y1[n]; + x1[n] = v; + y1[n] = y; + v = y; } - a->out[2 * i + 0] = a->y0[a->nstages - 1]; + out[2 * i + 0] = v; } } else if (a->out != a->in)