GNU Octaveで求める周波数特性

アクティブ・フィルタ

GNU Octaveで求める周波数特性

キルヒホッフの第一法則による連立方程式

多重帰還型(MFB)LPFを例に周波数特性の求め方を説明します。
(ただしここでは OPアンプは理想OPアンプとします。)

上図のように電流 i₁, i₂, i₃ を定めると、オームの法則と静電容量の定義式 Q ( = ∫ I dt ) = CV およびキルヒホッフの第一法則から、以下の連立方程式が得られます。

{\begin{cases} u = R_{1} (i_{1} + i_{2} + i_{3}) + \frac{1}{C_{1}} \int i_{1} ⅆ t \\ u = R_{1} (i_{1} + i_{2} + i_{3}) + R_{2} i_{2} \\ u = R_{1} (i_{1} + i_{2} + i_{3}) + R_{3} i_{3} + y \\ 0 = \frac{1}{C_{2}} \int i_{2} ⅆ t + y \end{cases}

上の連立方程式をラプラス変換すると、ラプラス変換の演算子 s により積分記号が消え、以下のように未知数 i₁(s), i₂(s), i₃(s), y(s) の連立一次方程式になります。

{\begin{cases} u (s) = R_{1} (i_{1} (s) + i_{2} (s) + i_{3} (s)) + \frac{1}{C_{1} s} i_{1} (s) \\ u (s) = R_{1} (i_{1} (s) + i_{2} (s) + i_{3} (s)) + R_{2} i_{2} (s) \\ u (s) = R_{1} (i_{1} (s) + i_{2} (s) + i_{3} (s)) + R_{3} i_{3} (s) + y (s) \\ 0 = \frac{1}{C_{2} s} i_{2} (s) + y (s) \end{cases}

この連立一次方程式を解けば伝達関数が得られます。

GNU Octaveでの演算手順

GNU Octave は Excelなどの表計算ソフトとは違い、直接、複素数の行列演算が可能です。

ラプラス変換後の連立一次方程式を行列表記すると

[\begin{matrix} u (s) \\ u (s) \\ u (s) \\ 0 \end{matrix}] = [\begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & 0 \\ R_{1} & R_{1} + R_{2} & R_{1} & 0 \\ R_{1} & R_{1} & R_{1} + R_{3} & 1 \\ 0 & \frac{1}{C_{2} s} & 0 & 1 \end{matrix}] [\begin{matrix} i_{1} (s) \\ i_{2} (s) \\ i_{3} (s) \\ y (s) \end{matrix}]

上の式中の正方行列をAとして

A = [\begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & 0 \\ R_{1} & R_{1} + R_{2} & R_{1} & 0 \\ R_{1} & R_{1} & R_{1} + R_{3} & 1 \\ 0 & \frac{1}{C_{2} s} & 0 & 1 \end{matrix}]

連立一次方程式は

[\begin{matrix} u (s) \\ u (s) \\ u (s) \\ 0 \end{matrix}] = A [\begin{matrix} i_{1} (s) \\ i_{2} (s) \\ i_{3} (s) \\ y (s) \end{matrix}]

これに左から A の逆行列 A^-1 をかけて

A^{- 1} [\begin{matrix} u (s) \\ u (s) \\ u (s) \\ 0 \end{matrix}] = [\begin{matrix} i_{1} (s) \\ i_{2} (s) \\ i_{3} (s) \\ y (s) \end{matrix}]

更にに左から行ベクトル [ 0 0 0 1 ] をかければ出力 y(s) が得られます。

[\begin{matrix} 0 & 0 & 0 & 1 \end{matrix}] A^{- 1} [\begin{matrix} u (s) \\ u (s) \\ u (s) \\ 0 \end{matrix}] = y (s)

伝達関数による方法

伝達関数の求め方は、VCVSフィルタとMFBフィルタの MFB(多重帰還型)フィルタで求めた方法以外にキルヒホッフの第一法則による連立方程式の最後に示した連立一次方程式を解くことによっても可能です。
むしろキルヒホッフの第一法則から連立方程式を立て、ラプラス変換で連立一次方程式化し、これを解く方が、手順として一般化が可能な分有効な手法かと思います。

連立一次方程式の解法は GNU Octaveでの演算手順の冒頭のように行列化すれば、クラメールの公式を適用することで手順を更に一般化できます。

$y (s) = \frac{| \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & u (s) \\ R_{1} & R_{1} + R_{2} & R_{1} & u (s) \\ R_{1} & R_{1} & R_{1} + R_{3} & u (s) \\ 0 & \frac{1}{C_{2} s} & 0 & 0 \end{matrix} |}{| \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & 0 \\ R_{1} & R_{1} + R_{2} & R_{1} & 0 \\ R_{1} & R_{1} & R_{1} + R_{3} & 1 \\ 0 & \frac{1}{C_{2} s} & 0 & 1 \end{matrix} |}$

あとは行列式をひたすら解いていきます。ここでの行列式の計算が苦にさえならなければ、とても便利な方法です。ではまず分子から

| \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & u (s) \\ R_{1} & R_{1} + R_{2} & R_{1} & u (s) \\ R_{1} & R_{1} & R_{1} + R_{3} & u (s) \\ 0 & \frac{1}{C_{2} s} & 0 & 0 \end{matrix} |

= \frac{u (s)}{C_{2} s} | \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & 1 \\ R_{1} & R_{1} & 1 \\ R_{1} & R_{1} + R_{3} & 1 \end{matrix} |

= \frac{u (s)}{C_{2} s} | \begin{matrix} \frac{1}{C_{1} s} & 0 & 0 \\ R_{1} & R_{1} & 1 \\ 0 & R_{3} & 0 \end{matrix} |

= - \frac{1}{C_{1} s} \frac{1}{C_{2} s} R_{3} u (s)

次は分母を

| \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} & 0 \\ R_{1} & R_{1} + R_{2} & R_{1} & 0 \\ R_{1} & R_{1} & R_{1} + R_{3} & 1 \\ 0 & \frac{1}{C_{2} s} & 0 & 1 \end{matrix} |

= | \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} \\ R_{1} & R_{1} + R_{2} & R_{1} \\ R_{1} & R_{1} - \frac{1}{C_{2} s} & R_{1} + R_{3} \end{matrix} |

