Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -54,4 +54,4 @@ TAGS$
semantic.cache
# Other text editors often create these
~.*

*.swp
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
Name: control
Version: 4.2.2+
Version: 4.2.99+
Date: XXXX-XX-XX
Author: Lukas Reichlin <lukas.reichlin@gmail.com>
Maintainer: Doug Stewart <doug.dastew@gmail.com>, Torsten Lilge <ttl-octave@mailbox.org>
Expand Down
6 changes: 6 additions & 0 deletions NEWS
Original file line number Diff line number Diff line change
@@ -1,6 +1,12 @@

Summary of important user-visible changes for releases of the control package

===============================================================================
control-4.2.99+ Upcoming 4.3.0, not yet released
===============================================================================

**

===============================================================================
control-4.2.2+ Upcoming 4.2.3, not yet released
===============================================================================
Expand Down
5 changes: 2 additions & 3 deletions inst/@tf/__c2d__.m
Original file line number Diff line number Diff line change
Expand Up @@ -35,8 +35,7 @@
sys=imp_invar(sys,1/tsam);

elseif (strncmpi (method, "m", 1)) # "matched"
## TODO: move this code to @zpk/__c2d__.m once ZPK models are implemented


if (! issiso (sys))
error ("tf: c2d: require SISO system for matched pole/zero method");
endif
Expand Down Expand Up @@ -68,7 +67,7 @@
w_d = exp (1j * w_c * tsam);
k_d = real (k_c * prod (1j*w_c - z_c) / prod (1j*w_c - p_c) * prod (w_d - p_d) / prod (w_d - z_d));

tmp = zpk (z_d, p_d, k_d, tsam);
tmp = tf (zpk (z_d, p_d, k_d, tsam));
sys.num = tmp.num;
sys.den = tmp.den;

Expand Down
7 changes: 3 additions & 4 deletions inst/@tf/__d2c__.m
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,6 @@
function sys = __d2c__ (sys, tsam, method = "zoh", w0 = 0)

if (strncmpi (method, "m", 1)) # "matched"
## TODO: move this code to @zpk/__d2c__.m once ZPK models are implemented

if (! issiso (sys))
error ("tf: d2c: require SISO system for matched pole/zero method");
Expand All @@ -48,11 +47,11 @@
tol = sqrt (eps);
while (any (abs ([p_d; z_d_orig] - w_d) < tol))
w_c += 0.1 / tsam;
w_d = exp (1j * w_c * tsam);
endwhile
w_d = exp (w_c * tsam);
k_c = real (k_d * prod (w_d - z_d_orig) / prod (w_d - p_d) * prod (w_c - p_c) / prod (w_c - z_c));
k_c = real (k_d * prod (w_d - z_d_orig) / prod (w_d - p_d) * prod (1j*w_c - p_c) / prod (1j*w_c - z_c));

tmp = zpk (z_c, p_c, k_c);
tmp = tf (zpk (z_c, p_c, k_c));
sys.num = tmp.num;
sys.den = tmp.den;

Expand Down
101 changes: 101 additions & 0 deletions inst/@zpk/__c2d__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Convert the continuous ZPK model into its discrete-time equivalent.
## The matched method maps each stored pole and zero directly via exp (s*tsam),
## avoiding the ill-conditioned polynomial round-trip of the TF representation.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: June 2026
## Version: 0.1

function sys = __c2d__ (sys, tsam, method = "zoh", w0 = 0)

if (strncmpi (method, "m", 1)) # "matched"

if (! issiso (sys))
error ("zpk: c2d: require SISO system for matched pole/zero method");
endif

z_c = sys.z{1};
p_c = sys.p{1};
k_c = sys.k;

p_d = exp (p_c * tsam);
z_d = exp (z_c * tsam);

if (any (! isfinite (p_d)) || any (! isfinite (z_d)))
error ("zpk: c2d: discrete-time poles and zeros are not finite");
endif

## continuous-time zeros at infinity are mapped to -1 in discrete-time
## except for one. for non-proper transfer functions, no zeros at -1 are added.
np = length (p_c); # number of poles
nz = length (z_c); # number of finite zeros, np-nz number of infinite zeros
z_d = vertcat (z_d, repmat (-1, np-nz-1, 1));

## the discrete-time gain k_d is matched at frequency w_c to continuous-time
## gain k_c. dc gain is taken (w_c=0) unless there are continuous-time
## poles/zeros near the imaginary axis at j*w_c. gain is evaluated on the
## imaginary axis (s=j*w_c) and unit circle (z=exp(j*w_c*tsam)) so that
## |H_d(exp(j*w_c*tsam))| = |H_c(j*w_c)| holds in the frequency domain.
w_c = 0; # start at dc
tol = sqrt (eps); # poles/zeros within tol of j*w_c are avoided
while (any (abs ([p_c; z_c] - 1j*w_c) < tol))
w_c += 0.1 / tsam;
endwhile
w_d = exp (1j * w_c * tsam);
k_d = real (k_c * prod (1j*w_c - z_c) / prod (1j*w_c - p_c) * prod (w_d - p_d) / prod (w_d - z_d));

