Compute 2N-point real FFT w/ a N-point complex FFT

fftre2n.py

import numpy as np

# Compute 2N-point real FFT w/ a N-point complex FFT, cf. TI SPRA 291

# random 2N-point real sequence
g = np.random.randn(16)

# size of sequence
n = len(g)//2
n2 = 2 * n

# form a N-point complex sequence and its DFT
x = g[::2] + 1j*g[1::2]
X = np.fft.fft(x)

# retain x0 for idft
X0 = X[0]

# some helper terms
W = lambda k, N: np.exp(-2j*np.pi*k/N)
A = lambda k: 0.5*(1 - 1j*W(k,n2))
B = lambda k: 0.5*(1 + 1j*W(k,n2))

# form the positive part of the DFT for the 2N-point real sequence
k = np.r_[:n]
G1 = X*A(k) + np.conj(X)[np.r_[0,n-k[1:]]]*B(k)

# check it's correct
G2 = np.fft.fft(g)[np.fft.fftfreq(n2)>=0]
np.testing.assert_allclose(G1, G2)

# back the other way w/ same fft func
k1 = np.r_[0,n-k[1:]]
X1 = G2*np.conj(A(k)) + np.conj(G2)[k1]*np.conj(B(k))
X1[0] = X0
np.testing.assert_allclose(X1, X)
x1 = np.conj(np.fft.fft(np.conj(X1)))/n
g1 = np.c_[x1.real, x1.imag].flat[:]
np.testing.assert_allclose(g1, g)