% Utilizzo della trasformata di Fourier per l'inversione della dinamica
clear all;
close all;

%======= Costante di tempo per sistema primo ordine ==========

tau = 33.5; % [s]         % Valore ricavato dall'analisi della risposta al gradino 


%=======   Inversione della dinamica   ==========

load gio_1_3_n2.txt;     % Lettura dei dati da FILE
                          % OSS: il FILE deve contenere esclusivamente  
                          % le tre colonne con i dati:
                          % 1a: tempo [s]
                          % 2a: Temperatura sonda [?C]
                          % 3a: Temperatura sonda prefiltrata (lasciarla stare)
                          % cancellare ogni eventuale aggiunta

TempReg = gio_1_3_n2;      % Cambio nome alla variabile

%***************************FILTRAGGIO*****************************************

[Num,Den] = Butter( ... , ...);     % Calcolo dei coefficienti di un filtro di Butterworth  

TempRegFil = filtfilt(Num,Den,TempReg(:,2)); % Filtraggio senza introduzione di ritardo
%***************************************************************************************

		figure
		plot(TempReg(:,1),TempReg(:,2));
		title('Andamento della temperatura con sistema a regime');
		xlabel('tempo [s]');
		ylabel('temperatura [?C]');

		hold on;
		plot(TempReg(:,1),TempRegFil,'m','LineWidth',2); %Visualiz dati filtrati
        legend('Dati originali','Dati Filtrati')
		grid on


%*****************************************************************************************
%======   Utilizzo della trasformata  per la compensazione dinamica =======
%*****************************************************************************************

% Estrazione dei periodi dal segnale non filtrato
TempCicli = '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

FreqCicli = fft(TempCicli(:,2));      % Trasformata di Fourier del segnale non filtrato


FreqCicliFil = fft(TempRegFil(:,2)); % Trasformata di Fourier del segnale filtrato

Tc = TempCicli(2,1) - TempCicli(1,1);  % estrazione delta tempo di campionamento

% Definizione del vettore delle frequenze:
f = '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!; 
%*****************************************************************************************

		figure
		plot(f,abs(FreqCicli),f,abs(FreqCicliFil));     % Visualizzazione del modulo della trasformata 
		
		title('Trasformata di Fourier del segnale originario e filtrato');
		xlabel('frequenza [Hz]');
		ylabel('modulo');
		set(gca,'YScale','log')

%*****************************************************************************************

FreqInvSonda = 1 + i*2*pi*f(1:end/2)' * tau;             % Calcolo della funzione di trasferimento inversa 
                                                                               % per frequenze positive

FreqInvSonda = '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   % aggiunta in coda della parte negativa     
%*****************************************************************************************

		figure
		subplot(2,1,1)
        Mod = abs(FreqInvSonda);
        Mod = Mod(1:end/2);
		plot(f(1:end/2) , Mod)
		set(gca,'YScale','log')
		title('Modulo della funzione inversa ');
		grid on
		subplot(2,1,2)
        Fase = angle(FreqInvSonda)*180/pi;
        Fase = Fase(1:end/2);
		plot(f(1:end/2) , Fase)
		title('Fase della funzione inversa ');
		xlabel('frequenza [Hz]');
		grid on

%*****************************************************************************************
%===========   Qu? avviene la compensazione della dinamica     ==============
%*****************************************************************************************
FreqIngSonda = FreqCicliFil.*FreqInvSonda;     % Calcolo della trasformata di Fourier 
                                                                          % del segnale di temperatura all'ingresso
                                                                          % della sonda (Inversione della dinamica)
                                             
TempIngSonda = ifft(FreqIngSonda); % Antitrasformata di Fourier del segnale di temperatura
%*****************************************************************************************
                                             
		figure 
		plot(TempCicli(:,1),TempCicli(:,2));         % Visualizzazione dei risultati
		hold on
		plot(TempCicli(:,1),sign(real(TempIngSonda)).*abs(TempIngSonda),'r',...
             'LineWidth',2);
		
		
		title('Inversione della dinamica');
		xlabel('tempo [s]');
		ylabel('temperatura [?C]');
		grid on

