傅里叶变换简单解释
傅里叶变换把无限长的三角函数作为基函数,这个基函数会伸缩,平移(其实是两个正交基的分解)。缩得窄,对应高频;伸得宽,对应低频。然后这个基函数不断和信号做相乘。某一个尺度(宽窄)下乘出来的结果,就可以理解成信号所包含的当前尺度对应频率成分有多少。于是,基函数会在某些尺度下,与信号相乘得到一个很大的值,因为此时二者有一种重合关系。那么我们就知道信号包含该频率的成分的多少。
使用傅立叶变换清理时间序列数据噪声代码展示
将干净的数据与噪声混合
创建两个正弦波并将它们合并为一个正弦波,然后故意用 np.random.randn(len(t)) 生成的数据污染干净的波。
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [16,10]
plt.rcParams.update({'font.size':18})
#Create a simple signal with two frequencies
data_step = 0.001
t = np.arange(start=0,stop=1,step=data_step)
f_clean = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
f_noise = f_clean + 2.5*np.random.randn(len(t))
plt.plot(t,f_noise,color='c',Linewidth=1.5,label='Noisy')
plt.plot(t,f_clean,color='k',Linewidth=2,label='Clean')
plt.legend()
(将两个信号组合成第三个信号也称为卷积或信号卷积。这是另一个有趣的话题,现在在神经网络模型中被大量使用)带有带噪音的波浪,黑色是我们想要的波浪,绿线是噪音。
如果我隐藏图表中的颜色,我们几乎无法将噪声从干净的数据中分离出来,但是 傅立叶变换在这里可以提供帮助。我们需要做的就是将数据转换到另一个角度,从时间视图(x 轴)到频率视图(x 轴将是波频率)。
从时域到频域的转换
这里可以使用 numpy.fft 或 scipy.fft(pytorch1.8以后也增加了torch.fft这里就不详细说了)。我发现 scipy.fft 非常方便且功能齐全,所以在本文中使用 scipy.fft,但是如果想使用其他模块或者根据公式构建自己的一个也是没问题的(代码见最后)。
from scipy.fft import rfft,rfftfreq
n = len(t)
yf = rfft(f_noise)
xf = rfftfreq(n,data_step)
plt.plot(xf,np.abs(yf))
在代码中,我使用 rfft 而不是 fft。r 意味着reduce(我认为)只计算正频率。所有负镜像频率将被省略。因为他的速度更快。rfft 函数的 yf 结果是一个复数,形式类似于 a+bj。np.abs() 函数将为复数计算 √(a² + b²)。
这是我们原始波的神奇的频域视图。x轴表示频率。
一些在时域看起来很复杂的东西现在被转换成非常简单的频域数据。这两个峰代表两个正弦波的频率。一种波是50Hz,另一种是120Hz。再回顾一下生成正弦波的代码。
f_clean = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
其他频率就是噪音,并且在下一个步骤中很容易去除。
去除噪声频率
在Numpy的帮助下,我们可以很容易地将这些频率数据设置为0,除了50Hz和120Hz。
yf_abs = np.abs(yf)
indices = yf_abs>300 # filter out those value under 300
yf_clean = indices * yf # noise frequency will be set to 0
plt.plot(xf,np.abs(yf_clean))
现在,所有的噪音都被清除了。
回到时域数据
最后一步就是逆转换了,将频域数据转换回时域数据
from scipy.fft import irfft
new_f_clean = irfft(yf_clean)
plt.plot(t,new_f_clean)
plt.ylim(-6,8)
结果表明,所有的噪声波都被去除了。