FPGAでFIRフィルタ

by K.I

2017/06/05〜2017/06/15

Index

概要
- ピン接続
- RaspberryPiのI2S出力を接続
- I2Sのタイムチャート
FIRフィルタの実装
- FIRフィルタの係数は
- FIRフィルタの特性を計算
- 係数を整数化
- FIRフィルタのプログラム
- シミュレーション
- インプリメント出来ない？
動作確認
- Excelでの計算と比較
- FIRフィルタの効果
- MaximaでFIRフィルタを計算
FIRフィルタプログラムの別解
- インプリメントしてみる
まとめ
- 算術右シフト
- 参考

概要

FPGAでFIRフィルタを作ってみよう！

ということで、 FPGAでデジタルフィルタを作ったが、乗算回路は避けていた。

その続きとして、今度は乗算回路を含むFIRフィルタを実装してみようと思う

FPGAは、やっぱりLatticeのMachxo2 Breakoutボードを使用する

I2S入力のデータにデジタルフィルタを掛けて、I2S入力のD級アンプとスピーカーに出力する

設定やモジュール呼出しは省略しているので、詳細は FPGAでデジタルフィルタを見てください。

ピン接続

MachXo2 Breakoutボードのヘッダピン接続一覧

以前にヘッダピンを半田付けした箇所を利用している

J2-1	接続	説明
113	lrclk_i	LRCLK入力
115	bclk_i	BCLK入力
119	data_i	DATA入力
GND
121
125
127

J5-2	接続	説明
61	lrclk_o	LRCLK出力（そのまま）
GND
58	bclk_o	BCLK出力（そのまま）
57	data_o	DATA出力（そのまま）
GND
54	datax	DATA出力（フィルタ処理後）
52

RaspberryPiのI2S出力を接続

Machxo2 Breakoutボードに、 RaspberryPiのI2S出力を入力として接続

出力として、I2S入力のD級アンプとスピーカーを接続する

フィルタを通した音を再生する場合は、出力をdata_o→dataxに接続し直す

I2Sのタイムチャート

RaspberryPiのI2S出力の、L/Rの１データあたりのBCLKは、16bitになっている

LRCLKの変化から、1BCLK分遅らせて、DATAのMSBからBCLKの立ち上がりで読み込む

[top]

FIRフィルタの実装

ここまで乗算器を避けてきたが、どうしても積算が必要になるFIRフィルタを作ってみる。

FIRフィルタの差分方程式は、

まぁ、FIRフィルタのダイアグラムのままプログラムするのが簡単だろう

FIRフィルタの係数は

Machxo2は規模が小さいし、モジュールとしての乗算器が無いので、係数は少ないものにしたい。

ということで、Maximaで FIRフィルタを計算した時に、使用した係数を使うことにする（fir_para.dat）

  4.036409219416838e-03
 -2.208366021329011e-18
 -2.346162858786038e-02
 -3.367761421639100e-02
  7.203267485172524e-02
  2.840739092673040e-01
  4.000000000000000e-01
  2.840739092673040e-01
  7.203267485172524e-02
 -3.367761421639100e-02
 -2.346162858786038e-02
 -2.208366021329011e-18
  4.036409219416838e-03

FIRフィルタの特性を計算

まずはCで記述した FIRフィルタのプログラムで、周波数特性を求めてみる

% pulse 44100 | firf fir_para.dat 16 | fftw | gp_fftl > firf_fft.png
% white 44100 | imul 50 | firf fir_para.dat 16 | fftw | gp_fftl > firf_fftw.png

インパルス応答による周波数特性と、擬似ホワイトノイズによる周波数特性を求めてみた

インパルス応答の方が、キレイに特性が出ている。
- 畳み込み積分なので、擬似ホワイトノイズの誤差が累積されるんじゃないかと思う

係数を整数化

この FIRフィルタのプログラムは内部的に係数を整数化して、整数演算しているので、

インパルス応答で、その整数化した係数を求めることが出来る

