function [A_f, B_f, C_f, D_f] = ...
    createInverseFeedthrough(A, B, C, Aa, Ba, L1, L2, K, T, ...
        useObservers)
    
% Full design model
if useObservers
    sys = designTSOB(A, B, C, Aa, Ba, L1, L2, K, T);
else    
    sys = designTS(A, B, C, Aa, Ba, L1, L2, T);
end
A_d = sys.A;
B_d = sys.B;
C_d = sys.C;
D_d = sys.D;

fprintf('Checking zeros of the design model\n');
tzeros = tzero(A_d, B_d, C_d, D_d);
for tz = tzeros'
    if abs(tz) > 1
        fprintf('WARNING: %.4f + %.4fi is outside the unit circle\n', ...
            real(tz), imag(tz)); 
    end
end

fprintf('Calculating inverse feedfoward filter\n');

if T == 0
    % Continuous
    for d = 1:10
        condition = cond(C_d*A_d^(d-1)*B_d);
        if (condition >0 && condition ~= Inf)
            fprintf('delay = %d cond = %f\n',d,condition);
            %break;
        end
        % Calculate inverse feedthrough
        D_f = inv(C_d*A_d^(d-1)*B_d);
        C_f = -D_f*C_d*A_d^d;
        B_f = B_d*D_f;
        A_f = A_d + B_d*C_f;
        stable = 1;
        for la = eig(A_f)'
            if abs(la) > 1
                stable = 0;
            end
        end
        if stable == 1
            break;
        end
    end
else
    for d = 2:100
        condition = cond(C_d*A_d^(d-1)*B_d);
        if (condition >0 && condition ~= Inf)
            fprintf('delay = %d cond = %f\n',d,condition);
        end
        % Calculate inverse feedthrough
        D_f = inv(C_d*A_d^(d-1)*B_d);
        C_f = -D_f*C_d*A_d^d;
        B_f = B_d*D_f;
        A_f = A_d + B_d*C_f;
        stable = true;
        eig(A_f);
        for e = eig(A_f)'
            if abs(e) >= 1
                stable = false;
            end
        end
        if stable == true
            break;
        end
    end
end