%% Problem 18.1 discretization of Butterworth filter
% transfer function and CT state-space model
num=[1];
den=[1/(2000*pi)^2 sqrt(2)/(2000*pi) 1];
[A,B,C,D]=tf2ss(num,den);

T=[0 1; 1 0]; % permutation matrix to get same form as in Chapter 11
A=inv(T)*A*T;
B=inv(T)*B;
C=C*T;
H=ss(A,B,C,D);

% bilinear transf
Ts=0.0001
Ab=inv(eye(2)-0.5*Ts*A)*(eye(2)+0.5*Ts*A);
Bb=Ts*inv(eye(2)-0.5*Ts*A)*inv(eye(2)-0.5*Ts*A)*B;
Cb=C;
Db=D+0.5*Ts*C*inv(eye(2)-0.5*Ts*A)*B;
Hb=ss(Ab,Bb,Cb,Db,Ts);

% Frequency response of CT Butterworth and its discretized version
w=logspace(1,log10(10000*pi),200);
w=w(1,1:199);
[MAG,PHASE] = bode(H,w);
[MAGbilin,PHASEbilin] = bode(Hb,w);

figure(1)
semilogx(w,20*log10(MAGbilin(:,:)),w,20*log10(MAG(:,:)))
figure(2)
semilogx(w,PHASEbilin(:,:),w,PHASE(:,:))

% Step responses
figure(3)
t=[0:0.00001:.00299];  % time vector to plot step resp of CT system
[y,ts,x]=lsim(H,ones(1,300),t);
plot(ts,y)
hold on
[yb,tsb,xb]=lsim(Hb,ones(1,30));
plot(tsb,yb,'o')
hold off

