MATLABでbodeコマンドを使わずにbode線図を描く方法

　ということがだいたいできたので今日は満足です。伝達関数を状態変数表示にして（tf→ss）
D+C(sI-A)^(-1)B で周波数応答が出てくる。今度暇なときにコード書き出そうかな。

☆追記(20120703)------

　上に書いてある事は間違い。tf形式のまま直接計算できる。

　bode線図は制御系の開ループ伝達関数をG(s)と表すと
ゲイン線図：x軸→周波数ω（対数目盛）　y軸→20×log_10(|G(jω)|) [dB]
位相線図：　x軸→周波数ω（対数目盛）　y軸→G(jω)の位相[deg]
となるようなグラフである。

○伝達関数を分子係数と分母係数にわける。

　伝達関数G(s)がtf形式で与えられているときG(s)→hとすると
[num,den] = tfdata(h,'v')
で分子係数num と分母係数den　に分けることができる。

　また状態空間モデルss形式で与えられているときはss2tfコマンドを使い以下のようにする。
[num,den] = ss2tf(A,B,C,D,iu);
これはiu番目の入力から出力への伝達関数をtf形式に変換し、分子係数と分母係数をわけている。
1入力系ならiu=1で問題はないはず。

○bode線図を描画

　例として、G(s)=10/(s^2+3s+2)　のときのbode線図を描くコードを下に書く。
(参考：http://www.mathworks.co.jp/support/solutions/ja/data/1-9GCLEM/index.html?solution=1-9GCLEM
http://infoshako.sk.tsukuba.ac.jp/ShakoDoc/MATLAB5/jhelp/techdoc/ref/angle.html
http://infoshako.sk.tsukuba.ac.jp/ShakoDoc/MATLAB5/jhelp/techdoc/ref/unwrap.html)

num=[0 0 10];den=[1 3 2]; %G(s)の分子係数はnum 分母係数はden
w = logspace(-1,2,100);　%ω→w　を10^(-1)から10^2まで対数的に等間隔にとる。wは1×100のベクトル。
j = sqrt(-1); %虚数 jを定義
Gjw = polyval(num,j*w)./polyval(den,j*w);　%　※1
Mag_dB = 20*log10(abs(Gjw));　%absで複素数の大きさを計算
Ph = unwrap(angle(Gjw))/(pi/180);　%※2
subplot(211)
semilogx(w,Mag_dB)
xlabel('Frequency (rad/s)'); ylabel('Magnitude (dB)');
subplot(212)
semilogx(w,Ph)
xlabel('Frequency (rad/s)'); ylabel('Phase (deg)');

※1

　polyval(B,x)は多項式Bにxを代入した値を計算するコマンド。このときBは次数の降順に多項式の係数をもつベクトルで、s^2+3s+2ならB=[1 3 2]。
　また、" ./ " は行列の要素ごとの除算である。wは1×100のベクトルなので、polyval(num,j*w)とpolyval(den,j*w)の計算結果をそれぞれ[A(jω1)　A(jω2)　　…　A(jω100)] [B(jω1)　B(jω2)　　…　B(jω100)]とすると、polyval(num,j*w)./polyval(den,j*w)　は
　[A(jω1)/B(jω1)　A(jω2)/B(jω2)　　…　A(jω100)/B(jω100)]
　となる。
　ここで注意するのは"./"は行列のサイズが同じでないと計算できないということ。たとえ分母と分子の次数が異なる場合はnum=[0 0 10]のように、少ない方に0を入れて調整する必要がある。

※2
　angle(G(jw))で伝達関数の位相を計算し、ラジアン単位で出力する。
　unwrap(P) は、連続する配列要素間の変化の絶対値が π より大きいとき、±2 の倍数を加えて配列 P 内の位相角を修正する。P が行列ならば、unwrap は列ごとに演算する。これによって位相曲線が連続的に出力される。角度単位で出力するので(pi/180)で除す。

　以上。何か間違えているところがあれば教えてください。

☆追記終わり----------------

　そして久しぶりのブックレビュー。