周波数応答とは
ゲインと位相
安定で線形なシステムに正弦波入力を加えると，定常状態では出力も正弦波となる．出力の正弦波の周波数は入力と同じになるが，振幅の変化や位相差が発生する．
次の図は，伝達関数$\displaystyle \frac{1}{s^2 + 0.5s + 1}$に$\omega= 0.1, 1, 10$ [rad/s]の3種類の正弦波を入力した時の応答をMATLAB（Octave）のlsim関数で計算した結果である．グラフから，過渡状態はおおよそ開始から10秒強で，その後は定常状態となっていることがわかる．これらのグラフから，この伝達関数の場合定常状態において，
- 周波数が低い場合は入出力にほとんど差がない
- 周波数が1 [rad/s]の時は入力より出力の振幅の方が大きく，また，入力に対し出力が1/4周期程度遅れている
- 周波数が高い場合は出力の振幅が非常に小さくなる
ことがわかる．このように周波数に応じて入出力の振幅の関係や時間差（位相ずれ）が変化する．
num = [1];
den = [1 0.5 1];
sys = tf(num, den);
tim = 0:0.01:20;
u1 = sin(0.1*tim);
u2 = sin(1*tim);
u3 = sin(10*tim);
y1 = lsim(sys, u1, tim);
y2 = lsim(sys, u2, tim);
y3 = lsim(sys, u3, tim);
subplot(3,1,1);
plot(tim, u1, tim, y1)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 0.1 (rad/s)');
legend('u', 'y');
grid( 'on');
subplot(3,1,2);
plot(tim, u2, tim, y2)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 1 (rad/s)');
legend('u', 'y');
grid( 'on');
subplot(3,1,3);
plot(tim, u3, tim, y3)
axis([0 20 -2 2]);
xlabel('time (s)');
ylabel('u, y');
title('\omega = 10 (rad/s)');
legend('u', 'y');
grid( 'on');
そこで入力の正弦波を$u(t)$，定常状態での出力を$y(t)$とし，それぞれ次式で表されるとする．
$$ \begin{array}{l l l} u(t) &=& A_i \sin \omega t \\ y(t) &=& A_o \sin(\omega t + \phi) \end{array} $$
定常状態での波形の例を図に示す．
このとき，周波数応答は次のゲインと位相の二つの情報の組で表現される．これらが周波数に応じて変化することになる．
- ゲイン: $g = \displaystyle \frac{A_o}{A_i}$
- 位相（位相差）: $\phi = -\dfrac{\tau}{T} \times 2\pi \, \mathrm{[rad] , or} -\dfrac{\tau}{T}\times 360^{\circ}$
周波数応答関数
システムのゲインと位相は，上記のように実際に正弦波信号を入力して調べなくても，システムの周波数応答関数（周波数伝達関数）から求めることができる．周波数応答関数とは，伝達関数のラプラス演算子$s$を$j\omega$に置き換えたものである．周波数応答関数は角周波数$\omega$を含んだ複素関数であり，その絶対値と偏角がそれぞれゲインと位相に対応する．
$$ \begin{array}{l l l} g &=& |G(j \omega)| \\ \phi &=& \angle G(j \omega) \end{array} $$
なお，ある複素数$z = a + b j$の絶対値と偏角はそれぞれ以下のように求められる．
$$ \begin{array}{l l l} |z| &=& \sqrt{a^2 + b^2} \\ \angle z &=& \tan^{-1} \dfrac{b}{a} \end{array} $$
MATLAB（Octave）での周波数応答の計算
MATLAB（Octave）での複素数の取扱い
MATLAB（Octave）では複素数は以下のように入力し，そのまま実数と同じように扱える．なお，入力する際は，虚数単位はiでもjでもよい．
1 + 2i
例えば，-1の平方根は
> sqrt(-1)
ans = 0 + 1i
と出力される．
しかし，変数としてiやjを使っていると，虚数単位と混同するので，複素数を表す関数complex()を使用した方が間違いがない．また，絶対値は関数abs()，偏角は関数angle()で求められる（単位はラジアン）．以下のこれらの使用例を示す．
> z = complex(1, -1)
z = 1 - 1i
> abs(z)
ans = 1.4142
> angle(z)
ans = -0.78540
MATLAB（Octave）での有理式の取扱い
伝達関数は一般に有理式（分子・分母が多項式）で表される（むだ時間など一部除く）．MATLAB（Octave）では多項式を係数を降べきの順に並べたベクトルで表し，有理式は分子・分母それぞれの係数ベクトルのペアで表現する．例えば，伝達関数$\displaystyle G(s) = \frac{s + 2}{s^2 + 3s + 2}$は，次のように二つのベクトルnumとdenのペアで表す．
> num = [1 2];
> den = [1 3 2];
また，MATLAB（Octave）では，多項式の変数に値を代入して計算した結果を返してくれる関数polyval()が用意されている．使用例を以下に示す．このように複素数にも対応するので，周波数応答の計算に利用できる．
> f = [1 2 1];
> polyval(f, 1)
ans = 4
> polyval(f, complex(1, 1))
ans = 3 + 4i
周波数応答（ボード線図）のプロット
周波数応答関数のゲインと位相をそれぞれ別のグラフで表したものがボード線図である．横軸は角周波数で対数目盛りとし，ゲインはdBに換算して表示する． ここではMATLAB（Octave）によるボード線図のプロットについて，以下の二つの方法を紹介する．
polyval関数を利用して原理に基づいて計算する方法- Control System Toolboxの
bode関数を利用する方法
polyval関数を利用した方法
以下にpolyval関数を利用してボード線図をプロットするスクリプト例を示す．
変数minpとmaxpにそれぞれ周波数軸の下限と上限を10のベキ数で指定する．この例では$10^{-2}=0.01$[rad/s]から$10^2=100$[rad/s]までを指定している．
グラフ表示のウィンドウをsubplot命令で2行1列に分割し，上側にゲイン線図，下側に位相線図をそれぞれ横軸対数表示であるsemilogx関数を使ってプロットしている．図2に実行結果のグラフを示す．
polyval関数を利用したボード線図の描画スクリプト}
num = [1]; %伝達関数の分子
den = [1 0.5 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
gain =[];
phase = [];
omega = [];
for k = 1 : n
w = 10^(minp+(maxp-minp)/(n-1)*(k-1));
z = polyval(num, complex(0, w))/polyval(den, complex(0, w));
gain(k) = 20*log10(abs(z));
phase(k) = angle(z)/3.14159*180;
omega(k) = w;
end
subplot(2,1,1);
semilogx(omega, gain);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid('on');
subplot(2,1,2);
semilogx(omega, phase);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');
以下のリストは，MATLAB（Octave）固有の要素ごとの演算を利用してforループを用いないスクリプト例である．ここで，logspace関数は等比数列を出力する関数である．
polyval関数を利用したボード線図の描画スクリプト（forループ無しバージョン）
num = [1]; %伝達関数の分子
den = [1 0.5 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n); %周波数ベクトルの生成
z = polyval(num, complex(0, omega))./polyval(den, complex(0, omega));
gain = 20*log10(abs(z));
phase = angle(z)/3.14159*180;
subplot(2,1,1);
semilogx(omega, gain);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');
subplot(2,1,2);
semilogx(omega, phase);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');
bode関数を利用した方法
MATLABのControl System
Toolboxに含まれるbode関数を利用すると，簡単にボード線図を描画できる．以下にスクリプト例と実行結果のグラフを示す．なお，この例では変数omegaで周波数ベクトルを指定しているが，省略することもできる．その場合はbode関数が自動的に周波数範囲を設定して描画する．
bode関数を利用したボード線図の描画スクリプト
num = [1]; %伝達関数の分子
den = [1 1 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n); %周波数ベクトルの生成
sys = tf(num, den); %伝達関数からシステムマトリックスへの変換
bode(sys, omega);
先のリストの例のようにbode関数単体でグラフまで描画するが，ゲイン線図だけを描画したい時や，実験結果と重ねてグラフ化したい時などは，次のように，bode関数の左辺に変数をおくと，グラフ表示は行わず周波数応答の計算結果だけが得られる（ゲインは絶対値，位相は度の単位）．
[gain phase omega] = bode(sys, omega);
これを利用したスクリプト例を以下のリストに示す．これはボード線図のゲイン線図だけをプロットする例である．実行結果を図に示す．
なお，bode関数が多入出力に対応しているため，MATLABの場合左辺の変数gainとphaseは3次元の行列として出力される（対象システム1入力1出力であっても）．よって，plot関数で使用する前にこのスクリプト例のようにsqueeze関数で1次元ベクトルに変換する必要があるので注意する．
bode関数を利用したゲイン線図のスクリプト
num = [1]; %伝達関数の分子
den = [1 0.5 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n); %周波数ベクトルの生成
sys = tf(num, den); %伝達関数からシステムマトリックスへの変換
[gain phase omega] = bode(sys, omega);
gain = squeeze(gain); % 1次元ベクトルへの変換 gain=gain(:); でもよい
gaindB = 20*log10(gain);
semilogx(omega, gaindB);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');
周波数応答（ベクトル軌跡）のプロット
ベクトル軌跡は周波数応答関数を$\omega = 0 \sim \infty$の範囲で複素平面上にそのままプロットしたものである．MATLAB（Octave）では複素数ベクトルをplot関数でプロットすると複素平面上で表されるので，それをそのまま利用できる．次のリストがスクリプトの例であり，実行結果は図vector_loci_2_6である．なお，MATLAB（Octave）にはnyquist関数が用意されており，それを利用してもベクトル軌跡を描くことができる．
二次遅れ系のベクトル軌跡を描くスクリプト
num = [1]; %伝達関数の分子
den = [1 0.5 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n);
z = polyval(num, complex(0, omega))./polyval(den, complex(0, omega));
plot(z)
axis([-1.5 1.5 -1.5 1.5]);
axis equal; %グラフの縦横比を1:1に指定
grid('on');
応用編
むだ時間要素のボード線図
むだ時間要素$G(s)=e^{-Ls}$は無理式であるので，分子分母にわけて係数ベクトルで表現することはできない．MATLABのControl SystemToolboxではシステムのむだ時間を指定する方法があるが，Octaveではまだ実装されていないようである．ここでは，むだ時間のボード線図を原理に基づいてプロットすることを考える．
むだ時間要素の周波数応答関数は$G(j\omega)=e^{-Lj\omega}$であり，その絶対値と偏角は以下となる．
$$ \begin{array}{l l l} g &=& |G(j \omega)| = 1 \\ \phi &=& \angle G(j \omega) = -L\omega \end{array} $$
よって，むだ時間要素のボード線図を描くスクリプト例は以下のリストとなる．実行結果も示す．
ここで，関数onesは指定したサイズで要素がすべて1の行列を生成する関数であり，関数sizeで周波数ベクトルの大きさを取り出して利用している．
むだ時間要素のボード線図を描くスクリプト
L = 1; %むだ時間 [s]
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n); %周波数ベクトルの生成
gain = 1*ones(size(omega)); %omegaと同サイズのベクトルを生成
gaindB = 20*log10(gain);
phase = -L*omega/3.14159*180;
subplot(2,1,1);
semilogx(omega, gaindB);
axis([0.01, 100, -80, 20]);
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid('on');
subplot(2,1,2);
semilogx(omega, phase);
axis([0.01, 100, -300, 0]);
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');
実験結果との重ね合わせの例
実験科目や卒業研究などで，実験装置での実測結果と理論式から求めた理論曲線を比較する場面がよく出てくる．次のリストは実験で求めた周波数応答と理論曲線を重ねて描くスクリプトの例である．
実験による測定結果の周波数，ゲイン，位相をそれぞれomega_e, gaindB_e,phase_eという変数に代入しておき，理論曲線と重ねている．この際，実験データの方はマーカーで示すべきであるので，semilog関数内で線類を指定している．また，関数legendでグラフに凡例を付け加えることができる．実行結果を図に示す．
実験結果と理論曲線を重ねて描くスクリプト例
omega_e = [0.01 0.03 0.05 0.1 0.3 0.5 1 3 5 10 30 50 100];
gaindB_e = [0 -2.1 3.1 0.1 0.4 1.0 -0.5 -20.1 -27.5 -45.5 -55.2 -69.1 -78.2];
phase_e = [-0.5 -2.1 -3.1 -5.5 -22.1 -32.1 -90.1 -150 -173 -185 -192 -205 -210];
num = [1]; %伝達関数の分子
den = [1 1 1]; %伝達関数の分母
minp = -2; %周波数下限のベキ数
maxp = 2; %周波数上限のベキ数
n = 100; %周波数分割数
omega = logspace(minp, maxp, n);
sys = tf(num, den);
[gain phase omega] = bode(sys, omega);
gain = squeeze(gain);
gaindB = 20*log10(gain);
phase = squeeze(phase);
subplot(2,1,1);
semilogx(omega, gaindB, '-', omega_e, gaindB_e, '*');
xlabel('Frequency (rad/s)');
ylabel('Gain (dB)');
grid( 'on');
legend('Theoretical', 'Experimental');
subplot(2,1,2);
semilogx(omega, phase, '-', omega_e, phase_e, '*');
xlabel('Frequency (rad/s)');
ylabel('Phase (deg)');
grid('on');
legend('Theoretical', 'Experimental');