傅立叶变换 (`scipy.fft` )¶

傅立叶分析是一种将函数表示为周期分量之和，并从这些分量恢复信号的方法。当函数及其傅里叶变换都用离散化的对应物代替时，称为离散傅立叶变换(DFT)。DFT已经成为数值计算的中流砥柱，部分原因是一种计算它的非常快的算法，称为快速傅立叶变换(FFT)，高斯(1805年)知道这种算法，并由Cooley和Tukey以目前的形式曝光 [CT65]. Press等人。 [NR07] 介绍傅立叶分析及其应用。

快速傅立叶变换 ¶

一维离散傅立叶变换 ¶

快速傅立叶变换 y[k] 长度的 $N$ 的长度 -$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
>>> 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的定义可以看出，

\[是的 [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
>>> # 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.ignal的Blackman窗口，并显示了窗口的效果(出于说明的目的，FFT的零分量已被截断)。

>>> from scipy.fft import fft, fftfreq
>>> # 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 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
>>> # 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的厄米对称性隐含 (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])

为了恢复原始的奇数长度信号，我们 must 将输出形状传递给 n 参数。

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

2-D和N-D离散傅立叶变换 ¶

功能 fft2 和 ifft2 分别提供二维FFT和IFFT。同样， fftn 和 ifftn 分别提供N-D FFT和IFFT。

对于实输入信号，类似于 rfft ，我们有功能 rfft2 和 irfft2 对于二维实变换； rfftn 和 irfftn 用于N-D实数变换。

下面的示例演示了2DIFFT，并绘制了结果(2D)时域信号。

>>> from scipy.fft import ifftn
>>> import matplotlib.pyplot as plt
>>> import matplotlib.cm as cm
>>> 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]; 但是，只有前4种类型在Scipy中实现。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 )：

\[是的 [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\]

请注意，仅当输入大小大于1时才支持DCT-I。

II型DCT ¶

SciPy使用非标准化DCT-II的以下定义 (norm=None )：

\[是的 [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 ：

\[F=\BEGIN{CASES}\sqrt{1/(4N)}，&\text{if$k=0$}\sqrt{1/(2N)}， &\text{否则}\end{case}\，。\]

在本例中，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 )：

\[是的 [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 )：

\[是的 [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' ：

\[是的 [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和两个重建 ($x_{{20}}$ 和 $x_{{15}}$ )从信号的DCT系数。该信号 $x_{{20}}$ 是从前20个DCT系数重建的， $x_{{15}}$ 是从前15个DCT系数重建的。可以看出，使用20个系数的相对误差仍然很小(~0.1%)，但提供了5倍的压缩率。

>>> 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 [Mak] 使用函数 dst 以及与该函数对应的IDST idst 。

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

第I类DST ¶

dst-i假设输入在n=-1和n=N附近是奇数，本网站使用非标准化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.\]

另请注意，仅当输入大小大于1时才支持DST-I。(未标准化的)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 )：

\[是的 [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 )：

\[是的 [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' ：

\[是的 [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_{\µ}(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: 詹姆斯·W·库利和约翰·W·图基，1965，“机器计算复傅立叶级数的算法”， 数学课。电脑。 19：297-301。
NR07: Press，W.，Teukolsky，S.，Vetterline，W.T.，和Flannery，B.P.，2007， 数字食谱：科学计算的艺术 、长春花属(Ch.)12-13。剑桥大学。PRESS，英国剑桥。
Mak(1,2): J.Makhoul，1980，“一维和二维的快速余弦变换”， IEEE Transactions on acoustics, speech and signal processing 卷。28(1)，第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

插值 (scipy.interpolate )

信号处理 (scipy.signal )

傅立叶变换 (scipy.fft )¶

傅立叶变换 (`scipy.fft` )¶