【Python】通过Fourier变换实现频域滤波

fftpack将被弃用，推荐使用 scipy.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}$$

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

• i 表示逆变换；
• r 表示实数域，如开头无r表示在复数域

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

实操

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

import numpy as np
from scipy import fftpack

PI = np.pi*2
fs = 60     #采样频率
T = 100     #采样周期数
N = fs*T    #采样点
t = np.linspace(0, T, N)
s = 2 * np.sin(PI * t) + 3 * np.sin(22 * PI * t) + 2 * np.random.randn(*t.shape)
F = fftpack.fft(s)
f = fftpack.fftfreq(N, 1.0/fs)
1
2
3
4
5
6
7
8
9
10
11

其中，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|")
ax = fig.add_subplot(2,2,3)
ax.stem(f, f_abs)
ax.set_xlim([0,2])
ax.set_title("fs vs |F|")
ax = fig.add_subplot(2,2,4)
ax.stem(f, f_abs)
ax.set_title("fs vs |F|")
ax.set_xlim([21,23])
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

结果为

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

频域滤波

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

fig = plt.figure(1)
f_filt = F * (np.abs(f) < 2)
s_filt = fftpack.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()
1
2
3
4
5
6
7
8
9

效果如下