Last active
February 4, 2021 16:32
-
-
Save endolith/6322721 to your computer and use it in GitHub Desktop.
Relationship between DCT, RFFT, FFT, DFT in SciPy
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
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 | |
""" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment