% comparación entre respuesta sistema discreto y sistema continuo


% sistema continuo
num=[1];
den=[3 1];
g=tf(num,den);
[y_cont, t_cont] = step(g);

% ecuación en diferencias
% y(n) = b0*x(n) + b1*x(n-1) + ... + bnb*x(n-nb)
%                          - a1*y(n-1) - ... - ana*y(n-na)
                          
T = 1;                      % periodo de muestreo
td = 0:T:t_cont(end);       % vector de tiempos

n=length(td);
u_k=ones(1,n);
y_k=zeros(1,n);
y_k(1)=0.1429*u_k(1);        % salida inicial
for cc=2:n
    y_k(cc)=0.7143*y_k(cc-1)+0.1429*u_k(cc)+0.1429*u_k(cc-1);
end

% implementación IFIR
a=[0 -0.7143] ;
b=[0.1429 0.1429];
% inicialización
y_k_mat=zeros(1,length(a));
u_k_mat=zeros(1,length(b));
y=zeros(1,n);
for k=1:n
    % desplazamiento entrada
    for cc=length(u_k_mat):-1:2
        u_k_mat(cc)=u_k_mat(cc-1);
    end
    % actualización entrada
    u_k_mat(1)=u_k(k);   
    % ec. diferencias
    suma=0;
    for cc=1:length(b)
        suma=suma+b(cc)*u_k_mat(cc)-a(cc)*y_k_mat(cc);
    end
    % desplazamiento salida
    for cc=length(u_k_mat):-1:2
        y_k_mat(cc)=y_k_mat(cc-1);
    end    
    % actualización salida
    y_k_mat(end)=suma;     
    % salida
    y(k)=suma;    
end


% transformación discreta
gd=c2d(g,T,'tustin');
[y_dis, t_dis] = step(gd);

% comparación de respuestas
figure(1);hold off;
plot(t_cont,y_cont,'b');grid on;hold on;
stairs(t_dis,y_dis,'r');
stairs(td,y_k,'g');
stairs(td,y,'m');
title('continous (b), c2d (r), for (g), IFIR (m)');


