「サイエンティフィックPythonのための」IDE、Spyder上にてScientific Python Lecturesの実習中。前回はFFT、今回は信号処理です。前回手抜きだったので今回はしっかりするつもりなのに、お手本はあまり力が入ってないご様子。そこで自主練習を企画。信号処理といったらフィルタだね、という短絡。
※「 ソフトな忘却力」投稿順 Index はこちら
※Scientific Python Lectures様のコースは例題だけでなく、エクササイズなども充実、それを全部順番に解いていったら必ずや立派な人になれるだろ～と思います。でも老い先短い年寄には量が多過ぎて多分死ぬまでに終わりません。適当な練習でお茶を濁してます。今回の「お勉強」は「5.9 Signal processing: scipy.signal」です。
Signal Processing (scipy.signal)
信号処理についてのSciPyのチュートリアルは以下にあります。
https://docs.scipy.org/doc/scipy/tutorial/signal.html
今回は信号処理ならまずフィルタだねという短絡思考。フィルタの中からありがちなFIRフィルタでローパス、ハイパス、バンドパスという具合にフィルタを作ってその特性を眺めてみたいと思います。
FIRフィルタ係数を生成するための関数 firwin と、その応答を求める関数freqz についてのAPIページは以下にあります。
firwin
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html
freqz
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html
上記のページを読んでいてちょっと引っかかったのが正規化周波数の表記方法です。大多数の信号処理の教科書（そしてソフトウエアではScilab）では、サンプリング周波数を1.0、ナイキスト周波数を0.5 とするのが「普通の」正規化周波数の表記法じゃないかと思います。しかし、SciPyはMatlab式みたいです。Matlab式はナイキスト周波数を1.0としている筈。世に隠れも無い偉大なソフトウエアMatlabの流儀なので、これはこれで一つの流派をなしているみたい。でも知らないと戸惑います。昔、Matlab使ったときに引っかかったのを思い出しました。本件については、「hark00」様が以下のページを書かれてます。ありがとうございます。
SciPyはMatlab式なので、教科書的な正規化周波数で考えている場合、その２倍の数値をAPI関数どもに与えないとならないみたい。ホントか？
FIRフィルタの周波数特性３つを表示するスクリプト
以下のスクリプトでは
-
- 内部ではSciPy（Matlab)方式で、ナイキスト周波数を1.0とする表現で計算
- グラフでは、一般的なナイキスト周波数を0.5とする表現で表示
という妥協案？にて処理してます。タップ数３１の
-
- ローパスフィルタ、カットオフ周波数０．１（一般的な意味の）
- ハイパスフィルタ、カットオフ周波数０．４（一般的な意味の）
- バンドパスフィルタ、カットオフ周波数０．１５、０．３５（一般的な意味の）
を生成してその周波数応答のグラフを描くもの。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thr Jul 3 2025
@author: jhalfmoon
SciPy filter sample
"""
import math
import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
def plotResponse(filt, mes):
wT, hT = sp.signal.freqz(filt)
plt.plot(wT/math.pi/2, 20*np.log10(np.abs(hT)), 'b')
plt.title(mes)
plt.ylabel('Amplitude [dB]')
plt.xlabel('Norm. Frequency')
plt.grid()
def main():
lpf1 = sp.signal.firwin(31, 0.1*2, pass_zero='lowpass')
plt.subplot(1,3,1)
plotResponse(lpf1, "LPF: 0.1")
hpf1 = sp.signal.firwin(31, 0.4*2, pass_zero='highpass')
plt.subplot(1,3,2)
plotResponse(hpf1, "HPF: 0.4")
bpf1 = sp.signal.firwin(31, [0.15*2, 0.35*2], pass_zero='bandpass')
plt.subplot(1,3,3)
plotResponse(bpf1, "BPF: 0.15~0.35")
plt.tight_layout()
plt.show()
if __name__ == '__main__':
main()
スクリプトの実行結果
実行結果のグラフが以下に。左からローパス、ハイパス、バンドパスです。
それらしい感じがするので、フィルタは作れるみたいね。次は実際に何かフィルタで処理してみないと。