傅立葉轉換 (scipy.fft)

傅立葉分析是一種將函數表示為週期性成分總和，並從這些成分中恢復訊號的方法。當函數及其傅立葉轉換都替換為離散化對應物時，它被稱為離散傅立葉轉換 (DFT)。DFT 已成為數值計算的主要支柱，部分原因在於一種非常快速的演算法，用於計算它，稱為快速傅立葉轉換 (FFT)，高斯 (1805) 已知此演算法，並由 Cooley 和 Tukey [CT65] 以目前的形式公諸於世。Press 等人 [NR07] 提供了傅立葉分析及其應用程式的易懂介紹。

快速傅立葉轉換

一維離散傅立葉轉換

長度為 \(N\) 的 FFT y[k]，對於長度為 \(N\) 的序列 x[n]，定義如下

\[y[k] = \sum_{n=0}^{N-1} e^{-2 \pi j \frac{k n}{N} } x[n] \, ,\]

且反轉換定義如下

\[x[n] = \frac{1}{N} \sum_{k=0}^{N-1} e^{2 \pi j \frac{k n}{N} } y[k] \, .\]

這些轉換可以分別透過 fft 和 ifft 來計算，如下例所示。

>>> from scipy.fft import fft, ifft
>>> import numpy as np
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
>>> y = fft(x)
>>> y
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
       -1.83155948+1.60822041j, -1.83155948-1.60822041j,
        2.08155948+1.65109876j])
>>> yinv = ifft(y)
>>> yinv
array([ 1.0+0.j,  2.0+0.j,  1.0+0.j, -1.0+0.j,  1.5+0.j])

從 FFT 的定義可以看出

\[y[0] = \sum_{n=0}^{N-1} x[n] \, .\]

在範例中

>>> np.sum(x)
4.5

這對應於 \(y[0]\)。對於偶數 N，元素 \(y[1]...y[N/2-1]\) 包含正頻率項，而元素 \(y[N/2]...y[N-1]\) 包含負頻率項，依負頻率遞減順序排列。對於奇數 N，元素 \(y[1]...y[(N-1)/2]\) 包含正頻率項，而元素 \(y[(N+1)/2]...y[N-1]\) 包含負頻率項，依負頻率遞減順序排列。

在序列 x 為實數值的情況下，正頻率的 \(y[n]\) 值是負頻率的 \(y[n]\) 值的共軛（因為頻譜是對稱的）。通常，只繪製對應於正頻率的 FFT。

此範例繪製了兩個正弦波之和的 FFT。

