% FIGURES 8.2  and 8.3 as well as 14.2 to 14.4

m0=poly2th(1,0.001*[0 0 10 7.4 0.924 0.1764],1,1,...
    [1 -2.14 1.553 -0.4387 0.042]);
u=idinput(1000,'prbs');
y=idsim(u,m0);
z=[y u];
ma=arx(z,[2 2 1]);
mo=oe(z,[2 2 1]);
g0=th2ff(m0,1,logspace(-2,pi,100)); [om,g0a]=getff(g0);
ga=th2ff(ma,1,logspace(-2,pi,100)); [om,gaa]=getff(ga);
go=th2ff(mo,1,logspace(-2,pi,100)); [om,goa]=getff(go);
figure(1) % FIGURE 8.2
subplot(121)
% Could also use bodeplot([g0 ghh],'a');
loglog(om,goa), line(om,g0a,'linewidth',2),...
    title('a'),xlabel('frequency (rad/s)')
axis([0.01 10 0.001 10])
set(gca,'ytick',[0.01 0.1  1]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[ 0.01 0.1 1])
set(gca,'xticklabel',[ 0.01 0.1 1])
subplot(122)
loglog(om,gaa), line(om,g0a,'linewidth',2),...
    title('b'),xlabel('frequency (rad/s)')
axis([0.01 10 0.001 10])
set(gca,'ytick',[0.01 0.1  1]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[ 0.01 0.1 1])
set(gca,'xticklabel',[ 0.01 0.1 1])

% The weighting function:
a=th2poly(ma);
we=poly2th(1,[],a);
gw=th2ff(we,-1,logspace(-2,pi,100)); [om,gwa]=getff(gw);
figure(2),subplot % FIGURE 8.3
loglog(om,gwa)
axis([0.01 10 0.0001 100])
set(gca,'ytick',[0.001  0.1  10 ]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[ 0.001 0.1 10  ])
set(gca,'xticklabel',[ 0.01 0.1 1])
xlabel('frequency (rad/s)')
 
% NOW THE FILTERED VERSIONS:

zhp=idfilt(z,5,0.5/pi,'hp');
mof=oe(zhp,[2 2 1]);
[zlp,L2]=idfilt(z,5,0.5/pi);
maf=arx(zlp,[2 2 1]);

gaf=th2ff(maf,1,logspace(-2,pi,100)); [om,gafa]=getff(gaf);
gof=th2ff(mof,1,logspace(-2,pi,100)); [om,gofa]=getff(gof);
figure(3)  % FIGURE 14.2
subplot(111)
%  Could also use bodeplot([g0 ghh],'a');
loglog(om,gofa), line(om,g0a,'linewidth',2),xlabel('frequency (rad/s)')
axis([0.01 10 0.001 10])
set(gca,'ytick',[0.01 0.1  1]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[ 0.01 0.1 1])
set(gca,'xticklabel',[ 0.01 0.1 1])
 
figure(5) % FIGURE !4.3
subplot(111)
loglog(om,gafa), line(om,g0a,'linewidth',2),xlabel('frequency (rad/s)')
axis([0.01 10 0.001 10])
set(gca,'ytick',[0.01 0.1  1]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[ 0.01 0.1 1])
set(gca,'xticklabel',[ 0.01 0.1 1])

% The weighting function |L_2 A|^2:
af=th2poly(maf);[l2n,l2d]=th2tf(L2);
wef=poly2th(l2d,[],conv(af,l2n)/l2n(1));
gwf=th2ff(wef,-1,logspace(-2,pi,100)); [om,gwfa]=getff(gwf);gwfa=gwfa/1000;
print -deps fig143new.eps

figure(4),subplot % FIGURE 14.4
loglog(om,gwfa)
axis([0.01 10 0.000001 200])
set(gca,'ytick',[0.00001 0.001  0.1  10 ]),set(gca,'xtick',[0.01 0.1 1])
set(gca,'yticklabel',[0.00001 0.001 0.1 10  ])
set(gca,'xticklabel',[ 0.01 0.1 1])
xlabel('frequency (rad/s)')
 