傅里叶变换 (scipy.fft)#

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

快速傅里叶变换#

一维离散傅里叶变换#

长度为 \(N\) 的序列 x[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 的厄米对称性（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\) 的项被限制为纯实数，因为根据厄米性质，它们是自身的复共轭。

相应的函数 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；以及用于 N 维实数变换的 rfftn 和 irfftn。

下面的示例演示了二维 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)。

DCT I 型#

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 II 型#

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 III 型#

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 IV 型#

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 种类型。

DST I 型#

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 II 型#

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 III 型#

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 IV 型#

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, “复数傅里叶级数机器计算算法,” Math. Comput. 19: 297-301.

[NR07]

Press, W., Teukolsky, S., Vetterline, W.T., and Flannery, B.P., 2007, 数值食谱：科学计算的艺术, ch. 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