Fast Hartley transform in Python
Computes the discrete Hartley transform in O(n log n) with only real arithmetic.
- File
fast_hartley_transform.py, 77 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/fast-hartley-transform/. The authors are listed in prov/JPR6KM.jsonl.
import math
def fht(x: list[float]) -> None:
"""Replace x by its discrete Hartley transform, H[k] = sum x[j] cas(2 pi j k / n).
Precondition: n = len(x) is a power of two and the entries are real.
Bracewell's radix-2 decimation-in-time algorithm, iterative and in place.
Cost: n lg n - 3n + 4 multiplications and (3/2) n lg n - (3/2) n + 2
additions for n >= 2, plus the n/2 cosine and sine evaluations.
"""
n = len(x) # @line L1
if n < 1 or n & (n - 1):
raise ValueError("the length must be a power of two")
t = n.bit_length() - 1
rev = [0] * n
for j in range(1, n):
rev[j] = (rev[j >> 1] >> 1) | ((j & 1) << (t - 1))
for j in range(n): # @line L2
if j < rev[j]:
x[j], x[rev[j]] = x[rev[j]], x[j]
h = 1 # @line L3
while h < n: # @line L4
for s in range(0, n, 2 * h): # @line L5
u, v = x[s], x[s + h]
x[s], x[s + h] = u + v, u - v # @line L6
if h >= 2: # @line L7
q = h // 2
u, v = x[s + q], x[s + h + q]
x[s + q], x[s + h + q] = u + v, u - v
for k in range(1, h // 2): # @line L8
c, sn = math.cos(math.pi * k / h), math.sin(math.pi * k / h) # @line L9
a, b = x[s + h + k], x[s + 2 * h - k]
t1, t2 = c * a + sn * b, sn * a - c * b # @line L10
e1, e2 = x[s + k], x[s + h - k] # @line L11
x[s + k], x[s + h + k] = e1 + t1, e1 - t1 # @line L12
x[s + h - k], x[s + 2 * h - k] = e2 + t2, e2 - t2 # @line L13
h *= 2 # @line L14
def inverse_fht(y: list[float]) -> None:
"""Replace y by the x with fht(x) == y. Cost: one fht and n divisions."""
fht(y) # @line L15
n = len(y)
for j in range(n): # @line L16
y[j] /= n
def dft_from_dht(big_h: list[float]) -> list[complex]:
"""Return the DFT X[k] = sum x[j] exp(-2 pi i j k / n) from H = DHT(x). Cost: O(n)."""
n = len(big_h)
return [
complex((big_h[k] + big_h[-k % n]) / 2, -(big_h[k] - big_h[-k % n]) / 2)
for k in range(n)
]
def real_cyclic_convolution(a: list[float], b: list[float]) -> list[float]:
"""Cyclic convolution of two real vectors of the same power-of-two length by FHTs.
Uses the Hartley convolution theorem. Cost: three FHTs and O(n) further operations.
"""
n = len(a)
if len(b) != n:
raise ValueError("the vectors must have the same length")
ha, hb = [float(v) for v in a], [float(v) for v in b]
fht(ha)
fht(hb)
hz = [
(ha[k] * (hb[k] + hb[-k % n]) + ha[-k % n] * (hb[k] - hb[-k % n])) / 2
for k in range(n)
]
inverse_fht(hz)
return hz
test_fast_hartley_transform.py 92 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/fast-hartley-transform/. The authors are listed in prov/JPR6KM.jsonl.
import cmath
import math
import random
import pytest
from fast_hartley_transform import dft_from_dht, fht, inverse_fht, real_cyclic_convolution
def naive_dht(x: list[float]) -> list[float]:
n = len(x)
return [
sum(x[j] * (math.cos(2 * math.pi * j * k / n) + math.sin(2 * math.pi * j * k / n)) for j in range(n))
for k in range(n)
]
def naive_dft(x: list[float]) -> list[complex]:
n = len(x)
return [sum(x[j] * cmath.exp(-2j * math.pi * j * k / n) for j in range(n)) for k in range(n)]
def close(a: list, b: list, scale: float) -> bool:
n = max(len(a), 1)
tol = 1e-13 * n * max(n.bit_length(), 1) * max(scale, 1.0)
return len(a) == len(b) and all(abs(u - v) <= tol for u, v in zip(a, b))
def transformed(x: list[float]) -> list[float]:
y = [float(v) for v in x]
fht(y)
return y
def test_examples() -> None:
assert close(transformed([1, 2, 3, 4]), [10, -4, -2, 0], 4)
r = 4 * math.sqrt(2)
assert close(transformed([1, 2, 3, 4, 5, 6, 7, 8]), [36, -8 - r, -8, -r, -4, -8 + r, 0, r], 8)
assert close(dft_from_dht([10, -4, -2, 0]), [10, -2 + 2j, -2, -2 - 2j], 4)
assert close(real_cyclic_convolution([1, 2, 3, 4], [1, 1, 0, 0]), [5, 3, 5, 7], 10)
def test_edge_cases() -> None:
assert transformed([2.5]) == [2.5]
assert transformed([3, 5]) == [8, -2]
for n in (1, 2, 4, 8, 16, 32, 64):
assert transformed([0.0] * n) == [0.0] * n
assert close(transformed([1.0] + [0.0] * (n - 1)), [1.0] * n, 1)
c = transformed([2.0] * n)
assert close(c, [2.0 * n] + [0.0] * (n - 1), 2 * n)
for p in range(n):
imp = [0.0] * n
imp[p] = 1.0
assert close(transformed(imp), naive_dht(imp), 1)
def test_random_against_oracle() -> None:
rng = random.Random(1984)
for case in range(600):
n = 1 << rng.randint(0, 6)
x = [rng.uniform(-10, 10) for _ in range(n)]
scale = 10 * math.sqrt(n)
y = transformed(x)
assert close(y, naive_dht(x), scale)
assert close(dft_from_dht(y), naive_dft(x), scale)
back = list(y)
inverse_fht(back)
assert close(back, x, 10)
if case % 3 == 0:
b = [rng.uniform(-10, 10) for _ in range(n)]
want = [sum(x[i] * b[(k - i) % n] for i in range(n)) for k in range(n)]
assert close(real_cyclic_convolution(x, b), want, 100 * n)
def test_involution_large() -> None:
rng = random.Random(7)
for n in (256, 1024, 4096):
x = [rng.uniform(-1, 1) for _ in range(n)]
y = transformed(x)
assert abs(sum(v * v for v in y) - n * sum(v * v for v in x)) < 1e-9 * n * n
inverse_fht(y)
assert close(y, x, 1)
def test_bad_length() -> None:
with pytest.raises(ValueError):
fht([1.0, 2.0, 3.0])
with pytest.raises(ValueError):
real_cyclic_convolution([1.0, 2.0], [1.0])
What the next level needs
- the tests did not run in this build