diff --git a/+LTV_signalmodelfault/LTV_SigModF_FOH.m b/+LTV_signalmodelfault/LTV_SigModF_FOH.m
new file mode 100644
index 0000000000000000000000000000000000000000..809839e2f06a8d0bc8bc0967a9188219220486db
--- /dev/null
+++ b/+LTV_signalmodelfault/LTV_SigModF_FOH.m
@@ -0,0 +1,280 @@
+%% Simulation for fault identification of a linear time variant system of the form
+% The fault and the disturbance are assumed to be described by a signal
+% model of the form
+% | x(t) | |A GP_d E_1P_f | | x(t) | | B |
+%d/dt | v_d(t) | = |0 S_d 0 | | v_d(t) | + | 0 |u(t)
+% | v_f(t) | |0 0 S_f | | v_f(t) | | 0 |
+%
+% | x(t) |
+% y(t) = |C 0 E_2P_f | | v_d(t) |
+% | v_f(t) |
+
+
+%% Define the simulation time t_end, the interval length T and the discretization step dt
+
+ T = 2;
+ dt = 1e-2;
+ t_end = 3*T;
+
+ t =(-T:dt:t_end)';
+ tau = (0:dt:T);
+
+%% Define System Dynamic
+% the system matrices A(t), B(t), C(t), G(t), E_1(t), E_2(t)
+% the matrices are defined as function handle matrices
+
+ A = @(t) [ -2 , 0 ;
+ 0 , -4];
+
+
+ B = @(t) [1;
+ 0];
+
+ C = @(t) [ -2 ,-3 ];
+
+ G = @(t) [ 0 ;
+ 1 ];
+
+ E_1 = @(t) [ -7 ;
+ 0 ];
+
+ E_2 = @(t) [0];
+ %% Define the fault signal model
+ % d/dt v_f(t) = S_f(t)v_f(t)
+ % f(t) = P_f(t)v_f(t)
+
+ S_f = @(t) 0;
+ P_f = @(t) 1;
+
+ %% Define the disturbance signal model
+ % d/dt v_d(t) = S_d(t)v_d(t)
+ % d(t) = P_d(t)v_d(t)
+ w0 = 2;
+
+ S_d = @(t) [0 , 1 ;
+ -(w0)^2 , 0];
+ P_d = @(t) [1 , 0 ];
+
+%% define the extended System
+
+n_x = size(A(1),1 );
+n_u = size(B(1),2 );
+n_y = size(C(1),1 );
+n_d = size(G(1),2 );
+n_f = size(E_1(1),2);
+n_vd = size(S_d(1),1);
+n_vf = size(S_f(1),1);
+
+A_E = @(t) [A(t) , G(t)*P_d(t) , E_1(t)*P_f(t);
+ zeros(n_vd,n_x) , S_d(t) , zeros(n_vd,n_vf);
+ zeros(n_vf,n_x) , zeros(n_vf,n_vd), S_f(t)];
+
+B_E = @(t) [B(t); zeros(n_vd,n_u); zeros(n_vf,n_u)];
+
+C_E = @(t) [C(t), zeros(n_y,n_vd), E_2(t)*P_f(t)];
+
+n_E = size(A_E(1),1);
+
+%% Determining the Output
+% here enter the intial values of the system
+t1 =(0:dt:t_end);
+opts = odeset('RelTol',1e-6,'AbsTol',1e-9);
+fac = 70;
+dt_dis = fac*dt;
+t_dis = 0:dt_dis:t_end;
+
+u = @(t) exp(-t).*sin(t);
+%u = @(t) 0;
+xE = ones(n_E,1);
+
+[t_dis, xE ] = ode45(@(t,xE) A_E(t)*xE + B_E(t)*u(t), t_dis,[zeros(n_x,1);0;10;6], opts);
+
+
+v_d = xE(:,(n_x+1):(n_vd+n_x)).';
+v_f = xE(:,(n_vd+n_x+1):n_E).';
+
+for k = 1:length(t_dis)
+ d(k,:) = P_d(t_dis(k)) * v_d(:,k);
+ y(k,:) = C_E(t_dis(k)) * xE(k,:).';
+ f(k,:) = P_f(t_dis(k)) * v_f(:,k);
+end
+
+% figure
+% hold on
+% plot(t_dis,f,'--');
+% plot(t_dis,y,':','LineWidth',3);
+% plot(t_dis,d);
+% hold off
+
+% t1 = t1.';
+
+%% Discretization of the -A_E.' and C_E matrices for the whole simualtion interval [-T, t_end]
+
+ for k = 1:size(t,1)
+
+ Aek(k,:,:) = -A_E(t(k)).';
+ Ce_Glob(:,:,k)= C_E(t(k)).';
+ end
+
+%% Simulating the Transition matrix for the whole simualtion interval [-T, t_end]
+
+Phi = misc.fundamentalMatrix.odeSolver_par(Aek, t);
+
+% n(tau,t) =
+% C_E.'(tau,t)Phi.'(tau,tau_0,t)W^(-1)(tau_0,T,t)Phi(tau_0,T,t)m_E,i1
+% m_E = Phi(tau,tau_0,t)int(Phi(tau_0,tau,t)C_E(tau,t)n(tau,t),tau,0,tau)
+% with the controlabillity Gramian
+% W(tau_0,T,t) =
+% int(Phi(tau_0,tau,t)C_E(tau,t)C_E.'(tau,t)Phi.'(tau_0,tau,t),tau,0,T)
+
+winlng = size(tau,2);
+I = eye(n_vf);
+N = zeros(winlng,size(t1,2),n_y,n_vf);
+m_E = zeros(n_E,length(tau),length(t1),n_vf);
+
+% The only a segment of the Phi matrix on the interval [t-T,t] is extracted
+% and the calculation for m_E and n is done numerically for the interval [0, t_end]
+for indtMinT = 1:size(t1,2)
+
+ indt = indtMinT + winlng -1;
+
+ Phitau = Phi(indtMinT:indt,indtMinT,:,:);
+
+ Ce_loc = Ce_Glob(:,:,indtMinT:indt);
+ %% Calculating W numerically
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+
+ Phit0t = inv(Phitt0(:,:,k));
+
+ w(k,:,:) = Phit0t*Ce_loc(:,:,k)*(Ce_loc(:,:,k).')*(Phit0t.');
+
+ end
+
+
+ tmp = reshape(w, length(tau), n_E^2);
+ tmp_I = trapz(tau, tmp);
+ W = reshape(tmp_I, n_E , n_E);
+
+
+%% Calculating n denoted as N because its a matrix :)
+
+
+ PhiTt0(:,:,1) = Phitau(end,1,:,:);
+ Phit0T = inv(PhiTt0);
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+ Phit0t = inv(Phitt0(:,:,k));
+
+ for ki = 1:n_vf
+ N(k,indtMinT,:,ki) = (Ce_loc(:,:,k).')*(Phit0t.')*inv(W)*(Phit0T*[zeros(n_E-n_vf,1);I(:,ki)]);
+ end
+ % N(tau,t,n_y,r_sum)
+ end
+ %% calculating the states m_E
+
+ for ki = 1:n_vf
+ m_E(:,:,indtMinT,ki) = calculateStates(dt,tau,Phitau,Ce_loc,N(:,indtMinT,:,ki));
+ end
+ % m_E(nk,tau,t,r_sum)
+end
+
+%% calculating < m_E,B_Eu >
+ u_dis = u(t_dis);
+ u_foh = interp1(t_dis,u_dis,t1);
+
+
+ for ki = 1:n_vf
+ mB(:,ki) = kernTransformTwo(t,t1,tau,m_E(:,:,:,ki),B_E,u_foh);
+ end
+
+%% calculating < n,y >
+% y is assumed to be zero on the interval from [-T,0] and a segment of y is
+% extracted on the interval from [t-T,t] and then multiplied with n for
+% every t, the result is stored in NY
+
+ y_foh = interp1(t_dis,y,t1).';
+
+
+ for ind2 = 1:length(t1)
+ Summe = zeros(n_vf,1);
+ yk = zeros(winlng,n_y);
+ diff = winlng - ind2;
+
+ if((t1(ind2)-T) < 0)
+ ind1 = 1;
+ yk((ind1+diff):winlng,:) = y_foh(ind1:ind2,:);
+ elseif((t1(ind2)-T)>=0)
+ ind1 = ind2-winlng+1;
+ yk = y_foh(ind1:ind2,:);
+ end
+
+ for j = 1:n_y
+ for ki = 1:n_vf
+ Summe(ki,1) = trapz(tau,N(:,ind2,j,ki).*yk(:,j)) + Summe(ki,1);
+ end
+ end
+ NY(ind2,:) = Summe;
+ end
+
+ % NY(t,r_sum)
+
+%% calculating ZOH
+% delta_t = dt_dis;
+% N_k = find(t_dis == T)-1;
+%
+%
+% for t_ind = 1:length(t1)
+% yk = zeros(winlng,n_y);
+% diff = winlng - t_ind;
+% nt = N(:,t_ind);
+% sum = 0;
+% if((t1(t_ind)-T) < 0)
+% tminT_ind = 1;
+% yk((tminT_ind+diff):winlng,:) = y_zoh(tminT_ind:t_ind,:);
+% ny(t_ind,1) = 0;
+% elseif((t1(ind2)-T)>=0)
+% tminT_ind = t_ind-winlng+1;
+% yk = y_zoh(tminT_ind:t_ind,:);
+%
+%
+% y_temp = unique(yk);
+% t_dis_temp = t_dis(t_dis > t1(tminT_ind) & t_dis < t1(t_ind));
+% limits = t_dis_temp - t1(tminT_ind);
+% limits = [limits;T];
+%
+% tau_1 = 1;
+% for i = 1:N_k
+% tau_2 = find(abs(tau - limits(i)) < 1e-6);
+% sum = sum + y_temp(i)*trapz(tau(tau_1:tau_2),nt(tau_1:tau_2));
+% tau_1 = tau_2;
+% end
+%
+% ny(t_ind,1) = sum;
+% end
+% end
+% figure
+% plot(t1,ny)
+
+ %% calculating and ploting the fault f_1, f_2,...., f_nf
+ % FHAT is the result of the equation for the coefficients
+ % vartheta = -<n,y>- <m,Bu>
+
+ FHAT = NY + mB;
+
+ for k = 1:length(t1)
+ f_Hat(k,:) = P_f(k)*FHAT(k,:);
+ end
+
+
+ for i = 1:n_f
+ subplot(n_f,1,i)
+ hold on;
+ plot(t_dis,f(:,i),'-.','LineWidth',3);
+ plot(t1,f_Hat(:,i),'LineWidth',3);
+ hold off;
+ end
\ No newline at end of file
diff --git a/+LTV_signalmodelfault/LTV_SigModF_ZOH.m b/+LTV_signalmodelfault/LTV_SigModF_ZOH.m
new file mode 100644
index 0000000000000000000000000000000000000000..06c65982593685cbe7fc6298d6f6f9a41f05077a
--- /dev/null
+++ b/+LTV_signalmodelfault/LTV_SigModF_ZOH.m
@@ -0,0 +1,280 @@
+%% Simulation for fault identification of a linear time variant system of the form
+% The fault and the disturbance are assumed to be described by a signal
+% model of the form
+% | x(t) | |A GP_d E_1P_f | | x(t) | | B |
+%d/dt | v_d(t) | = |0 S_d 0 | | v_d(t) | + | 0 |u(t)
+% | v_f(t) | |0 0 S_f | | v_f(t) | | 0 |
+%
+% | x(t) |
+% y(t) = |C 0 E_2P_f | | v_d(t) |
+% | v_f(t) |
+
+
+%% Define the simulation time t_end, the interval length T and the discretization step dt
+
+ T = 1;
+ dt = 1e-2;
+ t_end = 3*T;
+
+ t =(-T:dt:t_end)';
+ tau = (0:dt:T);
+
+%% Define System Dynamic
+% the system matrices A(t), B(t), C(t), G(t), E_1(t), E_2(t)
+% the matrices are defined as function handle matrices
+
+ A = @(t) [ -2 , 0 ;
+ 0 , -4];
+
+
+ B = @(t) [1;
+ 0];
+
+ C = @(t) [ -2 ,-3 ];
+
+ G = @(t) [ 0 ;
+ 1 ];
+
+ E_1 = @(t) [ -7 ;
+ 0 ];
+
+ E_2 = @(t) [0];
+ %% Define the fault signal model
+ % d/dt v_f(t) = S_f(t)v_f(t)
+ % f(t) = P_f(t)v_f(t)
+
+ S_f = @(t) 0;
+ P_f = @(t) 1;
+
+ %% Define the disturbance signal model
+ % d/dt v_d(t) = S_d(t)v_d(t)
+ % d(t) = P_d(t)v_d(t)
+ w0 = 2;
+
+ S_d = @(t) [0 , 1 ;
+ -(w0)^2 , 0];
+ P_d = @(t) [1 , 0 ];
+
+%% define the extended System
+
+n_x = size(A(1),1 );
+n_u = size(B(1),2 );
+n_y = size(C(1),1 );
+n_d = size(G(1),2 );
+n_f = size(E_1(1),2);
+n_vd = size(S_d(1),1);
+n_vf = size(S_f(1),1);
+
+A_E = @(t) [A(t) , G(t)*P_d(t) , E_1(t)*P_f(t);
+ zeros(n_vd,n_x) , S_d(t) , zeros(n_vd,n_vf);
+ zeros(n_vf,n_x) , zeros(n_vf,n_vd), S_f(t)];
+
+B_E = @(t) [B(t); zeros(n_vd,n_u); zeros(n_vf,n_u)];
+
+C_E = @(t) [C(t), zeros(n_y,n_vd), E_2(t)*P_f(t)];
+
+n_E = size(A_E(1),1);
+
+%% Determining the Output
+% here enter the intial values of the system
+t1 =(0:dt:t_end);
+opts = odeset('RelTol',1e-6,'AbsTol',1e-9);
+fac = 20;
+dt_dis = fac*dt;
+t_dis = 0:dt_dis:t_end;
+
+u = @(t) exp(-t).*sin(t);
+%u = @(t) 0;
+xE = ones(n_E,1);
+
+[t_dis, xE ] = ode45(@(t,xE) A_E(t)*xE + B_E(t)*u(t), t_dis,[zeros(n_x,1);0;10;6], opts);
+
+
+v_d = xE(:,(n_x+1):(n_vd+n_x)).';
+v_f = xE(:,(n_vd+n_x+1):n_E).';
+
+for k = 1:length(t_dis)
+ d(k,:) = P_d(t_dis(k)) * v_d(:,k);
+ y(k,:) = C_E(t_dis(k)) * xE(k,:).';
+ f(k,:) = P_f(t_dis(k)) * v_f(:,k);
+end
+
+% figure
+% hold on
+% plot(t_dis,f,'--');
+% plot(t_dis,y,':','LineWidth',3);
+% plot(t_dis,d);
+% hold off
+
+% t1 = t1.';
+
+%% Discretization of the -A_E.' and C_E matrices for the whole simualtion interval [-T, t_end]
+
+ for k = 1:size(t,1)
+
+ Aek(k,:,:) = -A_E(t(k)).';
+ Ce_Glob(:,:,k)= C_E(t(k)).';
+ end
+
+%% Simulating the Transition matrix for the whole simualtion interval [-T, t_end]
+
+Phi = misc.fundamentalMatrix.odeSolver_par(Aek, t);
+
+% n(tau,t) =
+% C_E.'(tau,t)Phi.'(tau,tau_0,t)W^(-1)(tau_0,T,t)Phi(tau_0,T,t)m_E,i1
+% m_E = Phi(tau,tau_0,t)int(Phi(tau_0,tau,t)C_E(tau,t)n(tau,t),tau,0,tau)
+% with the controlabillity Gramian
+% W(tau_0,T,t) =
+% int(Phi(tau_0,tau,t)C_E(tau,t)C_E.'(tau,t)Phi.'(tau_0,tau,t),tau,0,T)
+
+winlng = size(tau,2);
+I = eye(n_vf);
+N = zeros(winlng,size(t1,2),n_y,n_vf);
+m_E = zeros(n_E,length(tau),length(t1),n_vf);
+
+% The only a segment of the Phi matrix on the interval [t-T,t] is extracted
+% and the calculation for m_E and n is done numerically for the interval [0, t_end]
+for indtMinT = 1:size(t1,2)
+
+ indt = indtMinT + winlng -1;
+
+ Phitau = Phi(indtMinT:indt,indtMinT,:,:);
+
+ Ce_loc = Ce_Glob(:,:,indtMinT:indt);
+ %% Calculating W numerically
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+
+ Phit0t = inv(Phitt0(:,:,k));
+
+ w(k,:,:) = Phit0t*Ce_loc(:,:,k)*(Ce_loc(:,:,k).')*(Phit0t.');
+
+ end
+
+
+ tmp = reshape(w, length(tau), n_E^2);
+ tmp_I = trapz(tau, tmp);
+ W = reshape(tmp_I, n_E , n_E);
+
+
+%% Calculating n denoted as N because its a matrix :)
+
+
+ PhiTt0(:,:,1) = Phitau(end,1,:,:);
+ Phit0T = inv(PhiTt0);
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+ Phit0t = inv(Phitt0(:,:,k));
+
+ for ki = 1:n_vf
+ N(k,indtMinT,:,ki) = (Ce_loc(:,:,k).')*(Phit0t.')*inv(W)*(Phit0T*[zeros(n_E-n_vf,1);I(:,ki)]);
+ end
+ % N(tau,t,n_y,r_sum)
+ end
+ %% calculating the states m_E
+
+ for ki = 1:n_vf
+ m_E(:,:,indtMinT,ki) = calculateStates(dt,tau,Phitau,Ce_loc,N(:,indtMinT,:,ki));
+ end
+ % m_E(nk,tau,t,r_sum)
+end
+
+%% calculating < m_E,B_Eu >
+ u_dis = u(t_dis);
+ u_zoh = repelem(u_dis,fac);
+ u_zoh = u_zoh(1:length(t1));
+
+ for ki = 1:n_vf
+ mB(:,ki) = kernTransformTwo(t,t1,tau,m_E(:,:,:,ki),B_E,u_zoh);
+ end
+
+%% calculating < n,y >
+% y is assumed to be zero on the interval from [-T,0] and a segment of y is
+% extracted on the interval from [t-T,t] and then multiplied with n for
+% every t, the result is stored in NY
+
+ y_zoh = repelem(y,fac);
+ y_zoh = y_zoh(1:length(t1));
+
+ for ind2 = 1:length(t1)
+ Summe = zeros(n_vf,1);
+ yk = zeros(winlng,n_y);
+ diff = winlng - ind2;
+
+ if((t1(ind2)-T) < 0)
+ ind1 = 1;
+ yk((ind1+diff):winlng,:) = y_zoh(ind1:ind2,:);
+ elseif((t1(ind2)-T)>=0)
+ ind1 = ind2-winlng+1;
+ yk = y_zoh(ind1:ind2,:);
+ end
+
+ for j = 1:n_y
+ for ki = 1:n_vf
+ Summe(ki,1) = trapz(tau,N(:,ind2,j,ki).*yk(:,j)) + Summe(ki,1);
+ end
+ end
+ NY(ind2,:) = Summe;
+ end
+
+ % NY(t,r_sum)
+
+%% calculating ZOH
+% delta_t = dt_dis;
+% N_k = find(t_dis == T)-1;
+%
+%
+% for t_ind = 1:length(t1)
+% yk = zeros(winlng,n_y);
+% diff = winlng - t_ind;
+% nt = N(:,t_ind);
+% sum = 0;
+% if((t1(t_ind)-T) < 0)
+% tminT_ind = 1;
+% yk((tminT_ind+diff):winlng,:) = y_zoh(tminT_ind:t_ind,:);
+% ny(t_ind,1) = 0;
+% elseif((t1(ind2)-T)>=0)
+% tminT_ind = t_ind-winlng+1;
+% yk = y_zoh(tminT_ind:t_ind,:);
+%
+%
+% y_temp = unique(yk);
+% t_dis_temp = t_dis(t_dis > t1(tminT_ind) & t_dis < t1(t_ind));
+% limits = t_dis_temp - t1(tminT_ind);
+% limits = [limits;T];
+%
+% tau_1 = 1;
+% for i = 1:N_k
+% tau_2 = find(abs(tau - limits(i)) < 1e-6);
+% sum = sum + y_temp(i)*trapz(tau(tau_1:tau_2),nt(tau_1:tau_2));
+% tau_1 = tau_2;
+% end
+%
+% ny(t_ind,1) = sum;
+% end
+% end
+% figure
+% plot(t1,ny)
+
+ %% calculating and ploting the fault f_1, f_2,...., f_nf
+ % FHAT is the result of the equation for the coefficients
+ % vartheta = -<n,y>- <m,Bu>
+
+ FHAT = NY + mB;
+
+ for k = 1:length(t1)
+ f_Hat(k,:) = P_f(k)*FHAT(k,:);
+ end
+
+ figure
+ for i = 1:n_f
+ subplot(n_f,1,i)
+ hold on;
+ plot(t_dis,f(:,i),'-.','LineWidth',3);
+ plot(t1,f_Hat(:,i),'LineWidth',3);
+ hold off;
+ end
\ No newline at end of file
diff --git a/+LTV_signalmodelfault/LTV_SignalModelFaults.m b/+LTV_signalmodelfault/LTV_SignalModelFaults.m
new file mode 100644
index 0000000000000000000000000000000000000000..72b7a179e74fbebb483c956a035149979240d738
--- /dev/null
+++ b/+LTV_signalmodelfault/LTV_SignalModelFaults.m
@@ -0,0 +1,236 @@
+%% Simulation for fault identification of a linear time variant system of the form
+% The fault and the disturbance are assumed to be described by a signal
+% model of the form
+% | x(t) | |A GP_d E_1P_f | | x(t) | | B |
+%d/dt | v_d(t) | = |0 S_d 0 | | v_d(t) | + | 0 |u(t)
+% | v_f(t) | |0 0 S_f | | v_f(t) | | 0 |
+%
+% | x(t) |
+% y(t) = |C 0 E_2P_f | | v_d(t) |
+% | v_f(t) |
+
+
+%% Define the simulation time t_end, the interval length T and the discretization step dt
+
+ T = 1;
+ dt = 1e-2;
+ t_end = 3*T;
+
+ t =(-T:dt:t_end)';
+ tau = (0:dt:T);
+
+%% Define System Dynamic
+% the system matrices A(t), B(t), C(t), G(t), E_1(t), E_2(t)
+% the matrices are defined as function handle matrices
+
+ A = @(t) [ -2 , 0 ;
+ 0 , -4];
+
+
+ B = @(t) [1;
+ 0];
+
+ C = @(t) [ -2 ,-3 ];
+
+ G = @(t) [ 0 ;
+ 1 ];
+
+ E_1 = @(t) [ -7 ;
+ 0 ];
+
+ E_2 = @(t) [0];
+ %% Define the fault signal model
+ % d/dt v_f(t) = S_f(t)v_f(t)
+ % f(t) = P_f(t)v_f(t)
+
+ S_f = @(t) 0;
+ P_f = @(t) 1;
+
+ %% Define the disturbance signal model
+ % d/dt v_d(t) = S_d(t)v_d(t)
+ % d(t) = P_d(t)v_d(t)
+ w0 = 2;
+
+ S_d = @(t) [0 , 1 ;
+ -(w0)^2 , 0];
+ P_d = @(t) [1 , 0 ];
+
+%% define the extended System
+
+n_x = size(A(1),1 );
+n_u = size(B(1),2 );
+n_y = size(C(1),1 );
+n_d = size(G(1),2 );
+n_f = size(E_1(1),2);
+n_vd = size(S_d(1),1);
+n_vf = size(S_f(1),1);
+
+A_E = @(t) [A(t) , G(t)*P_d(t) , E_1(t)*P_f(t);
+ zeros(n_vd,n_x) , S_d(t) , zeros(n_vd,n_vf);
+ zeros(n_vf,n_x) , zeros(n_vf,n_vd), S_f(t)];
+
+B_E = @(t) [B(t); zeros(n_vd,n_u); zeros(n_vf,n_u)];
+
+C_E = @(t) [C(t), zeros(n_y,n_vd), E_2(t)*P_f(t)];
+
+n_E = size(A_E(1),1);
+
+%% Determining the Output
+% here enter the intial values of the system
+t1 =(0:dt:t_end);
+opts = odeset('RelTol',1e-6,'AbsTol',1e-9);
+
+
+u = @(t) exp(-t).*sin(t);
+%u = @(t) 0;
+xE = ones(n_E,1);
+
+[t_dis, xE ] = ode45(@(t,xE) A_E(t)*xE + B_E(t)*u(t), t1,[zeros(n_x,1);0;10;6], opts);
+
+
+v_d = xE(:,(n_x+1):(n_vd+n_x)).';
+v_f = xE(:,(n_vd+n_x+1):n_E).';
+
+for k = 1:length(t_dis)
+ d(k,:) = P_d(t_dis(k))*v_d(:,k);
+ y(k,:) = C_E(t_dis(k)) * xE(k,:).';
+ f(k,:) = P_f(t_dis(k))* v_f(:,k);
+end
+
+% figure
+% hold on
+% plot(t_dis,f,'--');
+% plot(t_dis,y,':','LineWidth',3);
+% plot(t_dis,d);
+% hold off
+
+% t1 = t1.';
+
+%% Discretization of the -A_E.' and C_E matrices for the whole simualtion interval [-T, t_end]
+
+ for k = 1:size(t,1)
+
+ Aek(k,:,:) = -A_E(t(k)).';
+ Ce_Glob(:,:,k)= C_E(t(k)).';
+ end
+
+%% Simulating the Transition matrix for the whole simualtion interval [-T, t_end]
+
+Phi = misc.fundamentalMatrix.odeSolver_par(Aek, t);
+
+% n(tau,t) =
+% C_E.'(tau,t)Phi.'(tau,tau_0,t)W^(-1)(tau_0,T,t)Phi(tau_0,T,t)m_E,i1
+% m_E = Phi(tau,tau_0,t)int(Phi(tau_0,tau,t)C_E(tau,t)n(tau,t),tau,0,tau)
+% with the controlabillity Gramian
+% W(tau_0,T,t) =
+% int(Phi(tau_0,tau,t)C_E(tau,t)C_E.'(tau,t)Phi.'(tau_0,tau,t),tau,0,T)
+
+winlng = size(tau,2);
+I = eye(n_vf);
+N = zeros(winlng,size(t1,2),n_y,n_vf);
+m_E = zeros(n_E,length(tau),length(t1),n_vf);
+
+% The only a segment of the Phi matrix on the interval [t-T,t] is extracted
+% and the calculation for m_E and n is done numerically for the interval [0, t_end]
+for indtMinT = 1:size(t1,2)
+
+ indt = indtMinT + winlng -1;
+
+ Phitau = Phi(indtMinT:indt,indtMinT,:,:);
+
+ Ce_loc = Ce_Glob(:,:,indtMinT:indt);
+ %% Calculating W numerically
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+
+ Phit0t = inv(Phitt0(:,:,k));
+
+ w(k,:,:) = Phit0t*Ce_loc(:,:,k)*(Ce_loc(:,:,k).')*(Phit0t.');
+
+ end
+
+
+ tmp = reshape(w, length(tau), n_E^2);
+ tmp_I = trapz(tau, tmp);
+ W = reshape(tmp_I, n_E , n_E);
+
+
+%% Calculating n denoted as N because its a matrix :)
+
+
+ PhiTt0(:,:,1) = Phitau(end,1,:,:);
+ Phit0T = inv(PhiTt0);
+
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+ Phit0t = inv(Phitt0(:,:,k));
+
+ for ki = 1:n_vf
+ N(k,indtMinT,:,ki) = (Ce_loc(:,:,k).')*(Phit0t.')*inv(W)*(Phit0T*[zeros(n_E-n_vf,1);I(:,ki)]);
+ end
+ % N(tau,t,n_y,r_sum)
+ end
+ %% calculating the states m_E
+
+ for ki = 1:n_vf
+ m_E(:,:,indtMinT,ki) = calculateStates(dt,tau,Phitau,Ce_loc,N(:,indtMinT,:,ki));
+ end
+ % m_E(nk,tau,t,r_sum)
+end
+
+%% calculating < m_E,B_Eu >
+
+ for ki = 1:n_vf
+ mB(:,ki) = kernTransform(t,t1,tau,m_E(:,:,:,ki),B_E,u);
+ end
+
+%% calculating < n,y >
+% y is assumed to be zero on the interval from [-T,0] and a segment of y is
+% extracted on the interval from [t-T,t] and then multiplied with n for
+% every t, the result is stored in NY
+ for ind2 = 1:length(t1)
+ Summe = zeros(n_vf,1);
+ yk = zeros(winlng,n_y);
+ diff = winlng - ind2;
+
+ if((t1(ind2)-T) < 0)
+ ind1 = 1;
+ yk((ind1+diff):winlng,:) = y(ind1:ind2,:);
+ elseif((t1(ind2)-T)>=0)
+ ind1 = ind2-winlng+1;
+ yk = y(ind1:ind2,:);
+ end
+
+ for j = 1:n_y
+ for ki = 1:n_vf
+ Summe(ki,1) = trapz(tau,N(:,ind2,j,ki).*yk(:,j)) + Summe(ki,1);
+ end
+ end
+ NY(ind2,:) = Summe;
+ end
+
+ % NY(t,r_sum)
+
+
+ %% calculating and ploting the fault f_1, f_2,...., f_nf
+ % FHAT is the result of the equation for the coefficients
+ % vartheta = -<n,y>- <m,Bu>
+
+FHAT = NY + mB;
+
+ for k = 1:length(t1)
+ f_Hat(k,:) = P_f(k)*FHAT(k,:);
+ end
+
+
+
+ for i = 1:n_f
+ subplot(n_f,1,i)
+ hold on;
+ plot(t1,f(:,i),'-.','LineWidth',3);
+ plot(t1,f_Hat(:,i),'LineWidth',3);
+ hold off;
+ end
\ No newline at end of file
diff --git a/+LTV_signalmodelfault/calculateStates.m b/+LTV_signalmodelfault/calculateStates.m
new file mode 100644
index 0000000000000000000000000000000000000000..4322ff44bb37ba561787385e4b6842ca5080a4b6
--- /dev/null
+++ b/+LTV_signalmodelfault/calculateStates.m
@@ -0,0 +1,39 @@
+function x = calculateStates(dt,tau,Phitau,B,u)
+%CALCULATESTATES Numeric solution of the state space equation
+% x(t) = Phi(t,t0)[x0 + int( Phi(t0,tau) * B(tau) * u(tau),tau,0,t)]
+
+u = squeeze(u).';
+ns = size(B,1);
+ for k = 1:length(tau)
+
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+
+ Phit0t = inv(Phitt0(:,:,k));
+
+ temp1(:,k) = Phit0t * B(:,:,k) * u(:,k);
+
+ end
+ temp2 = temp1.';
+
+
+ states = zeros(length(tau),ns);
+
+
+
+ for k = 2:length(tau)
+ Phitt0(:,:,k) = Phitau(k,1,:,:);
+
+ states(k,:) = trapz(0:dt:tau(k),temp2(1:k,:));
+
+ temp3 = states.';
+
+
+ x(:,k) = Phitt0(:,:,k) * temp3(:,k);
+
+
+ end
+
+
+
+end
+
diff --git a/+LTV_signalmodelfault/kernTransform.m b/+LTV_signalmodelfault/kernTransform.m
new file mode 100644
index 0000000000000000000000000000000000000000..1908621fbe87a9508f3218fcc5be0ea2bfa3f3fb
--- /dev/null
+++ b/+LTV_signalmodelfault/kernTransform.m
@@ -0,0 +1,29 @@
+function mB = kernTransform(t,t1,tau,m,B,u)
+%KERNTRANSFORM returns the result of < m,B*u >
+%
+winlng = size(tau,2);
+ for i = 1:length(t)
+
+ Bu(:,i) = B(t(i))*u(t(i));
+
+ end
+
+
+ for i = 1:length(t1)
+ ind1 = i;
+ ind2 = winlng + ind1 - 1;
+ Butemp = Bu(:,ind1:ind2);
+ for k = 1:length(tau)
+
+ mb(k,1) = (m(:,k,i).')*Butemp(:,k);
+ end
+
+
+
+ mB(i,:) = trapz(tau,mb);
+
+
+ end
+
+end
+
diff --git a/+LTV_signalmodelfault/kernTransformTwo.m b/+LTV_signalmodelfault/kernTransformTwo.m
new file mode 100644
index 0000000000000000000000000000000000000000..986e826572c464557e24d3552d2edc1a21be0869
--- /dev/null
+++ b/+LTV_signalmodelfault/kernTransformTwo.m
@@ -0,0 +1,36 @@
+function mB = kernTransformTwo(t,t1,tau,m,B,u)
+%KERNTRANSFORM returns the result of < m,B*u >
+%
+winlng = size(tau,2);
+n = size(B(1),1);
+T = tau(end);
+ for i = 1:length(t1)
+
+ Bu(:,i) = B(t1(i))*u(i);
+
+ end
+
+
+ for i = 1:length(t1)
+ ind2 = i;
+ Butemp = zeros(n,winlng);
+ diff = winlng - ind2;
+ if((t1(ind2)-T) < 0 )
+ ind1 = 1;
+ Butemp(:,(ind1+diff):winlng) = Bu(:,ind1:ind2);
+ elseif( (t1(ind2)-T) >= 0)
+ ind1 = ind2-winlng+1;
+ Butemp = Bu(:,ind1:ind2);
+ end
+
+ for k = 1:length(tau)
+
+ mb(k,1) = (m(:,k,i).')*Butemp(:,k);
+ end
+
+ mB(i,:) = trapz(tau,mb);
+
+ end
+
+end
+
diff --git a/+LTV_signalmodelfault/mBu.m b/+LTV_signalmodelfault/mBu.m
new file mode 100644
index 0000000000000000000000000000000000000000..e35b4b4df09a993c0dec3cc8b0c2c6c767878dbe
--- /dev/null
+++ b/+LTV_signalmodelfault/mBu.m
@@ -0,0 +1,35 @@
+%% test <m,Bu>
+function mB = mBu(t,t1,tau,MQP,B,u)
+ winlng = size(tau,2);
+ n = size(B(1),1);
+
+ d = size(MQP,4);
+ for i = 1:length(t1)
+
+ Bu(:,i) = B(t1(i))*u(t1(i));
+
+ end
+
+ Bu_E = cat(2,zeros(n,winlng-1),Bu);
+ Bu_k = zeros(n,winlng);
+ M = MQP(1:n,:,:,:);
+
+ for ind_t_1 = winlng:length(t)
+ ind_t_T = ind_t_1 - winlng + 1;
+
+ Bu_k = Bu_E(:,ind_t_T:ind_t_1);
+ Summe = zeros(d,1);
+ ind_t_T = ind_t_1 - winlng +1;
+
+ for j = 1:n
+
+ for ki = 1:d
+ m_temp = squeeze(M(j,:,:,ki));
+ Summe(ki,1) = trapz(tau,(m_temp(:,ind_t_T).').*Bu_k(j,:)) + Summe(ki,1);
+ end
+ end
+ mB(ind_t_T,:) = Summe.';
+
+
+ end
+end
\ No newline at end of file
diff --git a/+LTV_signalmodelfault/simulateLTV.m b/+LTV_signalmodelfault/simulateLTV.m
new file mode 100644
index 0000000000000000000000000000000000000000..24d881dc067300c45e0826ed6295f1c985ddb3d0
--- /dev/null
+++ b/+LTV_signalmodelfault/simulateLTV.m
@@ -0,0 +1,19 @@
+function [t, y, x] = simulateLTV(t, A, B, C, E,E2, G, u, f,dis)
+% simulates the output of linear time variant System
+ n = size(A(1),1);
+
+ x = ones(n,1);
+ opts = odeset('RelTol',1e-6,'AbsTol',1e-9);
+ [t, x] = ode45(@(t,x) A(t)*x + B(t)*u(t) + E(t)*f(t) + G(t)*dis(t), t, zeros(size(A(1),1),1), opts);
+
+
+ for k = 1:length(t)
+ y(k,:) = C(t(k))*(x(k,:)') + E2(t(k))*f(t(k));
+ end
+
+end
+
+
+
+%function dx = ODE(t, x, A, B, E, u, f)
+ %dx = A(t) * x + B(t) * u(t);