% Example ID3.3 Illustrate FIR model estimation
clc;echo on
% Copyright 1996, Babu Joseph, Washington University in St.Louis
%
% This program illustrates the use of FIR models to identify a process
%
% First create some data by exciting a process
% We will use a first order process 
% with gain 1 and time constant 5
%
num=[1];den=[5 1 ];
[numd,dend]=c2dm(num,den,1,'zoh');
% Define input as a random  binary sequence
N=100; % no of data points
u=sign(randn(N,1));
%
% Generate and plot response data with noise added
%
y=dlsim(numd,dend,u)+.2* randn(N,1) ;
idplot([y u]); pause
clg;clc;
% Now generate the matrices for least square estimation of 
% FIR ( finite impulse response coefficients) 5 coeffs are considered

nb=15  ;% no of impulse coeffs

% Generate matrices for identification, see Text

A=zeros(N-nb+1,nb);
for i=1:nb 
  A(:,i)=u(nb-i+1:N-i+1);
end;
b=y(nb:N);

%

% The estimated coefficients of the transfer function are
%
h=A\b;h'
e=b-A*h;   % error in prediction
norm(e)^2  % sum of the square of the errors

% The condition no of the matrix is
cond(A)

%
% The following are the actual coefficients
%
nu=zeros(nb,1);nu(1)=1;% This defines a pulse input
xactual=dlsim(numd,dend,nu);xactual'
%
% Computed coefficients can be compared to actual coefficients
% see graph
clg;plot(xactual);hold on; plot( h);
pause 
% 
% Next we compare how good the model is in predicting the actual response
%
% Next let us fit this data to a first order transfer function model 
%
A=[-y(1:N-1) u(1:N-1)];b=y(2:N);
%
% The estimated coefficients of the transfer function are
%
x=A\b
%
% The following are the actual coefficients
[numd,dend]
%
% The sum of the square of errors  is given by
%
e=b-A*x;norm(e)^2;
%
% The sum of the square of errors using FIR model is given by
%
%
% create some new data for testing
%
u=sign(randn(N,1));
y=dlsim(numd,dend,u)+.0* randn(N,1) ;
% estimate sum of squares of error in prediction
%
A=zeros(N-nb+1,nb);
for i=1:nb
  A(:,i)=u(nb-i+1:N-i+1);
end;
yp1=A*h;ei=y(nb:N)-yp1;norm(e)/sqrt((N-nb))
%
%Now estimate error in prediction using first order transfer fn model
%
A=[-y(1:N-1) u(1:N-1)];b=y(2:N);
ef=b-A*x;norm(e)/sqrt(N)
efs=ef(nb-1:N-1);
clg;hold on;plot(efs);plot(ei);
%