>>> from scipy.fft import fft, fftfreq
>>> import numpy as np
>>> # Number of sample points
>>> N = 600
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
>>> yf = fft(y)
>>> xf = fftfreq(N, T)[:N//2]
>>> import matplotlib.pyplot as plt
>>> plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
>>> plt.grid()
>>> plt.show()

FFT 輸入訊號本質上是被截斷的。這種截斷可以建模為無限訊號與矩形窗函數的乘法。在頻譜域中，這種乘法變成訊號頻譜與窗函數頻譜的卷積，形式為 \(\sin(x)/x\)。此卷積是稱為頻譜洩漏效應的原因（請參閱 [WPW]）。使用專用窗函數對訊號進行加窗有助於減輕頻譜洩漏。以下範例使用 scipy.signal 中的 Blackman 窗，並顯示加窗的效果（為了說明目的，FFT 的零組件已被截斷）。

>>> from scipy.fft import fft, fftfreq
>>> import numpy as np
>>> # Number of sample points
>>> N = 600
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
>>> yf = fft(y)
>>> from scipy.signal.windows import blackman
>>> w = blackman(N)
>>> ywf = fft(y*w)
>>> xf = fftfreq(N, T)[:N//2]
>>> import matplotlib.pyplot as plt
>>> plt.semilogy(xf[1:N//2], 2.0/N * np.abs(yf[1:N//2]), '-b')
>>> plt.semilogy(xf[1:N//2], 2.0/N * np.abs(ywf[1:N//2]), '-r')
>>> plt.legend(['FFT', 'FFT w. window'])
>>> plt.grid()
>>> plt.show()

在序列 x 為複數值的情況下，頻譜不再是對稱的。為了簡化 FFT 函數的使用，scipy 提供了以下兩個輔助函數。

函數 fftfreq 傳回 FFT 樣本頻率點。

>>> from scipy.fft import fftfreq
>>> freq = fftfreq(8, 0.125)
>>> freq
array([ 0., 1., 2., 3., -4., -3., -2., -1.])

以類似的精神，函數 fftshift 允許交換向量的下半部和上半部，使其適合顯示。

>>> from scipy.fft import fftshift
>>> x = np.arange(8)
>>> fftshift(x)
array([4, 5, 6, 7, 0, 1, 2, 3])

以下範例繪製了兩個複數指數的 FFT；請注意非對稱頻譜。

>>> from scipy.fft import fft, fftfreq, fftshift
>>> import numpy as np
>>> # number of signal points
>>> N = 400
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.exp(50.0 * 1.j * 2.0*np.pi*x) + 0.5*np.exp(-80.0 * 1.j * 2.0*np.pi*x)
>>> yf = fft(y)
>>> xf = fftfreq(N, T)
>>> xf = fftshift(xf)
>>> yplot = fftshift(yf)
>>> import matplotlib.pyplot as plt
>>> plt.plot(xf, 1.0/N * np.abs(yplot))
>>> plt.grid()
>>> plt.show()

函數 rfft 計算實數序列的 FFT，並僅輸出頻率範圍一半的複數 FFT 係數 \(y[n]\)。剩餘的負頻率成分由實數輸入的 FFT 的 Hermitian 對稱性暗示 (y[n] = conj(y[-n]))。在 N 為偶數的情況下：\([Re(y[0]) + 0j, y[1], ..., Re(y[N/2]) + 0j]\)；在 N 為奇數的情況下 \([Re(y[0]) + 0j, y[1], ..., y[N/2]\)。明確顯示為 \(Re(y[k]) + 0j\) 的項被限制為純實數，因為根據 Hermitian 屬性，它們是自身的複數共軛。

對應的函數 irfft 計算具有此特殊排序的 FFT 係數的 IFFT。

>>> from scipy.fft import fft, rfft, irfft
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5, 1.0])
>>> fft(x)
array([ 5.5 +0.j        ,  2.25-0.4330127j , -2.75-1.29903811j,
        1.5 +0.j        , -2.75+1.29903811j,  2.25+0.4330127j ])
>>> yr = rfft(x)
>>> yr
array([ 5.5 +0.j        ,  2.25-0.4330127j , -2.75-1.29903811j,
        1.5 +0.j        ])
>>> irfft(yr)
array([ 1. ,  2. ,  1. , -1. ,  1.5,  1. ])
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
>>> fft(x)
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
       -1.83155948+1.60822041j, -1.83155948-1.60822041j,
        2.08155948+1.65109876j])
>>> yr = rfft(x)
>>> yr
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
        -1.83155948+1.60822041j])

請注意，奇數和偶數長度訊號的 rfft 具有相同的形狀。預設情況下，irfft 假定輸出訊號應為偶數長度。因此，對於奇數訊號，它將給出錯誤的結果

>>> irfft(yr)
array([ 1.70788987,  2.40843925, -0.37366961,  0.75734049])

為了恢復原始的奇數長度訊號，我們必須透過 n 參數傳遞輸出形狀。

>>> irfft(yr, n=len(x))
array([ 1. ,  2. ,  1. , -1. ,  1.5])

二維與 N 維離散傅立葉轉換

函數 fft2 和 ifft2 分別提供二維 FFT 和 IFFT。同樣地，fftn 和 ifftn 分別提供 N 維 FFT 和 IFFT。

對於實數輸入訊號，與 rfft 類似，我們有函數 rfft2 和 irfft2 用於二維實數轉換；rfftn 和 irfftn 用於 N 維實數轉換。

以下範例示範了二維 IFFT，並繪製了產生的（二維）時域訊號。

>>> from scipy.fft import ifftn
>>> import matplotlib.pyplot as plt
>>> import matplotlib.cm as cm
>>> import numpy as np
>>> N = 30
>>> f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col', sharey='row')
>>> xf = np.zeros((N,N))
>>> xf[0, 5] = 1
>>> xf[0, N-5] = 1
>>> Z = ifftn(xf)
>>> ax1.imshow(xf, cmap=cm.Reds)
>>> ax4.imshow(np.real(Z), cmap=cm.gray)
>>> xf = np.zeros((N, N))
>>> xf[5, 0] = 1
>>> xf[N-5, 0] = 1
>>> Z = ifftn(xf)
>>> ax2.imshow(xf, cmap=cm.Reds)
>>> ax5.imshow(np.real(Z), cmap=cm.gray)
>>> xf = np.zeros((N, N))
>>> xf[5, 10] = 1
>>> xf[N-5, N-10] = 1
>>> Z = ifftn(xf)
>>> ax3.imshow(xf, cmap=cm.Reds)
>>> ax6.imshow(np.real(Z), cmap=cm.gray)
>>> plt.show()

離散餘弦轉換

SciPy 使用函數 dct 提供 DCT，並使用函數 idct 提供對應的 IDCT。DCT 有 8 種型別 [WPC]、[Mak]；但是，scipy 中僅實作了前 4 種型別。「The」DCT 通常指第二型 DCT，而「the」反 DCT 通常指第三型 DCT。此外，DCT 係數可以以不同的方式正規化（對於大多數型別，scipy 提供 None 和 ortho）。dct/idct 函數呼叫的兩個參數允許設定 DCT 型別和係數正規化。

對於單維陣列 x，dct(x, norm='ortho') 等於 MATLAB dct(x)。

第一型 DCT

SciPy 使用以下未正規化 DCT-I 的定義 (norm=None)

\[y[k] = x_0 + (-1)^k x_{N-1} + 2\sum_{n=1}^{N-2} x[n] \cos\left(\frac{\pi nk}{N-1}\right), \qquad 0 \le k < N.\]

請注意，DCT-I 僅支援輸入大小 > 1。

第二型 DCT

SciPy 使用以下未正規化 DCT-II 的定義 (norm=None)

\[y[k] = 2 \sum_{n=0}^{N-1} x[n] \cos \left({\pi(2n+1)k \over 2N} \right) \qquad 0 \le k < N.\]

在正規化 DCT (norm='ortho') 的情況下，DCT 係數 \(y[k]\) 乘以比例因子 f

\[\begin{split}f = \begin{cases} \sqrt{1/(4N)}, & \text{if $k = 0$} \\ \sqrt{1/(2N)}, & \text{otherwise} \end{cases} \, .\end{split}\]

在這種情況下，DCT「基底函數」\(\phi_k[n] = 2 f \cos \left({\pi(2n+1)k \over 2N} \right)\) 變為正交

\[\sum_{n=0}^{N-1} \phi_k[n] \phi_l[n] = \delta_{lk}.\]

第三型 DCT

SciPy 使用以下未正規化 DCT-III 的定義 (norm=None)

\[y[k] = x_0 + 2 \sum_{n=1}^{N-1} x[n] \cos\left({\pi n(2k+1) \over 2N}\right) \qquad 0 \le k < N,\]

或者，對於 norm='ortho'

\[y[k] = {x_0\over\sqrt{N}} + {2\over\sqrt{N}} \sum_{n=1}^{N-1} x[n] \cos\left({\pi n(2k+1) \over 2N}\right) \qquad 0 \le k < N.\]

第四型 DCT

SciPy 使用以下未正規化 DCT-IV 的定義 (norm=None)

\[y[k] = 2 \sum_{n=0}^{N-1} x[n] \cos\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,\]

或者，對於 norm='ortho'

\[y[k] = \sqrt{2\over N}\sum_{n=0}^{N-1} x[n] \cos\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N\]

DCT 與 IDCT

（未正規化）DCT-III 是（未正規化）DCT-II 的反函數，比例因子為 2N。正交正規化 DCT-III 正好是正交正規化 DCT-II 的反函數。函數 idct 執行 DCT 和 IDCT 型別之間的映射，以及正確的正規化。

以下範例顯示了不同型別和正規化的 DCT 和 IDCT 之間的關係。

>>> from scipy.fft import dct, idct
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

DCT-II 和 DCT-III 互為反函數，因此對於正交正規轉換，我們會返回原始訊號。

>>> dct(dct(x, type=2, norm='ortho'), type=3, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])

