1. 概要
今度はバンドリジェクションフィルター。係数は、前回までは中心周波数しか考慮していなかったが、今度は窓函数法を用いて帯域幅も考慮に入れて計算する。
2.係数の計算
Pythonでscipy.signal.firwin()函数を利用して下のように計算した。
import scipy.signal
import numpy as np
"""
構文: scipy.signal.firwin(numtaps, cutoff, width=None, window='hamming', pass_zero=True, scale=True, nyq=None, fs=None)
pass_zero : 'bandpass', 'lowpass', 'highpass', 'bandstop'
window : 'boxcar', 'blackman', 'hamming', 'hann', 'flattop', 'bohman', 'blackmanharris', その他
"""
COEFF_WIDTH = 16 #係数のビット数
NUM_COEFFS = 101 #係数の個数(タップ数)。奇数で指定する。
SAMPLE_FREQ = 40_000
CUT_OFF = [1_000, 4_000]
FILTER_TYPE = "bandstop"
WINDOW = "hamming"
# 係数を計算して、
coeffs = scipy.signal.firwin(NUM_COEFFS, CUT_OFF, window=WINDOW, pass_zero=FILTER_TYPE, fs=SAMPLE_FREQ)
# print(*coeffs, sep=',')
# 係数が範囲(-1≦ 係数 < 1)を超えていないことを念のため確かめて、
print("No Good") if (np.min(coeffs) < -1) or (np.max(coeffs) >= 1) else print("IS OK")
# 係数を固定小数点数化(整数化。符号に1ビット、整数部なし、残りを小数ビットに割り当て)して、
coeffs = np.asarray(np.round(coeffs * 2**(COEFF_WIDTH - 1)), dtype=np.int64)
# print(*coeffs, sep=",")
# それが負数であったら2の補数にずらして、
coeffs = np.asarray(np.where(coeffs < 0, coeffs + 2**COEFF_WIDTH, coeffs), dtype=np.uint64)
# print(*coeffs, sep=",")
# array(0 to TAPS - 1)(COEFF_WIDTH - 1 downto 0)の形に体裁を整えて出力する。
print(*[str(COEFF_WIDTH) + 'd"' + str(coeff) + '"'for coeff in coeffs], sep=",")
下のような周波数特性が得られるはずである。http://dsp.jpn.org/dfdesign/fir/mado.shtmlで描画した。
3. 窓函数法によるバンドリジェクションフィルターのためのVHDL
下のVHDLをコンパイルしてMAX10 FPGAに書き込む。
LEを大量に消費した。今回使ったMAX10 (10M08SAE144C8G)ではもう限界である。
係数(signed)は、トップレベルエンティティfir_filter_experiment.vhdの下の部分に記述する。
constant coeffs_list : coefficient_array := (
16d"17",16d"27",16d"34",16d"36",16d"30",16d"17",16d"0",16d"65521",16d"65514",16d"65519",16d"0",16d"21",16d"36",16d"31",16d"0",16d"65479",16d"65411",16d"65351",16d"65323",16d"65340",16d"65398",16d"65475",16d"0",16d"11",16d"65491",16d"65377",16d"65246",16d"65157",16d"65163",16d"65295",16d"0",16d"289",16d"535",16d"654",16d"601",16d"397",16d"137",16d"65495",16d"0",16d"338",16d"948",16d"1689",16d"2331",16d"2615",16d"2342",16d"1437",16d"0",16d"63828",16d"62211",16d"61057",16d"27756",16d"61057",16d"62211",16d"63828",16d"0",16d"1437",16d"2342",16d"2615",16d"2331",16d"1689",16d"948",16d"338",16d"0",16d"65495",16d"137",16d"397",16d"601",16d"654",16d"535",16d"289",16d"0",16d"65295",16d"65163",16d"65157",16d"65246",16d"65377",16d"65491",16d"11",16d"0",16d"65475",16d"65398",16d"65340",16d"65323",16d"65351",16d"65411",16d"65479",16d"0",16d"31",16d"36",16d"21",16d"0",16d"65519",16d"65514",16d"65521",16d"0",16d"17",16d"30",16d"36",16d"34",16d"27",16d"17"
);
constant BITS_NEEDED : natural := 25;
4. 全体の接続図
このときと同じ。
5. 実行結果
ファンクションジェネレーターから正弦波を入力して出力の周波数特性を見てみる。計算どおり阻止帯域1 kHz～4 kHzのバンドリジェクションフィルターになっている。ナイキスト周波数のところはグチャグチャになる。
6. その他の係数計算ツール
係数の計算は、冒頭に示したscipy.signal.firwin()函数のほかに以下のサイトが便利。