% Tema 1: Identificarea instalației ASTANK2

% **Încărcarea datelor**
load('Lab4_date_ASTANK2.mat'); % Încarcă datele din fișier (structura Did conține y (ieșiri) și u (intrări))

% **Extragerea datelor de intrare și ieșire**
y1_ident = Did.y(:,1); % Ieșire sistem canal 1
y2_ident = Did.y(:,2); % Ieșire sistem canal 2
u1_ident = Did.u(:,1); % Intrare sistem canal 1
u2_ident = Did.u(:,2); % Intrare sistem canal 2

% **Crearea obiectelor IDDATA pentru modelele ARX și ARMA**
D1 = iddata(Did.y(:,1), [Did.u(:,1) Did.u(:,2)]); % Obiect IDDATA pentru canalul 1
D2 = iddata(Did.y(:,2), [Did.u(:,1) Did.u(:,2)]); % Obiect IDDATA pentru canalul 2

% **Normalizarea ieșirilor (pentru comparații corecte)**
y1_norm = norm(detrend(Did.y(:,1), 'constant')); % Norma ieșirii canalului 1
y2_norm = norm(detrend(Did.y(:,2), 'constant')); % Norma ieșirii canalului 2

% **Generarea seturilor aleatoare de parametri pentru modelele ARX**
index1ARX = randi([1 30],20,3); % Parametrii ARX pentru canalul 1: aleatori între 1 și 30
index2ARX = randi([1 30],20,3); % Parametrii ARX pentru canalul 2

% **Generarea seturilor aleatoare de parametri pentru modelele ARMA**
index1ARMA = randi([1 60],10,2); % Parametrii ARMA pentru canalul 1: aleatori între 1 și 60
index2ARMA = randi([1 60],10,2); % Parametrii ARMA pentru canalul 2

% **Inițializarea variabilelor pentru modelele și rezultate**
matrice1ARX = cell(20, 1); % Modelele ARX canal 1
matrice2ARX = cell(20, 1); % Modelele ARX canal 2
matrice1ARMA = cell(10, 1); % Modelele ARMA canal 1
matrice2ARMA = cell(10, 1); % Modelele ARMA canal 2

v1 = cell(20,1); % Reziduurile ieșire reală - ieșire simulată pentru canal 1
v2 = cell(20,1); % Reziduurile ieșire reală - ieșire simulată pentru canal 2

e1 = cell(10,1); % Erorile pentru modelele ARMA canal 1
e2 = cell(10,1); % Erorile pentru modelele ARMA canal 2

% Variabile pentru calitatea modelelor (SNR)
snr1ARX = zeros(20,1); % SNR pentru ARX canal 1
snr2ARX = zeros(20,1); % SNR pentru ARX canal 2
snr1ARMA = zeros(10,1); % SNR pentru ARMA canal 1
snr2ARMA = zeros(10,1); % SNR pentru ARMA canal 2

ys1ARX = cell(20,1); % Ieșiri simulate ARX canal 1
ys2ARX = cell(20,1); % Ieșiri simulate ARX canal 2
ys1ARMA = cell(10,1); % Ieșiri simulate ARMA canal 1
ys2ARMA = cell(10,1); % Ieșiri simulate ARMA canal 2

