2016-05-01 5 views
3

私は解析したフーリエ解析を使用した場合のパワースペクトルを3時間ごとの温度データで示しています。Python:フーリエ解析後の時系列フィルタの設計

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

ありませんwは信号が1と2の周波数で最も強いことを知っています。これらの優位な周波数を維持するためにデータを平滑化できるフィルタ(非ボックスカー)を作成したいと思います。

これを行う特定のnumpyまたはscipy関数はありますか?これは主なパッケージの外で作らなければならないものでしょうか?

+0

乗算は、IFFT時間領域に、次に、関心の周波数を分離します。 – roadrunner66

答えて

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=fftfreq(n,dt) # a nice helper function to get the frequencies right 
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') 

enter image description here

関連する問題