signal.firwin関数によるFIRフィルタ係数の作成

scipyライブラリを使用して、signal.firwin関数によるFIRフィルタ係数を作成します。FIRフィルタ係数はファイル化して、「C#言語によるFIRディジタルフィルタの作成」でFIRフィルタ係数として入力し、wavファイルが作成できるようにします。

FIRフィルタ係数作成プログラムの作成

signal.firwin関数を使って、ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタのFIRフィルタ係数を作成し、signal.freqz関数を使って、それぞれのフィルタ係数の周波数領域を表示します。なお作成したFIRフィルタ係数はファイルに保存します。なお、カットオフ周波数は750Hz,バンド幅は750Hz-3kHz、タップ数は255、サンプリング周波数は44.1ｋHzとします。

firwintest.py

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

fs = 44100
nyq = fs / 2.0  # ナイキスト周波数

# フィルタの設計
# ナイキスト周波数が1になるように正規化
fe1 = 750.0 / nyq       # カットオフ周波数1
fe2 = 3000.0 / nyq      # カットオフ周波数2
numtaps = 255           # フィルタ係数（タップ）の数（要奇数）

# b = signal.firwin(numtaps, fe1)                         # Low-pass
#b = signal.firwin(numtaps, fe2, pass_zero=False)        # High-pass
#b = signal.firwin(numtaps, [fe1, fe2], pass_zero=False) # Band-pass
b = signal.firwin(numtaps, [fe1, fe2])                  # Band-stop

# ファイル係数の保存
np.savetxt("parmBS750-3000", b, fmt='% 7f')

w, h = signal.freqz(b)

plt.figure(figsize=(12, 5))  # (width, height)
plt.title('Digital filter frequency response')
plt.xscale("log")

a = 20 * np.log10(abs(h));
f = w / (2 * np.pi) * fs
plt.plot(f, a  )
plt.title('Digital filter frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.xlim(10, 44100)
plt.ylim([-15,15])
plt.grid(which='both', axis='both')
plt.show()

FIRフィルタ係数作成プログラムの実行

次のコマンドでFIRフィルタ係数作成プログラムを実行します。現状、ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタのsignal.firwin関数のコメントをプログラムの実行ごとにそれぞれ外し、作成したFIRフィルタ係数の保存ファイル名を「parmLP750」、「parmHP750」、「parmBP750-3000」、「parmBS750-3000」とします。

$ firwintest.py

ローパスフィルタ、ハイパスフィルタ、バンドパスフィルタ、バンドストップフィルタのFIRフィルタ係数の内容とフィルタ係数の周波数領域を次に示します。

ローパスフィルタ：フィルタ係数「parmLP750」

ハイパスフィルタ：フィルタ係数「parmHP750」

バンドパスフィルタ：フィルタ係数「parmBP750-3000」

バンドストップフィルタ：フィルタ係数「parmBS750-3000」