% fresp.m - generate a frequency response (by hand)
% RMM, 25 Aug 05
aminit; % initialize MATALB
%
% This file shows how to compute a frequency response by hand. This
% would almost never be done in practice, but is used in Chapter 2 to
% show how models can be used for more than just time responses.
%
% System definition: use the same one as we did in the intro
A = [-0.2 2 0; -0.5 -0.2 4; 0 0 -10];
B = [0; 0; 0.5]; C = [2.6 0 0]; D = [0];
% Scale to get the plots to look like I want
sys = 1.4*ss(A,B,C,D); % state space object
W1list = [0.1 0.4 1 3];
style = {'b-', 'b-', 'b-', 'b-'};
figure(1); clf; % clear the graphics window
subplot(221); % set size of window
clear G1list; % remove list in case len changes
for i = 1:length(W1list)
fprintf('Frequency: %g\n', W1list(i));
% Simulate 10 cycles of input
count = 1:10000;
time = 2*pi/10 * (count - 1) / max(W1list);
u = sin(W1list(i) * time);
y = lsim(sys, u, time);
% Plot the output on top of the others
xy = plot(time(1:1000), y(1:1000), char(style(i)));
hold on;
% Compute the magnitude of the response based on the last two cycles of data
start = floor(10000 - max(W1list)/W1list(i)*20);
G1list(i) = max(y(start:10000));
end;
% Finish up the labelling for the time plot
amaxis([0 50 -4 4]); grid on;
xlabel('Time [s]'); ylabel('Output y','Rotation', 90);
amprint('fresptime.eps');
%
% Now plot the frequency response
%
% We redo the plots with more points, marking those points that
% were in the first list
%
W2list = [0.1 0.2 0.4 1 1.6 3 8 10];
clear Glist;
for i = 1:length(W2list)
fprintf('Frequency: %g\n', W2list(i));
% Simulate 10 cycles of input
count = 1:10000;
time = 2*pi/10 * (count - 1) / max(W1list);
u = sin(W2list(i) * time);
y = lsim(sys, u, time);
% Compute the magnitude of the response based on the last two cycles of data
start = floor(10000 - max(W2list)/W2list(i)*20);
G2list(i) = max(y(start:10000));
end;
figure(2); clf; subplot(221);
loglog(W2list, G2list, 'o-');
hold on; xy = loglog(W1list, G1list, '.');
set(xy, 'MarkerSize', 20); % fill the circle
amaxis([0.1 10 1e-2 10]);
grid on; set(gca, 'YMinorGrid', 'off');
xlabel('Frequency [rad/s]');
ylabel('Gain','Rotation', 90);
amprint('frespfreq.eps');