$ pulse 13 | firf fir_para.dat 16
132
0
-768
-1103
2360
9308
13107
9308
2360
-1103
-768
0
132

係数をプロットしてみる

$ pulse 13 | firf fir_para.dat 16 | gp_pulse > fir_para.png

この係数で、FIRフィルタを作ってみようと思う

FIRフィルタのプログラム

求めた係数で、FIRのプログラムをVerilogで組んでみる

module fir_filter(
        input clk,
        input rst,
        input enable,
        input [15:0] idata,
        output [15:0] odata
);
        reg signed [15:0] data[0:12];
        reg signed [31:0] sum;
        reg signed [15:0] para[0:12];
        
        initial begin
                para[0] = 132;
                para[1] = 0;
                para[2] = -768;
                para[3] = -1103;
                para[4] = 2360;
                para[5] = 9308;
                para[6] = 13107;
                para[7] = 9308;
                para[8] = 2360;
                para[9] = -1103;
                para[10] = -768;
                para[11] = 0;
                para[12] = 132;
        end

        integer i;

        always @(posedge clk or posedge rst) begin
                if (rst) begin
                        
//                      for (i=0; i<13 ;i=i+1) data[i] <= 0;
                        data[0] <= 0;
                        data[1] <= 0;
                        data[2] <= 0;
                        data[3] <= 0;
                        data[4] <= 0;
                        data[5] <= 0;
                        data[6] <= 0;
                        data[7] <= 0;
                        data[8] <= 0;
                        data[9] <= 0;
                        data[10] <= 0;
                        data[11] <= 0;
                        data[12] <= 0;
                        
                end else if (enable) begin
                        
                        sum <=  data[0]*para[0] + data[1]*para[1] + data[2]*para[2] +
                                data[3]*para[3] + data[4]*para[4] + data[5]*para[5] +
                                data[6]*para[6] + data[7]*para[7] + data[8]*para[8] +
                                data[9]*para[9] + data[10]*para[10] + data[11]*para[11] +
                                data[12]*para[12];

                        data[0] <= idata;

//                      for (i=0; i<12 ;i=i+1) data[i+1] <= data[i];    
                        data[1] <= data[0];
                        data[2] <= data[1];
                        data[3] <= data[2];
                        data[4] <= data[3];
                        data[5] <= data[4];
                        data[6] <= data[5];
                        data[7] <= data[6];
                        data[8] <= data[7];
                        data[9] <= data[8];
                        data[10] <= data[9];
                        data[11] <= data[10];
                        data[12] <= data[11];
                end
        end

        assign odata = sum[31:16];

endmodule

乗算器はIPで作るべきだと思うけど、RTLで、FIRフィルタのブロック図をそのまま記述してみる

でも、足し算の繰り返しがちょっと気持ち悪い¹。
FIRフィルタの構成通りなんだけど、もっとうまく書けないかな

符号付演算は、ちょっと面倒になってきたので、signedを使ってみる

オーバーフローは、、えーと、考えないことにする。。。

畳み込み積分をforでうまく記述できなかったので、単純にベタに記述している

シミュレーションなら動いたんだけど、インプリメント出来なかった²

後から気がついたけど、出力が1LRCLK分遅れている。

最初の係数は、idataから掛けるようにした方が良かったかも。。
- まぁ、結果は変わらないから、良しとする³。。。

シミュレーション

ActiveHDLで、シミュレーションしてみる

FPGAでデジタルフィルタでの、テストベンチのシミュレーション時間を、10倍にして実行

テストベンチのdatalの元データは、こうなっている⁴ →fir_datal.datとする

CによるFIRプログラムに掛けて、RTLのsumの値を計算してみる

% cat fir_datal.dat | firf fir_para.dat 16
0
-3844896
1922316
21409080
17098904
-77290177
-207313618
-246314908
-225317216
-235948468
-247598900
-236195421
-239100806
-245851220
-234273105
-240062030
-245851220
-234273105
-240062030
-245851220
-234273105
-240062030
-245851220
-234273105
-240062030

