傅里叶变换 (`scipy.fft`)#

傅里叶分析是一种将函数表示为周期成分之和的方法，并且可以从这些成分中恢复信号。当函数及其傅里叶变换都被替换为离散的对应物时，称为离散傅里叶变换（DFT）。DFT 之所以成为数值计算的主要支柱，部分原因在于一种计算它的非常快速的算法，称为快速傅里叶变换（FFT），该算法由高斯（1805 年）所知，并在库利和图基 [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
>>> # 采样点数量
>>> N = 600
>>> # 采样间隔
>>> 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()

"此代码生成一个 X-Y 图，显示 Y 轴上的振幅与 X 轴上的频率。一条蓝色轨迹的振幅在整个范围内为零，除了两个峰值。较高的第一个峰值在 50 Hz，第二个峰值在 80 Hz。"

FFT 输入信号本质上是截断的。这种截断可以建模为无限信号与矩形窗函数的乘积。在频域中，这种乘积变成了信号频谱与窗函数频谱的卷积，形式为 $\sin(x)/x$。这种卷积是导致称为频谱泄漏效应的原因（参见 [WPW]）。使用专门的窗函数对信号进行加窗有助于减轻这种效应。频谱泄漏。下面的示例使用了来自 scipy.signal 的布莱克曼窗，并展示了加窗的效果（为了说明目的，FFT 的零分量已被截断）。

>>> from scipy.fft import fft, fftfreq
>>> import numpy as np
>>> # 采样点数
>>> N = 600
>>> # 采样间隔
>>> 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-Y 对数线性图，Y 轴为幅度，X 轴为频率。第一条轨迹是 FFT，在 50 和 80 Hz 处有两个峰值，噪声基底约为 1e-2 的幅度。第二条轨迹是加窗后的 FFT，具有相同的两个峰值，但由于窗函数的作用，噪声基底显著降低，约为 1e-7 的幅度。"

如果序列 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；注意频谱的不对称性。

alt:: “此代码生成一个X-Y图，Y轴为幅度，X轴为频率。图中除了在-80和50 Hz处有两个尖峰外，其余部分均为零值。右侧50 Hz处的峰值高度是左侧的两倍。”

>>> from scipy.fft import fft, fftfreq, fftshift
>>> import numpy as np
>>> # 信号点数
>>> N = 400
>>> # 采样间隔
>>> 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的厄米对称性隐含表示（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])

alt:: “此代码生成一个2x3网格排列的六个热图。第一行显示的图像大部分是空白的画布，除了每张图像上有两个微小的红色峰值。第二行显示了每个上方图像的逆FFT的实部。第一列的顶部图像中有两个水平排列的点，底部图像是一个平滑的灰度图，显示了5条黑色垂直条纹，代表二维时域信号。第二列的顶部图像中有两个垂直排列的点，底部图像是一个平滑的灰度图，显示了5条黑色水平条纹，代表二维时域信号。在最后一列中，顶部图像有两个对角线排列的点；对应的下方图像可能有20条黑色条纹，呈60度角。”

>>> 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 种类型。“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{如果 $k = 0$} \\ \sqrt{1/(2N)}, & \text{否则} \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，并使用相同的类型，以获得正确归一化的结果。

>>> # 归一化逆变换：无缩放因子
>>> 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])
>>> # 通过 DCT-I 的未归一化往返：缩放因子 2*(N-1) = 8
>>> dct(dct(x, type=1), type=1)
array([ 8. ,  16.,  8. , -8. ,  12.])
>>> # 归一化逆变换：无缩放因子
>>> 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])
>>> # 未归一化的DCT-IV往返变换：缩放因子2*N = 10
>>> dct(dct(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>> # 归一化的逆变换：无缩放因子
>>> 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()

"此代码生成一个X-Y图，显示Y轴上的幅度和X轴上的时间。第一条蓝色轨迹是原始信号，从幅度1开始并振荡到0幅度，类似于频率啁啾。第二条红色轨迹是使用DCT的x_20重建，在高幅度区域紧密跟随原始信号，但在图的右侧不清晰。第三条绿色轨迹是使用DCT的x_15重建，不如x_20重建精确，但仍然与x相似。"

离散正弦变换 #

SciPy 提供了 DST [Mak] 函数 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）之间的关系。

>>> 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])
>>>  # 缩放因子 2*(N+1) = 12
>>> dst(dst(x, type=1), type=1)
array([ 12.,  24.,  12., -12.,  18.])
>>>  # 无缩放因子
>>> 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])
>>>  # 缩放因子 2*N = 10
>>> dst(dst(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>>  # 无缩放因子
>>> 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, Numerical Recipes: The Art of Scientific Computing, ch. 12-13. Cambridge Univ. Press, Cambridge, UK.

[Mak] (1,2)

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)

1. 1. 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

傅里叶变换 (scipy.fft)#

傅里叶变换 (`scipy.fft`)#