% **Bucla principală pentru generarea modelelor ARX**
for i = 1:20
   % Încearcă să creeze modelul ARX pentru canalul 1
   try
       matrice1ARX{i} = armax(D1, [index1ARX(i, 1) [index1ARX(i, 2) index1ARX(i, 3)] 0 [1 1]]);
   catch
       matrice1ARX{i} = []; % Dacă apare eroare, salvează un model gol
   end

   % Încearcă să creeze modelul ARX pentru canalul 2
   try
       matrice2ARX{i} = armax(D2, [index2ARX(i, 1) [index2ARX(i, 2) index2ARX(i, 3)] 0 [1 1]]);
   catch
       matrice2ARX{i} = [];
   end

   % Determinăm lungimea maximă a coeficienților
   max1 = max(index1ARX(i,:));
   max2 = max(index2ARX(i,:));

   % Salvează primele valori simulate (folosite pentru inițializare)
   ys1ARX{i} = y1_ident(1:max1);
   ys2ARX{i} = y2_ident(1:max2);

   % Calculează ieșirea simulată pe baza coeficienților (dacă modelul e valid)
   if ~isempty(matrice1ARX{i})
       a_matrice1 = matrice1ARX{i}.a(2:end);
       b1_matrice1 = matrice1ARX{i}.b{1}(2:end);
       b2_matrice1 = matrice1ARX{i}.b{2}(2:end);

       ys = ys1ARX{i};
       for j = (max1 + 1):length(y1_ident)
           yst = b1_matrice1 * flip(u1_ident((j-length(b1_matrice1)):(j-1))) + ...
                 b2_matrice1 * flip(u2_ident((j-length(b2_matrice1)):(j-1))) - ...
                 a_matrice1 * flip(ys((end-length(a_matrice1)+1):end));
           ys = [ys; yst];
       end
       ys1ARX{i} = ys; % Salvează rezultatul simulării
   else
       ys1ARX{i} = [];
   end
   v1{i} = y1_ident - ys1ARX{i}; % Calculul reziduurilor pentru canalul 1

   % Proces similar pentru canalul 2
   if ~isempty(matrice2ARX{i})
       a_matrice2 = matrice2ARX{i}.a(2:end);
       b1_matrice2 = matrice2ARX{i}.b{1}(2:end);
       b2_matrice2 = matrice2ARX{i}.b{2}(2:end);

       ys = ys2ARX{i};
       for j = (max2 + 1):length(y2_ident)
           yst = b1_matrice2 * flip(u1_ident((j-length(b1_matrice2)):(j-1))) + ...
                 b2_matrice2 * flip(u2_ident((j-length(b2_matrice2)):(j-1))) - ...
                 a_matrice2 * flip(ys((end-length(a_matrice2)+1):end));
           ys = [ys; yst];
       end
       ys2ARX{i} = ys;
   else
       ys2ARX{i} = [];
   end
   v2{i} = y2_ident - ys2ARX{i}; % Calculul reziduurilor pentru canalul 2
end

% **Selectarea celui mai bun model ARX pe baza SNR**
for i = 1:20
    if ~isempty(ys1ARX{i})
        snr1ARX(i) = 20 * log10(y1_norm / norm(v1{i}));
    end
    if ~isempty(ys2ARX{i})
        snr2ARX(i) = 20 * log10(y2_norm / norm(v2{i}));
    end
end

[MSnr1, poz1ARX] = max(snr1ARX); % Modelul ARX optim pentru canalul 1
[MSnr2, poz2ARX] = max(snr2ARX); % Modelul ARX optim pentru canalul 2

% **Generarea modelelor ARMA pe reziduurile ARX**
D1ARMA = iddata(v1{poz1ARX}); % Datele reziduale pentru canalul 1
D2ARMA = iddata(v2{poz2ARX}); % Datele reziduale pentru canalul 2

% **Generarea modelelor ARMA pentru reziduurile ARX**
for i = 1:10
   % Încearcă să creeze un model ARMA pentru canalul 1
   try
       matrice1ARMA{i} = armax(D1ARMA, [index1ARMA(i, 1) index1ARMA(i, 2)]); % Model ARMA pentru canalul 1
   catch
       matrice1ARMA{i} = []; % Dacă apare o eroare, salvează model gol
   end

   % Încearcă să creeze un model ARMA pentru canalul 2
   try
       matrice2ARMA{i} = armax(D2ARMA, [index2ARMA(i, 1) index2ARMA(i, 2)]); % Model ARMA pentru canalul 2
   catch
       matrice2ARMA{i} = []; % Dacă apare o eroare, salvează model gol
   end

   % Simulează ieșirile pentru ARMA canal 1 și calculează erorile
   if ~isempty(matrice1ARMA{i})
       ys1ARMA{i} = predict(matrice1ARMA{i}, D1ARMA); % Prevede valorile pentru canalul 1
       e1{i} = v1{poz1ARX} - ys1ARMA{i}.y; % Calculul erorilor (reziduu - predicție)
   else
       e1{i} = []; % Dacă nu există model, lasă gol
   end

   % Simulează ieșirile pentru ARMA canal 2 și calculează erorile
   if ~isempty(matrice2ARMA{i})
       ys2ARMA{i} = predict(matrice2ARMA{i}, D2ARMA); % Prevede valorile pentru canalul 2
       e2{i} = v2{poz2ARX} - ys2ARMA{i}.y; % Calculul erorilor (reziduu - predicție)
   else
       e2{i} = [];
   end
end

