CQF预备知识：Python相关库 -- 傅里叶变换 scipy.fft

文中内容仅限技术学习与代码实践参考，市场存在不确定性，技术分析需谨慎验证，不构成任何投资建议。

傅里叶变换（scipy.fft）

傅里叶分析是一种将函数表示为周期分量之和，并从这些分量中恢复信号的方法。当函数及其傅里叶变换都被离散化时，就称为离散傅里叶变换（DFT）。由于存在一种非常快速的计算算法，即快速傅里叶变换（FFT），DFT 已成为数值计算的基石之一。FFT 在高斯（1805）时期就已知晓，并由库利和图基以当前形式提出¹。Press 等人²提供了傅里叶分析及其应用的通俗介绍。

快速傅里叶变换

一维离散傅里叶变换

长度为 $N$ 的序列 x[n] 的长度为 $-N$ 的 FFT y[k] 定义为

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

逆变换定义如下

$$x[n]=\frac{1}{N}\sum _{k=0}^{N-1}{e}^{2\pi j\frac{kn}{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()

FFT 输入信号本质上是截断的。这种截断可以被建模为无限信号与矩形窗函数的乘积。在频谱域中，这种乘积变成了信号频谱与窗函数频谱的卷积，窗函数频谱的形式为 $\sin (x)/x$。这种卷积是导致所谓频谱泄漏（参见³）的原因。使用专用窗函数对信号进行窗化有助于减轻频谱泄漏。以下示例使用 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 是复值的，那么频谱不再对称。为了简化与 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
>>> # 信号点数
>>> 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])

二维和 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()

离散余弦变换

SciPy 提供了函数 dct 来进行 DCT，以及对应的函数 idct 来进行逆 DCT。DCT 有 8 种类型^{4 5}，但 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),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(\frac{\pi (2n+1)k}{2N}\right)0\le k<N.$$

在归一化 DCT（norm='ortho'）的情况下，DCT 系数 $y[k]$ 被乘以一个缩放因子 f ：

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

在这种情况下，DCT“基函数”${\varphi }_k[n]=2f\cos \left(\frac{\pi (2n+1)k}{2N}\right)$ 变为正交归一：

$$\sum _{n=0}^{N-1}{\varphi }_k[n]{\varphi }_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(\frac{\pi n(2k+1)}{2N}\right)0\le k<N,$$

或者，对于 norm='ortho' ：

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

类型 IV DCT

SciPy 使用以下定义的未归一化 DCT-IV（norm=None）：

$$y[k]=2\sum _{n=0}^{N-1}x[n]\cos \left(\frac{\pi (2n+1)(2k+1)}{4N}\right)0\le k<N,$$

或者，对于 norm='ortho' ：

$$y[k]=\sqrt{\frac{2}{N}}\sum _{n=0}^{N-1}x[n]\cos \left(\frac{\pi (2n+1)(2k+1)}{4N}\right)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()

离散正弦变换

SciPy 提供了函数 dst 来进行 DST⁵，以及对应的函数 idst 来进行逆 DST。

理论上，对于不同的偶/奇边界条件和边界偏移的组合⁶，DST 有 8 种类型，但 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(\frac{\pi (n+1)(k+1)}{N+1}\right),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(\frac{\pi (n+1/2)(k+1)}{N}\right),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(\frac{\pi (n+1)(k+1/2)}{N}\right),0\le k<N.$$

类型 IV DST

SciPy 使用以下定义的未归一化 DST-IV（norm=None）：

$$y[k]=2\sum _{n=0}^{N-1}x[n]\sin \left(\frac{\pi (2n+1)(2k+1)}{4N}\right),0\le k<N,$$

或者，对于 norm='ortho' ：

$$y[k]=\sqrt{\frac{2}{N}}\sum _{n=0}^{N-1}x[n]\sin \left(\frac{\pi (2n+1)(2k+1)}{4N}\right),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])
>>> # 缩放因子 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 是连续汉克尔变换的离散版本，定义为⁷：

$$A(k)={\int }_0^{\infty }a(r)J_{\mu }(kr)kdr,$$

其中 $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 卷积的循环性质，必须尽量减少数值振铃。为了确保低振铃条件⁷成立，输出数组可以通过 fhtoffset 函数计算的偏移量进行微小的平移。

参考文献

风险提示与免责声明
本文内容基于公开信息研究整理，不构成任何形式的投资建议。历史表现不应作为未来收益保证，市场存在不可预见的波动风险。投资者需结合自身财务状况及风险承受能力独立决策，并自行承担交易结果。作者及发布方不对任何依据本文操作导致的损失承担法律责任。市场有风险，投资须谨慎。

---

1. Cooley, James W., and John W. Tukey, 1965, “An algorithm for the machine calculation of complex Fourier series,” Math. Comput. 19: 297-301. ↩︎

2. 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. ↩︎

3. https://en.wikipedia.org/wiki/Window_function ↩︎

4. https://en.wikipedia.org/wiki/Discrete_cosine_transform ↩︎

5. 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 ↩︎ ↩︎

6. https://en.wikipedia.org/wiki/Discrete_sine_transform ↩︎

7. 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 ↩︎ ↩︎