function [P, info] = build_tl_model(G, Fr, Fd, Fn, tau, args)
    %BUILD_PLL_MODEL Constructs a generalized auto-correlation transmission line plant [z, e] = P * [w, u] from input parameters.
    %   Uses default parameters if not overwritten by function arguments.
    
    arguments
        G (1,1) lti = zpk([], 0, db2mag(80))
        Fr (1,1) lti = tf([1, 10, 100e-3], [1, 0, 0], TimeUnit='milliseconds')  % build_powerlaw_tf(-60, [1e4, 1e1]);
        Fd (1,1) lti = tf([3.162, 316.2, 31.62], [1, 0, 0], TimeUnit='milliseconds')  % build_powerlaw_tf(-50, [1e5, 1e2]);
        Fn (1,1) lti = tf(db2mag(-70))
        tau (1,1) double = 1e-3
        args.Bandlimit (1,1) double = inf
        args.Rolloff (1,1) double = 0
        args.TimeUnit (1,1) string = ''
        args.Tolerance = 1e-9
    end
    
    tol = args.Tolerance;
    
    if tol > 0
        G = minreal(G, tol, false);
        Fr = minreal(Fr, tol, false);
        Fd = minreal(Fd, tol, false);
        Fn = minreal(Fn, tol, false);
    end
    
    time_unit = args.TimeUnit;
    if ~strlength(time_unit), time_unit = G.TimeUnit; end
    time_unit = ltipack.matchKey(time_unit,ltipack.getValidTimeUnits());
    if isempty(time_unit), error(message('Control:ltiobject:chgunit1')); end
    tunit_fctr = tunitconv('seconds', time_unit);
    
    bandlimit = 2*pi * args.Bandlimit / tunit_fctr;
    rolloff = 2*pi * args.Rolloff  / tunit_fctr;

    
    %% Models
    %- Helpers
    info.rename_ex_inputs = ss(eye(3), InputName='w', OutputName={'wr', 'wd', 'wn'});
    
    %- Voltage-controlled oscillator phase model
    G = chgTimeUnit(G, time_unit);
    G.InputName = 'u'; G.OutputName = 'yy';
    
    %- Noise colouring filters
    %-- Reference noise
    Fr = chgTimeUnit(Fr, time_unit);
    Fr = ensure_limits(Fr, bandlimit, rolloff);
    Fr.InputName = 'wr'; Fr.OutputName = 'r';
    
    %-- Local noise
    Fd = chgTimeUnit(Fd, time_unit);
    Fd = ensure_limits(Fd, bandlimit, rolloff);
    Fd.InputName = 'wd'; Fd.OutputName = 'd';
    
    %-- Measurement noise
    bandlimit = inf;  % Trick: We don't need to enforce bandlimit on Fn and don't need the variable hereon.
    Fn = chgTimeUnit(Fn, time_unit);
    Fn = ensure_limits(Fn, bandlimit, rolloff);
    Fn.InputName = 'wn'; Fn.OutputName = 'n';

    %-- Delay
    dly = tf(1, 1, InputDelay=tau);
    dly = chgTimeUnit(dly, time_unit);
    ddly = tf(1, 1, InputDelay=2*dly.InputDelay, TimeUnit=time_unit);
    tau_unc = ureal("tau", 0); % place holder uncertainty must be already normalized!
    % dly_ph = tf(tau_unc, 1); % placeholder for lft decomposition

    % build manual LFR
    warning('off', 'Control:combination:OutputNameClash');
    lfr.M = [0, 0, Fr, 0, 0, G;
             1, 0, 0, 2*Fd, 0, 0;
             -1, 0, Fr, -Fd, 0, 0;
             0, -1, Fr, 0, Fn, -G];
    lfr.M.InputName = ["", "", "w(1)", "w(2)", "w(3)", "u"];
    lfr.M.OutputName = ["", "", "z", "e"];
    lfr.Delta = blkdiag(tau_unc, tau_unc);
    lfr.BlkStruct = struct(Name='tau', Size=[1 1], Type='ureal', Occurrences=2, Simplify=2);

    %% Generalized plant
    %- Interconnection [z, e] = G * [wr, wd, wn, u];
    %P_ = connect(Fr, Fd, Fn, G, sumblk('y = yy + d'), sumblk('ym = y + n'), sumblk('e = r - ym'), sumblk('z = r - y'), rename_ex_inputs, {'w', 'u'}, {'z', 'e'});
    %P = [(1 - dly_ph) * Fr, -Fd, 0, -dly_ph * G;
    %     (1 - dly_ph*dly_ph) * Fr, -2 * dly_ph * Fd, Fn, -(1 + dly_ph*dly_ph) * G];
    %[lfr.M, lfr.Delta, lfr.BlkStruct] = lftdata(P, "tau");
    P = [(1 - dly) * Fr, -Fd, 0, -dly * G;
         (1 - ddly) * Fr, -2 * dly * Fd, Fn, -(1 + ddly) * G];  % constructing interconnection with tf/zpk models appears to work better than with ss wrt later numerical properties.
    P = lft(blkdiag(dly, dly), lfr.M);
    warning('on', 'Control:combination:OutputNameClash');
    P.InputName = {'w(1)', 'w(2)', 'w(3)', 'u'}; P.OutputName = {'z', 'e'};
    
    %- Improve
    P = ss(P);
    % P = minreal(P, tol, false);  % Remove redundant states. NOTE: take care of the tolerance parameter.
    % P = truncate_coeff(ss(P), tol);  % Get rid of numerical noise.
    P = prescale(P);  % Improve numerical conditioning.
    
    %- Extract helper variables
    info.nx = order(P);
    widx = startsWith(P.InputName, 'w'); info.widx = widx; info.nw = sum(widx);
    uidx = startsWith(P.InputName, 'u'); info.uidx = uidx; info.nu = sum(uidx);
    zidx = startsWith(P.OutputName, 'z'); info.zidx = zidx; info.nz = sum(zidx);
    yidx = startsWith(P.OutputName, 'e'); info.yidx = yidx; info.ny = sum(yidx);
    
    info.A = P.A;
    info.Bw = P.B(:,widx);  % Also known as B1
    info.Bu = P.B(:,uidx);  % Also known as B
    info.Cz = P.C(zidx,:);  % Also known as C1
    info.Cy = P.C(yidx,:);  % Also known as C
    info.Dzw = P.D(zidx, widx); assert(~any(info.Dzw, 'all'), "Non-zero feedthrough from z to w. Closed-loop has infinite H2-norm.")
    info.Dzu = P.D(zidx, uidx);  % Also known as E or D1 or D12
    info.Dyw = P.D(yidx, widx);  % Also known as F or D2 or D21
    info.Dyu = P.D(yidx, uidx); assert(~any(info.Dyu, 'all'), "Non-zero feedthrough from y to u cannot be non-zero. Required transforms not implemented.")
    
    info.G = G;
    info.Fr = Fr;
    info.Fd = Fd;
    info.Fn = Fn;
    
    info.Pyw = P(yidx, widx);
    info.Pzu = P(zidx, uidx);

    info.lfr = lfr;
    info.dly = dly;
end