%% Transformada de Fourier, ejemplos Matlab sencillos
% © 2019, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en  http://personales.upv.es/asala/YT/V/tfour.html
% Este código ejecutó correctamente en Matlab R2018b
% Tabla de Contenidos
% 1.-Ejemplos de aplicación sobre señales
% 2.-Respuesta ante escalón (como FdT)
%     2.1- Zero-Phase filtering
%     2.2- Implementación con filtfilt (signal processing toolbox)
%

syms t w real %usaremos la symbolic toolbox para las manipulaciones
s=1j*w; %analogía con Laplace (señales NO persistentes)

%% 1.-Ejemplos de aplicación sobre señales
% Exponencial causal:

FF1=1/(s+3)
rr=ifourier(FF1)
fplot(rr), grid on

% Exponencial anticausal:

FF2=1/(-s+3)
rr=ifourier(FF2)
fplot(rr), grid on

% Exponencial causal+anticausal

FF=simplify(FF1+FF2) %6/(w^2+9)
rr=ifourier(FF)
fplot(rr),grid on, axis([-4 4 0 1])

%% 2.-Respuesta ante escalón (como FdT)

fourier(heaviside(t)) %como el escalón es "persistente" cos(0t), tiene términos "delta Dirac"
escalon_t=heaviside(t)*exp(-0.0001*t); %aproximación no persistente para poder utilizar la analogía con Laplace

tablaLaplaceExponencial=1/(s+0.0001)

escalon_f=fourier(escalon_t)
%escalon_f=1/s NO ES CORRECTO... para que Laplace y Fourier sean
%equivalentes, no debe haber "polos" en eje imaginario, todo debe ser
%"exponencialmente estable".

% Caso causal estable:

FF=3/(s+3)*escalon_f
rr1=simplify(ifourier(FF))
fplot(rr1), grid on, hold on, fplot(escalon_t,[-4 4]), hold off

% Caso "inestable" en Laplace, que en Fourier cuando converge es "anticausal"

FF=3/(-s+3)*escalon_f
rr2=simplify(ifourier(FF))
fplot(rr2), grid on,hold on, fplot(escalon_t,[-4 4]), hold off

%% 2.1- Zero-Phase filtering
% Caso causal+anticausal, muy común en "suavizado" en la práctica.
% La functión de transferencia 0.5(3/(s+3)+3/(-s+3)) tiene un Bode  de fase CERO (real positivo). [Comprobadlo]

FF=9/(w^2+9)*escalon_f
rr=simplify(ifourier(FF))
fplot(rr,[-4 4]), grid on, hold on, fplot(escalon_t,[-4 4]), hold off
%todo junto:
fplot([rr rr1 rr2],[-4 4]), grid on, hold on, fplot(escalon_t,[-4 4]), hold off
legend('filtro bilateral','filtro causal','filtro anticausal','escalon','Location','best')

%% 2.2- Implementación con filtfilt (signal processing toolbox)

s=tf('s');
G=3/(s+3);
Ts=0.04;
Gd=c2d(G,Ts,'tustin'); %discretización aproximada
u=[zeros(1,100) ones(1,100)];
y=filtfilt(Gd.num{:},Gd.den{:},u);
plot((0:199)*Ts-100*Ts,[y;u]'), grid on, 
hold on, fplot(rr,[-4 4]),hold off, legend('filtfilt','entrada','Fourier continuo')