scipy通过快速傅里叶变换实现滤波

fft模块简介

scipy官网宣称，fftpack模块将不再更新，或许不久之后将被废弃，也就是说fft将是唯一的傅里叶变换模块。

Fourier变换极其逆变换在数学上的定义如下

$$\begin{gathered} F(\omega )={\int }_{?\infty }^{\infty }f(t)e^{?i\omega t}\text{d}t \\ f(t)=\frac{\pi }{2}{\int }_{?\infty }^{\infty }F(\omega )e^{i\omega t}\text{d}\omega \end{gathered}$$

下表整理出一部分与Fourier变换相关的函数，其中FFT为快速Fourier变换(Fast Fourier Transform)；DCT为离散余弦变换(Discrete Cosine Transform)；DST为离散正弦变换(discrete sine transform)，另外，函数的前缀和后缀有如下含义

• i表示逆变换；
• 2, n分别表示2维和n维

	正变换	逆变换
通用	fft, fft2, fftn	ifft, ifft2, ifftn
实数域	rfft, rfft2, rfftn	irfft, irfft2, irfftn
厄米对称	hfft, hfft2, hfftn	ihfft, ihfft2, ihfftn
离散余弦变换	dct, dctn	idct, idctn
离散正弦变换	dst, dstn	idst, idstn
汉克尔变换	fht	ifht
移动零频	fftshift	ifftshift
DFT采样频率	fftfreq	ifftfreq

fft函数示例

在数值计算中，一切输入输出均为离散值，所以实际上用到的是离散Fourier变换，即DFT，其功能是将离散的采样变换为一个离散的频谱分布。

下面将手动创建一组采样点，并添加一点噪声，然后通过FFT获取其频域信息。

import numpy as np
from scipy import fft

PI = np.pi*2
fs = 60     #采样频率
T = 100     #采样周期数
N = fs*T    #采样点
t = np.linspace(0, T, N)
noise = 2 * np.random.randn(*t.shape)
s = 2 * np.sin(PI * t) + 3 * np.sin(22 * PI * t) + noise
F = fft.fft(s)
f = fft.fftfreq(N, 1.0/fs)

其中，t为时间序列，s为模拟的采样点，F是Fourier变换之后的结果。但由于fft默认是在复数域上的，故而可以查看其实部、虚部、模和辐角的值。

下面对采样点以及Fourier变换的结果进行绘制

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(2,2,1)
ax.plot(t, s)
ax.set_title("t vs s")
f_abs = np.abs(F)
ax = fig.add_subplot(2,2,2)
ax.plot(f, f_abs)
ax.set_title("fs vs |F|")

xlims = [[0,2], [21,23]]
for i, xlim in enumerate(xlims):
    ax = fig.add_subplot(2,2,3+i)
    ax.stem(f, f_abs)
    ax.set_title("fs vs |F|")
    ax.set_xlim(xlim)

plt.show()

结果为

即$f=1$和$f=22$处被筛选了出来。

滤波

有了这个，就可以在频域上对数据进行滤波，其思路是，对f_abs中的值进行阈值分割，例如，只筛选出低频部分，然后看一下滤波效果

fig = plt.figure(1)
f_filt = F * (np.abs(f) < 2)
s_filt = fft.ifft(f_filt)
ax = fig.add_subplot()
ax.plot(t, s, lw=0.2)
ax.plot(t, s_filt.real, lw=2)
ax.set_title("threshold=2")
ax.set_xlim([0,10])
plt.show()

效果如下