sys.z{1} = z_d;
sys.p{1} = p_d;
sys.k = k_d;

else
## zoh/foh/tustin/prewarp/impulse are not per-root maps. their natural
## representation is state-space. convert back to zpk for type consistency.
sys = zpk (__c2d__ (ss (sys), tsam, method, w0));
endif

endfunction


%!test
%! ## single pole: p_d = exp(p_c * Ts)
%! sys = zpk ([], [-1], 1);
%! sys_d = c2d (sys, 0.1, 'matched');
%! assert (isa (sys_d, 'zpk'));
%! [~, p_d] = zpkdata (sys_d, 'v');
%! assert (p_d, exp (-0.1), 1e-14);

%!test
%! ## 25 poles clustered near the imaginary axis: matched c2d keeps them
%! ## stable and maps each one exactly to exp(p*Ts) with no rounding error
%! N = 25; Ts = 1/1000;
%! p_s = (-0.001 + 1i * linspace (1, 25, N))' * 2*pi*4;
%! sys_d = c2d (zpk ([], p_s, 1), Ts, 'matched');
%! [~, p_d] = zpkdata (sys_d, 'v');
%! assert (all (abs (p_d) < 1));
%! assert (max (abs (p_d - exp (p_s * Ts))), 0, 1e-12);

%!test
%! ## non-matched methods still work on zpk and return zpk
%! sys_d = c2d (zpk ([], [-1], 1), 0.1, 'zoh');
%! assert (isa (sys_d, 'zpk'));
%! assert (get (sys_d, 'tsam'), 0.1);
28 changes: 28 additions & 0 deletions inst/@zpk/__ctranspose__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Complex-conjugate transpose ZPK model by delegating to TF.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: July 2026
## Version: 0.2

function sys = __ctranspose__ (sys, ct)
sys = zpk (__ctranspose__ (tf (sys), ct));
endfunction
76 changes: 76 additions & 0 deletions inst/@zpk/__d2c__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,76 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Convert the discrete ZPK model into its continuous-time equivalent.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: June 2026
## Version: 0.1

function sys = __d2c__ (sys, tsam, method = "zoh", w0 = 0)

if (strncmpi (method, "m", 1)) # "matched"

if (! issiso (sys))
error ("zpk: d2c: require SISO system for matched pole/zero method");
endif

z_d = sys.z{1};
p_d = sys.p{1};
k_d = sys.k;

if (any (abs (p_d) < eps) || any (abs (z_d) < eps))
error ("zpk: d2c: discrete-time poles and zeros at 0 not supported because log(0) is -Inf");
endif

z_d_orig = z_d;
z_d(abs (z_d+1) < sqrt (eps)) = []; # remove zeros added at -1 by c2d

p_c = log (p_d) / tsam;
z_c = log (z_d) / tsam;

w_c = 0;
w_d = 1;
tol = sqrt (eps);
while (any (abs ([p_d; z_d_orig] - w_d) < tol))
w_c += 0.1 / tsam;
w_d = exp (1j * w_c * tsam);
endwhile
k_c = real (k_d * prod (w_d - z_d_orig) / prod (w_d - p_d) * prod (1j*w_c - p_c) / prod (1j*w_c - z_c));

sys.z{1} = z_c;
sys.p{1} = p_c;
sys.k = k_c;

else
sys = zpk (__d2c__ (ss (sys), tsam, method, w0));
endif

endfunction


%!test
%! ## matched round trip is the identity on poles/zeros
%! sys = zpk ([-3], [-1; -2], 4);
%! sysd = c2d (sys, 0.1, 'matched');
%! sysc = d2c (sysd, 'matched');
%! [z, p, k] = zpkdata (sysc, 'v');
%! assert (sort (real (p)), [-2; -1], 1e-10);
%! assert (real (z), -3, 1e-10);
%! assert (k, 4, 1e-10);
28 changes: 28 additions & 0 deletions inst/@zpk/__freqresp__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Compute frequency response of ZPK model by delegating to TF.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: June 2026
## Version: 0.1

function H = __freqresp__ (sys, w, cellflag = false)
H = __freqresp__ (tf (sys), w, cellflag);
endfunction
39 changes: 39 additions & 0 deletions inst/@zpk/__get__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Access key values of ZPK objects.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: June 2026
## Version: 0.1

function val = __get__ (sys, key)

switch (key)
case {"z", "zeros"}
val = sys.z;
case {"p", "poles"}
val = sys.p;
case {"k", "gain"}
val = sys.k;
otherwise
error ("zpk: get: invalid key name '%s'", key);
endswitch

endfunction
28 changes: 28 additions & 0 deletions inst/@zpk/__minreal__.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
## Copyright (C) 2026 Mitchell Thompkins <mitchell.thompkins@pm.me>
##
## This file is part of the control package for GNU Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING. If not,
## see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## Return minimal realization of ZPK model by delegating to TF.

## Author: Mitchell Thompkins <mitchell.thompkins@pm.me>
## Created: July 2026
## Version: 0.2

function sys = __minreal__ (sys, tol)
sys = zpk (__minreal__ (tf (sys), tol));
endfunction
Loading