Relationship between DCT, RFFT, FFT, DFT in SciPy

DCT Type I.png

DCT Type II.png

DCT_vs_cos.py

"""
Illustrate correspondence between DCT bins and cosine signals.
"""
from scipy.fftpack import idct
from numpy import cos, arange, zeros, pi
import matplotlib.pyplot as plt
from numpy.testing import assert_allclose

N = 8  # only even numbers
k = arange(N)


###################################### Type 1
plt.figure(1)

sig = cos(2*pi*k/N)
plt.plot(k,   sig, '.-')
plt.plot(k+N, sig, '.-')

spectrum = zeros(N//2+1)
spectrum[1] = 1

plt.plot(k[:N//2+1], idct(spectrum, 1)/2, 'o-')
plt.margins(0, 0.1)

assert_allclose(sig[:N//2+1], idct(spectrum, 1)/2, rtol=1e-05, atol=1e-08)


###################################### Type 2
plt.figure(2)

sig = cos(2*pi*(k/N + 0.5/N))  # half-sample shift left
plt.plot(k,   sig, '.-')
plt.plot(k+N, sig, '.-')

spectrum = zeros(N//2)
spectrum[1] = 1

plt.plot(k[:N//2], idct(spectrum, 2)/2, 'o-')
plt.margins(0, 0.1)

assert_allclose(sig[:N//2], idct(spectrum, 2)/2, rtol=1e-05, atol=1e-08)

plt.show()

DCT_vs_FFT.py

"""
Type I DCT produces the same output as RFFT, FFT if signal is real and
symmetrical but DCT is twice as fast as RFFT and RFFT is twice as fast as FFT.
"""

import numpy as np

# scipy's fft does rfft automatically, so don't use it for comparisons
# scipy's rfft outputs in a weird "packed complex" format, so use numpy's
# instead
from numpy.fft import fft, rfft
from scipy.fftpack import dct

# Prove they're the same
a = np.array([1.,  2.,  3.,  4.,  3.,  2.])
N = len(a)

print(fft(a).real)
print(rfft(a).real)
print(dct(a[:N//2+1], 1))

"""
FFT:  [ 15.  -4.   0.  -1.   0.  -4.]
RFFT: [ 15.  -4.   0.  -1.]
DCT:  [ 15.  -4.   0.  -1.]
"""

# Test speed on real, even-symmetric signal
N = 10000

# Create symmetrical real signal
a = np.random.rand(N//2+1)
a = np.concatenate((a, a[-2:0:-1]))

# Confirm it's symmetrical
assert np.allclose(fft(a).real, fft(a))

# Confirm that they produce the same output
assert np.allclose(fft(a)[:N//2+1], rfft(a))
assert np.allclose(rfft(a), dct(a[:N//2+1], 1))

"""
N = 2**20

In [87]: timeit fft(a)
1 loops, best of 3: 240 ms per loop

In [88]: timeit rfft(a)
10 loops, best of 3: 126 ms per loop

In [89]: timeit dct(a[:N/2+1], 1)
10 loops, best of 3: 48.6 ms per loop

DCT 2.6x RFFT
DCT 4.9x FFT
"""

"""
N = 65537

In [95]: timeit fft(a)
100 loops, best of 3: 10.2 ms per loop

In [96]: timeit rfft(a)
100 loops, best of 3: 5.18 ms per loop

In [97]: timeit dct(a[:N/2+1], 1)
1000 loops, best of 3: 1.62 ms per loop

DCT 3.2x RFFT
DCT 6.3x FFT
"""