1.5.4 快速傅立叶变换：scipy.fftpack

scipy.fftpack 模块允许计算快速傅立叶变换。例子，一个（有噪音）的信号输入是这样：

In [12]:

time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + \
       0.5 * np.random.randn(time_vec.size)

观察者并不知道信号的频率，只知道抽样时间步骤的信号sig。假设信号来自真实的函数，因此傅立叶变换将是对称的。scipy.fftpack.fftfreq() 函数将生成样本序列，而将计算快速傅立叶变换：

In [13]:

from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d=time_step)
sig_fft = fftpack.fft(sig)

因为生成的幂是对称的，寻找频率只需要使用频谱为正的部分：

In [14]:

pidxs = np.where(sample_freq > 0)
freqs = sample_freq[pidxs]
power = np.abs(sig_fft)[pidxs]

寻找信号频率：

In [15]:

freq = freqs[power.argmax()]
np.allclose(freq, 1./period)  # 检查是否找到了正确的频率

Out[15]:

True

现在高频噪音将从傅立叶转换过的信号移除：

In [16]:

sig_fft[np.abs(sample_freq) > freq] = 0

生成的过滤过的信号可以用scipy.fftpack.ifft()函数：

In [17]:

main_sig = fftpack.ifft(sig_fft)

查看结果：

In [18]:

import pylab as plt
plt.figure()
plt.plot(time_vec, sig)
plt.plot(time_vec, main_sig, linewidth=3)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')

/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/numpy/core/numeric.py:462: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Out[18]:

<matplotlib.text.Text at 0x107484b10>

numpy.fft

Numpy也有一个FFT(numpy.fft)实现。但是，通常scipy的实现更受欢迎，因为，他使用更高效的底层实现。

实例：寻找粗略周期

实例：高斯图片模糊

弯曲:

$f_1(t) = \int dt'\, K(t-t') f_0(t')$

$\tilde{f}_1(\omega) = \tilde{K}(\omega) \tilde{f}_0(\omega)$

练习：月球登陆图片降噪

• 检查提供的图片moonlanding.png，图片被周期噪音污染了。在这个练习中，我们的目的是用快速傅立叶变换清除噪音。

• 用pylab.imread()加载图片。

• 寻找并使用在scipy.fftpack中的2-D FFT函数，绘制图像的频谱（傅立叶变换）。在可视化频谱时是否遇到了麻烦？如果有的话，为什么？

• 频谱由高频和低频成分构成。噪音被包含在频谱的高频部分，因此将那些部分设置为0（使用数组切片）。

• 应用逆傅立叶变换来看一下结果图片。