% **Selectarea celui mai bun model ARMA pe baza SNR**
for i = 1:10
    % Calculul SNR pentru canalul 1
    if ~isempty(e1{i})
        snr1ARMA(i) = 20 * log10(norm(v1{poz1ARX}) / norm(e1{i}));
    end

    % Calculul SNR pentru canalul 2
    if ~isempty(e2{i})
        snr2ARMA(i) = 20 * log10(norm(v2{poz2ARX}) / norm(e2{i}));
    end
end

% Determinarea celor mai bune modele ARMA
[MSnr1ARMA, poz1ARMA] = max(snr1ARMA); % Cel mai bun model ARMA canal 1
[MSnr2ARMA, poz2ARMA] = max(snr2ARMA); % Cel mai bun model ARMA canal 2

% **Combinarea modelelor ARX și ARMA**
Yfinal1 = ys1ARX{poz1ARX} + ys1ARMA{poz1ARMA}.y; % Ieșirea finală simulată canal 1
Yfinal2 = ys2ARX{poz2ARX} + ys2ARMA{poz2ARMA}.y; % Ieșirea finală simulată canal 2

% **Vizualizări**

% a) Intrările și ieșirile măsurate pentru fiecare canal
figure;
subplot(2,1,1);
plot(u1_ident);
title('Date de intrare - Canal 1');
xlabel('Timp');
ylabel('Amplitudine');
subplot(2,1,2);
plot(y1_ident);
title('Date de ieșire - Canal 1');
xlabel('Timp');
ylabel('Amplitudine');

figure;
subplot(2,1,1);
plot(u2_ident);
title('Date de intrare - Canal 2');
xlabel('Timp');
ylabel('Amplitudine');
subplot(2,1,2);
plot(y2_ident);
title('Date de ieșire - Canal 2');
xlabel('Timp');
ylabel('Amplitudine');

% b) Ieșirea simulată a modelului ARX optim, comparată cu ieșirea măsurată
figure;
plot(ys1ARX{poz1ARX});
hold on;
plot(y1_ident);
legend('ARX simulată', 'Ieșire măsurată');
title('Comparare ARX Canal 1');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(y1_ident), ['SNR ARX1 = ', num2str(snr1ARX(poz1ARX))]);
hold off;

figure;
plot(ys2ARX{poz2ARX});
hold on;
plot(y2_ident);
legend('ARX simulată', 'Ieșire măsurată');
title('Comparare ARX Canal 2');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(y2_ident), ['SNR ARX2 = ', num2str(snr2ARX(poz2ARX))]);
hold off;

% c) Ieșirea simulată a modelului ARMA optim, comparată cu zgomotul rezidual
figure;
plot(ys1ARMA{poz1ARMA}.y);
hold on;
plot(v1{poz1ARX});
legend('ARMA simulată', 'Zgomot colorat');
title('Comparare ARMA Canal 1');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(v1{poz1ARX}), ['SNR ARMA1 = ', num2str(snr1ARMA(poz1ARMA))]);
hold off;

figure;
plot(ys2ARMA{poz2ARMA}.y);
hold on;
plot(v2{poz2ARX});
legend('ARMA simulată', 'Zgomot colorat');
title('Comparare ARMA Canal 2');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(v2{poz2ARX}), ['SNR ARMA2 = ', num2str(snr2ARMA(poz2ARMA))]);
hold off;

% d) Ieșirea totală simulată (ARX + ARMA), comparată cu ieșirea măsurată
figure;
plot(Yfinal1);
hold on;
plot(y1_ident);
legend('Simulare totală', 'Ieșire măsurată');
title('Simulare Totală ARX+ARMA Canal 1');
xlabel('Timp');
ylabel('Amplitudine');
hold off;

figure;
plot(Yfinal2);
hold on;
plot(y2_ident);
legend('Simulare totală', 'Ieșire măsurată');
title('Simulare Totală ARX+ARMA Canal 2');
xlabel('Timp');
ylabel('Amplitudine');
hold off;

% e) Zgomot alb rezidual
figure;
plot(e1{poz1ARMA});
title('Zgomot alb rezidual Canal 1');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(e1{poz1ARMA}), ['\lambdae1^2 = ', num2str(norm(e1{poz1ARMA}))]);

figure;
plot(e2{poz2ARMA});
title('Zgomot alb rezidual Canal 2');
xlabel('Timp');
ylabel('Amplitudine');
text(500, max(e2{poz2ARMA}), ['\lambdae2^2 = ', num2str(norm(e2{poz2ARMA}))]);