function G=ps63_6242_2004(n,d,e)
% function G=ps63_6242_2004(n,d,e)
%
% filter design for Problem 6.3(b), 6.242/Fall2004

% default parameters
if nargin<1, n=10; end
if nargin<2, d=0.1; end
if nargin<3, e=d; end

r=sqrt(1+d);

% generate A,B,CG,DG,CF,DF,CFi,DFi
A=toeplitz([-1;1;zeros(n-2,1)],[-1 zeros(1,n-1)]);
B=[1;zeros(n-1,1)];
CG=[zeros(1,n-1) 1;zeros(n)];
for k=1:n,
    CG(k+1,:)=CG(k,:)*A;
end
DG=[zeros(n,1);1];
DF=[1;zeros(n,1)];
CF=[zeros(1,n);eye(n)];
CFi=zeros(n+1,n);
for k=1:n,
    CFi(k,n-k+1)=r^(k-1);
end
DFi=[zeros(n,1);r^n];

% initialize LMI environment
abst_init_lmi;           

% define the LMI coefficients matrices
L1=abst([eye(n);zeros(1,n)]);
L2=abst([A';B']);
LG3=abst([CG';DG']);
LF3=abst([CF';DF']);
LF3i=abst([DFi';CFi']);  
EE=zeros(n+1);
EE(n+1,n+1)=1;
E=abst(EE);

% define decision variables
a=diagonal(n);       % diag(a_0,a_1,...,a_{n-1})
b=diagonal(n+1);     % diag(b_0,b_1,...,b_n)
Q1=symmetric(n); 
Q2=symmetric(n);
Q3=symmetric(n);
P1=symmetric(n);
P2=symmetric(n);
y=symmetric;

% define LMI's
af=L1*a*L1'+E;
L1*Q1*L2'+L2*Q1*L1'+LG3*af*LG3'+y*E>0;
L1*Q2*L2'+L2*Q2*L1'+LG3*b*LG3'+y*E>0;
L1*Q3*L2'+L2*Q3*L1'+LG3*((1+e)*af-b)*LG3'+y*E>0;
L1*P1*L2'+L2*P1*L1'+LF3*(b-(1-e)*af)*LF3'+y*E>0;
L1*P2*L2'+L2*P2*L1'+LF3i*(e*af-b)*LF3i'+y*E>0;

% run LMI optimization
[cost,xopt]=lmi_mincx_tbx(y);

disp(['ymin = ' num2str(cost)])

aa=[diag(value(a));1];
bb=diag(value(b));
aa=aa(n+1:-1:1);
bb=bb(n+1:-1:1);
w=logspace(-3,3,2000);
f=polyval(bb,w.^2)./polyval(aa,w.^2);

aaa=[zeros(1,2*n) aa(n+1)];
bbb=[zeros(1,2*n) bb(n+1)];
for k=1:n,
        aaa(2*k)=0;
        bbb(2*k)=0;
        aaa(2*k-1)=(-1)^(n-k+1)*aa(k);
        bbb(2*k-1)=(-1)^(n-k+1)*bb(k);
end
sp=roots(bbb);
sq=roots(aaa);
p=poly(sp(real(sp)<0));
q=poly(sq(real(sq)<0));
p=p*abs(polyval(q,j/2)/polyval(p,j/2))*sqrt(polyval(bb,1/4)/polyval(aa,1/4));
G=tf(p,q);

% plot the amplitude responce
close(gcf)
subplot(2,1,1);
semilogx(w,abs(polyval(p,j*w)./polyval(q,j*w)).^2,w,1+e,w,1-e,w,e)
subplot(2,1,2);
w1=linspace(0.5,2,1000);
plot(w1,abs(polyval(p,j*w1)./polyval(q,j*w1)).^2,w1,1+e,w1,1-e,w1,e);grid