% nominal.m % Purpose: to obtain nominal plant from frequency response % Input: data/spk1.dat, data/spk2.dat; % Output: data/nominal.mat close all; clear all; % z(V) <-----| Gzw Gzu |<----- w(V) % y(V) <-----| Gyw Gyu |<----- u(V) % % spk1.dat ... frequency data from w to z and y % spk2.dat ... frequency data from u to z and y % % 1st column: frequency (Hz) % 2nd and 3rd columns: gain and phase characteristic from w or u to y % 4th and 5th columns: gain and phase characteristic from w or u to z fl = 10; % the lowest frequency to make graph fh = 2000; % the highest frequency to make graph load data/spk1.dat; % load frequency response data load data/spk2.dat; % f = spk1(:,1); % frequency w = 2*pi*f; % angular frequency n = size(spk1,1); % number of test frequencies Gyw_jw = spk1(:,2).*exp(1i*pi/180.0*spk1(:,3)); Gzw_jw = spk1(:,4).*exp(1i*pi/180.0*spk1(:,5)); Gyu_jw = spk2(:,2).*exp(1i*pi/180.0*spk2(:,3)); Gzu_jw = spk2(:,4).*exp(1i*pi/180.0*spk2(:,5)); for k = 1:n resp(1,1,k) = Gzw_jw(k); resp(1,2,k) = Gzu_jw(k); resp(2,1,k) = Gyw_jw(k); resp(2,2,k) = Gyu_jw(k); end G_g = frd(resp, w); G_jw = []; for k=1:n tmp = [Gzw_jw(k), Gzu_jw(k); Gyw_jw(k), Gyu_jw(k)]; G_jw = [G_jw; tmp]; end G = vpck(G_jw, w); m = 2; p = 2; [A, B, C, D] = subspace(G, m, p, 24); G0 = ss(A, B, C, D); G0_g = frd(G0, w); figure(2); bode(G_g, 'b-', G0_g, 'r-'); legend('exp.', 'nominal plant'); save data/nominal.mat G0 G_g G0_g w