Add BDCT kernel: paired radix-4 half-shift transform + DCT-IV#41
Merged
Conversation
Adds a new header-only kernel (src/detail/bdct_kernel.hpp) alongside the
Bruun real kernel, computing the native half-bin-shifted transform
H[k] = sum_n x[n] exp(-2*pi*i*(k+1/2)*n/N)
for real input via a paired radix-4 decimation. Each output position k is
paired with its partner kp = M-1-k inside the radix-4 combine so the two
share the three child twiddle multiplies (t, t^2, t^3), halving the twiddle
multiplies of the naive radix-4 half-shift combine. Output is the contiguous
packed lower half H[0..N/2-1]; no DC/Nyquist special case, no final reorder.
The inverse is the exact algebraic inverse (machine-precision roundtrip).
A length-N DCT-IV is built on top via a length-2N half-shift of the
zero-padded input plus one output rotation.
Details:
- Templated for float and double (no long double; matches the real kernel).
- Reuses heap_array, complex_t/complex_f32_t and the SIMD-backend macros
from bruun_kernel.hpp; only the decimation arithmetic is new.
- Explicit 128-bit SoA SIMD specialization for the double forward combine,
running on NEON and SSE2 through the shared Bruun V2 macros; the scalar
path remains the exact reference for every backend, precision and tail.
Public API alongside the bfft headers:
- include/bfft/bdct.h (C ABI; reuses bfft_complex / bfft_status)
- include/bfft/bdct.hpp (RAII bfft::half_shift / bfft::dctiv)
- src/bdct.cpp implements the C ABI into libbfft.
- Makefile and CMakeLists wire the new source and headers into both builds.
Validation: examples/dctiv_benchmark.cpp checks correctness in float and
double via analytic exact tests (impulse half-shift, DCT-IV basis vector),
an independent same-precision direct-DFT cross-check, and the exact inverse
roundtrip, plus timing. Existing test suite still passes.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What
Adds BDCT, a new kernel alongside the Bruun real kernel that computes the native half-bin-shifted transform
for real input, plus a DCT-IV built on top of it. New files:
src/detail/bdct_kernel.hpp— the kernel (header-only)include/bfft/bdct.h,include/bfft/bdct.hpp— public C ABI + RAII C++ wrapper, alongside the bfft headerssrc/bdct.cpp— C ABI implementation compiled intolibbfftexamples/dctiv_benchmark.cpp— validation + benchmarkMakefile,CMakeLists.txt— build wiringWhy / how
Same Bruun-style decimation logic, specialized to the half-shifted spectrum. The decimation is radix-4: for
N = 4M, each child is itself a length-Mhalf-shift transform of real data, soA_r[M-1-k] = conj(A_r[k]). Pairing each output positionkwith its partnerkp = M-1-kinside the radix-4 combine lets the two share the three child twiddle multiplies (t, t^2, t^3) — halving the twiddle multiplies of the naive radix-4 half-shift combine, the same way the real kernel pairs conjugate residues.Because the input is real,
H[N-1-k] = conj(H[k]), so the kernel emits only the contiguous packed lower halfH[0..N/2-1]— no DC/Nyquist special case, no final reorder. The inverse is the exact algebraic inverse of the forward (machine-precision roundtrip, no1/N).A length-
NDCT-IV is recovered from a length-2Nhalf-shift of the zero-padded input followed by one output rotation.Notes for reviewers
long double— it is justdoubleon the NEON target). The kernel is templated;half_shift_plan/_f32anddctiv_plan/_f32.heap_array,complex_t/complex_f32_tand the SIMD-backend macros frombruun_kernel.hpp; only the decimation arithmetic is new.bdct_combine_fwd_simd_f64) that runs on NEON and SSE2 through the shared BruunV2_*macros — children/twiddles are transposed into real/imag lane registers and partner positions are written by address arithmetic (no vector reversal). The scalar path remains the exact reference for every backend, precision, and tail. Wider AVX / float / inverse vectorization are left as follow-ups rather than shipping intrinsics that can't be compiled or tested on the arm64 host.Validation
examples/dctiv_benchmark.cppchecks correctness in both precisions via:|H[k]| = 1, closed-form phase) and DCT-IV of a basis vector (exactN/2spike);Double: roundtrip ~1e-15, analytic/DFT agreement at the reference's own accumulation floor. Float: ~1e-7. The existing
make testsuite still passes.