时间序列学习经典案例（3）离散傅⾥叶变换DFT（案例：时序去噪）
1.傅⾥叶定理
法国科学家傅⾥叶提出，任何⼀条周期性曲线，⽆论多么跳跃或不规则，都能表⽰成⼀组光滑正弦曲线叠加之和。
2.离散傅⾥叶变换
离散傅⾥叶变换即是把 ⼀条周期性曲线 拆解成 ⼀组光滑正弦曲线 的过程。在离散傅⾥叶级数变换中，“频率”⽤“级数长度的分数”来表⽰。
【以下傅⾥叶变换均指离散傅⾥叶变换】
离散傅⾥叶变换的⽬的是可将时域（即时间域）上的信号转变为频域（即频率域）上的信号，随着域的不同， 对同⼀个事物的认知⾓度也随之改变，因此在时域中某些不好处理的地⽅，在频域就可以较为简单地处理。这就可以⼤量减少处理信号存储量。把信号从时间域变换到频率域，进⽽研究信号的频谱结构和变化规律。
假设有⼀时间域函数：y = f(x)，根据傅⾥叶的理论它可以被分解为⼀系列正弦函数的叠加，他们的振幅A，频率ω或初相位φ不同：
　所以离散傅⾥叶变换可以把⼀个⽐较复杂的函数转换为多个简单函数的叠加，看问题的⾓度也从时间域转到了频率域，有些问题处理起来就会⽐较简单。
np.fft和scipy.fftpack中的实现有所不同：NumPy 提供规范化标志，⽽fftpack没有。但请注意：归⼀化不会产⽣真实的振幅。正确的幅度可以通过乘以 1 /  . 但是，这种归⼀化⾜以⽐较不同长度的时间序列。
numpy中为了通过irfft对时域进⾏正确的重采样，rfft和irfft中的 归⼀化标志norm 必须相等(None或“ortho”)。
3.相关定理
（1）采样频率
采样频率，也称为采样速度或者采样率，定义了每秒从连续信号中提取并组成离散信号的采样个数，它⽤赫兹（Hz）来表⽰。采样频率的倒数是采样周期或者叫作采样时间，它是采样之间的时间间隔。通俗的讲采样频率是指计算机每秒钟采集多少个信号样本。
（2）采样定理
所谓采样定理 ，⼜称⾹农采样定理，奈奎斯特采样定理，是信息论，特别是通讯与信号处理学科中的⼀个重要基本结论。采样定理指出，如果信号是带限的，并且采样频率⾼于信号带宽的两倍，那么，原来的连续信号可以从采样样本中完全重建出来。
定理的具体表述为：在进⾏模拟/数字信号的转换过程中，当采样频率fs ⼤于 信号中最⾼频率fmax的2倍时，即 fs>2*fmax
Note：采样之后的数字信号完整地保留了原始信号中的信息，⼀般实际应⽤中保证采样频率为信号最⾼频率的2.56～4倍；采样定理⼜称奈奎斯特定理。
4.numpy进⾏傅⾥叶（变换）分析
过程简述：
（1）导⼊快速傅⾥叶变换所需模块
import numpy.fft as nf
（2）通过原函数值得序列经过快速傅⾥叶变换得到⼀个复数数组，复数的模代表的是振幅，复数的辐⾓代表初相位
nf.fft(原函数序列值)->⽬标函数序列值（复数数组）
计算⼀维离散傅⾥叶变换。此函数使⽤⾼效的快速傅⾥叶变换(FFT)算法[CT]计算⼀维n点离散傅⾥叶变换(DFT)。
（3）通过⼀个复数数组（复数的模代表振幅，复数的辐⾓代表初相位）经过逆傅⾥叶变换得到合成的函数值数组。nf.ifft(⽬标函数值序列（复数）)->原函数值序列
计算⼀维离散傅⾥叶逆变换。此函数计算由计算的⼀维n点离散傅⾥叶变换的逆fft。
（4）通过采样数域采样周期求得傅⾥叶变换分解所得曲线的频率序列
freqs = nf.fftfreq(采样熟练，采样周期)
（1）获取数据
import numpy as np
import matplotlib.pyplot as plt
import numpy.fft as nf
x = np.linspace(-2* np.pi,2* np.pi,1000)
y = np.zeros(x.size)
for i in range(1,1000):
y +=4* np.pi /(2* i -1)* np.sin((2* i -1)* x)
（2）获取复制数组、计算傅⾥叶变换获取复数实部
# 针对⽅波y做fft，获取复制数组
compy_arr = nf.fft(y)
# 计算傅⾥叶变换，并取复数数组的实部
y2 = nf.ifft(compy_arr).real
（3）采⽤获取频率序列
# 采样数域采样周期求得傅⾥叶变换分解所得曲线的频率序列
freqs = nf.fftfreq(y.size, x[1]- x[0])
linspace numpyfreqs
到此，freqs 就是所求的 傅⾥叶变换结果。
（4）完整代码
import matplotlib.pyplot as plt
import numpy.fft as nf
import matplotlib
<("font",family='YouYuan')
x = np.linspace(-2* np.pi,2* np.pi,1000)
y = np.zeros(x.size)
for i in range(1,1000):
y +=4* np.pi /(2* i -1)* np.sin((2* i -1)* x)
# 针对⽅波y做fft，获取复制数组
compy_arr = nf.fft(y)
# 计算傅⾥叶变换，并取复数数组的实部
y2 = nf.ifft(compy_arr).real
# 采样数域采样周期求得傅⾥叶变换分解所得曲线的频率序列
freqs = nf.fftfreq(y.size, x[1]- x[0])
# 可视化原始数据---------------------------
plt.figure('FFT', facecolor='lightgray')
plt.subplot(121)
plt.title('时域', fontsize=16)
plt.plot(x, y, label=r'$y$')
# 可视化傅⾥叶变换后的实部----------------
plt.plot(x, y2, color='orangered', linewidth=5, alpha=0.5, label=r'$y$')
# 可视化频域图形------------------------
plt.subplot(122)
pows = np.abs(comp_arr)# 复数的模
plt.title('频域', fontsize=16)
plt.plot(freqs[freqs >0], pows[freqs >0], color='orangered', label='frequency')
# 打印结果
plt.legend()
plt.savefig('fft.png')
plt.show()
5.scipy进⾏傅⾥叶(变换)分析
numpy中有⼀个fft的库，scipy中也有⼀个fftpack的库，各⾃都有fft函数，两者的⽤法基本是⼀致的（不完全，还是有⼀定区别的）。
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl
#采样点选择1400个，因为设置的信号频率分量最⾼为600赫兹，根据采样定理知采样频率要⼤于信号频率2倍，所以这⾥设置采样频率为1400赫兹（即⼀秒内有1400个采样点，⼀样意思的）
x=np.linspace(0,1,1400)
#设置需要采样的信号，频率分量有200，400和600
y=7*np.sin(2*np.pi*200*x)+5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x)
# np.linspace(-2 * np.pi, 2 * np.pi, 1400)
# y = np.zeros(x.size)
# for i in range(1, 1000):
#    y += 4 * np.pi / (2 * i - 1) * np.sin((2 * i - 1) * x)
fft_y=fft(y)#快速傅⾥叶变换
N=1400
x = np.arange(N)# 频率个数
half_x = x[range(int(N/2))]#取⼀半区间
abs_y=np.abs(fft_y)# 取复数的绝对值，即复数的模(双边频谱)
angle_y=np.angle(fft_y)#取复数的⾓度
normalization_y=abs_y/N            #归⼀化处理（双边频谱）
normalization_half_y = normalization_y[range(int(N/2))]#由于对称性，只取⼀半区间（单边频谱）
plt.subplot(231)
plt.plot(x,y)
plt.title('原始波形')
plt.subplot(232)
plt.plot(x,fft_y,'black')
plt.title('双边振幅谱(未求振幅绝对值)',fontsize=9,color='black')
plt.subplot(233)
plt.plot(x,abs_y,'r')
plt.title('双边振幅谱(未归⼀化)',fontsize=9,color='red')
plt.subplot(234)
plt.plot(x,angle_y,'violet')
plt.title('双边相位谱(未归⼀化)',fontsize=9,color='violet')
plt.subplot(235)
plt.plot(x,normalization_y,'g')
plt.title('双边振幅谱(归⼀化)',fontsize=9,color='green')
plt.subplot(236)
plt.plot(half_x,normalization_half_y,'blue')
plt.title('单边振幅谱(归⼀化)',fontsize=9,color='blue')
plt.show()
6.使⽤傅⽴叶变换清理时间序列数据噪声
将时域转成频域，然后处理频域中的数据，例如：去除噪声波。之后，可以使⽤逆向傅⾥叶 将频域数据转换回时域波。
（1）获取数据
创建两个正弦波并将它们合并为⼀个正弦波，然后故意⽤ np.random.randn(len(t)) ⽣成的数据污染⼲净的波。(将两个信号组合成第三个信号也称为卷积或信号卷积。)
import numpy as np
import matplotlib.pyplot as plt
#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 轴将是波频率）。
（2）从时域到频域的转换
可以使⽤ 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))

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。