rels_systemidentification(y,u,na,nb)

function rels_systemidentification(y,u,na,nb)
%--------------------------------------------------------------------------
%          Rekursive Erweiterte Fehlerquadrate-Identifikationsmethode                                                                                    %
%--------------------------------------------------------------------------

% -------------------------------------------------------------------------
%Prozessdifferenzengleichung:
% y(k)+a(1)*y(k-1)+...+a(na)*y(k-na) = b(0)*u(k)+b(1)*u(k-1)+b(nb)*u(k-nb)+e(k)
%
%   u: Eingang
%   y: Ausgang
%   na: Nullstellen-Anzahl
%   nb: Pole-Anzahl=Ordnung des Prozesses
%   
%   yest: Gescheatzter Ausgang
%   rela_fehler: Relativer Gleichungsfehler 
%   a = [a0 a1 ... ana]: Koeffizientenvektor vom Eingangssignal 
%   b = [b0 b1 ... bnb]: Koeffizientenvektor vom Ausgangsignal
%--------------------------------------------------------------------------

if nargin<1,   y=y;  end
if nargin<2,   u=u;   end
if nargin<3,   na=8;  end
if nargin<4,   nb=8;  end
if na<1,       na=8;  end
if nb<1,       nb=8;  end

N=1e3;
lamda0=0.1;
lamda= random('Normal',0,1,1,N);
% lamda= random('uniform',0,1,length(lamda),N);
% lamda=firls(N-1,[0 .4 .5 1],[1 1 0 0],[1 10]);

% for k=1:length(lamda),
%     lamda(k) = lamda(1)*lamda(k)+lamda(length(lamda))* (1-lamda(1));
% end;

n=na+nb;
delta = 1e3;
theta0=zeros(n);  % initialisierung von theta
P0=delta.*eye(n); %initialisierung von P
eps=zeros(n);

py=zeros(na);
pu=zeros(nb);
ymat=[y zeros(8,8);zeros(8,8) y]';

for i=n-1:-1:n-na
if i>0,
    py(n-i)=-y(i);
end
end
for i=n-1:-1:n-nb
if i>0, 
    pu(n-i)=u(i); 
end
end

phi = [py zeros(nb);zeros(na) pu];

for j=1:length(P0),
    
lamda=lamda0*lamda+1-lamda0;
x_vals=[lamda; lamda]; y_vals=[zeros(1,length(lamda)); lamda];

P =1/lamda(j)*P0*[eye(n) - P0*phi*phi'*[abs(lamda(j))+ phi'*P0*phi]^-1];
eps = [y  zeros(na);zeros(na) zeros(na)]- theta0'*phi;   
theta=theta0 + P*phi*eps;
theta0=theta;
yest=theta0'*phi;
P0=P;
end

figure,
stem(x_vals,y_vals,'fill','--');

if (P<1 | P>1e3), error('P soll passend sein');end
if (P<0 | P>100), disp ('Wahrnung:P ist im normalen Fall groesser als 1 und kleiner als unendlich');end

nu=length(u);
nmax=max(na,nb);
a1= [1; theta(1:na)'];
b1 = theta(na+1:n)';
% e1=filter(a1,1,y)-filter(b1,1,u);
% V1=(e1'*e1)/(nu-nmax);
% FPE1=V1*(1+n/nu)/(1-n/nu);
% rela_fehler = sqrt((e1'*e1)./(y'*y));

% for i=1:nmax e1(i)=0; end;    
disp('b1 Koeffizienten')
disp(b1)
disp('a1 Koeffizienten')
disp(a1)
% disp(['relativer Gleichungsfehler = ', num2str(100*rela_fehler),' %']);

% nullstell=roots(a1); % nullstelen von a1(1/z)
% 
% if all(abs(nullstell)<1)
%         sprintf('Das Prozess ist kausal')
% else
%         sprintf('Das Prozess ist nicht kausal')
% end


figure,
plot(y)
title('Systemausgang ') ;
xlabel('Abtastzeiten')
ylabel('echtes Ausgangsignal');
figure,
plot(yest,'r')
title('Schaetzerausgang') ;
xlabel('Abtastzeiten')
ylabel('geschaetztes Ausgangsignal');
figure,
semilogy((abs(eps))) ;
title('Gleichungsfehler') ;
xlabel('Abtastzeiten');
ylabel('Epsilon-Werte');
figure,
plot(theta)
title('Konvergenz de geschaetzten Koeffizienten') ;
xlabel('Abtastzeiten');
ylabel('Koeffizientenverlauf');
% axis([1  length(y)-500  min(min(theta)) max(max(theta)) ]);

end