Python：在傅立叶分析之后设计时间序列过滤器

Question

我已经分析了一个每小时3小时的温度数据的时间序列，并发现了使用傅里叶分析的功率谱。

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频率处最强，我想创建一个滤波器（非棚车），该滤波器可以使数据平滑以保持那些主导频率。

是否有特定的numpy或scipy函数可以做到这一点？这将是必须在主软件包之外创建的东西吗？

Answer 1

一些合成数据的例子：

# 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')