diff --git a/inst/@lti/__lti_group__.m b/inst/@lti/__lti_group__.m
index ec869ca6..1a273de4 100644
--- a/inst/@lti/__lti_group__.m
+++ b/inst/@lti/__lti_group__.m
@@ -54,6 +54,21 @@
## retlti.outgroup remains empty struct
endif
+ if (strcmpi (dim, "blkdiag"))
+ ## Block-diagonal append: I/O delays concatenate with zero cross-terms,
+ ## exactly like the state/b/c/d blocks. (Previously these were dropped,
+ ## which silently discarded delays through append/feedback/connect.)
+ retlti.indelay = [lti1.indelay(:); lti2.indelay(:)];
+ retlti.outdelay = [lti1.outdelay(:); lti2.outdelay(:)];
+ retlti.iodelay = blkdiag (lti1.iodelay, lti2.iodelay);
+ else
+ ## horzcat/vertcat/times share inputs and/or outputs; a general delay
+ ## merge is ambiguous, so keep the historical zeroing for those modes.
+ retlti.indelay = zeros (numel (retlti.inname), 1);
+ retlti.outdelay = zeros (numel (retlti.outname), 1);
+ retlti.iodelay = zeros (numel (retlti.outname), numel (retlti.inname));
+ endif
+
if (lti1.tsam == lti2.tsam)
retlti.tsam = lti1.tsam;
elseif (lti1.tsam == -1 && lti2.tsam > 0)
@@ -90,3 +105,21 @@
ret = cell2struct ([ca; cb], [fa; fb]);
endfunction
+
+
+%!test
+%! h = [tf(1,[1 1]); tf(1,[1 2])];
+%! assert (totaldelay (h), zeros (2, 1));
+
+%!test
+%! h1 = [tf(1,[1 1]); tf(1,[1 2])];
+%! assert (hasdelay (h1), false);
+%! h2 = [tf(1,[1 1]), tf(1,[1 2])];
+%! assert (hasdelay (h2), false);
+
+%!test
+%! h = [tf(1,[1 1]); tf(1,[1 2])];
+%! [p, m] = size (h);
+%! assert (size (h.InputDelay), [m, 1]);
+%! assert (size (h.OutputDelay), [p, 1]);
+%! assert (size (h.IODelay), [p, m]);
diff --git a/inst/@lti/__lti_keys__.m b/inst/@lti/__lti_keys__.m
index 5ffe42b3..8f0b9e79 100644
--- a/inst/@lti/__lti_keys__.m
+++ b/inst/@lti/__lti_keys__.m
@@ -28,6 +28,9 @@
## cell vector of lti-specific keys
keys = {"tsam";
+ "inputdelay";
+ "outputdelay";
+ "iodelay";
"inname";
"outname";
"ingroup";
@@ -38,6 +41,9 @@
## cell vector of lti-specific assignable values
vals = {"scalar (sample time in seconds)";
+ "m-by-1 real matrix (input delay: seconds for ct models, samples for dt models)";
+ "p-by-1 real matrix (output delay: seconds for ct models, samples for dt models)";
+ "p-by-m real matrix (I/O delay: seconds for ct models, samples for dt models)";
"m-by-1 cell vector of strings";
"p-by-1 cell vector of strings";
"struct with indices as fields";
diff --git a/inst/@lti/__lti_prune__.m b/inst/@lti/__lti_prune__.m
index 6fffdb4b..498e361d 100644
--- a/inst/@lti/__lti_prune__.m
+++ b/inst/@lti/__lti_prune__.m
@@ -60,6 +60,17 @@
lti.outname = lti.outname(out_idx);
lti.inname = lti.inname(in_idx);
+ ## delays follow the selected I/O (previously left at their old size,
+ ## harmless only while they were always zero)
+ if (! (ischar (out_idx) && strcmp (out_idx, ":")))
+ lti.outdelay = lti.outdelay(out_idx);
+ lti.iodelay = lti.iodelay(out_idx, :);
+ endif
+ if (! (ischar (in_idx) && strcmp (in_idx, ":")))
+ lti.indelay = lti.indelay(in_idx);
+ lti.iodelay = lti.iodelay(:, in_idx);
+ endif
+
endfunction
diff --git a/inst/@lti/__pade_order_vector__.m b/inst/@lti/__pade_order_vector__.m
new file mode 100644
index 00000000..256e5e7e
--- /dev/null
+++ b/inst/@lti/__pade_order_vector__.m
@@ -0,0 +1,87 @@
+## Copyright (C) 2026 Prateek Ganguli
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see .
+
+## Compute the per-nonzero-delay Pade approximation order vector for
+## pade(sys, n).
+##
+## Shared helper for @lti/pade.m (tf/zpk) and the future ss-specific pade
+## dispatches (Tasks 2/3) -- build the per-entry order list once here and
+## have every case index into it consistently.
+##
+## Ordering convention (this codebase's own, since MATLAB's internal
+## ordering isn't documented/verifiable): the total nonzero-delay count,
+## and the order in which "n" (if a vector) is consumed, walks:
+## 1. InputDelay nonzero entries, in input order.
+## 2. OutputDelay nonzero entries, in output order.
+## 3. IODelay nonzero entries, in column-major (Octave linear index) order.
+## 4. InternalDelay ports, in the existing internal port order returned
+## by get(sys, "internaldelay") (always empty for tf/zpk, since those
+## types can never carry InternalDelay).
+##
+## Inputs:
+## sys - an LTI model.
+## n - scalar or vector Pade order. If scalar, the same order is
+## used for every nonzero delay. If a vector, its length must
+## equal the total nonzero-delay count (computed per the ordering
+## above); anything else is an error.
+##
+## Outputs:
+## orders - column vector of per-delay-entry orders, length equal to the
+## total nonzero-delay count, in the documented order above.
+## idx - struct with fields "input", "output", "iod", "internal", each
+## a vector of linear/element indices (into InputDelay,
+## OutputDelay, IODelay(:), and internaldelay(:) respectively)
+## identifying which entries are nonzero, in the same order
+## "orders" was built in -- so callers can walk
+## idx.input/idx.output/idx.iod/idx.internal in lockstep with
+## orders(1:numel(idx.input)), orders(numel(idx.input)+1 : ...),
+## etc.
+function [orders, idx] = __pade_order_vector__ (sys, n)
+
+ if (nargin != 2 || ! isa (sys, "lti"))
+ print_usage ();
+ endif
+
+ [indelay, outdelay, iodelay] = get (sys, "inputdelay", "outputdelay", "iodelay");
+
+ ## InternalDelay is a ss-specific property; get (sys, "internaldelay")
+ ## errors for tf/zpk models, so only query it for ss (always empty for
+ ## tf/zpk, matching the documented ordering convention above).
+ if (isa (sys, "ss"))
+ internaldelay = get (sys, "internaldelay");
+ else
+ internaldelay = [];
+ endif
+
+ idx.input = find (indelay(:) != 0);
+ idx.output = find (outdelay(:) != 0);
+ idx.iod = find (iodelay(:) != 0);
+ idx.internal = find (internaldelay(:) != 0);
+
+ total = numel (idx.input) + numel (idx.output) + numel (idx.iod) + numel (idx.internal);
+
+ if (isscalar (n))
+ orders = repmat (n, total, 1);
+ else
+ if (numel (n) != total)
+ error ("pade: order vector length (%d) does not match the number of nonzero delays (%d)",
+ numel (n), total);
+ endif
+ orders = n(:);
+ endif
+
+endfunction
diff --git a/inst/@lti/c2d.m b/inst/@lti/c2d.m
index 2ccf17f4..826ff90d 100644
--- a/inst/@lti/c2d.m
+++ b/inst/@lti/c2d.m
@@ -76,6 +76,19 @@
error ("c2d: second argument is not a valid sample time");
endif
+ delay_modeling = "delay";
+
+ if (isstruct (method))
+ if (w0 != 0)
+ error ("c2d: cannot combine a c2dOptions struct with an explicit prewarp frequency argument");
+ endif
+
+ opt = method;
+ method = opt.Method;
+ w0 = opt.PrewarpFrequency;
+ delay_modeling = opt.DelayModeling;
+ endif
+
if (! ischar (method))
error ("c2d: third argument is not a string");
endif
@@ -84,8 +97,166 @@
error ("c2d: fourth argument is not a valid pre-warping frequency");
endif
- sys = __c2d__ (sys, tsam, lower (method), w0);
- sys.tsam = tsam;
+ origsys = sys;
+
+ if (hasinternaldelay (origsys))
+ ## Discretizing an InternalDelay system is exactly discretizing the
+ ## extended ordinary system (A, [B1 B2], [C1;C2], [[D11 D12];[D21 D22]])
+ ## via the unchanged __c2d__ (which only looks at matrix shapes), then
+ ## rounding tau to whole samples the same way InputDelay/OutputDelay are
+ ## rounded below.
+ ## __ss_ext_split__ copies InputDelay/OutputDelay/IODelay unchanged from
+ ## origsys (still in seconds), so if origsys also carries ordinary I/O
+ ## delay, round it to samples here, the same way the hasdelay() branch
+ ## below does for the InternalDelay-free case. Combining InternalDelay
+ ## with DelayModeling="state" or fractional (Thiran) ordinary delay is
+ ## out of scope; only the plain round-to-samples case is handled.
+ [ext_sys, nu, ny] = __ss_ext_build__ (origsys);
+ ext_sys = __c2d__ (ext_sys, tsam, lower (method), w0);
+ ext_sys.tsam = tsam;
+ tau = get (origsys, "internaldelay");
+ tau_samples = round (tau / tsam);
+
+ ## A nonzero InternalDelay that rounds to 0 samples would silently drop
+ ## that port's feedthrough during simulation instead of solving the
+ ## resulting algebraic loop (documented in __delay_lookup__/
+ ## __buffered_sim__'s tau>=1 assumption) -- convert this from a
+ ## silent-wrong-answer risk into a clear error, consistent with every
+ ## other guard on this branch.
+ if (any (tau != 0 & tau_samples == 0))
+ bad = tau (find (tau != 0 & tau_samples == 0, 1));
+ error (["c2d: InternalDelay of %g seconds is too small relative to ", ...
+ "the sampling time %g and would round to 0 samples"], ...
+ bad, tsam);
+ endif
+
+ if (strcmp (delay_modeling, "state"))
+ ## Close the discretized InternalDelay loop with an EXACT z^-k filter per
+ ## port (tau_samples is already a nonnegative integer -- no approximation
+ ## needed, unlike pade()'s Pade-filter closure for the continuous case),
+ ## instead of leaving InternalDelay as an explicit port. Mirrors
+ ## __pade_substitute_internal__'s ext_clean + lft() pattern in
+ ## @lti/pade.m: extract ext_sys's raw matrices via dssdata (this function
+ ## is an @lti method and cannot dot-access ext_sys's @ss-private fields
+ ## directly), rebuild a clean ordinary ss/dss, then close the loop.
+ [ea, eb, ec, ed, ee, etsam] = dssdata (ext_sys, []);
+ nports = columns (eb) - nu;
+
+ if (isempty (ee))
+ ext_clean = ss (ea, eb, ec, ed, etsam);
+ else
+ ext_clean = dss (ea, eb, ec, ed, ee, etsam);
+ endif
+
+ if (nports == 0)
+ ## Degenerate: internaldelay set but no actual delay ports -- nothing
+ ## to close, the clean plant already IS the delay-free equivalent.
+ sys = ext_clean;
+ else
+ G_close = __exact_discrete_delay_ss__ (tau_samples(1), tsam);
+ for kk = 2 : nports
+ G_close = append (G_close, __exact_discrete_delay_ss__ (tau_samples(kk), tsam));
+ endfor
+ sys = lft (ext_clean, G_close, nports, nports);
+ endif
+
+ if (hasdelay (origsys))
+ [indelay, outdelay, iodelay] = get (origsys, "inputdelay", "outputdelay", "iodelay");
+ sys = set (sys, "InputDelay", round (indelay / tsam), ...
+ "OutputDelay", round (outdelay / tsam), ...
+ "IODelay", round (iodelay / tsam));
+ sys = ss (absorbDelay (tf (sys)));
+ endif
+
+ return;
+ endif
+
+ sys = __ss_ext_split__ (origsys, ext_sys, nu, ny);
+ sys.tsam = tsam;
+
+ sys = set (sys, "internaldelay", tau_samples);
+
+ if (hasdelay (origsys))
+ [indelay, outdelay, iodelay] = get (origsys, "inputdelay", "outputdelay", "iodelay");
+ sys = set (sys, "InputDelay", round (indelay / tsam), ...
+ "OutputDelay", round (outdelay / tsam), ...
+ "IODelay", round (iodelay / tsam));
+ endif
+
+ return;
+ endif
+
+ if (hasdelay (origsys))
+ [indelay, outdelay, iodelay] = get (origsys, "inputdelay", "outputdelay", "iodelay");
+
+ indelay_samples = indelay / tsam;
+ outdelay_samples = outdelay / tsam;
+ iodelay_samples = iodelay / tsam;
+
+ all_samples = [indelay_samples(:); outdelay_samples(:); iodelay_samples(:)];
+
+ if (strcmp (delay_modeling, "state"))
+ % Delay must be zeroed before __c2d__ for the same recursive-re-entry
+ % reason documented below (ss() preserves delay fields), even though
+ % this branch never calls thiran -- otherwise a MIMO/recursive call
+ % into __c2d__ would re-enter this same "state" branch again.
+ sys_no_delay = set (origsys, "InputDelay", 0, "OutputDelay", 0, "IODelay", 0);
+
+ % absorbDelay errors on ss -- route ss inputs through tf, absorb
+ % there, then convert back, since absorbDelay's per-channel math
+ % already handles a general IODelay matrix correctly.
+ is_ss_input = isa (origsys, "ss");
+
+ if (is_ss_input)
+ sys = tf (sys_no_delay);
+ else
+ sys = sys_no_delay;
+ endif
+
+ sys = __c2d__ (sys, tsam, lower (method), w0);
+ sys.tsam = tsam;
+ sys = set (sys, "InputDelay", round (indelay_samples), ...
+ "OutputDelay", round (outdelay_samples), ...
+ "IODelay", round (iodelay_samples));
+ sys = absorbDelay (sys);
+
+ if (is_ss_input)
+ sys = ss (sys);
+ endif
+ elseif (any (abs (all_samples - round (all_samples)) > 1e-8))
+ [p, m] = size (origsys);
+
+ if (any (strcmpi (method, {"zoh", "std"})) && __c2d_frac_applicable__ (origsys, tsam))
+ sys = __c2d_frac_ss__ (origsys, tsam);
+ sys.tsam = tsam;
+ elseif (p != 1 || m != 1)
+ error ("c2d: fractional delays on MIMO systems are not yet supported (per-channel Thiran approximation is not yet implemented)");
+ else
+ % __c2d__ recurses into this same function via ss() for the general
+ % case, and ss() does not strip delay fields -- so discretizing
+ % origsys directly here would re-enter this fractional branch a
+ % second time and double-apply the Thiran filter below. Delay must
+ % be zeroed before the __c2d__ call, not just after.
+ sys_no_delay = set (origsys, "InputDelay", 0, "OutputDelay", 0, "IODelay", 0);
+ sys = __c2d__ (sys_no_delay, tsam, lower (method), w0);
+ sys.tsam = tsam;
+
+ total_delay = totaldelay (origsys);
+ filt = thiran (total_delay, tsam);
+ sys = sys * filt;
+ sys = set (sys, "InputDelay", 0, "OutputDelay", 0, "IODelay", 0);
+ endif
+ else
+ sys = __c2d__ (origsys, tsam, lower (method), w0);
+ sys.tsam = tsam;
+ sys = set (sys, "InputDelay", round (indelay_samples), ...
+ "OutputDelay", round (outdelay_samples), ...
+ "IODelay", round (iodelay_samples));
+ endif
+ else
+ sys = __c2d__ (origsys, tsam, lower (method), w0);
+ sys.tsam = tsam;
+ endif
endfunction
@@ -339,3 +510,501 @@
%!assert (Aexint, Aexint_exp, 1e-4);
+%!test # exact integer-sample InputDelay survives c2d
+%! sys = tf (1, [1 1], "InputDelay", 1.0);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (dsys.InputDelay, 2);
+%! assert (dsys.OutputDelay, 0);
+%! assert (dsys.IODelay, 0);
+
+%!test # no delay: unaffected (regression)
+%! sys = tf (1, [1 1]);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (hasdelay (dsys), false);
+
+
+%!test # SISO tf fractional InputDelay: exact split-ZOH via zoh (default)
+%! sys = tf (1, [1 1], "InputDelay", 1.33);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! ## d = floor(tau/tsam) = floor(1.33/0.5) = 2, plus 1 more because the
+%! ## fractional remainder (0.33 s) is nonzero -- same convention as
+%! ## __c2d_frac_zoh__/__c2d_frac_ss__.
+%! assert (hasdelay (dsys), true);
+%! assert (dsys.InputDelay, 3);
+%! [numd, dend] = tfdata (dsys, "v");
+%! ## independent ground truth: split-ZOH construction (Astrom & Wittenmark),
+%! ## computed directly via expm/integral (not via the code under test) --
+%! ## see the reproduction script in the task report for how these numbers
+%! ## were derived: A=-1, B=1, tsam=0.5, tau=1.33, d=floor(tau/tsam)=2,
+%! ## tau_frac=0.33.
+%! nume = [0.1563351834, 0.2371341569];
+%! dene = [1, -0.6065306597];
+%! assert (numd, nume, 1e-6);
+%! assert (dend, dene, 1e-6);
+%!
+%! ## Thiran approximation path is still reachable for a non-zoh method.
+%! dsys_tustin = c2d (sys, 0.5, "tustin");
+%! assert (hasdelay (dsys_tustin), false);
+
+%!test # SISO zpk fractional OutputDelay: exact split-ZOH via zoh (default)
+%! sys = zpk ([], -1, 1, "OutputDelay", 1.33);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (hasdelay (dsys), true);
+%! assert (dsys.InputDelay, 3);
+%! w = [0.1, 1, 5];
+%! ## independent ground truth: same split-ZOH construction as the tf
+%! ## InputDelay test above -- the dynamics (pole at -1, gain 1) are
+%! ## identical, since __c2d_frac_ss__ collapses all delay fields into a
+%! ## scalar totaldelay () before splitting, so the same reference numbers
+%! ## apply.
+%! nume = [0.1563351834, 0.2371341569];
+%! dene = [1, -0.6065306597];
+%! tsam = 0.5;
+%! expected = tf (nume, dene, tsam);
+%! expected = set (expected, "InputDelay", 3);
+%! assert (freqresp (dsys, w), freqresp (expected, w), 1e-6);
+%!
+%! ## Thiran approximation path is still reachable for a non-zoh method.
+%! dsys_tustin = c2d (sys, 0.5, "tustin");
+%! assert (hasdelay (dsys_tustin), false);
+
+%!test # SIMO (single input, multiple outputs) fractional InputDelay:
+%! ## Finding 1 fix-pass regression test. Before the fix, __c2d_frac_ss__'s
+%! ## trailing tf/zpk pole/zero-at-origin cancellation block was guarded
+%! ## only by "m == 1", which is trivially true here (m == 1) even though
+%! ## the system is not truly SISO (p == 2) -- the cancellation code calls
+%! ## tfdata(sys, "v") in a way that assumes a genuinely SISO system, and
+%! ## crashed with "ERROR: abs: not defined for cell". The guard is now
+%! ## "p == 1 && m == 1", so this case is correctly routed through the
+%! ## general per-column MIMO augmentation path instead (no cancellation
+%! ## attempted). Verify no crash, and cross-check each output channel's
+%! ## response against an independently, separately discretized SISO
+%! ## channel.
+%! sys = tf ({1;1}, {[1 1];[1 2]}, "InputDelay", 1.33);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! h1 = tf (1, [1 1], "InputDelay", 1.33);
+%! h2 = tf (1, [1 2], "InputDelay", 1.33);
+%! d1 = c2d (h1, 0.5, "zoh");
+%! d2 = c2d (h2, 0.5, "zoh");
+%! w = [0.1, 1, 5];
+%! resp = freqresp (dsys, w);
+%! assert (reshape (resp(1,1,:), 1, []), reshape (freqresp (d1, w), 1, []), 1e-10);
+%! assert (reshape (resp(2,1,:), 1, []), reshape (freqresp (d2, w), 1, []), 1e-10);
+
+%!test # MIMO fractional InputDelay, uniform per column: now handled exactly
+%! ## by the Task 2 per-column split-ZOH generalization (this used to error
+%! ## before MIMO support was added -- __c2d_frac_applicable__ now accepts
+%! ## it because each column's total delay is constant down that column).
+%! ## InputDelay is [2;0], not [3;0]: the fix-pass correction removed the
+%! ## blanket "+1" applied to every fractional-delay column, since the
+%! ## extra register state remains part of the augmented dynamics here (no
+%! ## tf/zpk pole-zero cancellation happens for p>1) and its own
+%! ## eigenvalue-zero pole already contributes exactly the missing sample.
+%! sys = tf ({1,1;1,1}, {[1 1],[1 2];[1 3],[1 4]}, "InputDelay", [1.33;0]);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (dsys.InputDelay, [2;0]);
+
+%!test # MIMO: column 1 exact-integer (nonzero) delay, column 2 fractional --
+%! ## only column 2 gets an extra register state; column 1's response must
+%! ## match a plain (no-extra-state) single-channel zoh discretization
+%! ## exactly, confirming the exact-integer column is untouched by the
+%! ## per-column augmentation applied to its fractional sibling.
+%! A = [-1, 0; 0, -3];
+%! B = [1, 0; 0, 1];
+%! C = [1, 0; 0, 1];
+%! D = [0, 0; 0, 0];
+%! sys = ss (A, B, C, D);
+%! sys = set (sys, "InputDelay", [0.2; 0.37]); # col1: exact 2 samples @ Ts=0.1
+%! tsam = 0.1;
+%! dsys = c2d (sys, tsam);
+%! assert (dsys.InputDelay(1), 2);
+%! assert (dsys.InputDelay(2), 3); # floor(0.37/0.1)=3; no +1, since the
+%! # register survives as a real state in
+%! # this ss/MIMO output (only cancelled-out
+%! # tf/zpk SISO conversions get the +1)
+%! sys1 = ss (A(1,1), B(1,1), C(1,1), D(1,1));
+%! sys1 = set (sys1, "InputDelay", 0.2);
+%! dsys1 = c2d (sys1, tsam);
+%! w = [0.1, 1, 5];
+%! resp_full = freqresp (dsys, w);
+%! resp_ref = freqresp (dsys1, w);
+%! assert (squeeze (resp_full(1,1,:)), squeeze (resp_ref), 1e-8);
+
+%!test # MIMO fractional OutputDelay, uniform across all outputs: in-scope
+%! ## (Finding 2 fix-pass): OutputDelay may be nonzero for MIMO as long as
+%! ## it is uniform across all outputs, since that is mathematically
+%! ## equivalent to attributing the same delay to InputDelay instead.
+%! ## Verify this succeeds (not "not yet supported") and matches the
+%! ## InputDelay-based formulation exactly (ssdata/freqresp).
+%! num = {1,1;1,1}; den = {[1 1],[1 2];[1 3],[1 4]};
+%! sys_out = tf (num, den, "OutputDelay", [1.33;1.33]);
+%! sys_in = tf (num, den, "InputDelay", [1.33;1.33]);
+%! d_out = c2d (sys_out, 0.5, "zoh");
+%! d_in = c2d (sys_in, 0.5, "zoh");
+%! assert (isdt (d_out), true);
+%! [Ao, Bo, Co, Do] = ssdata (d_out);
+%! [Ai, Bi, Ci, Di] = ssdata (d_in);
+%! assert (Ao, Ai, 1e-10);
+%! assert (Bo, Bi, 1e-10);
+%! assert (Co, Ci, 1e-10);
+%! assert (Do, Di, 1e-10);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (d_out, w), freqresp (d_in, w), 1e-10);
+
+%!error ...
+%! ## row-varying OutputDelay (differs across outputs) is genuinely out of
+%! ## scope -- it cannot be pushed through to the input side -- and must
+%! ## still fall through to the "not yet supported" error.
+%! c2d (tf ({1,1;1,1}, {[1 1],[1 2];[1 3],[1 4]}, "OutputDelay", [1.33;0.5]), 0.5, "zoh")
+
+%!error ...
+%! ## row-varying IODelay within a column is likewise genuinely out of
+%! ## scope and must still error.
+%! sys = ss ([-1, 0; 0, -2], [1, 0; 0, 1], [1, 0; 0, 1], [0, 0; 0, 0]);
+%! sys = set (sys, "IODelay", [1.33, 0; 0.5, 0]);
+%! c2d (sys, 0.5, "zoh")
+
+
+%!test # c2dOptions struct as 3rd argument, default DelayModeling
+%! opt = c2dOptions ();
+%! sys = tf (1, [1 1], "InputDelay", 1.0);
+%! dsys = c2d (sys, 0.5, opt);
+%! assert (dsys.InputDelay, 2);
+
+%!test # DelayModeling='state': fractional SISO delay rounds to integer samples, absorbed into extra dynamics
+%! opt = c2dOptions ("DelayModeling", "state");
+%! sys = tf (1, [1 1], "InputDelay", 1.33);
+%! dsys = c2d (sys, 0.5, opt);
+%! assert (isdt (dsys), true);
+%! assert (hasdelay (dsys), false);
+%! rational_dsys = c2d (tf (1, [1 1]), 0.5, "zoh");
+%! k = round (1.33 / 0.5);
+%! w = [0.1, 1, 5];
+%! expected = freqresp (rational_dsys, w) .* reshape (exp (-1i * w * 0.5 * k), 1, 1, []);
+%! assert (freqresp (dsys, w), expected, 1e-8);
+
+%!test # DelayModeling='state': MIMO fractional delay succeeds (unlike 'delay' mode), absorbed per channel
+%! opt = c2dOptions ("DelayModeling", "state");
+%! sys = tf ({1,1;1,1}, {[1 1],[1 2];[1 3],[1 4]}, "InputDelay", [1.33;0]);
+%! dsys = c2d (sys, 0.5, opt);
+%! assert (hasdelay (dsys), false);
+%! rational_dsys = c2d (tf ({1,1;1,1}, {[1 1],[1 2];[1 3],[1 4]}), 0.5, "zoh");
+%! total = [round(1.33/0.5), 0; round(1.33/0.5), 0];
+%! w = [0.1, 1, 5];
+%! resp = freqresp (rational_dsys, w);
+%! expected = zeros (2, 2, numel (w));
+%! for i = 1:2
+%! for j = 1:2
+%! expected(i,j,:) = reshape (resp(i,j,:), 1, []) .* exp (-1i * w * 0.5 * total(i,j));
+%! endfor
+%! endfor
+%! assert (freqresp (dsys, w), expected, 1e-8);
+
+%!test # DelayModeling='state': ss input absorbs InputDelay into extra dynamics
+%! sys = ss (-1, 1, 1, 0, "InputDelay", 0.65);
+%! opt = c2dOptions ("DelayModeling", "state");
+%! dsys = c2d (sys, 0.5, opt);
+%! assert (isa (dsys, "ss"), true);
+%! assert (isdt (dsys), true);
+%! assert (hasdelay (dsys), false);
+%! rational_dsys = c2d (ss (-1, 1, 1, 0), 0.5, "zoh");
+%! k = round (0.65 / 0.5);
+%! w = [0.1, 1, 5];
+%! expected = freqresp (rational_dsys, w) .* reshape (exp (-1i * w * 0.5 * k), 1, 1, []);
+%! assert (freqresp (dsys, w), expected, 1e-8);
+
+%!test # DelayModeling='state': ss input with per-channel IODelay matrix (MIMO)
+%! A = [-1, 0; 0, -2];
+%! B = [1, 0; 0, 1];
+%! C = [1, 0; 0, 1];
+%! D = [0, 0; 0, 0];
+%! sys = ss (A, B, C, D);
+%! sys = set (sys, "IODelay", [0.5, 0; 0, 1.0]);
+%! opt = c2dOptions ("DelayModeling", "state");
+%! dsys = c2d (sys, 0.5, opt);
+%! assert (hasdelay (dsys), false);
+%! rational_dsys = c2d (ss (A, B, C, D), 0.5, "zoh");
+%! total = [round(0.5/0.5), 0; 0, round(1.0/0.5)];
+%! w = [0.1, 1, 5];
+%! resp = freqresp (rational_dsys, w);
+%! expected = zeros (2, 2, numel (w));
+%! for i = 1:2
+%! for j = 1:2
+%! expected(i,j,:) = reshape (resp(i,j,:), 1, []) .* exp (-1i * w * 0.5 * total(i,j));
+%! endfor
+%! endfor
+%! assert (freqresp (dsys, w), expected, 1e-8);
+
+%!error c2d (tf (1, [1 1]), 0.5, c2dOptions (), 100)
+
+%!test # exact-integer-sample InternalDelay: feedback()-produced system, tau/tsam is a whole number
+%! T = 0.5;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! assert (hasinternaldelay (L), true);
+%! dsys = c2d (L, 0.25, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (hasinternaldelay (dsys), true);
+%! assert (dsys.internaldelay, T / 0.25, 1e-10);
+
+%!test # fractional InternalDelay: rounds to nearest sample
+%! T = 0.33;
+%! tsam = 0.1;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! dsys = c2d (L, tsam, "zoh");
+%! assert (hasinternaldelay (dsys), true);
+%! assert (dsys.internaldelay, round (T / tsam));
+
+%!test # zoh sanity check: discretized InternalDelay system tracks continuous freqresp at low frequency
+%! T = 0.3;
+%! G = ss (-2, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! tsam = 0.01;
+%! dsys = c2d (L, tsam, "zoh");
+%! w = [0.05, 0.2, 1];
+%! Hc = freqresp (L, w);
+%! Hd = freqresp (dsys, w);
+%! assert (Hd, Hc, 5e-2);
+
+%!test # InternalDelay with no delay ports (b2/c2/d12/d21/d22 empty): degenerate extended system still discretizes
+%! sys = set (ss (-1, 1, 1, 0), "internaldelay", 0.5);
+%! dsys = c2d (sys, 0.5, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (dsys.internaldelay, 1);
+
+%!test # combined InternalDelay AND ordinary InputDelay: both rounded to samples
+%! T = 0.5;
+%! tsam = 0.25;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! L = set (L, "InputDelay", 0.75);
+%! assert (hasinternaldelay (L), true);
+%! assert (hasdelay (L), true);
+%! dsys = c2d (L, tsam, "zoh");
+%! assert (isdt (dsys), true);
+%! assert (hasinternaldelay (dsys), true);
+%! assert (dsys.internaldelay, T / tsam, 1e-10);
+%! assert (dsys.InputDelay, round (0.75 / tsam));
+
+
+%!test # MIMO InternalDelay: two channels, DIFFERENT delays, dimension-correct multi-port c2d
+%! ## Genuine 2-in/2-out InternalDelay fixture built as append() of two
+%! ## independent SISO feedback loops (the isolated-SISO-block topology
+%! ## __sys_connect__ actually supports). Distinct delays (0.4 vs 0.8) so an
+%! ## index/port swap could not hide behind equal values.
+%! T1 = 0.4; T2 = 0.8; tsam = 0.2;
+%! G1 = ss (-1, 1, 1, 0, "IODelay", T1);
+%! G2 = ss (-2, 1, 1, 0, "IODelay", T2);
+%! sys = append (feedback (G1), feedback (G2));
+%! assert (get (sys, "internaldelay"), [T1; T2], 1e-12); # two channels
+%! dsys = c2d (sys, tsam, "zoh");
+%! assert (hasinternaldelay (dsys), true);
+%! assert (isdt (dsys), true);
+%! assert (get (dsys, "internaldelay"), [T1/tsam; T2/tsam], 1e-10); # 2 and 4 samples
+%! assert (size (dsys), [2, 2]);
+%! ## Independent reference: each diagonal channel must equal the SISO c2d of
+%! ## that loop alone; off-diagonals must be exactly zero (channels decoupled).
+%! d1 = c2d (feedback (G1), tsam, "zoh");
+%! d2 = c2d (feedback (G2), tsam, "zoh");
+%! w = [0.1, 0.7, 2];
+%! H = freqresp (dsys, w); H1 = freqresp (d1, w)(:); H2 = freqresp (d2, w)(:);
+%! for k = 1:numel (w)
+%! assert (H(1,1,k), H1(k), 1e-12);
+%! assert (H(2,2,k), H2(k), 1e-12);
+%! assert (H(1,2,k), 0, 1e-12);
+%! assert (H(2,1,k), 0, 1e-12);
+%! endfor
+
+## InternalDelay that rounds to 0 samples must error clearly, not silently
+## drop the port's feedthrough (see the tau>=1 assumption documented in
+## __delay_lookup__/__buffered_sim__).
+%!error c2d (feedback (ss (-1, 1, 1, 0, "IODelay", 0.01)), 1, "zoh")
+
+%!test # exact fractional IODelay via zoh: matches MATLAB's split-ZOH result
+%! h = tf (10, [1 3 10], "IODelay", 0.25);
+%! hd = c2d (h, 0.1); # method defaults to "zoh"
+%! [num, den] = tfdata (hd);
+%! assert (hd.InputDelay + hd.IODelay, 3); # z^-3 pure delay factored out
+%! assert (num{1}, [0.01187, 0.06408, 0.009721], 1e-4);
+%! assert (den{1}, [1, -1.655, 0.7408], 1e-3);
+
+%!test # DelayModeling="state" on InternalDelay: port removed, freqresp matches
+%! # the default ("delay") mode's own result exactly (both are exact, just
+%! # different representations -- explicit port vs. absorbed states).
+%! T = 0.3; tsam = 0.1;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! opt_state = c2dOptions ("DelayModeling", "state");
+%! dsys_state = c2d (L, tsam, opt_state);
+%! dsys_delay = c2d (L, tsam);
+%! assert (hasinternaldelay (dsys_state), false);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (dsys_state, w), freqresp (dsys_delay, w), 1e-8);
+
+%!test # DelayModeling="state" combined InternalDelay + ordinary InputDelay:
+%! # both absorbed, hasdelay and hasinternaldelay both false
+%! T = 0.3; tsam = 0.1;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = set (feedback (G), "InputDelay", 0.2);
+%! opt_state = c2dOptions ("DelayModeling", "state");
+%! dsys_state = c2d (L, tsam, opt_state);
+%! assert (hasinternaldelay (dsys_state), false);
+%! assert (hasdelay (dsys_state), false);
+%! w = [0.1, 1, 5];
+%! dsys_delay = c2d (L, tsam);
+%! assert (freqresp (dsys_state, w), freqresp (dsys_delay, w), 1e-6);
+
+%!test # DelayModeling="state" MIMO InternalDelay: two independent ports
+%! T1 = 0.3; T2 = 0.5; tsam = 0.1;
+%! G1 = ss (-1, 1, 1, 0, "IODelay", T1);
+%! G2 = ss (-2, 1, 1, 0, "IODelay", T2);
+%! sys = append (feedback (G1), feedback (G2));
+%! opt_state = c2dOptions ("DelayModeling", "state");
+%! dsys_state = c2d (sys, tsam, opt_state);
+%! assert (hasinternaldelay (dsys_state), false);
+%! w = [0.1, 1, 5];
+%! dsys_delay = c2d (sys, tsam);
+%! assert (freqresp (dsys_state, w), freqresp (dsys_delay, w), 1e-6);
+
+%!test # default ("delay") mode on InternalDelay: unaffected regression check
+%! T = 0.3; tsam = 0.1;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! dsys = c2d (L, tsam);
+%! assert (hasinternaldelay (dsys), true);
+%! assert (dsys.internaldelay, round (T/tsam));
+
+
+function tf_ok = __c2d_frac_applicable__ (sys, tsam)
+ [outdelay] = get (sys, "outputdelay");
+ total = totaldelay (sys); # p-by-m matrix, seconds
+ [p, m] = size (total);
+ ## A single-output system's OutputDelay commutes with its InputDelay(s)
+ ## (it is equivalent, for a scalar-output LTI system, to shifting every
+ ## input by the same amount), so it can always be folded into the
+ ## per-column totaldelay used below -- this is exactly what Task 1's
+ ## SISO OutputDelay regression test relies on. For p > 1 a nonzero,
+ ## per-row-varying OutputDelay would make the per-column delay
+ ## inconsistent across outputs; that case is already rejected by the
+ ## uniformity check below, so no separate early-return is needed here.
+ ##
+ ## Deliberate scope decision (not just an accepted side effect): for
+ ## MIMO (p > 1), OutputDelay may be nonzero as long as it is *uniform*
+ ## across all outputs -- that is mathematically equivalent to
+ ## attributing the same delay to InputDelay instead (every output being
+ ## delayed by the same amount commutes with delaying every input by that
+ ## amount), and is folded into the per-column totaldelay uniformity
+ ## check below exactly like any other InputDelay/IODelay contribution.
+ ## Only row-VARYING OutputDelay is rejected here, since that would
+ ## require genuine output-side treatment (delaying different outputs by
+ ## different amounts cannot be pushed through to the input side), which
+ ## this fix does not implement -- such systems fall through to the
+ ## "not yet supported" error below.
+ if (p > 1 && any (outdelay != outdelay(1)))
+ tf_ok = false;
+ return;
+ endif
+ for j = 1 : m
+ col = total (:, j);
+ if (any (abs (col - col(1)) > 1e-10))
+ tf_ok = false;
+ return;
+ endif
+ endfor
+ tf_ok = true;
+endfunction
+
+function sys = __c2d_frac_ss__ (origsys, tsam)
+ sys_ss = ss (origsys);
+ [A, B, C, D] = ssdata (sys_ss);
+ total = totaldelay (origsys);
+ [p, m] = size (total);
+ n = rows (A);
+
+ col_tau = total (1, :); # one delay value per column (uniform down
+ # each column, guaranteed by the
+ # applicability check)
+
+ extra_needed = false (1, m);
+ d_samples = zeros (1, m);
+ Bd_cols = cell (1, m);
+ g0_cols = cell (1, m);
+
+ for j = 1 : m
+ [Ad, Bd, extra] = __c2d_frac_zoh__ (A, B(:,j), tsam, col_tau(j));
+ d_samples(j) = extra.d;
+ Bd_cols{j} = Bd;
+ extra_needed(j) = extra.needed;
+ if (extra.needed)
+ g0_cols{j} = extra.g0;
+ ## NOTE: no "+1" here. The extra register state remains part of the
+ ## augmented dynamics (A_aug/B_aug below) and its own eigenvalue-zero
+ ## pole already contributes exactly one sample of delay to the
+ ## transfer function; adding 1 to the explicit InputDelay field here
+ ## as well would double-count that sample (verified against an
+ ## independent Astrom & Wittenmark split-ZOH ground truth: with the
+ ## register left in the dynamics, InputDelay must stay at d, not
+ ## d+1). The one case where the register's dynamics are explicitly
+ ## removed from the transfer function (the SISO tf/zpk pole/zero
+ ## cancellation below) re-adds the missing sample there, exactly
+ ## because the register's contribution is being deleted from the
+ ## dynamics at that point.
+ endif
+ endfor
+
+ n_extra = sum (extra_needed);
+ A_aug = zeros (n + n_extra, n + n_extra);
+ A_aug (1:n, 1:n) = Ad;
+ B_aug = zeros (n + n_extra, m);
+ C_aug = [C, zeros(rows(C), n_extra)];
+
+ row = n + 1;
+ for j = 1 : m
+ B_aug (1:n, j) = Bd_cols{j};
+ if (extra_needed(j))
+ A_aug (1:n, row) = g0_cols{j};
+ B_aug (row, j) = 1;
+ row += 1;
+ endif
+ endfor
+
+ sys = ss (A_aug, B_aug, C_aug, D, tsam);
+ sys = set (sys, "InputDelay", d_samples(:));
+
+ if (isa (origsys, "tf") || isa (origsys, "zpk"))
+ sys = tf (sys);
+
+ ## The extra register state added above (when extra.needed) has an
+ ## eigenvalue of exactly zero and contributes no real dynamics -- it
+ ## exists purely to carry the sub-sample remainder through one more
+ ## sample. Converting the augmented state-space realization to a
+ ## transfer function therefore leaves an exact pole/zero pair at the
+ ## origin (a leading zero numerator coefficient paired with a trailing
+ ## ~0 denominator coefficient): cancel it explicitly so num/den come
+ ## back at minimal (MATLAB-compatible) order instead of one degree too
+ ## high.
+ if (p == 1 && m == 1 && n_extra > 0)
+ [num, den] = tfdata (sys, "v");
+ if (abs (num(1)) < 1e-8 * max (abs (num)) ...
+ && abs (den(end)) < 1e-8 * max (abs (den)))
+ num = num(2:end);
+ den = den(1:end-1);
+ sys = tf (num, den, tsam);
+ ## The register's dynamics were just deleted from num/den above, so
+ ## (unlike the general case) the one sample of delay it used to
+ ## contribute via its own pole must now be re-added explicitly.
+ sys = set (sys, "InputDelay", d_samples + 1);
+ endif
+ endif
+
+ if (isa (origsys, "zpk"))
+ sys = zpk (sys);
+ endif
+ endif
+endfunction
diff --git a/inst/@lti/connect.m b/inst/@lti/connect.m
index 0550f212..cf8c6fd4 100644
--- a/inst/@lti/connect.m
+++ b/inst/@lti/connect.m
@@ -67,6 +67,11 @@
## inputs @var{in}.
## @end table
##
+## @strong{Remarks}
+## Input, output, and I/O delays on the operand systems are supported. If a delay
+## is trapped inside a closed loop, it will appear as an @code{InternalDelay}
+## property on the result.
+##
## @strong{Example}
##
## Consider the control loop with reference r, disturbances d1, d2
diff --git a/inst/@lti/ctranspose.m b/inst/@lti/ctranspose.m
index 30618a90..74225bdf 100644
--- a/inst/@lti/ctranspose.m
+++ b/inst/@lti/ctranspose.m
@@ -51,6 +51,15 @@
error ("lti: ctranspose: this is an unary operator");
endif
+ if (hasdelay (sys))
+ ## Pertransposition substitutes s -> -s (continuous) or z -> 1/z
+ ## (discrete), turning a delay e^(-tau s) into e^(+tau s) -- a pure
+ ## advance, which InputDelay/OutputDelay/IODelay (nonnegative-only)
+ ## cannot represent. Error clearly rather than silently return stale
+ ## pre-transpose delay data.
+ error ("lti: ctranspose: delay is not yet supported");
+ endif
+
[p, m] = size (sys);
ct = isct (sys);
@@ -62,3 +71,13 @@
sys.outgroup = struct ();
endfunction
+
+
+%!test # no delay: unaffected (regression)
+%! sys = tf (1, [1 1]);
+%! sys_ct = sys';
+%! assert (hasdelay (sys_ct), false);
+
+%!error tf (1, [1 1], "InputDelay", 0.5)'
+%!error ss (-1, 1, 1, 0, "OutputDelay", 0.5)'
+%!error zpk ([], -1, 1, "IODelay", 0.5)'
diff --git a/inst/@lti/d2c.m b/inst/@lti/d2c.m
index 84869611..a93ea8f8 100644
--- a/inst/@lti/d2c.m
+++ b/inst/@lti/d2c.m
@@ -73,9 +73,44 @@
error ("d2c: third argument is not a valid pre-warping frequency");
endif
- sys = __d2c__ (sys, sys.tsam, lower (method), w0);
+ origsys = sys;
+ tsam = sys.tsam;
+
+ if (hasinternaldelay (origsys))
+ ## Mirror c2d.m: discretize the extended ordinary system, then convert
+ ## samples back to seconds for internaldelay (samples are exact, so no
+ ## rounding needed here, unlike the tau/tsam rounding done in c2d.m).
+ ## __ss_ext_split__ copies InputDelay/OutputDelay/IODelay unchanged from
+ ## origsys (still in samples), so if origsys also carries ordinary I/O
+ ## delay, convert it to seconds here too, mirroring c2d.m.
+ [ext_sys, nu, ny] = __ss_ext_build__ (origsys);
+ ext_sys = __d2c__ (ext_sys, tsam, lower (method), w0);
+ sys = __ss_ext_split__ (origsys, ext_sys, nu, ny);
+ sys.tsam = 0;
+ tau = get (origsys, "internaldelay");
+ sys = set (sys, "internaldelay", tau * tsam);
+
+ if (hasdelay (origsys))
+ [indelay, outdelay, iodelay] = get (origsys, "inputdelay", "outputdelay", "iodelay");
+ sys = set (sys, "InputDelay", indelay * tsam, ...
+ "OutputDelay", outdelay * tsam, ...
+ "IODelay", iodelay * tsam);
+ endif
+
+ return;
+ endif
+
+ sys = __d2c__ (sys, tsam, lower (method), w0);
sys.tsam = 0;
+ if (hasdelay (origsys))
+ [indelay, outdelay, iodelay] = get (origsys, "inputdelay", "outputdelay", "iodelay");
+
+ sys = set (sys, "InputDelay", indelay * tsam, ...
+ "OutputDelay", outdelay * tsam, ...
+ "IODelay", iodelay * tsam);
+ endif
+
endfunction
@@ -146,3 +181,75 @@
%! Me = [A, B; C, D];
%!
%!assert (Mo, Me, 1e-4);
+
+
+%!test # integer sample delay converts to seconds correctly
+%! sys = tf (1, [1 -0.5], 0.5, "InputDelay", 2);
+%! csys = d2c (sys, "zoh");
+%! assert (isct (csys), true);
+%! assert (csys.InputDelay, 1.0, 1e-10);
+%! assert (csys.OutputDelay, 0, 1e-10);
+%! assert (csys.IODelay, 0, 1e-10);
+
+%!test # no delay: unaffected (regression)
+%! sys = tf (1, [1 -0.5], 0.5);
+%! csys = d2c (sys, "zoh");
+%! assert (hasdelay (csys), false);
+
+%!test # round-trip: delay survives c2d then d2c
+%! sys = tf (1, [1 1], "InputDelay", 1.0);
+%! sys2 = d2c (c2d (sys, 0.5, "zoh"), "zoh");
+%! assert (sys2.InputDelay, 1.0, 1e-10);
+
+%!test # InternalDelay round-trip: c2d then d2c recovers the original tau (samples convert back to seconds exactly)
+%! T = 0.3;
+%! tsam = 0.1;
+%! G = ss (-2, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! dsys = c2d (L, tsam, "zoh");
+%! csys = d2c (dsys, "zoh");
+%! assert (hasinternaldelay (csys), true);
+%! assert (csys.internaldelay, T, 1e-10);
+%! w = [0.05, 0.2, 1];
+%! assert (freqresp (csys, w), freqresp (L, w), 5e-2);
+
+%!test # InternalDelay with no delay ports: degenerate extended system still converts back
+%! sys = set (ss (-1, 1, 1, 0, 0.1), "internaldelay", 0.5);
+%! csys = d2c (sys, "zoh");
+%! assert (isct (csys), true);
+%! assert (csys.internaldelay, 0.05, 1e-10);
+
+%!test # combined InternalDelay AND ordinary InputDelay: both converted to seconds
+%! T = 0.5;
+%! tsam = 0.25;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! L = set (L, "InputDelay", 0.75);
+%! dsys = c2d (L, tsam, "zoh");
+%! csys = d2c (dsys, "zoh");
+%! assert (hasinternaldelay (csys), true);
+%! assert (csys.internaldelay, T, 1e-10);
+%! assert (csys.InputDelay, round (0.75 / tsam) * tsam, 1e-10);
+
+%!test # MIMO InternalDelay: two-channel round-trip c2d->d2c recovers both distinct taus
+%! ## Same isolated-SISO-block MIMO fixture as the c2d MIMO test; verify the
+%! ## extended-system split is dimension-correct for more than one delay port
+%! ## in BOTH directions (samples -> seconds per channel).
+%! T1 = 0.4; T2 = 0.8; tsam = 0.2;
+%! G1 = ss (-1, 1, 1, 0, "IODelay", T1);
+%! G2 = ss (-2, 1, 1, 0, "IODelay", T2);
+%! sys = append (feedback (G1), feedback (G2));
+%! dsys = c2d (sys, tsam, "zoh");
+%! csys = d2c (dsys, "zoh");
+%! assert (hasinternaldelay (csys), true);
+%! assert (isct (csys), true);
+%! assert (get (csys, "internaldelay"), [T1; T2], 1e-10); # both recovered
+%! ## freqresp of the round-tripped continuous system matches the original,
+%! ## per diagonal channel; off-diagonals stay zero.
+%! w = [0.05, 0.3, 1.2];
+%! H = freqresp (csys, w); H0 = freqresp (sys, w);
+%! assert (H, H0, 1e-8);
+%! for k = 1:numel (w)
+%! assert (H(1,2,k), 0, 1e-12);
+%! assert (H(2,1,k), 0, 1e-12);
+%! endfor
diff --git a/inst/@lti/feedback.m b/inst/@lti/feedback.m
index a6932253..ec808f04 100644
--- a/inst/@lti/feedback.m
+++ b/inst/@lti/feedback.m
@@ -51,6 +51,11 @@
## Resulting @acronym{LTI} model.
## @end table
##
+## @strong{Remarks}
+## Both operands may contain input, output, and I/O delays. If a delay is trapped
+## inside the closed loop, it will appear as an @code{InternalDelay} property on the
+## (always-@code{ss}) result.
+##
## @strong{Block Diagram}
## @example
## @group
@@ -178,6 +183,16 @@
in_idx = 1 : m1;
out_idx = 1 : p1;
+ ## A delay trapped inside the loop cannot be represented by tf/zpk, which
+ ## carry no InternalDelay. Force both operands to state-space so the delay
+ ## survives as an internal delay port (matches Matlab: such a result is ss).
+ d1 = isa (sys1, "lti") && (hasdelay (sys1) || hasinternaldelay (sys1));
+ d2 = isa (sys2, "lti") && (hasdelay (sys2) || hasinternaldelay (sys2));
+ if (d1 || d2)
+ sys1 = ss (sys1);
+ sys2 = ss (sys2);
+ endif
+
sys = __sys_group__ (sys1, sys2);
sys = __sys_connect__ (sys, M);
sys = __sys_prune__ (sys, out_idx, in_idx);
@@ -276,3 +291,76 @@
%!assert (S1.b, S2.b, 1e-4);
%!assert (S1.c, S2.c, 1e-4);
%!assert (S1.d, S2.d, 1e-4);
+
+
+%!test # no delay on either operand: existing behavior unaffected (regression)
+%! s1 = tf (1, [1 1]);
+%! s = feedback (s1);
+%! assert (hasdelay (s), false);
+
+## Internal delay: an input delay closing a feedback loop is no longer an
+## error; it is absorbed into the closed loop's InternalDelay (the delay
+## e^{-0.1 s} lives inside the loop 1/(s+1+e^{-0.1 s}) and cannot be
+## expressed as a pure I/O time-shift).
+%!test
+%! s = feedback (tf (1, [1 1], "InputDelay", 0.1));
+%! assert (isa (s, "ss")); # forced to ss to carry delay
+%! assert (get (s, "internaldelay"), 0.1, 1e-10);
+%! assert (hasinternaldelay (s), true);
+
+## Scalar loop with an IODelay in the forward path. A single scalar delay
+## closing a loop is irreducible, so the closed-loop InternalDelay is just T.
+## Hand derivation: G = e^{-sT} g(s), C = c(s), L = G C.
+## Closed loop feedback(G,C) = G/(1+GC); the e^{-sT} sits inside the loop,
+## so InternalDelay = T (0.3 here) and no residual I/O delay remains.
+%!test
+%! T = 0.3;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! C = ss (-2, 1, 1, 0);
+%! cl = feedback (G, C);
+%! assert (get (cl, "internaldelay"), T, 1e-10);
+%! assert (get (cl, "iodelay"), 0, 1e-10); # absorbed, no residual
+
+## feedback() and connect() must agree on InternalDelay for equivalent
+## topologies, since they share __sys_connect__.
+%!test
+%! T = 0.3;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! C = ss (-2, 1, 1, 0);
+%! cl_fb = feedback (G, C);
+%! ## build the same negative feedback loop with connect():
+%! ## inputs : [r ; u_G ; u_C] = grouped [G.in ; C.in] plus external r
+%! ## here use index-based cm on append (G, C):
+%! ## in 1 -> G, in 2 -> C ; out 1 -> G, out 2 -> C
+%! Gall = append (G, C);
+%! ## u_G = r - y_C (row: input 1 fed by -output 2)
+%! ## u_C = y_G (row: input 2 fed by +output 1)
+%! cm = [1, -2; 2, 1];
+%! cl_cn = connect (Gall, cm, 1, 1);
+%! assert (get (cl_cn, "internaldelay"), get (cl_fb, "internaldelay"), 1e-10);
+
+## MIMO plant whose IODelay column has more than one differently-valued
+## nonzero entry (one input feeding two differently-delayed outputs). A single
+## shared actuator column cannot carry two different delays, so the column is
+## decomposed onto an independent shadow-state block (one per decomposed
+## column) and each output taps it with its own delay. This used to error
+## ("not yet supported"); it is now realised exactly. Coupled A so both
+## outputs genuinely depend on the shared input.
+%!test
+%! a = [-1, 0.4; 0, -2];
+%! b = [1; 1]; # single input drives both states
+%! c = eye (2);
+%! d = [0; 0];
+%! iod = [0.2; 0.5]; # one input, two differently-delayed outputs
+%! G = ss (a, b, c, d, "IODelay", iod);
+%! cl = feedback (G, [1, 1]); # close the single loop (sum of outputs)
+%! assert (rows (ssdata (cl)), 4); # 2 original + 2 shadow (one decomposed col)
+%! assert (sort (get (cl, "internaldelay")), [0.2; 0.5], 1e-12);
+%! w = [0.5, 1.3, 3.0];
+%! for k = 1 : numel (w)
+%! s = 1i*w(k);
+%! T = c/(s*eye(2) - a)*b + d; # 2x1 open-loop transfer
+%! Gt = T .* exp (-1i*w(k)*iod); # per-entry delayed
+%! Href = Gt / (1 + [1, 1]*Gt); # feedback (G, [1 1]): e = u - [1 1]y
+%! assert (freqresp (cl, w(k)), Href, 1e-9);
+%! endfor
diff --git a/inst/@lti/freqresp.m b/inst/@lti/freqresp.m
index cc87ba3f..77f3d0d0 100644
--- a/inst/@lti/freqresp.m
+++ b/inst/@lti/freqresp.m
@@ -53,5 +53,44 @@
endif
H = __freqresp__ (sys, w);
+ H = __apply_freqresp_delay__ (sys, H, w, false);
endfunction
+
+%!test # no-delay system: unaffected (regression)
+%! sys = tf (1, [1, 1]);
+%! w = [0.1, 1, 5];
+%! H = freqresp (sys, w);
+%! rational_H = __freqresp__ (sys, w);
+%! assert (H, rational_H);
+
+%!test # SISO continuous tf with InputDelay: freqresp applies the delay
+%! sys = tf (1, [1, 1], "InputDelay", 0.3);
+%! rational = tf (1, [1, 1]);
+%! w = [0.1, 1, 5];
+%! H = freqresp (sys, w);
+%! expected = freqresp (rational, w) .* reshape (exp (-1i * w * 0.3), 1, 1, []);
+%! assert (H, expected, 1e-10);
+
+%!test # SISO discrete tf with integer InputDelay: freqresp applies the delay
+%! sys = tf (1, [1, -0.5], 0.1, "InputDelay", 2);
+%! rational = tf (1, [1, -0.5], 0.1);
+%! w = [0.1, 1, 5];
+%! H = freqresp (sys, w);
+%! expected = freqresp (rational, w) .* reshape (exp (-1i * w * 0.1 * 2), 1, 1, []);
+%! assert (H, expected, 1e-10);
+
+%!test # MIMO tf with per-channel IODelay: freqresp applies each channel's own delay
+%! sys = tf ({1, 1; 1, 1}, {[1, -0.5], [1, -0.6]; [1, -0.7], [1, -0.8]}, 0.1);
+%! sys = set (sys, "IODelay", [1, 0; 0, 2]);
+%! rational = tf ({1, 1; 1, 1}, {[1, -0.5], [1, -0.6]; [1, -0.7], [1, -0.8]}, 0.1);
+%! total = [1, 0; 0, 2];
+%! w = [0.1, 1, 5];
+%! resp = freqresp (rational, w);
+%! expected = zeros (2, 2, numel (w));
+%! for i = 1:2
+%! for j = 1:2
+%! expected(i,j,:) = reshape (resp(i,j,:), 1, []) .* exp (-1i * w * 0.1 * total(i,j));
+%! endfor
+%! endfor
+%! assert (freqresp (sys, w), expected, 1e-8);
diff --git a/inst/@lti/get.m b/inst/@lti/get.m
index ecf7eb97..4919e122 100644
--- a/inst/@lti/get.m
+++ b/inst/@lti/get.m
@@ -56,6 +56,12 @@
val = sys.outgroup;
case "tsam"
val = sys.tsam;
+ case "inputdelay"
+ val = sys.indelay;
+ case "outputdelay"
+ val = sys.outdelay;
+ case "iodelay"
+ val = sys.iodelay;
case "name"
val = sys.name;
case "notes"
diff --git a/inst/@lti/lti.m b/inst/@lti/lti.m
index 77133440..1025429a 100644
--- a/inst/@lti/lti.m
+++ b/inst/@lti/lti.m
@@ -38,7 +38,10 @@
"outgroup", struct (),
"name", "",
"notes", {{}},
- "userdata", []);
+ "userdata", [],
+ "indelay", zeros (m, 1),
+ "outdelay", zeros (p, 1),
+ "iodelay", zeros (p, m));
ltisys = class (ltisys, "lti");
diff --git a/inst/@lti/pade.m b/inst/@lti/pade.m
new file mode 100644
index 00000000..4e7b8c9a
--- /dev/null
+++ b/inst/@lti/pade.m
@@ -0,0 +1,578 @@
+## Copyright (C) 2026 Prateek Ganguli
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope. If not, see .
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{sys} =} pade (@var{sys}, @var{n})
+## Replace every delay (@var{InputDelay}, @var{OutputDelay}, @var{IODelay})
+## of the @acronym{LTI} model @var{sys} by an @var{n}th-order Pade
+## approximation, returning a delay-free equivalent model.
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## @var{tf}, @var{zpk}, or @var{ss} model. @var{ss} models carrying a
+## nonzero @var{InternalDelay} are supported: the delay loop is closed
+## through a rational Pade filter via a linear-fractional transformation
+## (@code{lft}).
+## @item n
+## Order of the Pade approximation. Either a scalar (applied to every
+## nonzero delay) or a vector with one entry per nonzero delay -- see
+## @code{__pade_order_vector__} for the exact ordering convention.
+## @end table
+##
+## @strong{Outputs}
+## @table @var
+## @item sys
+## Delay-free equivalent model with the Pade-approximated rational factors
+## absorbed into each entry's numerator/denominator (or zeros/poles, or
+## state-space realization). For @var{ss} inputs, the result is obtained
+## via a @code{tf} roundtrip and does @strong{not} preserve the original
+## state basis or state count (matching @acronym{MATLAB}'s own
+## @code{pade()}, which likewise only guarantees an equivalent model, not
+## a minimal one -- use @code{minreal} afterward if minimality is needed).
+## @end table
+## @end deftypefn
+
+## Author: Prateek Ganguli
+## Created: July 2026
+## Version: 0.1
+
+function sys = pade (sys, n)
+
+ if (nargin != 2 || ! isa (sys, "lti"))
+ print_usage ();
+ endif
+
+ if (isa (sys, "ss"))
+ if (hasinternaldelay (sys))
+ sys = __pade_substitute_internal__ (sys, n);
+ return;
+ endif
+ if (! hasdelay (sys))
+ return;
+ endif
+ sys_tf = __pade_substitute_ordinary__ (tf (sys), n);
+ sys = ss (sys_tf);
+ return;
+ endif
+
+ if (! (isa (sys, "tf") || isa (sys, "zpk")))
+ error ("pade: only tf, zpk, and ss models are supported");
+ endif
+
+ if (! hasdelay (sys))
+ return;
+ endif
+
+ sys = __pade_substitute_ordinary__ (sys, n);
+
+endfunction
+
+## Shared ordinary-delay (InputDelay/OutputDelay/IODelay) per-entry Pade
+## substitution logic for tf/zpk models. Assumes hasdelay(sys) is true.
+## Used directly for tf/zpk inputs, and via a tf-roundtrip for ss inputs
+## with no InternalDelay.
+function sys = __pade_substitute_ordinary__ (sys, n)
+
+ [indelay, outdelay, iodelay] = get (sys, "inputdelay", "outputdelay", "iodelay");
+ total = totaldelay (sys);
+ [pr, pc] = size (total);
+
+ if (isdt (sys))
+ ## n still validated (same path as continuous input) even though it has
+ ## no effect on the exact discrete result -- discarding the return
+ ## values is intentional. __pade_order_vector__ only checks that a
+ ## vector n has the right length; the actual per-entry order values
+ ## (integer, non-negative, causal) are validated by the low-level
+ ## pade() polynomial routine in the continuous branch below, so call
+ ## it here too (with a dummy dead time) to match that validation
+ ## exactly.
+ [orders, ~] = __pade_order_vector__ (sys, n);
+ for k = 1 : numel (orders)
+ ## nargout=2 (matching the continuous path's own call shape below)
+ ## is required here, not just cosmetic: pade() with nargout==0 runs
+ ## its plotting branch, which would pop open stray figure windows as
+ ## a side effect of what must be a pure validation call.
+ [~, ~] = pade (1, orders (k));
+ endfor
+ sys = set (sys, "InputDelay", indelay, "OutputDelay", outdelay, "IODelay", iodelay);
+ sys = absorbDelay (sys);
+ return;
+ endif
+
+ [orders, idx] = __pade_order_vector__ (sys, n);
+
+ n_in = numel (idx.input);
+ n_out = numel (idx.output);
+ n_iod = numel (idx.iod);
+
+ orders_input = orders (1 : n_in);
+ orders_output = orders (n_in + 1 : n_in + n_out);
+ orders_iod = orders (n_in + n_out + 1 : n_in + n_out + n_iod);
+
+ ## Build a per-entry order matrix. Priority (lowest to highest, later
+ ## assignments win on conflict): OutputDelay broadcast across a row,
+ ## InputDelay broadcast across a column, IODelay assigned directly to
+ ## its own (i,j) entry -- IODelay is the most entry-specific, so it
+ ## takes precedence if more than one delay kind is nonzero for the same
+ ## entry. Test fixtures in this task only ever set one delay kind at a
+ ## time per entry, so no conflict arises in practice.
+ order_matrix = nan (pr, pc);
+
+ for k = 1 : n_out
+ i = idx.output (k);
+ order_matrix (i, :) = orders_output (k);
+ endfor
+
+ for k = 1 : n_in
+ j = idx.input (k);
+ order_matrix (:, j) = orders_input (k);
+ endfor
+
+ for k = 1 : n_iod
+ [i, j] = ind2sub ([pr, pc], idx.iod (k));
+ order_matrix (i, j) = orders_iod (k);
+ endfor
+
+ origsys = sys;
+
+ if (isa (sys, "zpk"))
+
+ [z, p, k, tsam] = zpkdata (sys);
+
+ for i = 1 : pr
+ for j = 1 : pc
+ d = total (i, j);
+ if (d > 0)
+ order = order_matrix (i, j);
+ [num_p, den_p] = pade (d, order);
+ gain_factor = num_p(1) / den_p(1);
+ z{i,j} = [z{i,j}; roots(num_p)];
+ p{i,j} = [p{i,j}; roots(den_p)];
+ k(i,j) = k(i,j) * gain_factor;
+ endif
+ endfor
+ endfor
+
+ sys = zpk (z, p, k, tsam);
+
+ else
+
+ [num, den, tsam] = tfdata (sys);
+
+ for i = 1 : pr
+ for j = 1 : pc
+ d = total (i, j);
+ if (d > 0)
+ order = order_matrix (i, j);
+ [num_p, den_p] = pade (d, order);
+ num{i,j} = conv (num{i,j}, num_p);
+ den{i,j} = conv (den{i,j}, den_p);
+ endif
+ endfor
+ endfor
+
+ sys = tf (num, den, tsam);
+
+ endif
+
+ sys = set (sys, "lti", origsys);
+ sys = set (sys, "InputDelay", 0, "OutputDelay", 0, "IODelay", 0);
+
+endfunction
+
+## InternalDelay Pade substitution for ss models. Each internal delay
+## port is a signal z(t) fed back to the plant as w(t) = z(t-tau); the
+## exact freqresp closes this loop with Delta = diag(exp(-jw*tau)) via an
+## LFT (see @ss/__freqresp__.m). Here we replace that ideal-delay loop
+## closure by a rational one: build the extended ordinary plant (whose last
+## nports inputs are w and last nports outputs are z) with
+## __ss_ext_build__, then close the loop through a diagonal MIMO Pade
+## filter G_pade (one SISO Pade block per port) using the Redheffer star
+## product lft(). As the order grows, each SISO Pade block converges to
+## exp(-jw*tau), so the closed loop converges to the exact-delay system.
+##
+## Any ordinary delay (InputDelay/OutputDelay/IODelay) carried alongside
+## the InternalDelay is independent of the loop closure, so it is handled
+## afterwards by the same tf-roundtrip substitution used for the
+## InternalDelay-free ss case (__pade_substitute_ordinary__), fed the
+## ordinary slice of the SAME order assignment computed once up front.
+function sys = __pade_substitute_internal__ (sys, n)
+
+ ## Full per-delay order assignment, computed ONCE so the internal and
+ ## ordinary substitutions consume consistent slices (ordering per
+ ## __pade_order_vector__: InputDelay, OutputDelay, IODelay, InternalDelay).
+ [orders, idx] = __pade_order_vector__ (sys, n);
+ n_in = numel (idx.input);
+ n_out = numel (idx.output);
+ n_iod = numel (idx.iod);
+ n_int = numel (idx.internal);
+
+ orders_ordinary = orders (1 : n_in + n_out + n_iod);
+ orders_internal = orders (n_in + n_out + n_iod + 1 : end);
+
+ ## Original internal-delay values and any ordinary delay, extracted BEFORE
+ ## the extended plant is rebuilt (the rebuild deliberately drops both the
+ ## delay ports and the ordinary delay).
+ tau = get (sys, "internaldelay");
+ tau = tau(:);
+ had_ordinary = hasdelay (sys);
+ [oindelay, ooutdelay, oiodelay] = get (sys, "inputdelay", "outputdelay", "iodelay");
+
+ ## Build the extended ordinary plant. __ss_ext_build__ folds the delay
+ ## ports (b2/c2/d12/d21/d22) into ordinary extra inputs/outputs but leaves
+ ## the returned ext_sys's lti input/output NAMES (hence its reported
+ ## size()) and its residual b2/c2/.../internaldelay fields untouched, so
+ ## it is NOT a well-formed ordinary ss. Rebuild a clean ss from just the
+ ## extended (a,b,c,d[,e]) matrices: that gives correct input/output counts
+ ## for lft()'s size() calls AND clears every residual port field and the
+ ## internaldelay, so hasinternaldelay(ext_clean) is false.
+ ## __ss_ext_build__ lives in @ss and returns an ss whose a/b/c/d hold the
+ ## extended matrices; extract them with dssdata (the raw fields are not
+ ## dot-accessible from outside @ss). Passing [] keeps a descriptor E empty
+ ## rather than expanding it to the identity.
+ [ext_sys, nu, ny] = __ss_ext_build__ (sys);
+ [ea, eb, ec, ed, ee, etsam] = dssdata (ext_sys, []);
+ nports = columns (eb) - nu;
+ if (nports != rows (ec) - ny)
+ error ("pade: internal delay port dimension mismatch (%d inputs vs %d outputs)",
+ nports, rows (ec) - ny);
+ endif
+
+ ## Rebuild a clean ordinary ss from just the extended matrices. This gives
+ ## correct input/output counts for lft()'s size() calls (the raw ext_sys
+ ## keeps the ORIGINAL input/output names, so its reported size is stale)
+ ## AND clears every residual port field and the internaldelay, so
+ ## hasinternaldelay(ext_clean) is false.
+ if (isempty (ee))
+ ext_clean = ss (ea, eb, ec, ed, etsam);
+ else
+ ext_clean = dss (ea, eb, ec, ed, ee, etsam);
+ endif
+
+ if (nports == 0)
+ ## Degenerate: an internaldelay is set but there are no actual delay
+ ## ports (empty b2/c2/...), so the "delay" affects no signal. There is
+ ## no loop to close -- the delay-free equivalent is just the clean plant.
+ sys = ext_clean;
+ else
+ if (nports != numel (tau) || nports != n_int)
+ error ("pade: internal delay port count (%d) does not match delay count (%d/%d)",
+ nports, numel (tau), n_int);
+ endif
+
+ ## Diagonal MIMO Pade filter: nports inputs (each port's z, the delayed
+ ## signal's source) and nports outputs (each port's w, the destination).
+ ## For discrete sys, tau(k) is already a nonnegative-integer sample
+ ## count, so the loop is closed EXACTLY via __exact_discrete_delay_ss__
+ ## instead of a rational Pade approximation.
+ if (isdt (sys))
+ G_pade = __exact_discrete_delay_ss__ (tau(1), etsam);
+ for k = 2 : nports
+ G_pade = append (G_pade, __exact_discrete_delay_ss__ (tau(k), etsam));
+ endfor
+ else
+ G_pade = pade (tau(1), orders_internal(1));
+ for k = 2 : nports
+ G_pade = append (G_pade, pade (tau(k), orders_internal(k)));
+ endfor
+ endif
+
+ ## Close the loop. lft(sys1, sys2, nu, ny) connects the LAST nu inputs of
+ ## sys1 to the FIRST nu outputs of sys2, and the LAST ny outputs of sys1
+ ## to the FIRST ny inputs of sys2. With sys1 = ext_clean, sys2 = G_pade,
+ ## nu = ny = nports: ext_clean's last nports inputs (w) are driven by
+ ## G_pade's first nports outputs (w), and ext_clean's last nports outputs
+ ## (z) drive G_pade's first nports inputs (z) -- exactly the loop the ideal
+ ## delay used to close (w = Delta z, now w = G_pade z).
+ sys = lft (ext_clean, G_pade, nports, nports);
+ endif
+
+ ## Ordinary delay, if any, is untouched by the loop closure; reattach the
+ ## ORIGINAL ordinary delay to the delay-free closure result and
+ ## Pade-substitute it via the same tf-roundtrip path used for
+ ## InternalDelay-free ss inputs, using the ordinary slice of the order
+ ## assignment computed above. lft() preserves the ordinary I/O dimensions
+ ## (the first nu inputs / first ny outputs), so the ordinary delay vectors
+ ## still line up.
+ if (had_ordinary)
+ sys = set (sys, "InputDelay", oindelay, ...
+ "OutputDelay", ooutdelay, ...
+ "IODelay", oiodelay);
+ sys = ss (__pade_substitute_ordinary__ (tf (sys), orders_ordinary));
+ endif
+
+endfunction
+
+
+%!test # SISO tf with InputDelay: freqresp matches hand-multiplied reference
+%! T = 0.5;
+%! sys = tf (1, [1 2], "InputDelay", T);
+%! sys2 = pade (sys, 3);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! [num_p, den_p] = pade (T, 3);
+%! expected = freqresp (tf (1, [1 2]), w) .* freqresp (tf (num_p, den_p), w);
+%! assert (freqresp (sys2, w), expected, 1e-8);
+
+%!test # SISO tf with InputDelay, a different order
+%! T = 0.5;
+%! sys = tf (1, [1 2], "InputDelay", T);
+%! sys2 = pade (sys, 5);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! [num_p, den_p] = pade (T, 5);
+%! expected = freqresp (tf (1, [1 2]), w) .* freqresp (tf (num_p, den_p), w);
+%! assert (freqresp (sys2, w), expected, 1e-8);
+
+%!test # dense MIMO tf with differing per-channel IODelay
+%! sys = tf ({1, 1; 1, 1}, {[1 2], [1 3]; [1 4], [1 5]});
+%! sys = set (sys, "IODelay", [0.2, 0.3; 0.4, 0.5]);
+%! sys2 = pade (sys, 3);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! rational = tf ({1, 1; 1, 1}, {[1 2], [1 3]; [1 4], [1 5]});
+%! total = [0.2, 0.3; 0.4, 0.5];
+%! resp = freqresp (rational, w);
+%! expected = zeros (2, 2, numel (w));
+%! for i = 1:2
+%! for j = 1:2
+%! [num_p, den_p] = pade (total(i,j), 3);
+%! expected(i,j,:) = reshape (resp(i,j,:), 1, []) .* reshape (freqresp (tf (num_p, den_p), w), 1, []);
+%! endfor
+%! endfor
+%! assert (freqresp (sys2, w), expected, 1e-6);
+
+%!test # equivalent zpk case, SISO with OutputDelay
+%! T = 0.3;
+%! sys = zpk ([], -2, 3, "OutputDelay", T);
+%! sys2 = pade (sys, 4);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! [num_p, den_p] = pade (T, 4);
+%! expected = freqresp (zpk ([], -2, 3), w) .* freqresp (tf (num_p, den_p), w);
+%! assert (freqresp (sys2, w), expected, 1e-8);
+
+%!test # equivalent zpk case, dense MIMO with IODelay
+%! z = {[], []; [], []};
+%! p = {-2, -3; -4, -5};
+%! k = [1, 1; 1, 1];
+%! sys = zpk (z, p, k);
+%! sys = set (sys, "IODelay", [0.2, 0.3; 0.4, 0.5]);
+%! sys2 = pade (sys, 3);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! rational = zpk (z, p, k);
+%! total = [0.2, 0.3; 0.4, 0.5];
+%! resp = freqresp (rational, w);
+%! expected = zeros (2, 2, numel (w));
+%! for i = 1:2
+%! for j = 1:2
+%! [num_p, den_p] = pade (total(i,j), 3);
+%! expected(i,j,:) = reshape (resp(i,j,:), 1, []) .* reshape (freqresp (tf (num_p, den_p), w), 1, []);
+%! endfor
+%! endfor
+%! assert (freqresp (sys2, w), expected, 1e-6);
+
+%!test # per-delay order vector: each entry keeps its own assigned order
+%! # 1x2 tf, two nonzero IODelay entries with distinct Pade orders 2 and 5.
+%! # Each substituted entry's denominator degree must equal its own
+%! # rational-part degree (1, since den = [1 2] or [1 3]) plus its own
+%! # assigned Pade order -- checked individually, not in aggregate.
+%! sys = tf ({1, 1}, {[1 2], [1 3]});
+%! sys = set (sys, "IODelay", [0.2, 0.4]);
+%! sys2 = pade (sys, [2, 5]);
+%! [num, den] = tfdata (sys2);
+%! assert (length (den{1,1}) - 1, 1 + 2);
+%! assert (length (den{1,2}) - 1, 1 + 5);
+
+%!error pade (tf (1, [1 2], "InputDelay", 0.5), [1, 2])
+
+%!test # ss with dense per-entry IODelay, cross-checked against the tf path
+%! a = [-2, 0; 0, -3];
+%! b = [1, 0; 0, 1];
+%! c = [1, 1; 1, 1];
+%! d = [0, 0; 0, 0];
+%! sys_ss = ss (a, b, c, d);
+%! sys_ss = set (sys_ss, "IODelay", [0.2, 0.3; 0.4, 0.5]);
+%! sys_tf = tf ({1, 1; 1, 1}, {[1 2], [1 3]; [1 2], [1 3]});
+%! sys_tf = set (sys_tf, "IODelay", [0.2, 0.3; 0.4, 0.5]);
+%! ss2 = pade (sys_ss, 3);
+%! tf2 = pade (sys_tf, 3);
+%! assert (hasdelay (ss2), false);
+%! assert (hasdelay (tf2), false);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (ss2, w), freqresp (tf2, w), 1e-6);
+
+%!test # ss InternalDelay: feedback loop, freqresp matches Pade LFT closed form
+%! ## L = feedback (G), G = 1/(s+1) with loop delay T = 0.3 s, traps the delay
+%! ## internally. Exact closed loop: Gd/(1+Gd), Gd = G e^{-jwT}. The Pade
+%! ## approximation replaces e^{-jwT} by the scalar Pade transfer function's
+%! ## own freqresp Pd, so the reference is Gp/(1+Gp) with Gp = G*Pd -- an
+%! ## INDEPENDENT computation (scalar pade path), not a second call into the
+%! ## code under test.
+%! T = 0.3;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! assert (hasinternaldelay (L), true);
+%! order = 4;
+%! sys2 = pade (L, order);
+%! assert (hasinternaldelay (sys2), false);
+%! assert (hasdelay (sys2), false);
+%! w = [0.05, 0.5, 1.7, 4, 11];
+%! H = reshape (freqresp (sys2, w), 1, []);
+%! [np, dp] = pade (T, order);
+%! Pd = reshape (freqresp (tf (np, dp), w), 1, []);
+%! Gp = (1 ./ (1i*w + 1)) .* Pd;
+%! expected = Gp ./ (1 + Gp);
+%! assert (H, expected, 1e-9);
+
+%!test # ss InternalDelay: freqresp error against TRUE exact delay shrinks with order
+%! ## Strongest evidence the substitution converges to the exact-delay system:
+%! ## as the Pade order grows, the approximation error against the untouched
+%! ## exact-delay freqresp must strictly decrease.
+%! T = 0.3;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = feedback (G);
+%! w = [0.05, 0.5, 1.7, 4, 11];
+%! Hexact = reshape (freqresp (L, w), 1, []);
+%! err = zeros (1, 3);
+%! orders = [2, 4, 8];
+%! for k = 1:3
+%! sysk = pade (L, orders(k));
+%! Hk = reshape (freqresp (sysk, w), 1, []);
+%! err(k) = max (abs (Hk - Hexact));
+%! endfor
+%! assert (err(2) < err(1));
+%! assert (err(3) < err(2));
+
+%!test # ss MIMO InternalDelay: two independent ports, distinct delays and orders
+%! ## append() of two SISO feedback loops -> genuine 2-port InternalDelay.
+%! ## A per-port order VECTOR (3 for port 1, 6 for port 2) verifies each port
+%! ## consumes its own assigned order. Reference: decoupled per-channel Pade
+%! ## closed form; off-diagonals must be identically zero.
+%! T1 = 0.3; a1 = 1; o1 = 3;
+%! T2 = 0.7; a2 = 2; o2 = 6;
+%! G1 = ss (-a1, 1, 1, 0, "IODelay", T1);
+%! G2 = ss (-a2, 1, 1, 0, "IODelay", T2);
+%! sys = append (feedback (G1), feedback (G2));
+%! assert (get (sys, "internaldelay"), [T1; T2], 1e-12);
+%! sys2 = pade (sys, [o1, o2]);
+%! assert (hasinternaldelay (sys2), false);
+%! assert (hasdelay (sys2), false);
+%! w = [0.05, 0.5, 1.7, 4, 11];
+%! [n1, d1] = pade (T1, o1); P1 = reshape (freqresp (tf (n1, d1), w), 1, []);
+%! [n2, d2] = pade (T2, o2); P2 = reshape (freqresp (tf (n2, d2), w), 1, []);
+%! Gp1 = (1 ./ (1i*w + a1)) .* P1; e1 = Gp1 ./ (1 + Gp1);
+%! Gp2 = (1 ./ (1i*w + a2)) .* P2; e2 = Gp2 ./ (1 + Gp2);
+%! H = freqresp (sys2, w);
+%! for k = 1:numel (w)
+%! assert (H(1,1,k), e1(k), 1e-9);
+%! assert (H(2,2,k), e2(k), 1e-9);
+%! assert (H(1,2,k), 0, 1e-9);
+%! assert (H(2,1,k), 0, 1e-9);
+%! endfor
+
+%!test # ss combined InternalDelay AND ordinary OutputDelay: both approximated
+%! ## The loop-trapped internal delay and an untouched OutputDelay are
+%! ## independent; the result must equal the internal-delay Pade closure
+%! ## multiplied by the OutputDelay's own scalar Pade transfer function.
+%! T = 0.3; Tout = 0.2; order = 5;
+%! G = ss (-1, 1, 1, 0, "IODelay", T);
+%! L = set (feedback (G), "OutputDelay", Tout);
+%! assert (hasinternaldelay (L), true);
+%! assert (hasdelay (L), true);
+%! sys2 = pade (L, order);
+%! assert (hasinternaldelay (sys2), false);
+%! assert (hasdelay (sys2), false);
+%! w = [0.05, 0.5, 1.7, 4];
+%! [np, dp] = pade (T, order);
+%! Pd = reshape (freqresp (tf (np, dp), w), 1, []);
+%! Gp = (1 ./ (1i*w + 1)) .* Pd; internal_ref = Gp ./ (1 + Gp);
+%! [no, do_] = pade (Tout, order);
+%! Pout = reshape (freqresp (tf (no, do_), w), 1, []);
+%! expected = internal_ref .* Pout;
+%! H = reshape (freqresp (sys2, w), 1, []);
+%! assert (H, expected, 1e-8);
+
+%!test # ss InternalDelay with no delay ports (degenerate b2/c2 empty): still closes
+%! ## set() an internaldelay directly with empty ports -- the extended plant is
+%! ## degenerate (no extra I/O) but the substitution must not error and must
+%! ## drop the internal delay.
+%! sys = set (ss (-1, 1, 1, 0), "internaldelay", 0.5);
+%! sys2 = pade (sys, 3);
+%! assert (hasinternaldelay (sys2), false);
+
+%!test # discrete tf ordinary delay: exact absorption, matches absorbDelay directly
+%! sys = tf (1, [1, -0.5], 0.1, "InputDelay", 3);
+%! sys2 = pade (sys, 4);
+%! expected = absorbDelay (sys);
+%! assert (hasdelay (sys2), false);
+%! [num2, den2] = tfdata (sys2);
+%! [nume, dene] = tfdata (expected);
+%! assert (num2, nume, 1e-12);
+%! assert (den2, dene, 1e-12);
+
+%!test # discrete tf: pade order n has NO effect on the exact discrete result
+%! sys = tf (1, [1, -0.5], 0.1, "InputDelay", 2);
+%! sys_a = pade (sys, 2);
+%! sys_b = pade (sys, 9);
+%! [numa, dena] = tfdata (sys_a);
+%! [numb, denb] = tfdata (sys_b);
+%! assert (numa, numb, 1e-14);
+%! assert (dena, denb, 1e-14);
+
+%!test # discrete ss ordinary delay: tf-roundtrip through absorbDelay, cross-check freqresp
+%! sys = ss (0.5, 1, 1, 0, 0.1, "InputDelay", 2);
+%! sys2 = pade (sys, 3);
+%! assert (isa (sys2, "ss"));
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! expected = absorbDelay (tf (set (sys, "InputDelay", 2)));
+%! assert (freqresp (sys2, w), freqresp (expected, w), 1e-8);
+
+%!test # discrete ss InternalDelay: exact loop closure, matches the untouched
+%! # exact-delay discrete freqresp (no approximation error to expect, unlike
+%! # the continuous case -- this closure is exact by construction)
+%! G = ss (0.5, 1, 1, 0, 0.1, "IODelay", 3);
+%! L = feedback (G);
+%! assert (hasinternaldelay (L), true);
+%! sys2 = pade (L, 4);
+%! assert (hasinternaldelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (sys2, w), freqresp (L, w), 1e-9);
+
+%!test # discrete ss InternalDelay: order n has no effect on the exact result
+%! G = ss (0.5, 1, 1, 0, 0.1, "IODelay", 2);
+%! L = feedback (G);
+%! sys_a = pade (L, 2);
+%! sys_b = pade (L, 7);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (sys_a, w), freqresp (sys_b, w), 1e-12);
+
+%!test # discrete ss: combined InternalDelay AND ordinary delay, both exact
+%! G = ss (0.5, 1, 1, 0, 0.1, "IODelay", 2);
+%! L = set (feedback (G), "InputDelay", 3);
+%! assert (hasinternaldelay (L), true);
+%! assert (hasdelay (L), true);
+%! sys2 = pade (L, 4);
+%! assert (hasinternaldelay (sys2), false);
+%! assert (hasdelay (sys2), false);
+%! w = [0.1, 1, 5];
+%! assert (freqresp (sys2, w), freqresp (L, w), 1e-8);
+
+%!error pade (tf (1, [1, -0.5], 0.1, "InputDelay", 2), [1, 2])
diff --git a/inst/@lti/parallel.m b/inst/@lti/parallel.m
index 60d77cfc..c60a9f42 100644
--- a/inst/@lti/parallel.m
+++ b/inst/@lti/parallel.m
@@ -45,6 +45,17 @@
if (nargin == 2)
sys = sys1 + sys2;
+
+ if ((isa (sys1, "lti") && hasdelay (sys1)) || (isa (sys2, "lti") && hasdelay (sys2)))
+ d1 = totaldelay (sys1);
+ d2 = totaldelay (sys2);
+
+ if (! isequal (d1, d2))
+ error ("parallel: mismatched delays require internal delay support (not yet implemented)");
+ endif
+
+ sys = set (sys, "InputDelay", 0, "OutputDelay", 0, "IODelay", d1);
+ endif
## elseif (nargin == 6)
## TODO: implement "complicated" case sys = parallel (sys1, sys2, in1, in2, out1, out2)
@@ -54,3 +65,25 @@
endif
endfunction
+
+
+%!test # matching delays: exact composition
+%! s1 = tf (1, [1 1], "InputDelay", 0.1, "OutputDelay", 0.2);
+%! s2 = tf (1, [1 2], "InputDelay", 0.1, "OutputDelay", 0.2);
+%! s = parallel (s1, s2);
+%! assert (s.InputDelay, 0, 1e-10);
+%! assert (s.OutputDelay, 0, 1e-10);
+%! assert (s.IODelay, 0.3, 1e-10);
+
+%!test # no delay on either operand: result has no delay (regression)
+%! s1 = tf (1, [1 1]);
+%! s2 = tf (1, [1 2]);
+%! s = parallel (s1, s2);
+%! assert (hasdelay (s), false);
+
+%!error parallel (tf (1, [1 1], "InputDelay", 0.1, "OutputDelay", 0.2), tf (1, [1 2], "InputDelay", 0.3, "OutputDelay", 0.4))
+
+%!test # one operand is a plain numeric gain (no delay): should not error (regression)
+%! s1 = tf (1, [1 1]);
+%! s = parallel (s1, 2);
+%! assert (hasdelay (s), false);
diff --git a/inst/@lti/pole.m b/inst/@lti/pole.m
index 75dcab37..9827c3f4 100644
--- a/inst/@lti/pole.m
+++ b/inst/@lti/pole.m
@@ -49,6 +49,11 @@
function pol = pole (sys)
+ if (hasinternaldelay (sys))
+ error (["pole: InternalDelay is not yet supported for exact pole ", ...
+ "computation; use pade(sys, n) to approximate the delay first"]);
+ endif
+
if(nargin == 1) # pole(sys)
if(!(isa(sys, "lti")) && issquare(sys))
pol = eig(sys);
@@ -128,3 +133,7 @@
%!assert (infp, infp_exp);
%!assert (kronr, kronr_exp);
%!assert (kronl, kronl_exp);
+
+
+## Test for InternalDelay guard - should error with function-specific message
+%!error