但是，在預設正規化下執行相同的操作時，由於正向轉換未正規化，因此我們獲得了額外的比例因子 \(2N=10\)。

>>> dct(dct(x, type=2), type=3)
array([ 10.,  20.,  10., -10.,  15.])

因此，我們應該對兩者都使用相同型別的函數 idct，以獲得正確正規化的結果。

>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=2), type=2)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

對於 DCT-I，也可以看到類似的結果，DCT-I 是自身的反函數，比例因子為 \(2(N-1)\)。

>>> dct(dct(x, type=1, norm='ortho'), type=1, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])
>>> # Unnormalized round-trip via DCT-I: scaling factor 2*(N-1) = 8
>>> dct(dct(x, type=1), type=1)
array([ 8. ,  16.,  8. , -8. ,  12.])
>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=1), type=1)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

對於 DCT-IV，也是自身的反函數，比例因子為 \(2N\)。

>>> dct(dct(x, type=4, norm='ortho'), type=4, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])
>>> # Unnormalized round-trip via DCT-IV: scaling factor 2*N = 10
>>> dct(dct(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=4), type=4)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

範例

DCT 展現「能量集中特性」，這表示對於許多訊號，只有前幾個 DCT 係數具有顯著的幅度。將其他係數歸零會導致較小的重建誤差，這一事實在有損訊號壓縮（例如 JPEG 壓縮）中得到利用。

