scipyを使ったFIRフィルタの仕方について

Question

お世話になります。

scipyを使ったFIRフィルタの仕方について質問させていただきます。teratailでもいくつかFIRフィルタについて質問がありましたので、質問内容と以下のホームページ（https://momonoki2017.blogspot.com/2018/03/pythonfft-1-fft.html）を参考にしてプログラムを組んでみました。

以下ソースコードです。

これは、60Hzと120Hzを合計した信号をFIRフィルタ（バンドパスフィルタ。カットオフ周波数は100～130Hz）を通す処理を行っております。私自身の理解として、フィルタを通した後は120Hzの周波数特性だけが抜き出されるのかと思いましたが、私のイメージとは異なります。この波形は正解なのでしょうか。仕方が間違っている、もしくは理解が異なっておりましたらご指摘願います。

また、本サイトの質問の主旨から少し離れるかもしれませんが、タップ数についてあまり理解できていないので、それにつきましても教えていただきたく思います。今回はデータ数分のタップ数を設けておりますが、それで問題ないでしょうか。タップ数の決め方も教えていただければ幸いです。

よろしくお願いいたします。

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt

N=10.00                                       #データ最大値
samplingF=0.0001                              #サンプリング周波数
NN=N/samplingF
data=np.arange(0,N,samplingF) 
y=np.sin(2*np.pi*60*data)+np.sin(2*np.pi*120*data)   #60Hz+120Hz信号
fq = np.linspace(0, 1/samplingF,NN) 
set=y                                         #フーリエ変換する要素
numf=abs(np.fft.fft(set)*2/NN)                #フィルタ前のフーリエ変換
nyq=1/(samplingF*2)                           #ナイキスト周波数
fe1=100/nyq                                    #カットオフ周波数１
fe2=130/nyq                                   #カットオフ周波数２
numtaps =int(NN)                              #タップ数
bandpass=scipy.signal.firwin(numtaps,[fe1,fe2],pass_zero=False) #バンドパスフィルタ
passdata=scipy.signal.lfilter(bandpass,1,set) #フィルタ処理
numff=abs(np.fft.fft(passdata)*2/NN)          #フィルタ後のフーリエ変換


plt.plot(fq,numf)  #フィルタ前の波形
plt.plot(fq,numff) #フィルタ後の波形
plt.xlim(50,200)
plt.xlabel('freqency(Hz)', fontsize=14)
plt.ylabel('signal amplitude', fontsize=14)
plt.show()

Answer 1

こんにちは。
ご要望が多いため、完全な答えをご提供できませんが、下記が解決のヒントになれば幸いです。

結論から申し上げますと、フィルタ設計は適切にできていると思われます。

設計されたFIRフィルタを scipy.signal.freqz で分析すると下図のようになります。
正規化後のカットオフ周波数は fe1, fe2 = 0.05, 0.065 ですので、たしかにバンドパスフィルタとして機能しているように見えます。

下記のコードでこの内容が再現できます。

plt.figure()
plt.plot(w/np.pi, 20*np.log10(abs(h)))
plt.vlines(x=[fe1, fe2], ymin=-20, ymax=0, linestyle="dashed", colors='grey')
plt.xlim(0, fe2*1.25)
plt.ylim(-20,0)
plt.show()

フィルタは正しく設計できているので、プロットする元のデータやFFTの内容を見直してみてはいかがでしょうか。

Comments

[sham0909]

回答ありがとうございました。
データを見直して再検討したいと思います。