Python:在傅立叶分析后设计时间序列滤波器 [英] Python: Designing a time-series filter after Fourier analysis
本文介绍了Python:在傅立叶分析后设计时间序列滤波器的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
我有一个时间序列的 3 小时温度数据,我分析并找到了使用傅立叶分析的功率谱.
I have a time series of 3-hourly temperature data that I have analyzed and found the power spectrum for using Fourier analysis.
data = np.genfromtxt('H:/RData/3hr_obs.txt',
skip_header=3)
step = data[:,0]
t = data[:,1]
y = data[:,2]
freq = 0.125
yps = np.abs(np.fft.fft(y))**2
yfreqs = np.fft.fftfreq(y.size, freq)
y_idx = np.argsort(yfreqs)
fig = plt.figure(figsize=(14,10))
ax = fig.add_subplot(111)
ax.semilogy(yfreqs[y_idx],yps[y_idx])
ax.set_ylim(1e-3,1e8)
原始数据:
频谱:
功率谱:
既然我知道信号在频率 1 和 2 处最强,我想创建一个滤波器(非 Boxcar)来平滑数据以保留这些主要频率.
Now that I know that the signal is strongest at frequencies of 1 and 2, I want to create a filter (non-boxcar) that can smooth out the data to keep those dominant frequencies.
是否有特定的 numpy 或 scipy 函数可以做到这一点?这是否必须在主包之外创建?
Is there a specific numpy or scipy function that can do this? Will this be something that will have to be created outside the main packages?
推荐答案
一些合成数据的例子:
# fourier filter example (1D)
%matplotlib inline
import matplotlib.pyplot as p
import numpy as np
# make up a noisy signal
dt=0.01
t= np.arange(0,5,dt)
f1,f2= 5, 20 #Hz
n=t.size
s0= 0.2*np.sin(2*np.pi*f1*t)+ 0.15 * np.sin(2*np.pi*f2*t)
sr= np.random.rand(np.size(t))
s=s0+sr
#fft
s-= s.mean() # remove DC (spectrum easier to look at)
fr=np.fft.fftfreq(n,dt) # a nice helper function to get the frequencies
fou=np.fft.fft(s)
#make up a narrow bandpass with a Gaussian
df=0.1
gpl= np.exp(- ((fr-f1)/(2*df))**2)+ np.exp(- ((fr-f2)/(2*df))**2) # pos. frequencies
gmn= np.exp(- ((fr+f1)/(2*df))**2)+ np.exp(- ((fr+f2)/(2*df))**2) # neg. frequencies
g=gpl+gmn
filt=fou*g #filtered spectrum = spectrum * bandpass
#ifft
s2=np.fft.ifft(filt)
p.figure(figsize=(12,8))
p.subplot(511)
p.plot(t,s0)
p.title('data w/o noise')
p.subplot(512)
p.plot(t,s)
p.title('data w/ noise')
p.subplot(513)
p.plot(np.fft.fftshift(fr) ,np.fft.fftshift(np.abs(fou) ) )
p.title('spectrum of noisy data')
p.subplot(514)
p.plot(fr,g*50, 'r')
p.plot(fr,np.abs(filt))
p.title('filter (red) + filtered spectrum')
p.subplot(515)
p.plot(t,np.real(s2))
p.title('filtered time data')
这篇关于Python:在傅立叶分析后设计时间序列滤波器的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文