* Enhancement-306: ngspice — the E-241 twin in the fft expression function

Enhancement-241 fixed an amplitude normalization that divided by the ZERO-PADDED
transform size instead of the number of input samples -- in the `fft` COMMAND
(frontend/com_fft.c). The identical mistake survived in maths/cmaths/cmath4.c, the
vector-expression function reached by `let F = fft(v)`: a separate implementation of
the same computation.

Found by continuing the oracle campaign that produced E-241, and located with E-241's
own discriminator, the DC bin. One signal, 4001 samples padded to 4096, DC offset 2.0:

  fft s          ; the COMMAND    ->  mag(s)[0] = 2.000000    correct
  let F = fft(s) ; the FUNCTION   ->  mag(F)[0] = 1.953613    = 2.0 * 4001/4096

X[0] is the sum of the samples, D*length for a DC offset D, so dividing by the padded N
reads back D*length/N. As E-241 put it: a DC value cannot depend on how many samples
were taken.

This is a contradiction, not a matter of convention. cx_fft holds TWO complete
implementations -- one for complex input, one for real -- and each has an FFTW branch
and a Green's radix-2 branch. In BOTH, the FFTW branch already used the input length
while Green's used the padded size:

  real branch     FFTW: scale = ((double)length)/2.0    Green: ((double)N)/2   <- wrong
  complex branch  FFTW: scale = (double) fpts           Green: (double) N      <- wrong

the correct version sitting a few lines from the wrong one inside the same function --
the same shape as `avg` disagreeing with `integ` in E-302. HAVE_LIBFFTW3 is undefined in
this build, so Green's is the live path and the defect was reachable.

  oracle                              before      after       closed form
  real-input, DC bin                  1.953613    2.000000    2.0
  real-input, ifft(fft(x))            2.3e-02     1.1e-16     0
  complex-input, bin 0                0.9766780   0.9998683   0.9998683

The round trip is an INDEPENDENT confirmation: nothing here touches ifft, so its going
from 2.3% error to machine precision is evidence from a direction the fix did not aim
at -- the pair only inverts when the forward normalization is right. cx_ifft was
audited and deliberately left alone for that reason.

Every caller of the Green's radix-2 kernel was audited, not only the one that failed:
com_fft.c's fft (x2) and spec/PSD (x2) are correct from E-241; cx_ifft is correct;
trannoise/1-f-code.c cannot pad at all, since n_pts is grown to 2^n_exp by construction;
and fft/ifft are the only transform functions in the expression table, so there is no
spec twin to miss.

Verification: examples/fftexpr_examples/verify_fftexpr.py -- 6 checks under both
solvers, all against closed form (DC bin on both paths, round trip padded and unpadded,
complex-input bin 0 against the analytic mean of an RC response). It scores 3/6 on the
pre-fix binary, so it is a real regression guard. E-241's own suite (fftnorm_examples)
and ifftreal_examples pass unchanged. Full sweep 241/241 OK.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Meisam Bahadori 2026-08-02 19:56:26 +02:00 committed by Holger Vogt
parent 9ee0028132
commit 79ba4be561
1 changed files with 22 additions and 2 deletions

View File

@ -745,7 +745,17 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
*newlength = N;
outdata = alloc_c(N);
scale = (double) N;
/* Enhancement-306: normalise by the number of INPUT samples, not by the
zero-padded transform size. Green's radix-2 kernel pads `length` samples
up to the next power of two N, but zero padding only interpolates the
spectrum -- it adds no signal, so it must not change an amplitude. The
unambiguous case is bin 0: X[0] is the sum of the samples = D*length for
a DC offset D, so dividing by N reads back D*length/N instead of D, and
"a DC value cannot depend on how many samples were taken". This is the
twin of the Enhancement-241 bug, which fixed the identical mistake in the
`fft` COMMAND (frontend/com_fft.c, scale = length/2) but not here in the
vector-expression function reached by `let F = fft(v)`. */
scale = (double) length;
for (i = 0; i < N; i++) {
outdata[i].cx_real = datax[2*i]/scale;
outdata[i].cx_imag = datax[2*i+1]/scale;
@ -803,7 +813,17 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
rffts(datax, M, 1);
fftFree();
scale = ((double)N)/2;
/* Enhancement-306: normalise by the number of INPUT samples, not the
zero-padded transform size -- the FFTW twin four lines up already
uses ((double)length)/2.0. Green's radix-2 kernel pads `length`
samples to the next power of two N, but zero padding only
interpolates the spectrum; it adds no signal and must not change an
amplitude. Bin 0 is the unambiguous case: X[0] is the sum of the
samples = D*length for a DC offset D, so dividing by N reads back
D*length/N. This is the twin of Enhancement-241, which fixed the
identical mistake in the `fft` COMMAND (frontend/com_fft.c) but not
here in the vector-expression function reached by `let F = fft(v)`. */
scale = ((double)length)/2;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
outdata[0].cx_real = datax[0]/scale/2.0;
outdata[0].cx_imag = 0.0;