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