4byte整数で表せる範囲は、 -2147483648〜2147483647 なので、

とりあえずオーバーフローはしていないみたいだ

RTLのプログラムと同様に、下位16ビットを削る

% cat fir_datal.dat | firf fir_para.dat 16 | rshift 16
0
-59
29
326
260
-1180
-3164
-3759
-3439
-3601
-3779
-3605
-3649
-3752
-3575
-3664
-3752
-3575
-3664
-3752
-3575
-3664
-3752
-3575
-3664

FIRフィルタ後のシミュレーション

とりあえず、計算通りになっていることは確認できた
乗算による遅延を考慮したシミュレーションじゃないから、あまり意味は無いけど。。

インプリメント出来ない？

インプリメントしてみると、以下のエラー

ERROR - Design doesn't fit into device specified, refer to the Map report for more details.

ReportsタブのMapを参照すると、100%以上の項目がある

やはり乗算器はLogicを消費する。1200ZEでは入りきらなかったようだ

        ：
Target Device:  LCMXO2-1200ZETQFP144
        ：
Design Summary
   Number of registers:    525 out of  1604 (33%)
      PFU registers:          525 out of  1280 (41%)
      PIO registers:            0 out of   324 (0%)
   Number of SLICEs:      1402 out of   640 (219%)
      SLICEs as Logic/ROM:   1402 out of   640 (219%)
      SLICEs as RAM:            0 out of   480 (0%)
      SLICEs as Carry:       1352 out of   640 (211%)
   Number of LUT4s:        2804 out of  1280 (219%)
      Number used as logic LUTs:        100
      Number used as distributed RAM:     0
      Number used as ripple logic:      2704
      Number used as shift registers:     0
   Number of PIO sites used: 19 + 4(JTAG) out of 108 (21%)
   Number of block RAMs:  0 out of 7 (0%)
   Number of GSRs:  1 out of 1 (100%)
   EFB used :       No
   JTAG used :      No
   Readback used :  No
   Oscillator used :  No
   Startup used :   No
   POR :            On
   Bandgap :        On
   Number of Power Controller:  0 out of 1 (0%)
   Number of Dynamic Bank Controller (BCINRD):  0 out of 4 (0%)
   Number of Dynamic Bank Controller (BCLVDSO):  0 out of 1 (0%)
   Number of DCCA:  0 out of 8 (0%)
   Number of DCMA:  0 out of 2 (0%)
   Number of PLLs:  0 out of 1 (0%)
   Number of DQSDLLs:  0 out of 2 (0%)
   Number of CLKDIVC:  0 out of 4 (0%)
   Number of ECLKSYNCA:  0 out of 4 (0%)
   Number of ECLKBRIDGECS:  0 out of 2 (0%)
        ：

Project→Device...で、ターゲットデバイスを7000HEに変更してやり直すと、今度はインプリメントできた

ProgrammerのDeviceも変更するのを忘れないように

Target Device:  LCMXO2-7000HETQFP144
        ：
Design Summary
   Number of registers:    664 out of  7209 (9%)
      PFU registers:          664 out of  6864 (10%)
      PIO registers:            0 out of   345 (0%)
   Number of SLICEs:      1403 out of  3432 (41%)
      SLICEs as Logic/ROM:   1403 out of  3432 (41%)
      SLICEs as RAM:            0 out of  2574 (0%)
      SLICEs as Carry:       1349 out of  3432 (39%)
   Number of LUT4s:        2805 out of  6864 (41%)
      Number used as logic LUTs:        107
      Number used as distributed RAM:     0
      Number used as ripple logic:      2698
      Number used as shift registers:     0
   Number of PIO sites used: 19 + 4(JTAG) out of 115 (20%)
   Number of block RAMs:  0 out of 26 (0%)
   Number of GSRs:  1 out of 1 (100%)
   EFB used :       No
   JTAG used :      No
   Readback used :  No
   Oscillator used :  No
   Startup used :   No
   POR :            On
   Bandgap :        On
   Number of Power Controller:  0 out of 1 (0%)
   Number of Dynamic Bank Controller (BCINRD):  0 out of 6 (0%)
   Number of Dynamic Bank Controller (BCLVDSO):  0 out of 1 (0%)
   Number of DCCA:  0 out of 8 (0%)
   Number of DCMA:  0 out of 2 (0%)
   Number of PLLs:  0 out of 2 (0%)
   Number of DQSDLLs:  0 out of 2 (0%)
   Number of CLKDIVC:  0 out of 4 (0%)
   Number of ECLKSYNCA:  0 out of 4 (0%)
   Number of ECLKBRIDGECS:  0 out of 2 (0%)
        ：

