% confronto tra processi MA(1) ed AR(1), nel tempo, con la funzione di
% autocorrelazione e con gli spettri
%
% 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
% i processi AR(1) sono quelli analizzati negli esempi in L7-pp. 71-74 del
% materiale del corso

clc
close all
clear all


% in ingresso: rumore bianco di varianza lambda^2
lambda = 1;

% numero di campioni da simulare
N = 50000; % provare poi con 1500 <<-- notare come si modificano i grafici
          % dell'autocorrelazione e dello spettro
eta = lambda .* randn(1,N); 
% N campioni di rumore bianco, con pdf gaussiana a varianza lambda^2

% primo esempio: L7-pp. 57-60
c = 0.5;
eta1 = [0 eta(1:end-1)]; % sequenza d'ingresso ritardata di 1 passo
v = eta + c .* eta1; % v.a. in uscita

[gamma_v lags_v] = xcorr(v,31,'unbiased'); % autocorrelazione della v.a.

hf = figure;plot(v,'db-','MarkerSize',6, 'MarkerFaceColor','b','LineWidth',1.0);grid on;
title('realizzazione di un processo stazionario MA(1)');
xlabel('campioni');ylabel('v.a.');zoom on
legend('MA(1) c = +0.5');

hst = figure;stem(lags_v, gamma_v,'LineWidth',2);grid on;
xlabel('\tau');ylabel('\gamma(\tau)');legend('MA(1) c = +0.5');
title('realizzazione di un processo stazionario MA(1)');


spettro_1 = fftshift(fft(gamma_v));
L = length(spettro_1);
pulsazioni = linspace(-pi,pi,L);
hgf = figure;plot(pulsazioni, abs(spettro_1),'LineWidth',2);
xlabel('pulsazione \omega');ylabel('| \Gamma(j \omega) |');
legend('MA(1) c = +0.5');
title('realizzazione di un processo stazionario MA(1)');

disp('...tasto'); pause



%====================================================================
% secondo esempio: L7-pp. 57-60
c2 = -0.5;
eta2 = [0 eta(1:end-1)]; % sequenza d'ingresso ritardata di 1 passo
v2 = eta + c2 .* eta2; % v.a. in uscita

[gamma_v2 lags_v2] = xcorr(v2,31,'unbiased'); % autocorrelazione della v.a.

hf2 = figure;plot(v2,'sr-','MarkerSize',6, 'MarkerFaceColor','r','LineWidth',1.0);grid on;
title('realizzazione di un processo stazionario MA(1)');
xlabel('campioni');ylabel('v.a.');legend('MA(1) c = -0.5');


figure;stem(lags_v2, gamma_v2,'r','LineWidth',2);grid on;
xlabel('\tau');ylabel('\gamma(\tau)');legend('MA(1) c = -0.5');

spettro_2 = fftshift(fft(gamma_v2));
L2 = length(spettro_2);
pulsazioni2 = linspace(-pi,pi,L2);
hgf2 = figure;plot(pulsazioni2, abs(spettro_2),'r','LineWidth',2);
xlabel('pulsazione \omega');ylabel('| \Gamma(j \omega) |');
title('realizzazione di un processo stazionario MA(1)');
legend('MA(1) c = -0.5');

clc
disp('...tasto');
pause
figure(hf);hold on;
plot(v2, 'sr-','MarkerSize',6, 'MarkerFaceColor','r','LineWidth',1.0);
legend('MA(1) c = +0.5', 'MA(1) c = -0.5');

figure(hst);hold on; stem(lags_v2, gamma_v2,'r','LineWidth',2)
figure(hgf);hold on;plot(pulsazioni2, abs(spettro_2),'r','LineWidth',2);
legend('MA(1) c = +0.5', 'MA(1) c = -0.5');

disp('...tasto');pause;

