フィルタの周波数応答関数をプロットしてみる。
以下の例は、京都大学の講義「MATLABで時系列解析」の線形システムの伝達関数を求めるを翻案したものである。
想定としては、1ヶ月ずつのデータがあり、シグナルとしては、48ヶ月周期と6ヶ月周期とホワイトノイズが重なったシグナルに対して13ヶ月の移動平均をかけている。
理論的な周波数応答関数を求める関数はscipyのfreqz関数である(matlabのfreqz関数に対応する)。応答関数は入力 x のパワースペクトル密度関数 Pxx と出力のクロススペクトル密度関数 Pyx の比としても推定できる。スペクトル密度関数を求める関数psdに関しては何度か紹介した。
#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as plab
import matplotlib.mlab as mlab
from scipy.signal import freqz,lfilter
#--- Signal
Fs = 1 # frequency
T = 1000. # final time
t = np.arange(0.,T+1./Fs,1./Fs) # time
x = np.sin(2.*np.pi*t/6.) + np.sin(2.*np.pi*t/48.)+np.random.normal(0,0.1,len(t)) # input signal
#---- power spectral by Welch's method
plt.figure()
nfft = 256
window = np.hanning(256)
noverlap = 128
[Pxx,f] = mlab.psd(x,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # Power spectral
plt.subplot(2,2,1)
plt.plot(f,10*np.log10(Pxx)) # plot in decibel
plt.legend(('Pxx',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))
#--- Moving average
h = np.ones(13)/13.
#--- Calculate theoritical response
fs,Res = freqz(h)
plt.subplot(2,2,2)
plt.plot(fs/2./np.pi,np.abs(Res))
plt.legend(('theoretical response',))
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim((0,1.))
#--- Apply Filter
y = lfilter(h,1,x) # filter
[Pxy,f] = mlab.csd(x,y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # cross scectral density
[Pyy,f] = mlab.psd(y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # power spectral density
plt.subplot(2,2,3)
plt.plot(f,10.*np.log10(abs(Pyy)))
plt.legend(('Pyy',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))
#--- Calculate estimated response
Resc = Pxy/Pxx
plt.subplot(2,2,4)
plt.plot(f,abs(Resc));
plt.legend(('estimated response',))
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim([0,1])
#--- Save fig
plt.savefig("test_response.png")
plt.show()
以下がフィルタを12ヶ月移動平均にした場合。
参考） Scipy Cookbook FIRFilter
