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)
の形で結果を求めることができます。
