Skip to content

Latest commit

 

History

History
218 lines (167 loc) · 8.46 KB

File metadata and controls

218 lines (167 loc) · 8.46 KB

Python bindings

BFFT ships a small ctypes-based Python package that exposes numpy-friendly drop-in transforms. No prebuilt binaries are distributed: installation compiles the native library from source on your machine.

Install

From a clone of the repository:

pip install .

This compiles src/bfft.cpp and src/bodft.cpp with your C++ compiler and bundles the resulting shared library inside the installed package. The only requirements are a C++17 compiler and NumPy.

Optimization flags

Because the library is compiled on your own machine (no prebuilt binaries are distributed), the build tunes for the local CPU by default. It selects, when the compiler accepts them:

  • -O3
  • -march=native (or -mcpu=native on Apple-silicon clang) — emit AVX2 and other supported host-specific instructions.
  • -ffast-math — relaxed floating-point for faster math.

Each flag is probed against your compiler first, so the build degrades gracefully on toolchains that lack them. Control the defaults with environment variables:

Variable Effect
BFFT_NO_NATIVE=1 Skip -march=native / -mcpu=native (portable codegen).
BFFT_NO_FAST_MATH=1 Keep strict IEEE math (drop -ffast-math).
CXX=... Choose the compiler.
BFFT_CXXFLAGS="..." Append extra flags to the compile.
CXX=clang++ BFFT_NO_FAST_MATH=1 pip install .

-ffast-math assumes no NaNs/infinities and reorders operations, so results may differ in the last bits from a strict-IEEE build (still accurate to floating-point precision for these transforms). It also enables flush-to-zero/denormals-are-zero for the process when the library loads, which can affect denormal handling elsewhere. Set BFFT_NO_FAST_MATH=1 if you need bit-reproducible or strict denormal behavior.

Using a system-installed library

If you prefer to build the native library separately, install it first and the Python loader will discover it automatically:

make && sudo make install PREFIX=/usr/local
pip install .

You can also point the loader at a specific shared object:

export BFFT_LIBRARY=/path/to/libbfft.so

The loader searches, in order: $BFFT_LIBRARY, the library bundled in the package, a sibling build/ directory in a source checkout, and finally the system library search path.

Usage

import numpy as np
import bfft

x = np.random.randn(1024)          # power-of-two length

X = bfft.rfft(x)                   # == numpy.fft.rfft(x)        -> N/2 + 1 bins
x_back = bfft.irfft(X)             # == numpy.fft.irfft(X)       -> N samples

H = bfft.odft(x)                   # half-bin-shifted transform  -> N/2 bins
x_back2 = bfft.iodft(H)            # inverse of odft             -> N samples

The module-level functions cache the plan, the transform sizes, and reusable scratch buffers per length internally, so repeated calls at the same size avoid re-creating that state. They stay safe to call from multiple threads: a concurrent call at the same size that cannot reuse the shared scratch falls back to a private buffer.

Planned objects (hot loops)

For the lowest per-call overhead -- transforming the same size repeatedly in a tight loop -- use a planned object. It caches everything except the unavoidable output allocation and input-pointer fetch:

plan = bfft.Plan(N)                # standard real FFT at fixed size N
X = plan.rfft(x)                   # == numpy.fft.rfft(x)
x_back = plan.irfft(X)             # == numpy.fft.irfft(X)

oplan = bfft.OdftPlan(N)           # half-bin transform at fixed size N
H = oplan.odft(x)
x_back2 = oplan.iodft(H)

A planned object owns shared scratch and is not thread-safe: create one plan per thread (or use the module-level functions, which guard against concurrent use).

Pass a caller-owned output buffer with out= to avoid the per-call output allocation entirely (the pyfftw-style zero-allocation loop):

plan = bfft.Plan(N)
out = np.empty(plan.bins, np.complex128)   # allocate once
for chunk in stream:                        # chunk has length N
    plan.rfft(chunk, out=out)               # writes into out, no allocation

out= is accepted by Plan.rfft/irfft and OdftPlan.odft/iodft; it must be a C-contiguous array of the right dtype and length, and is returned as-is.

Calling BFFT from Numba (@njit)

