% Identificazione e Analisi dei Dati: Esercitazione #1

t=[0:100]';
N=size(t,1);
s=cos(0.04*pi*t); % segnale informazione
e=randn(N,1);     % rumore bianco, varianza=1

% il meccanismo di generazione del rumore:
Nv=[1];
Dv=[1 -0.7];
% il sensore:
Nm=[1];
Dm=[1 0.5 -0.36];

v=filter(Nv,Dv,e);     % rumore additivo che si somma al segnale
m=filter(Nm,Dm,e);     % segnale misurato dal sensore
z=s+v;                 % informazione corrotta da rumore

figure,
plot(t,v,'m',t,m,'g');
title('v (mag) e m (ver)');

%%%%%%%%%%%%
% Filtro FIR
%%%%%%%%%%%%

% ordine del filtro FIR
F=4;

% Dalla teoria: la covarianza incrociata di v ed m coincide 
% con quella di z ed m: R_vm(tau)=R_zm(tau)

% calcolo delle covarianze
R=covf([m z],F+1);
Rm=R(1,:)';     % covarianza di m
Rzm=R(3,:)';    % covarianza incrociata tra z e m

% i campioni della risposta impulsiva del filtro FIR 
g=inv(toeplitz(Rm))*Rzm;

vhat=filter(g',1,m);  % stima del rumore...
shat=z-vhat;          % ... e segnale ottenuto dopo la compensazione

figure,
plot(t,v,'m',t,vhat,'c');
title('v (mag) e vhat (cia)');

figure,
plot(t,s,'b',t,shat,'k',t,z,'r');
title('s (blu), z (ros) e shat (ner)');

% calcolo dell'EQM sperimentale
covv=cov(v);
EQMs=norm(v-vhat)^2/N;  % EQM sperimentale

% miglioramento ottenuto nella ricezione del segnale (in dB)
raz=10*log10(norm(z-s)^2/norm(s-shat)^2);

% -----------------------------------------------------
% Esercizio: Provare con diversi ordini del filtro FIR!
% -----------------------------------------------------