%====================================================================
% ora AR(1) - L7-pp. 71-74
% uso la medesima sequenza di rumore bianco
a1 = 0.5;
% devo calcolare ricorsivamente l'uscita del sistema AR(1)
v_ar1 = zeros(size(eta));
for ii = 2:N
    v_ar1(ii) = a1 * v_ar1(ii-1)+eta(ii);
end % for ii

hfa1 = figure;plot(v_ar1,'*g-','MarkerSize',6, 'MarkerFaceColor','g','LineWidth',1 );grid on;
title('realizzazione di un processo stazionario AR(1)');
xlabel('campioni');ylabel('v.a.');legend('AR(1) a = +0.5');

[gamma_va1 lags_va1] = xcorr(v_ar1,31,'unbiased'); % autocorrelazione della v.a.
hsta1 = figure;stem(lags_va1, gamma_va1, 'g', 'LineWidth', 2);grid on;
xlabel('\tau');ylabel('\gamma(\tau)');legend('AR(1) a = +0.5');

spettro_a1 = fftshift(fft(gamma_va1));
La1 = length(spettro_a1);
pulsazionia1 = linspace(-pi,pi,La1);
hgfa1 = figure;plot(pulsazionia1, abs(spettro_a1),'g', 'LineWidth',2);
xlabel('pulsazione \omega');ylabel('| \Gamma(j \omega) |');legend('AR(1) a = +0.5');

disp('...tasto');pause;

%====================================================================
% ora secondo AR(1) - L7-pp. 71-74
% uso la medesima sequenza di rumore bianco

a2 = -0.5;
% devo calcolare ricorsivamente l'uscita del sistema AR(1)
v_ar2 = zeros(size(eta));
for ii = 2:N
    v_ar2(ii) = a2 * v_ar2(ii-1)+eta(ii);
end % for ii

hfa2 = figure;plot(v_ar2,'vc-','MarkerSize',6, 'MarkerFaceColor','c','LineWidth',1  );grid on;
title('realizzazione di un processo stazionario AR(1)');
xlabel('campioni');ylabel('v.a.');legend('AR(1) a = -0.5');

[gamma_va2 lags_va2] = xcorr(v_ar2,31,'unbiased'); % autocorrelazione della v.a.
hsta2 = figure;stem(lags_va2, gamma_va2, 'c', 'LineWidth', 2);grid on;
xlabel('\tau');ylabel('\gamma(\tau)');legend('AR(1) a = -0.5');

spettro_a2 = fftshift(fft(gamma_va2));
La2 = length(spettro_a2);
pulsazionia2 = linspace(-pi,pi,La2);
hgfa2 = figure;plot(pulsazionia2, abs(spettro_a2),'c','LineWidth',2);
xlabel('pulsazione \omega');ylabel('| \Gamma(j \omega) |');legend('AR(1) a = -0.5');

disp('confronti finali');pause
clc
figure(hf);hold on;
plot(v_ar1,'*g-','MarkerSize',6, 'MarkerFaceColor','g','LineWidth',1  );
plot(v_ar2,'vc-','MarkerSize',6, 'MarkerFaceColor','c','LineWidth',1  );
title('confronti finali: realizzazioni');
legend('MA(1) c = 1/2',...
    'MA(1) c = -1/2', 'AR(1) a = 1/2',...
    'AR(1) a = -1/2');

figure(hst);hold on;stem(lags_va1, gamma_va1,'g','LineWidth',2);
stem(lags_va2, gamma_va2,'c','LineWidth',2);
title('confronti finali: autocorrelazione');
legend('MA(1) c = 1/2',...
    'MA(1) c = -1/2', 'AR(1) a = 1/2',...
    'AR(1) a = -1/2');

figure(hgf);hold on; plot(pulsazionia1, abs(spettro_a1),'g','LineWidth',2);
plot(pulsazionia2, abs(spettro_a2),'c','LineWidth',2);
title('confronti finali: spettro');
legend('MA(1) c = 1/2',...
    'MA(1) c = -1/2', 'AR(1) a = 1/2',...
    'AR(1) a = -1/2');
grid on;