今度は、余裕でインプリメント出来ている

¹係数が少ないからまだ良いけど、まともなFIRフィルタ記述したら、とんでもなく長くなってしまう。
²これは、あまり追求していないので、やりようによっては出来るのかも。
³やり直すのが、面倒なので。。
⁴あまり現実的な波形とは言えないが。。

[top]

動作確認

7000HEのボードも持ってるけど、半田付けしないとダメだから使いたくなかったんだけど、

まぁ仕方ないので、ヘッダを同様に半田付けして、配線する

FPGAでデジタルフィルタと同様に、

DOVA-SYNDROMEというサイトの、 Kyaaiさんの「さあでかけよう」という曲を使わせていただく。

ちゃんとFIRフィルタ後の音が出ることも確認出来た。

でも、ちょっと音が小さくなった程度の変化しか感じられないなぁ。。

Excelでの計算と比較

ロジックアナライザのプロトコル解析機能で読んだFIRフィルタの出力データを、

ExcelでFIRの計算をした結果と比較して、グラフにしてみた

Ch	dout	datax	計算値
1	0	0
1	-941	0	0
1	-1009	0	0
1	-865	-2	-2
1	-818	-3	-2
1	-685	9	9
1	-654	26	26
1	-596	-9	-8
1	-658	-148	-147
1	-674	-343	-342
1	-768	-470	-470
1	-814	-475	-474
1	-885	-412	-412
1	-872	-354	-353
1	-856	-321	-321
1	-765	-312	-312
1	-670	-322	-322
1	-524	-348	-347
1	-405	-382	-381
1	-272	-415	-414
1	-118	-438	-437
1	65	-443	-442
1	266	-425	-425
1	462	-387	-387
1	686	-333	-332
1	874	-270	-270
：	：	：	：

FIRフィルタ出力dataxは、ちゃんと計算値とほぼ一致している
乗算器は、あれでちゃんと論理合成出来たんだな〜

FIRフィルタの効果

スピーカーで再生するだけでは、周波数特性の変化が分からないので、

ロジックアナライザのプロトコル解析機能でI2Sデータを取得して、WAVファイルを作成する
- CR削除、3列目抽出、ヘッダ削除、CH1のみ抽出して、WAVファイル作成
```
% tr -d \\r < dfil_fir.csv | cut -d ',' -f 3 | sed -e '1d' | idws 2 > dfil_fir.dat
% txt2wav -44100 dfil_fir.dat
% lame dfil_fir.wav dfil_fir.mp3
```

それから、Windowsのサウンド編集アプリ Audacityでスペクトラム解析して、周波数成分の変化を見てみる

本当は50秒ほどの曲だが、ロジックアナライザであまり長時間データ取れないので、中間の10秒だけとしている

フィルタを通す前の曲データ

FIRフィルタを通した後の曲データ

音を聞いても、音量が小さい以上の違いはちょっと分かりにくいが

スペクトラム解析では、高音が削られているのが分かる

MaximaでFIRフィルタを計算

Maximaでも、 FIRフィルタの特性を計算してみよう

係数を整数化しているので、一応やり直してみた

kill(all)$
a: [
132,
0,
-768,
-1103,
2360,
9308,
13107,
9308,
2360,
-1103,
-768,
0,
132];

