function est = kest_rec(N,A,C,Q,R,x0_,x0,P0,T,u,sys,t,gamma)

% RECURSIVE
fprintf('Recursive Kalman filtering\n')
fprintf('Compute off-line paramenter...')

% Allocate memory
[outs,states] = size(C);
P_prior = zeros(states,states,N); 
P_post = P_prior;
G = zeros(outs,outs,N);
iG = G;
L = zeros(states,outs,N);
e_prior=zeros(states,N);
e_post=zeros(states,N);

% Initialize
P_prior(:,:,1) = P0;

% Compute off-line parameters
for j = 1:N    
    G(:,:,j) = C*P_prior(:,:,j)*C' + R;    
    iG(:,:,j) = pinv( G(:,:,j) );
    L(:,:,j) = gamma(j)*P_prior(:,:,j)*C'*iG(:,:,j); 
    P_post(:,:,j) = P_prior(:,:,j) - L(:,:,j)*C*P_prior(:,:,j);
    P_prior(:,:,j+1) = A * P_post(:,:,j) * A' + Q;
end

fprintf('done\n')

fprintf('Compute estimate...')
[y,t,x] = lsim(sys,squeeze(u(:,:)),t,x0); y = y'; x = x';
x_post = zeros(states,N);
x_prior = zeros(states,N+1);
x_prior(:,1) = x0_;
% Compute estimate
for j = 1:N
    x_post(:,j) = x_prior(:,j) + L(:,:,j) * ( y(:,j) - C*x_prior(:,j) );
    x_prior(:,j+1) = A*x_post(:,j);
end
y_ = C * x_post;
fprintf('done\n')

est.x=x;
est.y=y;
est.P_prior = P_prior;
est.P_post = P_post;
est.L = L;
est.x_prior = x_prior;
est.x_post = x_post;


% Calcolo della traccia delle matrici di varianza d'errore

est.NP=zeros(N,1);
est.NP_post=zeros(N,1);

for j = 1:N
    est.NP(j)=trace(est.P_prior(:,:,j));
    est.NP_post(j)=trace(est.P_post(:,:,j));
end

est.P=P_prior(:,:,1);


