
%
% Esercizio : costruire la funzione che implementa il generico filtro numerico:
%             function [y]=simula_DLTI(b,a,u) ...  e confrontarla con quella data:  "simula_DLTI_ARMA.m"

%%
%%  Calcoliamo l'uscita di un sistema DLTI, rappresentato mediante un'equazione alle differenze ARMA,
%%  in risposta ad un segnale sinusoidale discreto:
%%
r = 0.95;       % modulo delle radici
phi = 0.1*pi;   % fase delle radici
a = [1 -2*r*cos(phi) r^2];
b = [1 -0.1];

%
%
% Esercizio_1_a : dato il sistema DLTI-ARMA qui sopra (caratterizzato dai coefficienti "a" e "b"),
%   - calcolare e disegnare la risposta in frequenza usando "freqresp.m";
%   - scelta una sequenza di ingresso sinusoidale pura, calcolare ricorsivamente nel tempo la sequenza di uscita del sistema e verificare la corrispondenza 
%     tra amplificazione indicata nel modulo della risposta in frequenza e rapporto tra le ampiezze, rispettivamente, della sequenza sinusoidale in uscita
%     e della sequenza in ingresso;
%   - conoscendo la risposta in frequenza del sistema, determinare la frequenza del segnale in ingresso tale che 
%     l'amplificazione in uscita sia massima.
% 





% %%
% %%  calcoliamo la DFT di alcuni segnali tipici (es. una sinusoide pura) mediante algoritmo FFT
% %%
N = 256;
omega_nu = (2*pi/N) * 0.5*(N/2);
x = sin(omega_nu * [0:N-1]);
X = fft(x);
figure(3); plot(0:N-1, abs(X), 'b-');
% %
% % Esercizio : fare la FFT di una sequenza composta dalla somma di alcune sinusoidi ed una sequenza
% %             casuale e verificare la posizione delle righe dello spettro. Dove e' finito il
% %             contributo della sequenza casuale ?
% %


%%%
%%%  Calcoliamo in vari modi la risposta in frequenza del sistema meccanico con le tre masse:
%%%
N = 1024*10;
%
% costruisco il modello continuo del sistema meccanico:
[A,B,C,D]=build_sistema_meccanico_3gdl();
%
% definisco il periodo di campionamento:
Ts = 0.0001;
%
%
% calcolo analitico della risposta in frequenza (mediante trasformata zeta):
%
if 1  % possiamo eseguire questa sezione di codice solo se abbiamo il "control system toolbox"
  ss_d_3gdl = c2d(ss(A,B(:,1),C,D(:,1)), Ts);
  tf_d_3gdl = tf(ss_d_3gdl);
  [b,a] = tfdata(tf_d_3gdl,'v');
  figure(1); clf, freqz(b,a);
end
pause
%
%
%  Stima della risposta in frequenza da una sequenza di dati "sperimentali":
%
% 1. risposta al campione unitario:
u=zeros(1,N);  u(1)=1;
h = simula_sistema_meccanico_3gdl(A,B,C,D,u,Ts);
H = fft(h);
subplot(2,1,1); hold on, plot(0:2/N:1, 20*log10( abs(H(1:end/2+1)) ), 'm--');
subplot(2,1,2); hold on, plot(0:2/N:1, angle(H(1:end/2+1))*360/(2*pi), 'm--');
pause
U = fft(u)
Y = fft(h);
H = Y ./ U;
subplot(2,1,1); hold on, plot(0:2/N:1, 20*log10( abs(H(1:end/2+1)) ), 'r--');
subplot(2,1,2); hold on, plot(0:2/N:1, angle(H(1:end/2+1))*360/(2*pi), 'r--');
pause
figure(2); clf, title('DFT della sequenza delta di Kronecker')
           xlabel('Normalized Frequency')
subplot(2,1,1); hold on, plot(0:2/N:1, 20*log10( abs(U(1:end/2+1)) ), 'b-');
                ylabel('Magnitude (dB)')
subplot(2,1,2); hold on, plot(0:2/N:1, angle(U(1:end/2+1))*360/(2*pi), 'b-');
                ylabel('Phase (degrees)')
pause
%
% 2. risposta ad un ingresso sinusoidale:
omega_nu = (2*pi/N) * 0.1*(N/2);
u = cos(omega_nu * (1:N));
y = simula_sistema_meccanico_3gdl(A,B,C,D,u,Ts);
%figure(2); clf, hold on, plot(1:N/4,u(1:N/4),'g-', 1:N/4,y(1:N/4)*1e7,'b-'); title('"y" è amplificato !'); pause,
U = fft(u);           
Y = fft(y);
H = Y ./ U;
figure(1); subplot(2,1,1); hold on, plot(0:2/N:1, 20*log10( abs(H(1:end/2+1)) ), 'r--'); 
           subplot(2,1,2); hold on, plot(0:2/N:1, angle(H(1:end/2+1))*360/(2*pi), 'r--');
           title('stima (non parametrica) della risposta in frequenza'); 
subplot(2,1,1); axis([0 0.2 -190 140]);
subplot(2,1,2); axis([0 0.2 -250 0]);
pause
% --> notare che la stima ha validità nello spettro di frequenze del segnale di ingresso, che qui è una sinusoide pura
%

%%
%%  calcolo l'uscita del sistema meccanico nel dominio delle frequenze ed anti-trasformo: confronto
%%  con simulazione nel dominio del tempo
%%
U = fft(u);
Y = H .* U ;
yfft = ifft(Y);
figure(5); clf, 
           plot(N/2:N, y(N/2:N), 'b-'); hold on,
           plot(N/2:N, yfft(N/2:N), 'r--'); hold on,
           title('risposta calc nel tempo (blu) e tramite fft (red)');

%
% stima della sequenza in ingresso dalla misura della sequenza in uscita e dalla conoscenza della
% risposta in frequenza del modello (deconvoluzione in frequenza)
%
u=zeros(1,N);  u(1)=1;
h = lsim(tf_d_3gdl, u, []);
H = fft(h);
%
% ingresso sinusoidale:
omega_nu = (2*pi/N) * 0.3*(N/2);
u = (sin(omega_nu * (1:N)) + sin(1.5*omega_nu * (1:N)) + sin(2*omega_nu * (1:N)))';
y = simula_DLTI_ARMA(b, a, u);  % questa routine risolve l'equazione alle differenze vista alla "ESE_1"
figure(5); clf, hold on, plot(1:N/4,u(1:N/4),'g-', 1:N/4,y(1:N/4)*1e7,'b-'); title('"y" è amplificato !'); pause,
Y = fft(y);           
U_est = (Y) ./ (H);
u_est = real(ifft(U_est));
figure(6); clf, hold on, plot(1:N/4,u(1:N/4),'g-', 1:N/4,u_est(1:N/4),'r--'); title('sequenza di ingresso (verde) e stima (rosso)');  pause,
%
% ingresso a gradino:
u=[zeros(100,1); ones(200,1); zeros(N-300,1)]; 
y = simula_DLTI_ARMA(b, a, u);  % questa routine risolve l'equazione alle differenze vista alla "ESE_1"
figure(7); clf, hold on, plot(1:N,u,'g-', 1:N,y*1e7,'b-'); title('"y" è amplificato !'); pause,
Y = fft(y);           
U_est = (Y) ./ (H);
u_est = real(ifft(U_est));
figure(8); clf, hold on, plot(1:N,u,'g-', 1:N,u_est,'r--'); title('sequenza di ingresso (verde) e stima (rosso)');  pause,


