ew8bak/wdsp
1
0
Fork 0
mirror of https://git.vladimir.cc/vladimir/wdsp.git synced 2026-08-25 17:27:33 +00:00
Code Issues Packages Projects Releases Wiki Activity

phrot, cfcomp: cascade in a register; drop idiv from cfcomp's ring walks

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Uladzimir Karpenka
2026-07-10 00:31:32 +03:00
co-authored by Claude Opus 4.8
parent 91df5b1f2d
commit f39eea6b01
2 changed files with 75 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
cfcomp.c
+56 -27
View File
@@ -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,30 +374,41 @@ 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];
for (j = 0; j < m1; j++) oa[k + j] = sv[j];
for (; j < incr; j++) oa[j - m1] = sv[j];
}
else
a->outaccum[k] += a->save[sbuff][j];
}
}
a->saveidx = (a->saveidx + 1) % a->ovrlp;
a->oainidx = (a->oainidx + a->incr) % a->oasize;
}
for (i = 0; i < a->bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[a->oaoutidx];
a->out[2 * i + 1] = 0.0;
a->oaoutidx = (a->oaoutidx + 1) % a->oasize;
for (j = 0; j < m1; j++) oa[k + j] += sv[j];
for (; j < incr; j++) oa[j - m1] += sv[j];
}
}
if (++a->saveidx == ovrlp) a->saveidx = 0;
if ((a->oainidx += incr) >= oasize) a->oainidx -= oasize;
}
for (i = 0, k = a->oaoutidx; i < bsize; i++)
{
a->out[2 * i + 0] = a->outaccum[k];
a->out[2 * i + 1] = 0.0;
if (++k == oasize) k = 0;
}
a->oaoutidx = k;
}
else if (a->out != a->in)
memcpy (a->out, a->in, a->bsize * sizeof (complex));
}
iir.c
+22 -10
View File
@@ -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)