3次以上のデジタルフィルタの計算精度

以下の仕様で1次の連続時間系のHPFを求めます。

• カットオフ周波数：0.03Hz
• 通過域のゲインは１倍

2次と3次のHPFは、1次のHPFを直列接続して構成します。

この場合、サンプリング周波数5KHzでc2dを用いて離散時間系のフィルタを求めると、2次までは連続時間系と同じ周波数応答が得られますが、3次になるとカットオフ周波数以下で望む特性が得られません。

MATLABコードと結果のプロットを以下に示します。

clear ; format short e ; s=tf('s') ;

% frequency

f=(logspace(-2,3,501))' ; w = 2*pi*f ;

fs = 5e3 ; dT = 1/fs ;

% condition

fc = 0.03 ; % Hz

Tc = 1/(2*pi*fc) ;

% continuous

TF_HPF1 = Tc*s/(Tc*s+1) ; % 1st order HPF

TF_HPF2 = TF_HPF1^2 ; % 2nd order HPF

TF_HPF3 = TF_HPF1^3 ; % 3rd order HPF

% discrete

DF_HPF1 = c2d(TF_HPF1, dT) ;

DF_HPF2 = c2d(TF_HPF2, dT) ;

DF_HPF3 = c2d(TF_HPF3, dT) ;

% plot

figure(1) ; clf ;

bode(TF_HPF1, DF_HPF1) ; title('1st order') ; grid on ;

figure(2) ; clf ;

bode(TF_HPF2, DF_HPF2) ; title('2nd order') ; grid on ;

figure(3) ; clf ;

bode(TF_HPF3, DF_HPF3) ; title('3rd order') ; grid on ;

青：連続時間系HPF、赤：離散時間系HPF

ご覧のように1次と2次は連続時間と離散時間で完全に一致していますが、3次では低周波での特性が連続時間と離散時間で大きく異なります。4次以上でも異なります。

原因と解決方法がお分かりになる方はいらっしゃいますか？

Answer 1

[Naoya]

1次のディジタルフィルタは "DF_HPF1 = c2d(TF_HPF1, dT);" で表現される部分になると思います。

このフィルタ係数は、

>> [num,den] = tfdata(DF_HPF1,'v')

num =

1.0000e+00  -1.0000e+00


den =

1.0000e+00  -9.9996e-01

で取得できます。

もう少し小数点以下の精度を上げて確認してみますと、以下のように確認できます。

>> sprintf('%20.19f\n',num)

ans =

'1.0000000000000000000
-0.9999999999999998890
'

>> sprintf('%20.19f\n',den)

ans =

'1.0000000000000000000
-0.9999623015987595398

ご指摘の 3次のフィルタについては、この num, den を畳み込み演算していることになります。

>> num3 = conv(conv(num,num),num)

num3 =

1.0000e+00  -3.0000e+00   3.0000e+00  -1.0000e+00

>> den3 = conv(conv(den,den),den)

den3 =

1.0000e+00  -2.9999e+00   2.9998e+00  -9.9989e-01

今回3次のフィルタの低周波数側の振幅が落ち切らない要因としては、この倍精度浮動小数点レベルの畳み込みの内部演算過程で丸め誤差が生じているものと推測されます。

Signal Processing Toolbox をお持ちであれば、dfilt オブジェクト等を使って、3次の1つのフィルタではなく、 1次のフィルタを3つカスケード接続した多段フィルタというそのままの形で実現することができます。 これにより、3次の単一フィルタに変換した際の数値演算の丸め誤差の影響を抑えることができます。

% 1次のディジタルフィルタ
>> Hd = dfilt.df2t(num,den)

% 1次のフィルタを3段接続する
>> Hd3 = cascade(Hd,Hd,Hd)

% 周波数応答表示 (横軸単位は Hz)
freqz(Hd3,logspace(-3,1,1000),1/dT);
ax = findobj(gcf,'Type','Axes');
for n = 1:length(ax)
    set(ax(n),'XScale','log')
end

なお、このフィルタを用いて、フィルタリングする場合は、

>> y = filter(Hd3, x)

の形で結果を求めることができます。

[Takashi]

回答ありがとうございます。

dfiltは知りませんでしたので、よく調べてみます。