pg-fft: Fast Fourier Transforms in PostgreSQL

pg-fft is a PostgreSQL extension to compute Fast Fourier Transforms in SQL.

it's based on https://github.com/pme/pgfft and Mark Borgerding's kissfft. This post documents bringing it up to date and getting it passing tests against PostgreSQL 18.

The C code compiled cleanly against PostgreSQL 18 on the first try.

Fourteen years of Postgres releases and the C API kept every promise it made. Respect to the postgres community for not breaking an extension API in over 14 years.

So, the fork:

  • Re-implements the final function: proper Datum construction, runtime ereport checks instead of compiled-in asserts, {} for empty input, a real error for NULL samples. It reads the float array directly instead of round-tripping every element through deconstruct_array.
  • A new extension script (kissfft--0.1.0.sql) that works with CREATE EXTENSION, marked IMMUTABLE STRICT PARALLEL SAFE.
  • A Makefile you can read in one breath.
  • Actual regression tests: a known 8-point spectrum, peak detection on a two-sine signal, empty input, NULL rejection. make installcheck is green.
  • CI across PostgreSQL 13–18.
  • A README that opens with an example, and a META.json that parses.

Most interesting detail:

The 2012 code wrote the FFT results back like this:

((float *)data)[i] = power;   /* data is a Datum[] */

On a 32-bit build, a Datum is 4 bytes — same as a float — so this worked fine on the machine it was written on. On a 64-bit machine, a Datum is 8 bytes. The results land at the wrong offsets, and the back half of the output array is just... your input, handed back to you.

The fix is one line: data[i] = Float4GetDatum(power).

Does it work now?

Verified against numpy: an 8-point transform of cos(2πk/8) returns {0, 2, 0, 0, 0, 0, 0, 2}, matching np.abs(np.fft.fft(x))**2 / n exactly:

-- 1000 samples of 0.7·sin(2π·50t) + sin(2π·120t) at 1 kHz
SELECT bin - 1 AS freq_hz, round(power::numeric, 1) AS power
  FROM (SELECT fft_agg(v ORDER BY t) AS spectrum FROM samples) s,
       unnest(spectrum) WITH ORDINALITY AS u(power, bin)
 WHERE power > 10 AND bin - 1 < 500;

 freq_hz | power
---------+-------
      50 | 122.5
     120 | 250.0

Two sines in, two peaks out. Right frequencies, right relative power (0.7² : 1²).

Is it fast?

Here's the part that surprised me. A 2²⁰-point (1,048,576-sample) FFT through the aggregate, on an M-series MacBook against PostgreSQL 18:

operation time
array_agg(v ORDER BY i) — accumulate only ~102 ms
fft_agg(v ORDER BY i) — accumulate + FFT ~120 ms

The FFT costs ~20 ms per million points. Eighty-five percent of the wall time is Postgres plumbing — sorting rows and building the transition array — not Fourier math.

What's next

  • kiss_fftr for real-input transforms (2× on the FFT itself, and the natural N/2+1-bin output real signals want)
  • A non-aggregate fft(real[]) → real[] form
  • Complex/magnitude/phase output variants and a float8 version
  • Upgrade the bundled 2010-era KISS FFT to current upstream
  • Maybe a PGXN release

Thanks to Claude, Grok, and Gemini for doing the work and writing this essay.