Scale the Kaiser-Bessel window to unit peak.

[jenskeiner]

This PR is to improve the numerical properties of Kaiser-Bessel window function evaluations.

In contrast to other windows, the peak of the Kaiser-Bessel window function as used in the code grows with m. For example, in frequency space PHI_HUT(n,k,d) is I0(m·b) for k = 0, where b = π(2 − 1/σ) and σ is the oversampling factor, typically somewhere between 2 and 4. This grows quickly with m and applies per dimension. The time domain PHI(n,x,d) has similar properties where x = 0. For d = 3 and m >= 8, this overflows in single precision, making the transform unusable.

Since the NFFT algorithm decomposes into CONV ∘ FFT ∘ DECONV, where the DECONV step divides by PHI_HUT and CONV multiplies by PHI, a scaling by 1 / I0(m·b) of both PHI_HUT and PHI keeps the result identical, but removes the undesired scaling.

It is important to never evaluate I0(m·b) directly, but to use the representation exp(−L) with L = log I0(m·b). The exponent -L can then be directly integrated into the existing calculations that use exponentials already.

The scaling approach itself adds to another issue since the calculation relies on taking the difference of exponents of nearly equal magnitude in three different places. Done naively, this destroys accuracy. Cancellations are avoided by rewriting the respective results into a form without the problematic subtractions.

For a more detailed explanation, see kaiser-bessel-window-numerics.html.

[github-actions]

Accuracy report

2254 unchanged · 103 improved · 3 regressed

module	improved	regressed
nfct	18	2
nfft	69	0
nfst	16	1

📊 Full accuracy report

[jenskeiner]

Leaving as draft until benchmarks show no regressions.