Ocean Science Hack: フィルタの応答関数

フィルタの周波数応答関数をプロットしてみる。

以下の例は、京都大学の講義「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()