= | \begin{matrix} R_{1} + \frac{1}{C_{1} s} & R_{1} & R_{1} \\ R_{1} & R_{1} + R_{2} & R_{1} \\ R_{1} & R_{1} - \frac{1}{C_{2} s} & R_{1} + R_{3} \end{matrix} |

= | \begin{matrix} R_{1} + \frac{1}{C_{1} s} & - \frac{1}{C_{1} s} & - \frac{1}{C_{1} s} \\ R_{1} & R_{2} & 0 \\ R_{1} & - \frac{1}{C_{2} s} & R_{3} \end{matrix} |

= R_{1} \frac{1}{C_{1} s} (R_{2} + \frac{1}{C_{2} s}) + R_{3} [(R_{1} + \frac{1}{C_{1} s}) R_{2} + R_{1} \frac{1}{C_{1} s}]

したがって多重帰還型(MFB)LPF の伝達関数は

\frac{y (s)}{u (s)} = \frac{- \frac{1}{C_{1} s} \frac{1}{C_{2} s} R_{3}}{R_{1} \frac{1}{C_{1} s} (R_{2} + \frac{1}{C_{2} s}) + R_{3} [(R_{1} + \frac{1}{C_{1} s}) R_{2} + R_{1} \frac{1}{C_{1} s}]} = \frac{- R_{3}}{R_{1} [R_{2} R_{3} C_{1} C_{2} s^{2} + (R_{2} + R_{3} + \frac{R_{2} R_{3}}{R_{1}}) C_{2} s + 1]}

GNU Octaveスクリプトとgnuplot

実際に GNU Octave で周波数特性を計算し gnuplotを呼びだしグラフ表示させた例を示します。GNU Octaveでの演算手順と伝達関数による方法の両方で計算し出力しています。
(回路定数はTexas Instruments 8-channel Audio DAC PCM1691 のデータシート記載された Typical Post-LPF 値を使用しました。)

以下 GNU Octave のスクリプトです。

clear -all
R1=5.1e3;
R2=470;
R3=7.5e3;
C1=4.7e-9;
C2=4.7e-10;
f1=logspace(1,6,64);
for k=1:64
  s=2*pi*f1(k)*j;
  y=-R3/R1/(R2*R3*C1*C2*s*s+(R2+R3+R2*R3/R1)*C2*s+1);
  y1dB(k)=20*log10(abs(y));
end
f2=logspace(1,6,256);
for k=1:256
  s=2*pi*f2(k)*j;
  A=[R1+1/C1/s R1 R1 0; R1 R1+R2 R1 0; R1 R1 R1+R3 1; 0 1/C2/s 0 1];
  b=[1;1;1;0];
  y=[0 0 0 1]*inv(A)*b;
  y2dB(k)=20*log10(abs(y));
end
semilogx(f1,y1dB,'xb;by Transfer Function;',f2,y2dB,'c;by Matrix;')
axis([10,1000000,-50,10])
grid on

平衡入力MFBフィルタ

←

→

2次LPFのステップ応答