Cosine, sine and Hartley transforms

Real-to-real transforms map real data to real data. They are the transforms behind JPEG (DCT-II), spectral Poisson solvers with Dirichlet or Neumann boundary conditions (DST-I, DCT-I), and the Hartley transform. FFTL computes all eleven kinds, along any dimensions, each dimension with its own kind.

The kinds

The formulas are unnormalised. is the length along the transformed dimension.

constant	kind
REDFT00	DCT-I
REDFT10	DCT-II
REDFT01	DCT-III
REDFT11	DCT-IV
RODFT00	DST-I
RODFT10	DST-II
RODFT01	DST-III
RODFT11	DST-IV
DHT	Hartley
R2HC	real DFT, halfcomplex order	for  , then
HC2R	its inverse	from halfcomplex input

R2HC is the DFT   of real data, stored in the same reals as the input: the real parts for  , the imaginary parts for    at position  .

Pairs and scale factors

The inverses are the same kind, or a partner kind, times a factor. For each kind, applying it twice gives the input times the factor below, so the normalised inverse divides by it:

kind	applied twice (or with its partner) gives
DCT-I (REDFT00)
DST-I (RODFT00)
DCT-II / DCT-III, DST-II / DST-III, DCT-IV, DST-IV
DHT, R2HC with HC2R

The example checks the self-inverse kinds directly:

using FFTL, LinearAlgebra

n = 9
x = [sin(0.7 * j) + 0.1 * j for j in 0:n-1]
for (kind, factor) in ((FFTL.REDFT00, 2(n - 1)), (FFTL.REDFT11, 2n),
                       (FFTL.RODFT00, 2(n + 1)), (FFTL.RODFT11, 2n), (FFTL.DHT, n))
  y = FFTL.r2r(FFTL.r2r(x, kind), kind) ./ factor
  println(rpad(string(kind), 10), maximum(abs, y - x) < 1e-12)
end

REDFT00   true
REDFT11   true
RODFT00   true
RODFT11   true
DHT       true

The orthonormal DCT-II

dct is the DCT-II scaled to an orthonormal map, and idct its exact inverse. It is the transform most users want for compression and for energy-preserving spectra:

c = randn(16)
D = dct(c)
(norm(D) ≈ norm(c), idct(D) ≈ c)

(true, true)

Its values are the unnormalised DCT-II divided by for the zero frequency and by for the others, which is the usual orthonormal scaling.

R = FFTL.r2r(c, FFTL.REDFT10)
D[1] ≈ R[1] / (2 * sqrt(16)), D[2] ≈ R[2] / sqrt(2 * 16)

(true, true)

Several dimensions with different kinds

A region takes one kind per dimension. The example uses DCT-II along the rows and the real DFT along the columns:

M = randn(8, 6)
Y = FFTL.r2r(M, (FFTL.REDFT10, FFTL.R2HC), (1, 2))
size(Y)

(8, 6)

Plans work as for the complex transforms: plan_r2r returns a plan, inv of a plan gives its normalised inverse for the pairs above, and plan_r2r_pass transforms one dimension of a strided array.