%--------------------------------------------------------------------------
%
% Kuramoto Sivashinsky equation
%
% External description.
%
% An alternative description of the system, based on the Input/Output
% behavior of the system is presented. In particular, the response of the
% system is calculated via convolution sums and compared with the LTI
% system description.
%
% Reference figure(s): Fig. 7, Fig. 10
% 
% Nicolo' Fabbiane, March 2014
% KTH Mechanics
% nicolo@mech.kth.se
%
%--------------------------------------------------------------------------
clc, clear all
fprintf('\nExternal description.\n\n')




%% Initialization of the plant --------------------------------------------

% LTI system matrix A and x grid
nq = 400; % number of grid points

[A,x,I] = KS_init(nq);


% Inputs matrix B

% disturbance d (Gaussian shape at x_d with sigma_d variance)
x_d = 35; sigma_d = 4; 
Bd  = exp(-((x-x_d)/sigma_d).^2)/sigma_d;

% actuator u (Gaussian shape at x_u with sigma_u variance)
x_u = 400; sigma_u = 4; 
Bu  = exp(-((x-x_u)/sigma_u).^2)/sigma_u;


% Outputs matrix C

% measurement y (Gaussian shape at x_y with sigma_y variance)
x_y = 300; sigma_y = 4; 
Cy  = (exp(-((x-x_y)/sigma_y).^2)/sigma_y).' * I;

% output z (Gaussian shape at x_z with sigma_z variance)
x_z = 700; sigma_z = 4; 
Cz  = (exp(-((x-x_z)/sigma_z).^2)/sigma_z).' * I;




%% Convolution kernels Pyu, Pzu -------------------------------------------

% Parameters
dt   = 1;
Tker = 2000;

tker = 0:dt:Tker;
nker = length(tker);


% Time-discrete system (Eqn. 37-38)
Ad  = expm(A*dt); Bud = dt*Bu;
Cyd = Cy;
Czd = Cz;


% Convolution kernel Pqu (Eqn. 45)
Pvu = zeros(nq,nker);

Pvu(:,1) = Bud;
for i = 2:nker
    Pvu(:,i) = Ad*Pvu(:,i-1);
end


% Convolution kernels Pyu, Pzu (Eqn. 48)
Pyu = Cyd * Pvu;
Pzu = Czd * Pvu;

        

        
%% FIR reduction ----------------------------------------------------------
toll = 1e-2; % kernel reduction tollerance

% select the interval where |Pzu| > max(|Pzu|) * toll
ifir = (find( abs(Pzu) > max(abs(Pzu))*toll ,1,'first')) : ...
       (find( abs(Pzu) > max(abs(Pzu))*toll ,1,'last'));
   
Pzu_fir = Pzu(ifir);


        

        
%% Time integration -------------------------------------------------------

% Parameters
dt   = 1;
Tend = 50000;

t  = 0:dt:Tend;
nt = length(t);


% Time-stepper matrices (Crank-Nicholson)
CNi = (speye(nq) - dt/2*A);
CNe = (speye(nq) + dt/2*A);


% Noise inizialization
u = randn(1,nt);     % unitary variance Gaussian white noise
u = u - mean(u,2);   % enforce zero-mean
u = u / std(u,[],2); % enforce unitary variance


% Time-loop (LTI system)
v = zeros(nq,nt);
z = zeros(1 ,nt);

for i = 1:nt-1
    
    % outputs
    z(:,i) = Cz*v(:,i);
    
    % timestep
    v(:,i+1) = CNi\(CNe*v(:,i) + dt * Bu*u(:,i));
    
    % runtime output
    runtime_output

end


% Time-loop (External description - full kernel)
z_ker = zeros(1,nt);

u_buf = zeros(1,nker); % input signal buffer

for i = 1:nt-1
    
    % outputs (Eqn. 49)
    z_ker(:,i) = Pzu * u_buf';
    
    % update input buffer
    u_buf = [u(:,i) u_buf(:,1:end-1)];

end


% Time-loop (External description - FIR)
z_fir = zeros(1,nt);

u_buf = zeros(1,ifir(end)); % input signal buffer

for i = 1:nt-1
    
    % outputs (Eqn. 50)
    z_fir(:,i) = Pzu_fir * u_buf(ifir)';
    
    % update input buffer
    u_buf = [u(:,i) u_buf(:,1:end-1)];

end




%% Figure 7 ---------------------------------------------------------------
figure(7); clf

plot(tker      ,Pzu    ,'or',...
     tker(ifir),Pzu_fir,'sb','MarkerSize',5)
axis([0 1200 -0.3 0.4]);
xlabel('i \Deltat'); ylabel('P_{zu}(i)'), grid on

legend('P_{zu}(i)','P_{zu}(i) (FIR)','Location','NW')




%% Figure 10 --------------------------------------------------------------
figure(10); clf

subplot(4,1,1:2); imagesc(x,tker,Pvu'); hold on
                  plot([1 1]*x_y,[0 Tker],'-k', ...
                       [1 1]*x_z,[0 Tker],'-k','Linewidth',2); hold off
                  axis xy, caxis([-1 1]*.3)
                  xlabel('x'); ylabel('t'); title('P_{vu}(t)')
subplot(4,1,3);   plot(tker,Pyu,'-k','Linewidth',1.5)
                  axis([0 2000 -.01 .01]);
                  xlabel('t'); ylabel('P_{yu}(t)'), grid on
subplot(4,1,4);   plot(tker,Pzu,'-k','Linewidth',1.5)
                  axis([0 2000 -.3 .4]);
                  xlabel('t'); ylabel('P_{zu}(t)'), grid on




%% z signal comparison ----------------------------------------------------
figure(201); clf

subplot(2,1,1); plot(t,z,'-k','Linewidth',1.5); hold on
                plot(t,z_ker,'or',...
                     t,z_fir,'sb','MarkerSize',5); hold off
                axis([Tend-2500 Tend -6 6]);
                xlabel('t'); ylabel('z(t)'), grid on

                legend('LTI','full kernel','FIR','Location','NW')

subplot(2,1,2); plot(t,z-z_ker,'or',...
                     t,z-z_fir,'sb','MarkerSize',5); hold off
                ax = axis; axis([Tend-2500 Tend -0.1 0.1]);
                xlabel('t'); ylabel('error'), grid on