% confronto tra campioni di una realizzazione di processo bianco (WN) e di 
% un processo MA(1), nel tempo e tramite funzione di correlazione
%
% corso di Sistemi Dinamici
% a.a. 2021/2022 G. Fenu - UniTS
%
% i processi MA(1) sono quelli analizzati negli esempi in L7-pp. 57-60 del
% materiale del corso


close all
clear
clc

noise_var = 1e-1; % varianza del rumore bianco
N = 1200; % numero di campioni per ciascun processo

WN = zeros(N, 1); % prealloco un vettore per i campioni del processo bianco
MA1 = zeros(N, 1); % prealloco per il processo MA(1)
MA2 = zeros(N, 1); % prealloco per il processo MA(1)

c1 = 0.98; c11 = -0.98; 
% per avere gli esempi di L7 pp. 57-60 porre
% c1 = +0.5; c11 = -0.5;

WN0 = sqrt(noise_var)*randn(1,1); % un primo valore che uso come valore 
% all'istante "-1" nel'inizializzazione della equazione per il processo MA

WN = sqrt(noise_var)*randn(size(WN));

% MA1(1)=WN0; MA2(1)=MA1(1); % inizializzo
%  for ii=2:N
%      MA1(ii) = WN(ii)+c1*WN(ii-1);
%      MA2(ii) = WN(ii)+c11*WN(ii-1);
%  end % ii
 
 % in forma vettoriale
 WN1 = [WN0; WN(1:end-1)]; % memorizzo i campioni del processo bianco con un passo di ritardo
 % NB 
 
 MA1 = WN + c1*WN1;
 MA2 = WN +c11*WN1;

 
 figure('Name', 'realizzazioni di WN e processi MA(1)');
 % campioni di rumore bianco
 subplot(4,1,1);
 plot(WN,'bo-', 'Linewidth', 2);
 grid on;h1=gca;legend('WN');hz1 = zoom(h1);
 setAxesZoomMotion(hz1,h1,'horizontal');hz1.Enable = 'On';
 title('realizzazione di WN');

 % ora campioni del processo MA(1) con c = +0.5
 subplot(4,1,2);
 plot(MA1,'rd-', 'Linewidth', 2);
 grid on;h2=gca;legend('MA(1) c>0');hz2 = zoom(h2);
 setAxesZoomMotion(hz2,h2,'horizontal');hz2.Enable = 'On';
 title(['realizzazione di MA(1) con c_1=', num2str(c1)]);


 % ed ora invece il processo MA(1) con c = -0.5
 subplot(4,1,3);
 plot(MA2,'ks-', 'Linewidth', 2);
 grid on;h3=gca;legend('MA(1) c<0');hz3 = zoom(h3);
 setAxesZoomMotion(hz3,h3,'horizontal');hz3.Enable = 'On';
 title(['realizzazione di MA(1) con c_1=', num2str(c11)]);


 % ora tutti sovrapposti, per cercare di "vedere le differenze"
 subplot(4,1,4);h4=gca;grid on;hold on;
 plot(WN,'bo-', 'LineWidth', 1.0, 'MarkerSize', 6);
 plot(MA1,'rd-', 'LineWidth', 1.0, 'MarkerSize', 6);
 plot(MA2,'ks-', 'LineWidth', 1.0, 'MarkerSize', 6);
 legend('White Noise', ['MA(1) c = ', num2str(c1)], ...
     ['MA(1) c = ',num2str(c11)]);
 title('tutte le realizzazzioni')
 % li associo, per gestire lo zoom su tutti i grafici assieme
 hz4 = zoom(h4);
 setAxesZoomMotion(hz4,h4,'horizontal');hz4.Enable = 'On';
 linkaxes([h1,h2,h3,h4],'x');

 
 % le autocorrelazioni dei 3 processi
figure; subplot(3,1,1);
[c_wn,wnlags] = xcorr(WN,10,'unbiased');
stem(wnlags,c_wn, 'LineWidth', 2);grid on;h1=gca;
title('autocorrelazione del processo WN')

 subplot(3,1,2);
 [c_ma1,ma1lags] = xcorr(MA1,10,'unbiased');
stem(ma1lags,c_ma1, 'LineWidth', 2);grid on;h2=gca;
title(['autocorrelazione del processo MA(1) con c=', num2str(c1)])
  
subplot(3,1,3);
   [c_ma2,ma2lags] = xcorr(MA2,10,'unbiased');
stem(ma2lags,c_ma2, 'LineWidth', 2);grid on;h3=gca;
title(['autocorrelazione del processo MA(1) con c=', num2str(c11)])

 linkaxes([h1,h2,h3],'x');