% Identificazione e Analisi dei Dati: Esercitazione #2
% Identificazione parametrica 

% Sono necessari i sistemi simulink contenuti nei files:
% system1.mdl
% system1n.mdl

% PRIMA PARTE
% Si consideri il sistema system1.mdl
% 1) generare una sequenza di dati ingresso-uscita
% 2) applicare algoritmi di identificazione parametrica

%%%%%%%%%%%%%%%%%%%%%%%%
% Esempio di svolgimento
%%%%%%%%%%%%%%%%%%%%%%%%

% 1) GENERAZIONE DATI

% dati
m1=1;
b1=0.4;
k1=1;
m2=1;
b2=0.2;
k2=4;
sigma=0.1; % varianza del rumore additivo in system1n

% la vera funzione di trasferimento (tempo-continuo)
num=[b2 k2];
kt=k1+k2;
bt=b1+b2;
den=[m1*m2 bt*m2+b2*m1 b1*b2+k2*m1+kt*m2 k2*b1+k1*b2 k1*k2];
G=tf(num,den);

% discretizzazione
Tc=0.1;  % tempo di campionamento
Gd=c2d(G,Tc,'zoh');
[numd,dend]=tfdata(Gd);

figure, bode(G);

% generazione dell'ingresso
Tf=50;
T=[0:Tc:Tf]';
B=1;  % banda (normalizzata) del segnale di ingresso
U=idinput(size(T,1),'PRBS',[0 B]);
U=U(1:size(T,1));  % necessario se B e' minore di 1

% lancia la simulazione sul sistema simulink 
tout=sim('system1',Tf);        

Y=Y(1:size(T,1));  % uscita
z=iddata(Y,U,Tc);     % dati i/o (in formato iddata)
%z=[Y U];   % dati i/o
% plot dei dati
figure, idplot(z)

% 2) IDENTIFICAZIONE PARAMETRICA

% modelli arx
th1=arx(z,[2 2 1]);
figure, compare(z,th1,1);
%% oppure, per mettere l'uscita in un vettore...
%[yp1,fit1]=compare(z,th1,1);
%yp1=yp1{1};
%figure, plot(T,Y,'b',T,yp1.y,'g')

th2=armax(z,[4 4 3 1]);
figure, compare(z,th2,1);
[yp2,fit2]=compare(z,th2,1);

th3=oe(z,[4 4 1],'trace');
figure, compare(z,th3,1);
[yp3,fit3]=compare(z,th3,1);

% residui
figure, resid(z,th1);
figure, resid(z,th2);
figure, resid(z,th3);

% confronto in frequenza
ww=logspace(-2,log10(3.14/Tc),256);
[magd,phad]=bode(Gd,ww);  magd=squeeze(magd); phad=squeeze(phad);
figure, semilogx(ww,20*log10(magd),'b'), grid    % la vera fdt
hold on

g1=idfrd(th1,ww);
g1v=squeeze(g1.ResponseData);
semilogx(ww,20*log10(abs(g1v)),'r')
g2=idfrd(th2,ww);
g2v=squeeze(g2.ResponseData);
semilogx(ww,20*log10(abs(g2v)),'g')
g3=idfrd(th3,ww);
g3v=squeeze(g3.ResponseData);
semilogx(g3.Frequency,20*log10(abs(g3v)),'y')

% validazione su set di dati diverso
zid=z(1:400,:);
%zva=z(401:size(T,1),:);
thnew=bj(zid,[4 2 2 4 1]);
figure, compare(z,thnew,inf,401:size(T,1))
[ypnew,fitnew]=compare(z,thnew,1,401:size(T,1));


% SECONDA PARTE
% Ripetere tutto l'esercizio, utilizzando il sistema simulink 
% 'system1n', in cui e' presente un rumore sui dati in uscita.

% Provare a considerare diverse classi di modelli e a variare per ciascuna classe:
% - il numero dei dati
% - la banda dell'ingresso
% - il tipo di ingresso