sum(a[k+1]*exp(-k*j*\omega*T),k,0,12)/65536;  
subst(%i,j,%);
subst(2*%pi*f,\omega,%);
subst(1/44100,T,%);
wxplot2d(cabs(%),[f,0,22050],[logy],[y,0.01,1]);

10kHzぐらいまで、あまり落ちていないので、変化を感じづらいんだと思う

上のグラフは、Y方向のレンジを、 RCフィルタや移動平均フィルタと合わせているが、

レンジを自動にすると、こんな特性になってることがわかる
```
wxplot2d(cabs(%),[f,0,22050],[logy]);
```

DCでゲインが１になるように、調整した方が良いのかもしれない。

[top]

FIRフィルタプログラムの別解

固定小数点FIRフィルタかいてみたという記事を見たら、ちょっと違う記述方法だった。

気になったので、どんなもんか同様に記述して、試してみることにした

module fir_filter2(
        input clk,
        input rst,
        input enable,
        input signed [15:0] idata,
        output signed [15:0] odata
);
        reg signed [31:0] data[0:12];
        reg signed [15:0] para[0:12];
        
        initial begin
                para[0] = 132;
                para[1] = 0;
                para[2] = -768;
                para[3] = -1103;
                para[4] = 2360;
                para[5] = 9308;
                para[6] = 13107;
                para[7] = 9308;
                para[8] = 2360;
                para[9] = -1103;
                para[10] = -768;
                para[11] = 0;
                para[12] = 132;
        end

        integer i;

        always @(posedge clk or posedge rst) begin
                if (rst) begin
                        
//                      for (i=0; i<13 ;i=i+1) data[i] <= 0;
                        data[0] <= 0;
                        data[1] <= 0;
                        data[2] <= 0;
                        data[3] <= 0;
                        data[4] <= 0;
                        data[5] <= 0;
                        data[6] <= 0;
                        data[7] <= 0;
                        data[8] <= 0;
                        data[9] <= 0;
                        data[10] <= 0;
                        data[11] <= 0;
                        data[12] <= 0;
                        
                end else if (enable) begin
                        
                        data[12] <= idata * para[12];

//                      for (n=0; n<12 ;n=n+1) data[n] <= data[n+1] + idata * para[n];
                        data[11] <= data[12] + idata * para[11];
                        data[10] <= data[11] + idata * para[10];
                        data[9]  <= data[10] + idata * para[9];
                        data[8]  <= data[9]  + idata * para[8];
                        data[7]  <= data[8]  + idata * para[7];
                        data[6]  <= data[7]  + idata * para[6];
                        data[5]  <= data[6]  + idata * para[5];
                        data[4]  <= data[5]  + idata * para[4];
                        data[3]  <= data[4]  + idata * para[3];
                        data[2]  <= data[3]  + idata * para[2];
                        data[1]  <= data[2]  + idata * para[1];
                        data[0]  <= data[1]  + idata * para[0];
                end
        end

        assign odata = data[0] >> 16;

endmodule

プログラムを理解するためにブロックダイアグラムを描いてみると、こんな構成になっているようだ

なるほど。確かにFIRフィルタになってるみたいだ。この構成は思いつかなかった。

この方が、Verilogのコードは感覚的にスッキリしている感じがする

乗算した結果を遅延するので、保持するレジスタ幅は広くなるが、

これなら配線的には、かなり単純になるはずだ

インプリメントしてみる

実際、どっちの方が良いんだろうか。

Target Device:  LCMXO2-7000HETQFP144
        ：
