Matlabでフィルタの通過特性を計算する

Caution

この方法は¹次系のフィルタでは一致しましたが、2次系では会いませんでした。
間違った方法なので使わないでください。

目的

フィルタの通貨特性は理論式で追うのが理想だけど、Matlab上でRC回路ポチポチつないで計算できる方法も知りたく、調査した。その結果を忘備録として記します。

対象

今回は理論式もわかっている　一次のRCフィルタを対象とする。
理論式は（https://analogista.jp/transfer-function1/) を参考としました。

下図も上記サイトの引用です。

上記の回路で　R=1Ω, C =0.001F とすると　カットオフ周波数$f_c$は下式になります。

$f_c = \frac{1}{2\pi RC}=159.2 (Hz) $

このフィルタは１次のLowPassFilterなので、この159Hzで-3dB下がって、位相が45度回転する伝達特性が得られれば、やりたいことができている状態になります。

作成モデルの考え方

Vinは①の経路のインピーダンスZinで決まり、　Voutは②の経路のインピーダンスZoutで決まります。
なのでMatlab上でZinとZoutのインピーダンスを算出し、Zout/Zinを算出すれば　Vout/Vin すなわち通過特性が求まると考えます。（この考え方あってますかね？ちゃんと理論武装できず、保証できませんが・・・）

実行コードとソフト

%入力側インピーダンスZinと　出力側インピーダンスZoutを計算し
% Zout / Zin を伝達特性とする。

%% 理論式
R = 1;
C = 1e-3;
sys = tf([1],[R*C, 1]);

[mag,phase,wout] = bode(sys,Freq*2*pi); %入力はradian
mag = squeeze(mag);
phase = squeeze(phase);

Cutoff_Freq  = 1/(2*pi*R*C); %159Hz 

%% Matlabモデルからインピーダンス読み込み
Z_data = power_zmeter('NoiseFilter',logspace(0,5,10000)); %1Hz ー　10^5　Hzまで

%Zin 
Freq  = Z_data.Freq;
Gain_in = abs( Z_data.Z(:,1) );
Phase_in = angle( Z_data.Z(:,1) )/pi*180;

%Zout
Gain_out = abs( Z_data.Z(:,2) );
Phase_out = angle( Z_data.Z(:,2) )/pi*180;

%%　描画して比較
close all
figure()
subplot 211
semilogx(Freq,20*log10([Gain_out./Gain_in]),wout/2/pi, 20*log10(mag), ...
    [Cutoff_Freq Cutoff_Freq],[0 -60],... %カットオフ周波数
    [10 10000],[-3 -3]... % -3dBライン
    );

subplot 212
semilogx(Freq,Phase_out - Phase_in,wout/2/pi,phase,...
    [Cutoff_Freq Cutoff_Freq],[0 -90],... %カットオフ周波数
    [10 10000],[-45 -45]) % -3dBライン)
ylim([-90 90])

Simlinkでは、インピーダンスを計算するブロックがあり、それを入出力部につけることでZout,Zinを計算します。作成したSimulinkファイルの絵は以下です。

実行結果

上がゲイン特性　下が位相特性で、横軸は周波数(Hz) 縦軸はパワー(dB)と　位相(deg)です。
理論式とSimより算出した式が一致しているため、青いラインにオレンジ色が完全に一致しています。

また想定通り159Hz（黄色の縦ライン）の位置で　-3dB, 45degの位相遅れ（紫色の横線）が発生しています。想定通りのようです。

おわりに

ここでは簡単な1次のフィルタについて考えましたが、回路が複雑になってくるとこの方式が
力を発揮すると思われます。誰かの役に立てば幸いです。