% IDS_ARX  Identifikation mittels ARX für ein zeitdiskretes System
%          Vorlesung Identifikation dynamischer Systeme (IDS)

% M. Grobe  06/2007
% Technische Universitaet Braunschweig, IfR

%% Definition des Systems in ID-Form und als TF
% Abtastzeit
T=1;
% System
A = [1  -1.5  0.7];  % A(q)
% VL: [b2  b1   b0]
% M:  [b1  b2   b3]
B=[0 1 0.5];         % B(q)
% VL: [a2  a1   a0]
% M:  [a1  a2   a3]
Gid=idpoly(A,B);    % IDMODEL
num=B;
den=A;
Gtf=tf(num,den,T);  % TF

%% Signale generieren
% PRBS generieren
u=idinput(300,'prbs');  % als Signal
uid=iddata([],u);       % als Objekt für ID
% normalverteilte Störung generieren
r=randn(300,1);         % als Signal
rid=iddata([],r);       % als Objekt für ID

%% Simulation
% Zeitvektor
t=linspace(0,299,300)';
% Ausgang ohne Rauschen
y=lsim(Gtf,u,t);
% Ausgang mit Rauschen gefiltert durch 1/A(q) am Ausgang
yid=sim(Gid,[uid rid]);

%% Identifikation
% Messwerte / Anregung
z=[yid uid];
% Modellordnungen für Identifikation
na=2;  % Ordnung A(q)
nb=2;  % Ordnung A(q)
nk=1;  % zusätzliche Verzögerungen (default=1)
% ARX (LS)
Ghatid=arx(z,[na nb nk]);
% ARMAX
%Ghatid=arx(z,[na nb 2 nk]);
% BJ
%Ghatid=bj(z,[na nb 2 2 nk]);

%% Vergleich Sprungantwort reales System mit Identifikationsergebnis
yidstep=sim(Gid,ones(35,1));
yhatidstep=sim(Ghatid,ones(35,1));


%% Plots
figure(1)
clf
subplot(3,2,1)
stairs(yidstep)
ylabel('Sprungantwort reales System')
subplot(3,2,2)
plot(u)
ylabel('PRBS')
axis([0 300 -1.2 1.2])
subplot(3,2,3)
plot(y)
ylabel('Antwort auf PRBS')
subplot(3,2,4)
plot(r)
ylabel('Rauschen')
subplot(3,2,5)
plot(yid)
ylabel('Antwort auf PRBS + Rauschen')
title(''); xlabel('');
subplot(3,2,6)
stairs(yidstep,'b');
hold on;
stairs(yhatidstep,'r');
ylabel('Vergleich (Schätzung rot)')
