Matlabによる信号処理【フィルタの作り方】

DSP(Digital Signal Processing)でフィルタを使う訳ですが、今の知識だとFIRとIIRどちらを選べばいいのかわからないのでお勉強します。

【記事1】のNI様のところに分かりやすく書いてあります。

【記事1】
IIRフィルタとFIRフィルタ - DIAdem 2012ヘルプ - National Instruments

とりあえず、今後の基準は以下の通りと考えておきます。

■IIR
メリット: 係数が少ない/計算ボリュームが小さい
デメリット: 非線形位相応答
得意とする情報: 振幅

■FIR
メリット: 非線形位相応答
デメリット: 係数が多い
得意とする情報: 位相

こんなところです。

Matlabのドキュメンテーションに【記事2】があるので、基本的には記事にそって進めていきます。

【記事2】
jp.mathworks.com

まずは信号を設定します。信号はMatlabのサンプルデータとして用意されているchirpを使います。信号は以下の通りです。

f:id:cocosuzu:20170304124209p:plain

さて、フィルタを作ります。

blo = fir1(34,0.48,chebwin(35,30));

がフィルタの係数を生成する関数です。fir1はウィンドウを指定してフィルタの係数を生成する関数で、引数の最後にはウィンドウを定義するchebwinがいます。この辺はひとまずwindow関数が選べますよという程度にとどめておきます。

続いてフィルタを適用するためにfilter関数に処理したいデータとフィルタの係数セットを渡します。

すると時間波形は以下の通りとなります。

f:id:cocosuzu:20170305133126p:plain

大きな信号が消えました。ぱっと見で変化は分かりますが、より分かりやすくするために周波数スペクトルを見てみましょう。

f:id:cocosuzu:20170305164639p:plain

上段がフィルタ無しの生データ、下段がフィルタを適用した後の信号です。2000Hz以上の主ピークたちがごっそりいなくなっています。

これでフィルタが適用できていることがよく分かります。

ここで一歩戻って今回どんなフィルタを作ったのか見てみましょう。freqz関数でフィルタ係数の特性を見ることが出来ます。先ほど作ったbloを見てみると…

freqz(blo,1)

f:id:cocosuzu:20170305164351p:plain

どこにも2000Hzとは書いてません。正規化周波数で書かれています。ただし、ここで言う正規化周波数は書籍などでよく見るサンプリング周波数を1とする正規化周波数ではなく、ナイキスト周波数を1としています。それで正規化周波数の単位には2 $\pi$ rad/Samplingではなく $\pi$ rad/Samplingと書かれているのですね。

サンプリング周波数は8192Hzであることから、その半分である4096Hzを1とするとfir関数で引数としていた0.48は1966Hzに相当します。この1966Hzがローパスフィルタの周波数となっているわけですね。

さて、大体【記事2】の内容が分かってきたところでフィルタを自由自在に作る方法に移っていきます。

ここで使うのがFDA(Filter Design and Analisit) toolです。これはフィルタの形状を見ながらちょちょいとフィルタを作るには便利なツールです。

先ずは、ここまでで扱ってきたchirpデータの周波数スペクトルを見直してみます。

f:id:cocosuzu:20170305170509p:plain

バンドパスフィルタで主ピークの両端を削るようなフィルタを作ってみましょう。となるとバンドパスフィルタの周波数を2700～3500Hzにします。

f:id:cocosuzu:20170305170750p:plain

こんな感じの画面が出てきますので、以下のように設定します。細かい数値は適当です。

f:id:cocosuzu:20170305171614p:plain

エクスポートしてワークスペースに係数ベクトルをNumとして出力します。これを先のbloの代わりに代入してあげればバンドパスフィルタとして処理されるはずです。

ただ、メインのスクリプトがclear allで変数を都度リセットするようにしているので、都度Numを生成してくれるように作成したフィルタをmファイルとして出力しておきます。Matlabコードを生成を選ぶとmファイルを作ってくれます。

% All frequency values are in Hz.
Fs = 8192;  % Sampling Frequency

Fstop1 = 2500;            % First Stopband Frequency
Fpass1 = 2700;            % First Passband Frequency
Fpass2 = 3500;            % Second Passband Frequency
Fstop2 = 3700;            % Second Stopband Frequency
Dstop1 = 0.001;           % First Stopband Attenuation
Dpass  = 0.057501127785;  % Passband Ripple
Dstop2 = 0.0001;          % Second Stopband Attenuation
dens   = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 1 ...
                          0], [Dstop1 Dpass Dstop2]);

% Calculate the coefficients using the FIRPM function.
bpf  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);

ここで、Fsはスクリプト内で共通なので削除しておいて、生成されたフィルタ係数をfreqzで見てみましょう。

f:id:cocosuzu:20170305172558p:plain