load data/spk1.dat; load data/spk2.dat; f = spk1(:,1); % 周波数(Hz)。 w = 2*pi*f; % 角周波数(rad/s) n = size(spk1,1); % 周波数応答のデータ点数 Gzw_jw = spk1(:,2).*exp(i*pi/180.0*spk1(:,3)); % w から z までの周波数応答 Gyw_jw = spk1(:,4).*exp(i*pi/180.0*spk1(:,5)); % w から y までの周波数応答 Gzu_jw = spk2(:,2).*exp(i*pi/180.0*spk2(:,3)); % u から z までの周波数応答 Gyu_jw = spk2(:,4).*exp(i*pi/180.0*spk2(:,5)); % u から y までの周波数応答 Gzw_g = frd(Gzw_jw, w); % 角周波数とともに frd(frequency response data)モデルにまとめる。(以降のデータの取り扱いの容易のため) Gzu_g = frd(Gzu_jw, w); % 同上 Gyw_g = frd(Gyw_jw, w); % 同上 Gyu_g = frd(Gyu_jw, w); % 同上 fprintf(1, '----- menu -----\n'); fprintf(1, 'g,G ... gain for phase-delay controller\n'); fprintf(1, 'd,D ... delay for phase-delay controller\n'); fprintf(1, 'q ... quit\n'); GAIN_MAX = 2.0; GAIN_STEP = 0.1; gain = 0; delay = 0; ND = 200; while 1 fprintf(1, 'gain = %f, delay = %d (*0.2ms)\n', gain, delay); s = input('command?:', 's'); switch(s) case 'q' break; case 'g' gain = gain - GAIN_STEP; if gain < 0 gain = 0; end case 'G' gain = gain + GAIN_STEP; if gain > GAIN_MAX gain = GAIN_MAX; end case 'd' if delay > 0 delay = delay - 1; end case 'D' if delay < ND delay = delay + 1; end end tau = delay*0.2e-3; % むだ時間 K = gain*exp(-1i*w*tau); % phase-delay controller の周波数応答 K_g = frd(K, w); % frd モデルにまとめる figure(1); nyquist(-Gyu_g*K_g, 'r'); % ナイキスト軌跡のプロット axis([-1.5, 0.5, -1, 1]); figure(2); Gcl_g = Gzw_g + Gzu_g*K_g/(1-Gyu_g*K_g)*Gyw_g; % 閉ループ系の周波数応答 bodemag(Gzw_g, 'b', Gcl_g, 'r'); % そのゲイン特性の表示。 legend('open loop', 'closed loop'); end save data/phasedelay.mat gain delay; % 最後の gain と delay を調整された値として保存(別のプログラムで読み込んで使うために)