傅里叶变换 (scipy.fft)#

傅里叶分析是一种将函数表示为周期分量之和，并从这些分量中恢复信号的方法。当函数及其傅里叶变换都被离散化对应物取代时，它被称为离散傅里叶变换 (DFT)。DFT 已成为数值计算的主要内容之一，部分原因是它有一种计算速度非常快的算法，称为快速傅里叶变换 (FFT)，高斯 (1805) 已知该算法，并由 Cooley 和 Tukey [CT65] 以其当前形式呈现。Press 等人 [NR07] 提供了傅里叶分析及其应用的易懂介绍。

快速傅里叶变换#

一维离散傅里叶变换#

长度为 \(N\) 的序列 x[n] 的长度为 \(N\) 的 FFT y[k] 定义为

\[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()
"This code generates an X-Y plot showing amplitude on the Y axis vs frequency on the X axis. A single blue trace has an amplitude of zero all the way across with the exception of two peaks. The taller first peak is at 50 Hz with a second peak at 80 Hz."

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()
"This code generates an X-Y log-linear plot with amplitude on the Y axis vs frequency on the X axis. The first trace is the FFT with two peaks at 50 and 80 Hz and a noise floor around an amplitude of 1e-2. The second trace is the windowed FFT and has the same two peaks but the noise floor is much lower around an amplitude of 1e-7 due to the window function."

如果序列 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()
"This code generates an X-Y plot with amplitude on the Y axis vs frequency on the X axis. The trace is zero-valued across the plot except for two sharp peaks at -80 and 50 Hz. The 50 Hz peak on the right is twice as tall."

函数 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()
"This code generates six heatmaps arranged in a 2x3 grid. The top row shows mostly blank canvases with the exception of two tiny red peaks on each image. The bottom row shows the real-part of the inverse FFT of each image above it. The first column has two dots arranged horizontally in the top image and in the bottom image a smooth grayscale plot of 5 black vertical stripes representing the 2-D time domain signal. The second column has two dots arranged vertically in the top image and in the bottom image a smooth grayscale plot of 5 horizontal black stripes representing the 2-D time domain signal. In the last column the top image has two dots diagonally located; the corresponding image below has perhaps 20 black stripes at a 60 degree angle."

离散余弦变换#

SciPy 提供了带有函数 dct 的 DCT 和带有函数 idct 的相应 IDCT。DCT 有 8 种类型 [WPC], [Mak]；但是，scipy 中只实现了前 4 种类型。“DCT”通常指 DCT 类型 2，“逆 DCT”通常指 DCT 类型 3。此外，DCT 系数可以不同地归一化（对于大多数类型，scipy 提供 None 和 ortho）。dct/idct 函数调用的两个参数允许设置 DCT 类型和系数归一化。

对于一维数组 x，dct(x, norm='ortho') 等于 MATLAB dct(x)。

I 型 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。

II 型 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}.\]

III 型 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.\]

IV 型 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，其自身是逆的，相差一个因子 \(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()
"This code generates an X-Y plot showing amplitude on the Y axis and time on the X axis. The first blue trace is the original signal and starts at amplitude 1 and oscillates down to 0 amplitude over the duration of the plot resembling a frequency chirp. The second red trace is the x_20 reconstruction using the DCT and closely follows the original signal in the high amplitude region but it is unclear to the right side of the plot. The third green trace is the x_15 reconstruction using the DCT and is less precise than the x_20 reconstruction but still similar to x."

离散正弦变换#

SciPy 提供了带有函数 [Mak] 的 DST dst 和带有函数 idst 的相应 IDST。

理论上，DST 有 8 种类型，用于偶数/奇数边界条件和边界偏移的不同组合 [WPS]，scipy 中只实现了前 4 种类型。

I 型 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)。

II 型 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.\]

III 型 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.\]

IV 型 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，其自身是逆的，相差一个因子 \(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., and John W. Tukey, 1965, “An algorithm for the machine calculation of complex Fourier series,” Math. Comput. 19: 297-301。

[NR07]

Press, W., Teukolsky, S., Vetterline, W.T., and Flannery, B.P., 2007, Numerical Recipes: The Art of Scientific Computing, ch. 12-13. Cambridge Univ. Press, Cambridge, UK。

[Mak] (1,2)

J. Makhoul, 1980, ‘A Fast Cosine Transform in One and Two Dimensions’, 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, “Uncorrelated modes of the non-linear power spectrum”, MNRAS, 312, 257. DOI:10.1046/j.1365-8711.2000.03071.x