% FIGURE 9.1:
% fig 1 to 6 correspond to Figure 9.1 a to 9.1 f

for fig=[1:6]
switch fig
 case {1,5,6}, b=[0 1 0.5];
         a=[1 -1.5 0.7];
 case 2, b=[0 1 -0.9];
         a=[1 -0.8 0.2];
 case 3, b=[0 1 0];
         a=[1 -1.8 0.81];
 case 4, b=[0 1 0.5];
         a=[1 -1.4 0.98];
end
figext=['a','b','c','d','e','f'];

m0=poly2th(a,b);
m1=poly2th(a,1);

if fig==6, m1=poly2th(a,[1 -0.5]/0.5);end
gm1=th2ff(m1,1,logspace(-2,pi,1000));
[om1,amp1,dum,sdamp]=getff(gm1);

figure(fig)
loglog(1,1);
if fig~=5,line(om1,amp1.^2,'linestyle','-','linewidth',2),
       else line(om,baseres/2,'linestyle','-','linewidth',2)
end
if fig==5 ns=2;else ns=[2 10];end
for n=ns
 
   a1=[-a(2) -b(2) -a(3) -b(3), zeros(1,2*n-4)];
   a2=kron(eye(2),eye(n-1));
   A=[[a1;zeros(1,2*n)];[a2,zeros(2*n-2,2)]]; 
   B=zeros(2*n,2);B(1,2)=1;B(2,1)=1;
   if fig>4
      a1=[-a(2) -b(2) 0 -a(3) -b(3) 0,zeros(1,3*n-6)]; % ARMAX
      A=[[a1;zeros(2,3*n)];[kron(eye(3),eye(n-1)),zeros(3*n-3,3)]];
      B=zeros(3*n,2);B(1,2)=1;B(2,1)=1;B(3,2)=1;
      if fig==6, B(2,1)=0.5;A(2,2)=0.5;end
   end
   R=dlyap(A,B*B');
 
   om=logspace(-2,pi,50);
   for kk=1:length(om)
       omk=om(kk);
       if fig>4,krons=3;else krons=2;end
       W=kron(exp(-i*omk*[1:n]),eye(krons));
       ei=exp(-i*omk);
       M=[-[b(2)*ei+b(3)*ei*ei]/((1+a(2)*ei+a(3)*ei*ei).^2);...
	     1/(1+a(2)*ei+a(3)*ei*ei)];
       if fig>4 M=[M;0];end
 
        Varg(kk)=real(M.'*W*pinv(R)*W'*conj(M));
   end
   if fig==1&n==2
        baseres=Varg;
   end
   if n==2, line(om,(Varg/n)), else line(om,(Varg/n),'linestyle','--'),end

end
switch fig
   case {1,5,6}, maxa=1000;tick=2;
   case 2, maxa=10,tick=0;
   case {3,4}, maxa = 10^6, tick=5;
end
axis([0.01 pi 0.1 maxa])
set(gca,'xtick',[0.01 0.1 1],'xticklabel',[0.01 0.1 1])
if fig==2
   set(gca,'ytick',[sqrt(0.1) 1 sqrt(10) ],'yticklabel',[-0.5 0 0.5])
else
   set(gca,'ytick',[1 maxa/10 ],'yticklabel',[0 tick])
end
 


end
