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

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)

Original Data:

Frequency Spectrum:

Power Spectrum:

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.

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?

解决方案

An example with some synthetic data:

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