numpy.fft cannot be called from @njit(nopython=True) code -- it is a Python C-extension that only exists in object mode. BFFT can, because it is a plain C ABI taking raw pointers, which Numba lowers through its cffi support. cffi is a dependency of BFFT, so only numba itself needs to be installed alongside it:

import numpy as np
from numba import njit
import bfft.numba_support as bn
from bfft.numba_support import bfft_forward, ffi

N = 4096
plan, bins, work_n, scratch_n = bn.make_plan(N)   # plan is an int address

@njit(cache=True)
def rfft_into(plan, x, out_f64, work, scratch_f64):
    bfft_forward(plan,
                 ffi.from_buffer(x), ffi.from_buffer(out_f64),
                 ffi.from_buffer(work), ffi.from_buffer(scratch_f64))

x = np.random.randn(N)
out = np.empty(bins, np.complex128)
work = np.empty(work_n, np.float64)
scratch = np.empty(scratch_n, np.complex128)
rfft_into(plan, x, out.view(np.float64), work, scratch.view(np.float64))
# out == numpy.fft.rfft(x)

Two rules make it work with Numba: pass the plan as the integer address from make_plan (Numba can type an int but not a raw cffi pointer), and pass complex buffers as their real view (buf.view(np.float64), or buf.view(np.float32) for single precision) so ffi.from_buffer yields the pointer type the C function expects. A JIT-compiled loop then performs each transform with no Python-object interaction -- in practice at the bare C transform speed.

All four transforms are available from @njit in both precisions:

double (float64 / complex128) single (float32 / complex64) transform
bfft_forward bfft_forward_f32 real FFT (rfft)
bfft_inverse bfft_inverse_f32 inverse real FFT (irfft)
bodft_forward bodft_forward_f32 half-bin ODFT (odft)
bodft_inverse bodft_inverse_f32 inverse ODFT (iodft)

Create plans with make_plan(N) for the standard real FFT (pass dtype=np.float32 to size the single-precision work buffer) and make_odft_plan(N) for the ODFT. Both plan helpers return the same (plan, bins, work_n, scratch_n) tuple shape. The ODFT helpers return zero for work_n and scratch_n, and the exported bodft_forward / bodft_forward_f32 callables accept the same work and native_scratch arguments as bfft_forward / bfft_forward_f32 while ignoring them. This lets the same jitted call site switch between rfft and odft by changing only the plan factory and transform function.

API

Function Equivalent Notes
bfft.rfft(x) numpy.fft.rfft(x) Power-of-two N >= 4. Returns N/2 + 1 complex bins.
bfft.irfft(X, n=None) numpy.fft.irfft(X, n) n defaults to 2 * (len(X) - 1). Returns N real samples.
bfft.odft(x) half-bin phase shift + rfft H[k] = sum_n x[n] exp(-2j*pi*(k+1/2)*n/N), N >= 2. Returns N/2 complex bins.
bfft.iodft(H, n=None) inverse of bfft.odft n defaults to 2 * len(H). Returns N real samples.

Python real FFT transforms operate on power-of-two N >= 4 in double precision. The forward and inverse pairs round-trip to floating-point precision, and rfft and irfft match numpy.fft to within floating-point error.

Short-time transforms

bfft.STFTPlan is a reusable native plan for streaming short-time Fourier transforms. It uses the same BFFT real FFT or BODFT half-bin transform for every frame, returns a two-dimensional NumPy complex128 spectrogram with shape (n_bins, n_segs), and stores the inverse overlap-add buffer inside the plan. Call reset_buffer() before starting a fresh inverse stream.

import numpy as np
import bfft

tf = bfft.STFTPlan(n=24576, n_fft=512, hop_length=128)
x = np.random.randn(tf.n)
Zx = tf.stft(x)       # complex128, shape (tf.n_bins, tf.n_segs)
y = tf.istft(Zx)     # float64, length tf.n_segs * tf.hop_length
tf.reset_buffer()    # clear streaming overlap state

Pass transform="odft" to use the half-bin ODFT path. Pass a 1-D float64 window of length n_fft to override the default Hann window; the native plan derives the matching MSE-optimal synthesis window. bfft.hann_window(n_fft) returns the exact default window used by native STFT plans.