% program ex_id_3.m to solve problem 4 in the Test 2 of Che532 
% programmed by Kartik Arora April 28, 2001
%
% get prbs input signal covering freq from .1 to 10

load ch19pr1_t_d_x_y.dat
data=ch19pr1_t_d_x_y;
time=data(:,1);u=data(:,2);x=data(:,3);y=data(:,4);
% check the frequency content of the input signal

len=fix(length(y)/2);

u1=u(1:len);
x1=x(1:len);
y1=y(1:len);
ut=u(len+1:2*len);
xt=x(len+1:2*len);
yt=y(len+1:2*len);


z=fft(u1);
% compute the frequency axis
%
samp_time=.1;
len=length(z);
wlen=ceil(len/2); % use only first half of transform;
w=(1:wlen)*2*pi/samp_time/len; % compute actual frequencies
zm=abs(z)/len; % compute the magnitude of the fourier transform
% capture the first half of the fourier transform
zm1=zm(2:wlen+1);
% display the signal for viewing
stairs(u); ax=axis;ax(3)=2*ax(3);ax(4)=2*ax(4);axis(ax);
disp( 'the prbs generated is shown in graph');
disp( ' hit any key when done viewing')
semilogx(w,zm1);

pause

delay=11;



% now fit some models to this data

%
% create z vector for identification
% the following requires Identification Toolbox in Matlab



z1=[x1 u1];% display input output data
th1=arx(z1,[1 1 delay]);
th2=arx(z1,[2 2 delay]);
th3=arx(z1,[3 3 delay]);
th4=arx(z1,[4 4 delay]);
th5=arx(z1,[5 5 delay]);
y1sim=idsim(u1,th1);
y2sim=idsim(u1,th2);
y3sim=idsim(u1,th3);
y4sim=idsim(u1,th4);
y5sim=idsim(u1,th5);
pe1=y1-y1sim;pe1norm=sqrt(pe1'*pe1)
pe2=y1-y2sim;pe2norm=sqrt(pe2'*pe2)
pe3=y1-y3sim;pe3norm=sqrt(pe3'*pe3)
pe4=y1-y4sim;pe4norm=sqrt(pe4'*pe4)
pe5=y1-y5sim;pe5norm=sqrt(pe5'*pe5)
y1sim=idsim(ut,th1);
y2sim=idsim(ut,th2);
y3sim=idsim(ut,th3);
y4sim=idsim(ut,th4);
y5sim=idsim(ut,th5);
pe1=yt-y1sim;pe1norm=sqrt(pe1'*pe1)
pe2=yt-y2sim;pe2norm=sqrt(pe2'*pe2)
pe3=yt-y3sim;pe3norm=sqrt(pe3'*pe3)
pe4=yt-y4sim;pe4norm=sqrt(pe4'*pe4)
pe5=yt-y5sim;pe5norm=sqrt(pe5'*pe5)

% building an estimator for y using x
th1=arx([y x],[0 1 0]);
present (th1);
[num,den]=th2tf(th1);sys=tf(num,den,.1)
sysc=d2c(sys)