* E-241: fft/spec — fix amplitude normalization for non-power-of-2 records
First numerical-correctness (wrong-number, not crash) find of a correctness-audit campaign that checked ngspice analyses against analytic and numpy oracles. com_fft.c has two paths: FFTW3 (exact length-point transform, scale=length/2, correct) and a Green's radix-2 FFT (no FFTW3) that zero-pads the `length` input samples up to the next power of two N. The non-FFTW path normalized the single-sided amplitude by the PADDED size (scale=N/2) instead of the true sample count (length/2), so every FFT whose sample count is not a power of two reported amplitudes too small by length/N -- up to 2x. A .tran essentially never yields exactly 2^k points, so this bit the common case silently. The DC bin is the clean discriminator (always exactly bin 0, no scalloping): a 2.0 V DC offset read back as 2.0*length/N, e.g. 1.0009 for a 1025-sample record padded to 2048. spec had the same error in its power normalization (intres = N*N instead of length*length; the FFTW path already used length*length). Fix: scale = length/2 and intres = length*length in the non-FFTW path, matching FFTW. The frequency axis was already padding-aware, so only the magnitude scale needed changing; power-of-2 records (length == N) and FFTW-linked builds are unchanged. Validated against numpy: after the fix ngspice's fft matches numpy.fft.rfft of the same zero-padded record (normalized by length) bin-by-bin to ~1e-9. fft/spec amplitudes now change for non-power-of-2 records (they become correct). Regression 199/199. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
parent
4ed89173f2
commit
6a9610efd4
|
|
@ -80,7 +80,14 @@ com_fft(wordlist *wl)
|
|||
M++;
|
||||
}
|
||||
fpts = N/2 + 1;
|
||||
scale = ((double)N)/2;
|
||||
/* E-241: normalize the amplitude by the ACTUAL sample count `length`, not the
|
||||
* zero-padded FFT size N. The Green's radix-2 FFT pads `length` samples up to
|
||||
* the next power of two N; the single-sided amplitude of a bin is 2*|X|/length
|
||||
* (independent of padding -- padding only interpolates the spectrum), so the
|
||||
* scale must be length/2 to match the FFTW path above. Using N/2 made every
|
||||
* FFT whose sample count is not a power of two report amplitudes too small by
|
||||
* length/N (up to 2x): e.g. a pure DC offset read half its value. */
|
||||
scale = ((double)length)/2.0;
|
||||
#endif
|
||||
|
||||
win = TMALLOC(double, length);
|
||||
|
|
@ -460,7 +467,10 @@ com_psd(wordlist *wl)
|
|||
fftFree();
|
||||
|
||||
/* 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]). */
|
||||
intres = (double)N * (double)N;
|
||||
/* E-241: normalize the PSD by length^2 (the actual sample count), not the
|
||||
* zero-padded FFT size N^2 -- matching the FFTW path above. Using N^2 made
|
||||
* the PSD of a non-power-of-two-length record too small by (length/N)^2. */
|
||||
intres = (double)length * (double)length;
|
||||
fdvec[i][0].cx_real = reald[0]*reald[0]/intres;
|
||||
fdvec[i][0].cx_imag = 0;
|
||||
noipower = fdvec[i][0].cx_real;
|
||||
|
|
|
|||
Loading…
Reference in New Issue