双一次z変換でバターワース特性4次のハイパスフィルタを作る。
デジタルフィルタには色々ある
(1)IIRフィルタ
(2)FIRフィルタ
(3)FFTフィルタ
主にどんなものか(間違っている可能性有り)
IIRフィルタは主に実際の機器に使われている(と思う)。抵抗やコンデンサの組み合わせで実現し、プログラムで実現する場合はそれらの特性を計算して行う。少ない計算量でFIRより減衰を得られるが、フィードバックするから誤差が蓄積するとか、線形位相では無いという欠点があるものの、少なくとも計測の世界では一昔前の分析器やフィルタ機はIIRフィルタだった。近年の分析器はデジタル化してプログラムで行うので、何のフィルターを使っているか知らない。Matlabの関数のfiltfitのような手法を使えば、IIRフィルタでも位相ひずみを補正できる。一旦フィルタを通した波形のデータを逆順にして再度フィルタを掛けるという手法。
騒音の分析ではA特性感覚補正をIIRフィルタの組み合わせで出来る。また、時定数をかけるのもIIRフィルタ(RCフィルタ)である。
FIRフィルタは、IIRに比べて安定とか線形位相という利点はあるが、計算に時間がかかるという欠点はある。FFTを使用した高速化処理の方法もある。しかし何より使い辛いのが、減衰の鋭さがサンプリング周波数でも変わってくるので、バンドパスフィルターのように性能が規定されている場合は毎回確認が必要となる。
FFTフィルタの概念はフーリエ変換で周波数情報にしてから通したい帯域のゲインを弄って逆フーリエ変換で波形に戻すという代物。位相の遅れは生じないし一見、理想的なフィルタに思えるが、フーリエ変換時の周波数分解能に影響されると思う。例えばの話分解能が10Hz刻みだとしたら逆フーリエ変換で出力される波形は随分大雑把な正弦波の合成となってしまう。
取りあえず一番実用なIIRフィルタのプログラムを作成してみる。作るのは4次のハイパスフィルタだ。6dB*4/Octの遮断特性を持つ。
IIRフィルタのプログラムの作り方は色々ある。これから行うのは双一次z変換というアナログで設計してからデジタルに変換するという方法。
この方法で面倒なのはフィルタ係数を求める部分で、遅延回路の部分はシンプルだ。フィルタ係数を求めるのは次数や遮断特性でも随分違うので、プログラムの作成自体はFIRの方が楽だと思う。
まあ、取りあえず開始。
バターワース特性の4次の伝達関数
1/(s^2 + 0.7654s + 1)(s^2 + 1.8478s + 1)
HPFの場合
s→2πfc/s
2πfc/s→A/sとして2πfcを変数Aとしてしまう。同様に0.7654と1.8478も変数P1,P2に変えてしまう。
以下、分母の((A/s)^2 + P1*(A/s) + 1)((A/s)^2 + P2*(A/s) + 1)の展開。
sが(分母の)分母にあるので、分母と分子にs^4をかけてしまう。こういう感じ。1*s^4/((A/s)^2 + P1*(A/s) + 1)((A/s)^2 + P2*(A/s) + 1)*s^4
s^4/(s^2*A^2*P1*P2+s*A^3*P2+s^3*A*P2+s*A^3*P1+s^3*A*P1+A^4+2*s^2*A^2+s^4)
=s^4/(A^4+s*(A^3*P1+A^3*P2)+s^2*(2*A^2+A^2*P1*P2)+s^3*(A*P1+A*P2)+s^4)
更にA^3*P1+A^3*P2等も変数にしてしまう。
=s^4
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
A^4+s*T1+s^2*T2+s^3*T3+s^4
又も面倒な作業が。sを(1-z^-1)/(1+z^-1)に置き換える。面倒なので(1-z)/(1+z)で表す。
s→(1-z)/(1+z)
*1^4/(A^4+((1-z)^3*T3)/(z+1)^3+((1-z)^2*T2)/(z+1)^2+((1-z)*T1)/(z+1)+(1-z)^4/(z+1)^4)
今度は分母と分子に(1+z)^4をかける。
(1-z)^4
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
- -
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
(A^4)*(1+z)^4+((1+z)^3)*(1-z)*T1+((1+z)^2)*((1-z)^2)*T2+(1+z)*((1-z)^3)*T3+(1-z)^4
分子の(1-z)^4を展開するとこうなる。
z^4-4*z^3+6*z^2-4*z+1
フィルタ計数bnが、
z^4→b(5)=1
- 4*z^3→b(4)=-4
- 6*z^2→b(3)=6
- 4*z→b(2)=-4
- 1→b(1)=1
となる。
分母はちょっと大変。フィルタ計数anが、
(A^4-T1-T3+T2+1)*z^4→a(5)=A^4-T1-T3+T2+1
(4*A^4-2*T1+2*T3-4)*z^3→a(4)=4*A^4-2*T1+2*T3-4
(6*A^4-2*T2+6)*z^2→a(3)=6*A^4-2*T2+6
(4*A^4-2*T3+2*T1-4)*z→a(2)=4*A^4-2*T3+2*T1-4
- A^4+T1+T2+T3+1→a(1)=A^4+T1+T2+T3+1
此処からはScilabでの作成
//遮断周波数を変換する。
Lcf=**;遮断周波数Hz
fs=**;サンプリング周波数Hz
fc=tan(%pi*Lcf/fs)/(2*%pi);
//作った変数を設定
A=2*%pi*fc;
P1=0.76545;
P2=1.84785;
T1=A^3*P1+A^3*P2;
T2=2*A^2+A^2*P1*P2;
T3=A*P1+A*P2;
//フィルタ計数
a=[5];
b=[5];
a(1)=A^4+T1+T2+T3+1;
a(2)=-1*(4*A^4-2*T3+2*T1-4)/a(1);
a(3)=-1*(6*A^4-2*T2+6)/a(1);
a(4)=-1*(4*A^4-2*T1+2*T3-4)/a(1);
a(5)=-1*(A^4-T1-T3+T2+1)/a(1);
b(1)=1/a(1);
b(2)=-4/a(1);
b(3)=6/a(1);
b(4)=-4/a(1);
b(5)=1/a(1);
a(1)=1;
//出力配列
out=[];
out(1)=b(1)*Implus(1);
out(2)=b(1)*Implus(2)+b(2)*Implus(1)+a(2)*out(1);
out(3)=b(1)*Implus(3)+b(2)*Implus(2)+b(3)*Implus(1)+a(2)*out(2)+a(3)*out(1);
out(4)=b(1)*Implus(4)+b(2)*Implus(3)+b(3)*Implus(2)+b(4)*Implus(1)+a(2)*out(3)+a(3)*out(2)+a(4)*out(1);
for i=5:length(Implus)
out(i)=b(1)*Implus(i)+b(2)*Implus(i-1)+b(3)*Implus(i-2)+b(4)*Implus(i-3)+b(5)*Implus(i-4)+a(2)*out(i-1)+a(3)*out(i-2)+a(4)*out(i-3)+a(5)*out(i-4);
end
///////////////////////////////////////////////////////////////////////////////////////////////////
ScilabのIIR関数と結果を比較する。Scilabは数行でフィルタを掛ける事が出来る。
hz=iir(4,'hp','butt',[Lcf,0]/fs,[1,1]);
出力=flts(入力,hz)';
インパルス応答で確認する。結果が重なっているが、青線のScilabの関数による曲線の下にプログラムの曲線がある。遮断周波数は300Hz。
*1:1-z)/(1+z
