% Programma per la visualizzazione grafica di un generico segnale x(t)
% campionato in una finestra temporale di durata finita.
%
% Il programma chiede in ingresso:
% 
% 1 ) Estremi dell' Intervallo Temporale su cui si campiona il segnale
% 2 ) Espressione Analitica del Segnale
% 3 ) Frequenza di Campionamento
% 4 ) Risoluzione Spettrale DeltaF Minima
%
% Il programma è in grado di gestire le espressioni analitiche note a MATLAB 
%
%  sin( ) , cos( ) , exp( ) , log( ) , tan( ) , atan( ) , sqrt( ) , etc
%  
%  N.B. In Matlab sinc(x) = sin( pi*x)/(pi*x) ;
%
%
% e funzioni definite esternamente come 
%
% - rect( t ,T )  Rettangolo centrato nell'origine di Base T ed altezza 1
% - tri( t , T)   Triangolo centrato nell'origine di Semibase T ed altezza 1
%
% LA SINTASSI di immissione deve prevedere sempre come variabile indipendente il
% tempo 't' e le moltiplicazioni (*) tra diverse funzioni o elevamento a potenza (^)
% devono sempre essere in forma vettoriale ( Precedute dal punto '.*' '.^' )
%
% Esempio:     3 .* cos( 2*pi*10*t ) .* rect( t - 5 , 4 )  + t.^2
%
% Il programma rappresenta graficamente:
%
%  1) Andamento temporale e Campioni nel tempo 
%  2) Modulo e Fase dello Spettro stimato campionando lo Spettro tramite DFT
%
% --------------------------------------------------------------------------------
%   Esempio        x(t)     = 100*sinc( 100*t ) + 500*tri(t, 0.005 ) ;
%                  [t1 ,t2] = [ -0.4 , 0.4 ]
%                  Fc       = 1000
%                  DeltaF   = 1.25
% 
%  Paolo Banelli - 15/11/2001
%---------------------------------------------------------------------------------

clear all ;

% -------------  Inserimento Parametri --------------------------------------------------

xText  = input( 'Inserire espressione del segnale  x(t) = ', 's') ;
RangeT = input( 'Inserire Estremi Temporali di Osservazione  in sec.  (Es: [-1,5]): ' ); 
Fc     = input( 'Inserire Frequenza di Campionamento (in Hz):  ');
% ---------------------------------------------------------------------------------------

Tc         = 1/Fc        ; % Intervallo di Campionamento

VettIndici = round( RangeT / Tc ) ;
N1         = VettIndici(1);
N2         = VettIndici(2);


% DeltaT  = (N2-N1) * Tc  ; % Intervallo di Osservazione a multipli del periodo di campionamento
N       = N2-N1 ;     % Numero di Campioni prelevati nel tempo


% ---------------    Inserimento Risoluzione Spettrale ----------------------------------

DeltaF =  inf ; % Inizializzazione Variabile per Controllo 'while'
while DeltaF > (Fc/N) 
  textInput = ['Inserire Risoluzione Spettrale \DeltaF [ =< ', num2str(Fc/N) , ' Hz]  '];
  DeltaF    = input( textInput );
end
%----------------------------------------------------------------------------------------



t      = [N1 : N2-1] * Tc ;         % Vettore Istanti di Campionamento
eval( [ 'xc =' , xText ';' ] );     % Segnale Campionato

tplot  = [N1:0.1:N2] * Tc ;         % Vettore Istanti di Sovra-Campionamento per il Grafico
t      = tplot ;          ;     
eval( [ 'xplot =' , xText ';' ] );  % Vettore del Coseno per il Grafico ("Tempo Continuo")

t      = [N1 : N2-1] * Tc ;


Nfft   = max( N , fix(Fc/DeltaF) ) ; % Numero di Campioni su cui si calcola lo Spettro
                                     % Nfft >= Numero Campioni nel tempo

DeltaF = Fc / Nfft ;  % Risoluzione Spettrale 
                      % (Intervallo di Campionamento dell'asse delle Frequenze)


Xf     = fft( xc , Nfft )/Fc;         % Spettro Campionato del Segnale
f      = [ 0 : Nfft-1 ] * DeltaF ; % Asse delle Frequenze Campionato


PlotType = 0 ;
while ( PlotType ~= 1 & PlotType ~= 2 ) 
   PlotType = input('Figura Unica (1) o Separata (2) ? : ');
end


% -------------- Grafico nel tempo del Coseno Campionato

figure ;

if PlotType == 1
   subplot(211); % Ripartizione del Grafico su un'unica figura
   TitleText = ['Segnale x(t)= ' , xText  ,'  Campionato con Fc = ' , num2str(Fc) , ' (Hz)'  ];
else
   TitleText = ['Segnale x(t)= ' , xText  ];
end

plot( tplot , xplot ,'-b') ; grid on ; title( TitleText);

% hold on ; stem( t , xc ,'k' ); hold off ;

xlabel('Tempo  (sec)'); ylabel( 'x(t)' );

% --------------------------------------------------------------------------------------------------
                              

% -------------- Grafico del Modulo dello Spettro  ---------------------------------------------------------

if PlotType == 1
   subplot( 223 ); % Ripartizione dei Grafici su un'unica Finestra 
else
   figure ; subplot( 211 );
   TitleText = [' Spettro del segnale x(t)= ' , xText , '  Campionato con Fc = ' , num2str(Fc) , ' (Hz)'  ];
end

stem( f , abs(Xf) , 'r' ); grid on ; % Divido per N = DeltaT/Tc per compensare
% Il campionamento (Tc) ed il Finestramento (DeltaT) ylabel( '| X(f) |') ; xlabel( 'Frequenza (Hz)');

ylabel( 'Modulo[ X(f) ]' ) ; xlabel( 'Frequenza (Hz)'); 

% ---------------------------------------------------------------------------------------------------------



% -------------- Grafico della Fase dello Spettro ---------------------------------------------------------

if PlotType == 1
   subplot( 224 ); % Ripartizione dei Grafici su un'unica Finestra 
else
   title( TitleText);
   subplot( 212 );
end


to         = N1*Tc          ;                 % Traslazione del 1° campione nell'origine dell'asse dei tempi
Phase      = angle( Xf .* exp(-j*2*pi*f*to)); % Stima della Fase con Compensazione della Rotazione

stem( f , Phase *180/pi  ,'r'); grid on ; % Fase Disegnata in Gradi

 
% hold on ; stem( f(Ifo) , angle( Xf(Ifo) )*180/pi , 'k' ); hold off ; 

ylabel( 'Fase[ X(f) ]   (Gradi)') ; xlabel( 'Frequenza (Hz)'); axis([0,Fc, -180 , 180 ]);

% ---------------------------------------------------------------------------------------------------------