Repository navigation
Several improvements: const for *Setup, zconvolve w/o scaling, zero-phase convolution #101
Description
I've already proposed these to Julien of PFFFT. Repeating here as these apply to DOUBLE version as well.
-
As I wanted to optimize memory use in a multi-threaded application, I asked myself - can the pffft_new_setup() called once for a specific N, and then used for concurrent calls of pffft_transform(), only with differing "work" buffers? However, the setup argument of pffft_transform() is non-const. It would be great if a const PFFFT_Setup * is specified in the function's declaration - if it actually uses setup in read-only manner.
-
It would be great if you could add a pffft_zconvolve() function that puts the result into dft_ab (may be equal to dft_b) without accumulation and scaling. A scaling factor is unnecessary, because it's more efficient to rescale the dft_a because it's a static filter used numerous times, usually.
-
Another suggestion is to provide a function (e.g. pffft_zconvert_zp) that permits convolving a fft with "zero-phase" fft (impulse responses). Linear-phase filters all have zero-phase response if the time-domain FIR filter's center is at t=0 and left-hand taps are wrapped around the FFT block (placed right-to-left at the end of the input transform block). The FFT transform of such a filter yields zero complex components. I do that already in r8brain-free-src and this is very efficient. I transform filter's FFT into "zero-phase" FFT by copying real components to complex components (except first and last bin of course as they have no complex components). However, unordered pffft requires knowing the ordering of real and complex components. So, a function is necessary for the unordered PFFFT result, to obtain an "zero-phase" FFT first (a scaling factor can be specified in this function). Then a simple sequential multiplication would work - for both unordered and ordered fft; compilers usually vectorise that well. It seems that engineers probably do not know this simple optimization trick. Unordered PFFFT is substantially faster so I wanted to utilize it, but currently can't - if I use a usual complex-domain convolution, the efficiencies of zero-phase convolution are lost.
-
PFFFTD defines _USE_MATH_DEFINES unconditionally while sometimes _USE_MATH_DEFINES is defined as compilator-level option. It would be useful to add #ifndef _USE_MATH_DEFINES before defining it.
Activity
marton78 commentedon Aug 24, 2026
Hi @avaneev, all good suggestions, thanks! I would be grateful if you could review my attempt at implementing them: #102
avaneev commentedon Aug 25, 2026
Thank you. The issue with zconvolve is that it still applies inverse N gain. I can't get correct convolution gain with scale=1. If I change the call to pffftd_zconvert_zp( s -> setup, p, p, 1.0/N ), it works. Beside that, pffftd_zconvolve loses precision - with ordered FFT and hand-written complex multiplication in r8brain I can get -200 dB null-test while with pffftd_zconvolve I get only -167 dB.
The reason you need a special zconvolve_zp is because you are applying an inverse gain which should not be there. If you remove it, zconvolve_zp won't be necessary, it's a redundant function. For backwards compatibility, it's better to apply this gain inside the pffftd_zconvolve_no_accum() function.
avaneev commentedon Aug 25, 2026
Also, the performance of pffftd_zconvolve_zp() is not too great - note the code SSE2 in r8brain CDSPRealFFT multiplyBlocksZP() function - it runs much faster. With this function and unordered FFT I get 5% overall performance improvement, while with pffftd_zconvolve_zp() I get only 2% (on Ryzen 9950X with GCC 13). But currently I can't use pffftd_zconvert_zp() in production, because it lowers precision - looks like double-to-float conversion sneaks in somewhere.
marton78 commentedon Aug 26, 2026
Hi @avaneev,
Thank you for your insights.
1. Inverse-N gain
Yes, makes sense. Will do.
BTW, your claim that
zconvert_zpbecomes redundant is only true for ordered spectra. When unordered, whole vectors alternate real/imag, soY_im = X_im · H_reneeds cross-vector indexing2. Double→float precision loss
I don't think this comes from us. I asked the Ox Alpha model (which is free for a few days on Openrouter), it gave me:
pffft_double.cdoes#define float double; every double backend union isdouble f[4](pf_avx_double.h:60, sse2/neon/scalar variants too). No float-typed temporary exists in the shared impl.pffftd_z*convolve*/zconvert*functions on this build — zerofcvt*instructions.pffftd_zconvolveoutput matches hand-written complex multiply on ordered spectra to −321 dB (pure double rounding); full zp round-trip vs time-domain circular convolution: −308 dB.His −167 dB cannot originate in these kernels. Likely harness-side — during this verification I stepped twice into the same trap he may have: feeding an unordered spectrum to
transform_ordered(.., BACKWARD)or using the wrongzreorderdirection silently produces plausible-but-wrong output (I measured garbage at "−3.7 dB" that way). His −200 dB baseline also indicates other error sources dominate his null test either way.3.
zconvolve_zpperformanceAgain, tested by Ox Alpha:
Kernel-only on M2/NEON, ns/call (
.perf/bench101.c):No pathology — zp delivers its ~2-multiplies-vs-6-per-pair advantage, degrading toward memory-bound at large N. Since the multiply stage is small next to the FFTs, a ~1.5× kernel win correctly shrinks to a few percent end-to-end — his 2% here vs 5% with r8brain SSE2 is plausible without any bug. Side observation: plain
zconvolvealready beatszconvolve_scaleby 10–20% because the runtime scaling multiply can't be constant-folded — the unscaled entry point captures part of the optimization on its own.marton78 commentedon Aug 26, 2026
Re: unit-gain zero-phase convolution
I think it's best to fold the 1/N into
zconvert_zponly, because this would make the zp pipeline unit-gain:=> circular convolution, like r8brain's convention.
Why conversion-side and not in the kernel:
zconvolve_zp, it would cost one extraVMULper re/im vector pair in every call. The kernel advantage overzconvolve_scaleis already small (see above), and this would lose roughly a third of it.zconvolve_scalewould breakpffastconv, which feeds pre-scaled filters tozconvolve_scale.zconvolve_accumulatekeeps upstream's spectral-scaling semantics;zconvolve_scale/zconvolvestay spectral-domain operations.New
zconvert_zpsemantics:output = zero-phase spectrum * scaling/N, i.e. thescalingbecomes an additional gain on top of unit-gain normalization.Unit-gain zero-phase convolution: normalize zconvert_zp by 1/N
marton78 commentedon Aug 26, 2026
Please have a look at the new commit, @avaneev, and my comments above.
As for "the performance of pffftd_zconvolve_zp() is not too great", I have a performance branch in preparation, which I hope will improve things slightly.
avaneev commentedon Aug 28, 2026
Thank you.
avaneev commentedon Aug 28, 2026
You may be right that zconvolve_zp may be needed. I did assume ordered FFT where sequential multiplication works correctly.
avaneev commentedon Aug 28, 2026
I stand corrected about my critique. In r8brain-free I do spectrum mirroring which requires ordered FFT, so my null-test evaluations were incorrect and your latest commit is in fact completely correct.
marton78 commentedon Aug 28, 2026
amazing! then i'll merge
marton78 commentedon Sep 2, 2026
Hi @avaneev you mentioned earlier that "the performance of pffftd_zconvolve_zp() is not too great". Could you give it another try with the new
perfbranch? I hope it should become faster now.avaneev commentedon Sep 3, 2026
Hi @marton78 After a more thorough testing the performance was actually good in commit 102. I've utilized the "zp" functions already in r8brain-free-src. Thanks again for your effort.
The
pffftd_zconvolve_zp()in theperfbranch is actually 4% slower in my SSE2 tests compared to commit 102.avaneev commentedon Sep 3, 2026
I would also suggest to remove SCALAR1 branch completely along with #if ( SIMD_SZ == 4 ) branching - this would make SIMD_SZ=4 uniform throughout. SCALAR1 is adds no value from performance standpoint. And I guess scaling to e.g. SIMD_SZ=8 is technically impossible.
avaneev commentedon Sep 3, 2026
I'd like to also point out to
restrictrelated contract violations (which were inherited from the original PFFFT). E.g. in zconvolve_accumulate thevabmay point to bothvaorvbyet all pointers were defined with therestrictspecifier. The issues do not practically happen, but it's a formal violation. I have not checked other instances.restrictis tricky - it's a contract that a pointer does not alias any other pointer in the code context.marton78 commentedon Sep 3, 2026
@avaneev
Perfect question for AI. I am only mildly ashamed of being a meat proxy in this case and will paste you the relevant parts, edited for brevity:
The function is byte-identical between
masterandperf, it has no hunk touchingFUNC_ZCONVOLVE_ZP/FUNC_ZCONVERT_ZP, andpf_sse2_double.honly gained aVMSUBmacro that these kernels don't use. I cross-compiledpffft_double.cfor x86-64/SSE2 at both commits (clang -O3) and diffed the disassembly:_pffftd_zconvolve_zpis the same 45 instructions. Same forzconvert_zp,zconvolve,zconvolve_scale,zconvolve_accumulate. No build-flag changes for x86 either (target_optimizations.cmakediff is FMA/ARMv7 messages + WASM only). So an isolated-kernel 4 % delta is noise/thermal/code-placement, a few percent swing between two builds is normal.Did you measure the whole pipeline (transform → zp → transform) or just kernel only? What did change on SSE2 double:
passf2/3/5(SVMADD restructure, fused multiply-add instead of mul+add ops)passf8(N=32 only)radf4/radb4and inlinedradf3/radb3gotido==3fast paths (only active for N = 3 * 2^k)So on SSE2, which is a non-FMA target, the FFT itself is mostly a code-layout-level change.
If you want, you could try to pin it down with
perf stat -r 20or-falign-functions=64.SCALAR1 exists not for performance, but for sizes below/not multiple of 32. The cost of maintaining it is small, I'd keep it.
Correct, unfortunately. I looked into supporting AVX for floats, but that would be a HUGE rewrite for questionable benefit.
Holy shit, thats a GREAT catch, thanks! Fixed on
master.Fortunately, the generated code with/without
RESTRICTis instruction-identical for all six kernels on x86-64 and arm64, float and double (with clang -O3), so there is no performance cost to this fix, and that's also why it never happened in practice so far.marton78 commentedon Sep 3, 2026
Also added on
masterin the branch just merged: expanded floating point constants for the twiddle factors (sin&cosofpi*p/q) to double precision. Before, they had only floating point precision.Drop restrict from zconvolve kernels whose pointers may alias by cont…
avaneev commentedon Sep 3, 2026
That's okay about LLM reply. We'll all soon be senior developers riding AI spines. The performance difference I've got may be due to code layout differences, nothing serious. I've simply measured the
__rdtsc()tick count of thepffftd_zconvolve_zp()function call.