以下範例顯示訊號 x 和從訊號的 DCT 係數進行的兩個重建 (\(x_{20}\) 和 \(x_{15}\))。訊號 \(x_{20}\) 是從前 20 個 DCT 係數重建的，\(x_{15}\) 是從前 15 個 DCT 係數重建的。可以看出，使用 20 個係數的相對誤差仍然非常小（約 0.1%），但提供了五倍的壓縮率。

>>> from scipy.fft import dct, idct
>>> import matplotlib.pyplot as plt
>>> N = 100
>>> t = np.linspace(0,20,N, endpoint=False)
>>> x = np.exp(-t/3)*np.cos(2*t)
>>> y = dct(x, norm='ortho')
>>> window = np.zeros(N)
>>> window[:20] = 1
>>> yr = idct(y*window, norm='ortho')
>>> sum(abs(x-yr)**2) / sum(abs(x)**2)
0.0009872817275276098
>>> plt.plot(t, x, '-bx')
>>> plt.plot(t, yr, 'ro')
>>> window = np.zeros(N)
>>> window[:15] = 1
>>> yr = idct(y*window, norm='ortho')
>>> sum(abs(x-yr)**2) / sum(abs(x)**2)
0.06196643004256714
>>> plt.plot(t, yr, 'g+')
>>> plt.legend(['x', '$x_{20}$', '$x_{15}$'])
>>> plt.grid()
>>> plt.show()

離散正弦轉換

SciPy 使用函數 dst 提供 DST [Mak]，並使用函數 idst 提供對應的 IDST。

理論上，DST 有 8 種型別，適用於偶數/奇數邊界條件和邊界偏移的不同組合 [WPS]，scipy 中僅實作了前 4 種型別。

第一型 DST

DST-I 假設輸入在 n=-1 和 n=N 附近為奇數。SciPy 使用以下未正規化 DST-I 的定義 (norm=None)

\[y[k] = 2\sum_{n=0}^{N-1} x[n] \sin\left( \pi {(n+1) (k+1)}\over{N+1} \right), \qquad 0 \le k < N.\]

另請注意，DST-I 僅支援輸入大小 > 1。（未正規化）DST-I 是自身的反函數，比例因子為 2(N+1)。

第二型 DST

DST-II 假設輸入在 n=-1/2 附近為奇數，在 n=N 附近為偶數。SciPy 使用以下未正規化 DST-II 的定義 (norm=None)

\[y[k] = 2 \sum_{n=0}^{N-1} x[n] \sin\left( {\pi (n+1/2)(k+1)} \over N \right), \qquad 0 \le k < N.\]

第三型 DST

DST-III 假設輸入在 n=-1 附近為奇數，在 n=N-1 附近為偶數。SciPy 使用以下未正規化 DST-III 的定義 (norm=None)

\[y[k] = (-1)^k x[N-1] + 2 \sum_{n=0}^{N-2} x[n] \sin \left( {\pi (n+1)(k+1/2)} \over N \right), \qquad 0 \le k < N.\]