Design Summary
   Number of registers:    768 out of  7209 (11%)
      PFU registers:          768 out of  6864 (11%)
      PIO registers:            0 out of   345 (0%)
   Number of SLICEs:       628 out of  3432 (18%)
      SLICEs as Logic/ROM:    628 out of  3432 (18%)
      SLICEs as RAM:            0 out of  2574 (0%)
      SLICEs as Carry:        583 out of  3432 (17%)
   Number of LUT4s:        1243 out of  6864 (18%)
      Number used as logic LUTs:         77
      Number used as distributed RAM:     0
      Number used as ripple logic:      1166
      Number used as shift registers:     0
   Number of PIO sites used: 19 + 4(JTAG) out of 115 (20%)
   Number of block RAMs:  0 out of 26 (0%)
   Number of GSRs:  1 out of 1 (100%)
   EFB used :       No
   JTAG used :      No
   Readback used :  No
   Oscillator used :  No
   Startup used :   No
   POR :            On
   Bandgap :        On
   Number of Power Controller:  0 out of 1 (0%)
   Number of Dynamic Bank Controller (BCINRD):  0 out of 6 (0%)
   Number of Dynamic Bank Controller (BCLVDSO):  0 out of 1 (0%)
   Number of DCCA:  0 out of 8 (0%)
   Number of DCMA:  0 out of 2 (0%)
   Number of PLLs:  0 out of 2 (0%)
   Number of DQSDLLs:  0 out of 2 (0%)
   Number of CLKDIVC:  0 out of 4 (0%)
   Number of ECLKSYNCA:  0 out of 4 (0%)
   Number of ECLKBRIDGECS:  0 out of 2 (0%)
        ：

インプリメントしてみると、registersが微増しているものの、

SLICEsとLUT4sが、半分以下になっている！

スッキリしているので、少しは減ると思っていたが、こんなに差があるとは思わなかった。

小さな独立した回路の繰り返しみたいな構成にすると、FPGAでは論理合成しやすいのかもしれない

というか、配線が少なくなるようにすると良いのかな

インプリメントしたものを、実際にProgramしてみたけど、ちゃんと動作しました。

結果は同じなんだから、こっちの方が良いのか。。

[top]

まとめ

FIRフィルタはかなり自在な特性を持たせられるが、乗算を繰り返すのでどうしても回路が複雑になる

今回作成したものは係数が少なかったが、乗算回路でロジックをかなり消費している

Machxo2は乗算器が無いこともあり、今回はIPを使わずにRTLで乗算器を記述⁵したが、

乗算器は少ないし、速度的には特に厳しい回路でもないので、問題なく論理合成できた

まぁ乗算器を時分割で使いまわすようにしないと、まともなFIRフィルタは入りそうもない⁶

またFPGAマガジンNo.11の記事のように回路を工夫して乗算を減らす必要があると思う

それにしても、回路の構成による論理合成後のサイズの差は、予想以上にあると感じた。

そこらへんはいろいろとノウハウがありそうだ

FPGAでデジタルフィルタで考えた、クロックを同期せず通過させる、いい加減な回路だけど、

とりあえず、計算した通りにFIRフィルタとしての動作を確認できたので、
まぁ、これはこれで良かったかな。

算術右シフト

PCで計算する時に、掛け算でビット長を短くするのは面倒なので、

算術右シフトで下位ビットを削るためのCプログラム（rshift.c）

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[])
{
        int in, out, n;

        n = 0;
        if (argc>1)
          n = atoi(argv[1]);

        while (scanf("%d", &in) != EOF) {
          out = (n>0)?(in >> n):(in << -n);
          printf("%f\n",out);
        }

        return 0;
}

シフト値をマイナスにすると、算術左シフトになる

gccだと、これで算術シフトになるが、コンパイラによって異なるかもしれないので注意

参考

FPGAでデジタルフィルタ

この記事と同じ構成で、RCフィルタと移動平均フィルタを作ってみた

FIRフィルタでサンプルレートコンバータ

Cで作成したFIRフィルタのプログラム

MaximaでFIRフィルタを計算

Maximaでこの記事と同じ係数のFIRフィルタの計算を試した

固定小数点FIRフィルタかいてみた

ちょっと違う構成でFIRフィルタを記述しているので参考になった

⁵RTLで記述すると、やっぱり楽で良いなぁ。
⁶そもそもこのFPGAボードでは大きな回路は無理なんだけど、大規模なFPGAボードでも、そんなに余裕があるわけじゃない。

[top] [電子工作関連に戻る]