-
Notifications
You must be signed in to change notification settings - Fork 11
implement native @zpk class #31
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
base: dev
Are you sure you want to change the base?
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -54,4 +54,4 @@ TAGS$ | |
| semantic.cache | ||
| # Other text editors often create these | ||
| ~.* | ||
|
|
||
| *.swp | ||
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @ttl-octave I applied the same fix from #29; let me know if you think that's incorrect/inappropriate. |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,100 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @ttl-octave one other note here, I took the notice from https://github.com/MitchellThompkins/pkg-control/blob/037fc5f50fb9666518667de018fe84865bbbac29/inst/%40tf/__c2d__.m#L1-L23 but that doesn't actually match https://github.com/gnu-octave/pkg-control/blob/eaba61e11b22de7e9b48871f939423e618d12f0e/CONTRIBUTING.md#license-and-documentation. I'm not sure which style is right. (This applies to both the GPL notice and the
Member
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Regarding the texinfo section, the help text for internal function The GPL notice of existing files still contains the hint to the name "LTI Syncope", which was introduced by Lukas Reichlin as far as I know. As this name was and is not used officially for the control package it is not referred to in the example comment in the guidelines for contributing (where I just saw the copy & paste error, "statistics package"). |
||
| ## | ||
| ## 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 <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); | ||
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,27 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <http://www.gnu.org/licenses/>. | ||
|
|
||
| ## -*- texinfo -*- | ||
| ## Complex-conjugate transpose ZPK model by delegating to TF. | ||
|
|
||
| ## Author: Mitchell Thompkins <mitchell.thompkins@pm.me> | ||
| ## Created: June 2026 | ||
| ## Version: 0.1 | ||
|
|
||
| function sys = __ctranspose__ (sys, ct) | ||
| sys = __ctranspose__ (tf (sys), ct); | ||
| endfunction |
|
Contributor
Author
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. @ttl-octave I applied the same fix from #29; let me know if you think that's incorrect/inappropriate. |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,75 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <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); |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,27 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <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 |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,38 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <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 |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,27 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <http://www.gnu.org/licenses/>. | ||
|
|
||
| ## -*- texinfo -*- | ||
| ## Return minimal realization of ZPK model by delegating to TF. | ||
|
|
||
| ## Author: Mitchell Thompkins <mitchell.thompkins@pm.me> | ||
| ## Created: June 2026 | ||
| ## Version: 0.1 | ||
|
|
||
| function sys = __minreal__ (sys, tol) | ||
| sys = __minreal__ (tf (sys), tol); | ||
| endfunction |
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,35 @@ | ||
| ## Copyright (C) 2026 Mitchell Thompkins | ||
| ## | ||
| ## 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 <http://www.gnu.org/licenses/>. | ||
|
|
||
| ## -*- texinfo -*- | ||
| ## Return poles of ZPK model directly from stored data (no polynomial round-trip). | ||
|
|
||
| ## Author: Mitchell Thompkins <mitchell.thompkins@pm.me> | ||
| ## Created: June 2026 | ||
| ## Version: 0.1 | ||
|
|
||
| function pol = __pole__ (sys) | ||
|
|
||
| if (issiso (sys)) | ||
| pol = sys.p{1}; | ||
| else | ||
| warning ("Control:convert-to-state-space", | ||
| "zpk: pole: converting to minimal state-space for pole computation of mimo zpk\n"); | ||
| pol = pole (ss (sys)); | ||
| endif | ||
|
|
||
| endfunction |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@ttl-octave I'm not sure how appropriate or necessary it is to add my own name here. I was just trying to follow suit.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This depends on the amount of changes or extension that you have contributed to the file. Im this case (and maybe in some other of this commit) I think there is too little originality. Moreover, your contribution is recorded in the git history anyways.