周波数応答

周波数応答とは

ゲインと位相

安定な線形なシステムに正弦波入力を加えると，定常状態では出力も正弦波となる．出力の正弦波の周波数は入力と同じになるが，振幅の変化や位相差が発生する． 入力の正弦波を$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$を位相（位相差）と呼ぶ．ゲインと位相は入力の周波数に応じて変化する．

図1は，伝達関数$\displaystyle {\frac {1}{s^{2}+0.5s+1}}$に$\omega =$0.1, 1, 10 [rad/s]の3種類の正弦波を入力した時の応答をMATLAB（Octave）のlsim関数で計算した結果である（リストsample2_1）．

グラフから，過渡状態はおおよそ開始から10秒強で，その後は定常状態となっていることがわかる．定常状態でのゲインと位相は周波数によって変化していることがわかる．

図1 二次遅れ系の正弦波応答

リストsample2_1 二次遅れ系の正弦波応答比較スクリプト

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');

周波数応答関数

システムのゲインと位相は，上記のように実際に正弦波信号を入力して調べなくても，システムの周波数応答関数（周波数伝達関数）から求めることができる．周波数応答関数とは，伝達関数のラプラス演算子$s$を$j\omega$に置き換えたものである．周波数応答関数は角周波数$\omega$を含んだ複素関数であり，その絶対値と偏角がそれぞれゲインと位相に対応する．

${\begin{array}{l l l}g&=&|G(j\omega )|\\\phi &=&\angle G(j\omega )\end{array}}$

なお，ある複素数$z=a+bj$の絶対値と偏角はそれぞれ以下のように求められる．

${\begin{array}{l l l}|z|&=&{\sqrt {a^{2}+b^{2}}}\\\angle z&=&\tan ^{-1}{\frac {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関数を利用した方法

リストsample2_2にpolyval関数を利用してボード線図をプロットするスクリプト例を示す．

変数minpとmaxpにそれぞれ周波数軸の下限と上限を10のベキ数で指定する．この例では$10^{-2}=0.01$[rad/s]から$10^{2}=100$[rad/s]までを指定している．

グラフ表示のウィンドウをsubplot命令で2行1列に分割し，上側にゲイン線図，下側に位相線図をそれぞれ横軸対数表示であるsemilogx関数を使ってプロットしている．図2に実行結果のグラフを示す．

図2 $\displaystyle {\frac {1}{s^{2}+0.5s+1}}$のBode線図（polyval関数利用）

sample2_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');

リストsample2_3は，MATLAB（Octave）固有の要素ごとの演算を利用してforループを用いないスクリプト例である．ここで，logspace関数は等比数列を出力する関数である．

sample2_3 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関数を利用すると，簡単にボード線図を描画できる．リストsample2_4にスクリプト例を，図3に実行結果のグラフを示す．なお，この例では変数omegaで周波数ベクトルを指定しているが，省略することもできる．その場合はbode関数が自動的に周波数範囲を設定して描画する．

図3 $\displaystyle {\frac {1}{s^{2}+0.5s+1}}$のBode線図（bode関数利用）

sample2_4 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);

リストsample2_4の例のようにbode関数単体でグラフまで描画するが，ゲイン線図だけを描画したい時や，実験結果と重ねてグラフ化したい時などは，次のように，bode関数の左辺に変数をおくと，グラフ表示は行わず周波数応答の計算結果だけが得られる（ゲインは絶対値，位相は度の単位）．

[gain phase omega] = bode(sys, omega);

これを利用したスクリプト例をリストsample2_5に示す．これはボード線図のゲイン線図だけをプロットする例である．実行結果を図4に示す．

なお，bode関数が多入出力に対応しているため，MATLABの場合左辺の変数gainとphaseは3次元の行列として出力される（対象システム1入力1出力であっても）．よって，plot関数で使用する前にこのスクリプト例のようにsqueeze関数で1次元ベクトルに変換する必要があるので注意する．

図4 $\displaystyle {\frac {1}{s^{2}+0.5s+1}}$のBode線図（bode関数利用，ゲイン線図のみ）

sample2_5 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関数でプロットすると複素平面上で表されるので，それをそのまま利用できる．リストsample2_6がスクリプトの例であり，実行結果は図vector_loci_2_6である．なお，MATLAB（Octave）にはnyquist関数が用意されており，それを利用してもベクトル軌跡を描くことができる．

図5 $\displaystyle {\frac {1}{s^{2}+0.5s+1}}$のベクトル軌跡

sample2_6 二次遅れ系のベクトル軌跡を描くスクリプト

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 system Toolboxではシステムのむだ時間を指定する方法があるが，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}}$

よって，むだ時間要素のボード線図を描くスクリプト例はリストsample2_7となる．実行結果は図6である． ここで，関数onesは指定したサイズで要素がすべて1の行列を生成する関数であり，関数sizeで周波数ベクトルの大きさを取り出して利用している．

図6 むだ時間要素のボード線図（$L=1$）

sample2_7 むだ時間要素のボード線図を描くスクリプト

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');

実験結果との重ね合わせの例

実験科目や卒業研究などで，実験装置での実測結果と理論式から求めた理論曲線を比較する場面がよく出てくる．リストsample2_8は実験で求めた周波数応答と理論曲線を重ねて描くスクリプトの例である．

実験による測定結果の周波数，ゲイン，位相をそれぞれomega_e, gaindB_e, phase_eという変数に代入しておき，理論曲線と重ねている．この際，実験データの方はマーカーで示すべきであるので，semilog関数内で線類を指定している．また，関数legendでグラフに凡例を付け加えることができる．実行結果を図7に示す．

図7 実験結果と理論曲線を重ねる例

sample2_8 実験結果と理論曲線を重ねて描くスクリプト例

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');