第四型 DST

SciPy 使用以下未正規化 DST-IV 的定義 (norm=None)

\[y[k] = 2 \sum_{n=0}^{N-1} x[n] \sin\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,\]

或者，對於 norm='ortho'

\[y[k] = \sqrt{2\over N}\sum_{n=0}^{N-1} x[n] \sin\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,\]

DST 與 IDST

以下範例顯示了不同型別和正規化的 DST 和 IDST 之間的關係。

>>> from scipy.fft import dst, idst
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

DST-II 和 DST-III 互為反函數，因此對於正交正規轉換，我們會返回原始訊號。

>>> dst(dst(x, type=2, norm='ortho'), type=3, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])

但是，在預設正規化下執行相同的操作時，由於正向轉換未正規化，因此我們獲得了額外的比例因子 \(2N=10\)。

>>> dst(dst(x, type=2), type=3)
array([ 10.,  20.,  10., -10.,  15.])

因此，我們應該對兩者都使用相同型別的函數 idst，以獲得正確正規化的結果。

>>> idst(dst(x, type=2), type=2)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

對於 DST-I，也可以看到類似的結果，DST-I 是自身的反函數，比例因子為 \(2(N-1)\)。

>>> dst(dst(x, type=1, norm='ortho'), type=1, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])
>>>  # scaling factor 2*(N+1) = 12
>>> dst(dst(x, type=1), type=1)
array([ 12.,  24.,  12., -12.,  18.])
>>>  # no scaling factor
>>> idst(dst(x, type=1), type=1)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

對於 DST-IV，也是自身的反函數，比例因子為 \(2N\)。

>>> dst(dst(x, type=4, norm='ortho'), type=4, norm='ortho')
array([ 1. ,  2. ,  1. , -1. ,  1.5])
>>>  # scaling factor 2*N = 10
>>> dst(dst(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>>  # no scaling factor
>>> idst(dst(x, type=4), type=4)
array([ 1. ,  2. ,  1. , -1. ,  1.5])

快速漢克爾轉換

SciPy 提供了函數 fht 和 ifht，用於對數間隔的輸入陣列執行快速漢克爾轉換 (FHT) 及其反轉換 (IFHT)。

FHT 是由 [Ham00] 定義的連續漢克爾轉換的離散化版本

\[A(k) = \int_{0}^{\infty} \! a(r) \, J_{\mu}(kr) \, k \, dr \;,\]

其中 \(J_{\mu}\) 是 \(\mu\) 階的貝索函數。在變數變更 \(r \to \log r\)、\(k \to \log k\) 下，這變為

\[A(e^{\log k}) = \int_{0}^{\infty} \! a(e^{\log r}) \, J_{\mu}(e^{\log k + \log r}) \, e^{\log k + \log r} \, d{\log r}\]

這是對數空間中的卷積。FHT 演算法使用 FFT 對離散輸入資料執行此卷積。

必須小心以最大程度地減少由於 FFT 卷積的循環性質而引起的數值振鈴。為了確保低振鈴條件 [Ham00] 成立，可以使用 fhtoffset 函數計算的偏移量稍微移動輸出陣列。

參考文獻

[CT65]

Cooley, James W. 和 John W. Tukey，1965 年，「複雜傅立葉級數的機器計算演算法」，Math. Comput. 19: 297-301。

[NR07]

Press, W.、Teukolsky, S.、Vetterline, W.T. 和 Flannery, B.P.，2007 年，《數值食譜：科學計算的藝術》，第 12-13 章。Cambridge Univ. Press, Cambridge, UK。

[Mak] (1,2)

J. Makhoul，1980 年，「一維和二維快速餘弦轉換」，IEEE Transactions on acoustics, speech and signal processing vol. 28(1), pp. 27-34, DOI:10.1109/TASSP.1980.1163351

[Ham00] (1,2)

A. J. S. Hamilton，2000 年，「非線性功率譜的不相關模式」，MNRAS, 312, 257。 DOI:10.1046/j.1365-8711.2000.03071.x

[WPW]

https://en.wikipedia.org/wiki/Window_function

[WPC]

https://en.wikipedia.org/wiki/Discrete_cosine_transform

[WPS]

https://en.wikipedia.org/wiki/Discrete_sine_transform