VHDLによる簡単なディジタルフィルターの実験 / その7 of 7 / FIRフィルター / 窓函数法によるバンドリジェクションフィルター



#1. 概要
今度はバンドリジェクションフィルター。係数は、前回までは中心周波数しか考慮していなかったが、今度は窓函数法を用いて帯域幅も考慮に入れて計算する。

#2.係数の計算
Pythonで`scipy.signal.firwin()`函数を利用して下のように計算した。

```python
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](http://dsp.jpn.org/dfdesign/fir/mado.shtml)で描画した。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193936/d2c8cebb-fd66-5da7-824a-9d8e9609a06a.png)


#3. 窓函数法によるバンドリジェクションフィルターのためのVHDL

下のVHDLをコンパイルしてMAX10 FPGAに書き込む。

ファイル一式: https://github.com/ti-nspire/VHDL_for_Quartus_Prime/tree/main/fir_filter_experiment_band_rejection_hamming

LEを大量に消費した。今回使ったMAX10 \(10M08SAE144C8G\)ではもう限界である。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193936/f4bc82e7-aff6-3283-9de0-12824fb38389.png)

係数(`signed`)は、トップレベルエンティティ`fir_filter_experiment.vhd`の下の部分に記述する。

```vhdl
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. 全体の接続図
[このとき](https://qiita.com/cyq04000/items/c6b6f871ecd8e679bbb6#3-%E5%85%A8%E4%BD%93%E3%81%AE%E6%8E%A5%E7%B6%9A%E5%9B%B3)と同じ。


#5. 実行結果
ファンクションジェネレーターから正弦波を入力して出力の周波数特性を見てみる。計算どおり阻止帯域1 kHz～4 kHzのバンドリジェクションフィルターになっている。ナイキスト周波数のところはグチャグチャになる。
![scope_66.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/193936/84b8ca5d-00ff-f402-b164-c2119a937f35.png)



#6. その他の係数計算ツール
係数の計算は、冒頭に示した`scipy.signal.firwin()`函数のほかに以下のサイトが便利。

- [石川高専 山田洋士 研究室ホームページ「FIRフィルタの窓関数法による設計」](http://dsp.jpn.org/dfdesign/fir/mado.shtml)
- [HS SOFT「ディジタルフィルタ設計サイト」](http://hs-soft.exout.net/fir.php)