Improve numerical accuracy#47
Conversation
…y gain" This reverts commit ed93c83.
Butterworth and Elliptic now carry their analytically-computed poles and zeros as a FactoredTF through to the discretization step. ZOH skips the QR eigenvalue search on the CCF matrix; Matched-Z skips both companion matrix eigendecompositions entirely. Tustin is unchanged. Adds FactoredTF struct, matrix_exp_with_evals, zoh_discretize_with_evals, matched_z_discretize_factored, and discretize_with_factored dispatch overloads to discretize.hpp.
CCF eigenvectors are Vandermonde columns, so V is built directly from the analytic poles rather than via inverse iteration on the ill-conditioned CCF matrix. Fixes Butterworth ZOH accuracy to N=12.
📝 WalkthroughWalkthroughAdds factored pole/zero support for analog filters, new factored ZOH and Matched-Z discretization paths, and constructor wiring so Butterworth and Elliptic pass factored transfer functions into ChangesFactored analog filter and discretization
Estimated code review effort🎯 4 (Complex) | ⏱️ ~60 minutes 🚥 Pre-merge checks | ✅ 3 | ❌ 2❌ Failed checks (1 warning, 1 inconclusive)
✅ Passed checks (3 passed)
✨ Finishing Touches📝 Generate docstrings
🧪 Generate unit tests (beta)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
There was a problem hiding this comment.
Actionable comments posted: 2
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@include/constfilt/analog_filter.hpp`:
- Around line 40-49: Update the AnalogFilter documentation comment in
analog_filter.hpp so the Vandermonde eigenvector matrix wording applies only to
the ZOH path, not MatchedZ. In the comment near AnalogFilter(continuous_tf,
factored_tf, sample_rate_hz, method_tag), revise the description of factored_tf
to reflect that MatchedZ uses the supplied poles/zeros directly, while only ZOH
uses the known poles to build the Vandermonde matrix analytically. Keep the note
that both continuous_tf and factored_tf are required for direct callers.
In `@include/constfilt/discretize.hpp`:
- Line 112: The new factored helper templates in this file use consteig::Size
for the order parameter N, but the project guideline requires unsigned int for
filter order template parameters. Update FactoredTF and the related helper
template declarations identified in this file to use unsigned int N instead,
while leaving consteig::Size for loop indices and matrix-size calculations where
needed; keep the existing behavior and only adjust the template parameter type
in the affected helper templates.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: db351af4-41ae-4f5b-94d4-26feef90a7cc
⛔ Files ignored due to path filters (12)
docs/profiling/results/accuracy_gcc_15.2.0.csvis excluded by!**/*.csvdocs/profiling/results/accuracy_gcc_15.2.0.pngis excluded by!**/*.pngdocs/profiling/results/accuracy_gcc_15.2.0_butterworth.pngis excluded by!**/*.pngdocs/profiling/results/accuracy_gcc_15.2.0_constfilt.pngis excluded by!**/*.pngdocs/profiling/results/accuracy_gcc_15.2.0_elliptic.pngis excluded by!**/*.pngdocs/profiling/results/compile_times_gcc_15.2.0.csvis excluded by!**/*.csvdocs/profiling/results/compile_times_gcc_15.2.0.pngis excluded by!**/*.pngdocs/profiling/results/runtime_gcc_15.2.0.csvis excluded by!**/*.csvdocs/profiling/results/runtime_gcc_15.2.0.pngis excluded by!**/*.pngdocs/profiling/results/runtime_gcc_15.2.0_butterworth.pngis excluded by!**/*.pngdocs/profiling/results/runtime_gcc_15.2.0_constfilt.pngis excluded by!**/*.pngdocs/profiling/results/runtime_gcc_15.2.0_elliptic.pngis excluded by!**/*.png
📒 Files selected for processing (4)
include/constfilt/analog_filter.hppinclude/constfilt/butterworth.hppinclude/constfilt/discretize.hppinclude/constfilt/elliptic.hpp
📜 Review details
🧰 Additional context used
📓 Path-based instructions (5)
**/*.{hpp,h,cpp}
📄 CodeRabbit inference engine (AGENTS.md)
**/*.{hpp,h,cpp}: Use C++17 only; do not use C++20 features (no std::is_constant_evaluated, no concepts)
Enforce clang-format code formatting via CI; run make container.make.format before committing
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hppinclude/constfilt/butterworth.hppinclude/constfilt/elliptic.hpp
include/constfilt/**
📄 CodeRabbit inference engine (AGENTS.md)
Maintain header-only library structure; do not add .cpp translation units to the library itself
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hppinclude/constfilt/butterworth.hppinclude/constfilt/elliptic.hpp
**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
Implement constexpr everywhere possible; coefficient computation must work at compile time; Filter::operator() batch path must remain constexpr
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hppinclude/constfilt/butterworth.hppinclude/constfilt/elliptic.hpp
include/constfilt/**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
include/constfilt/**/*.{hpp,h}: Do not introduce hard-coded filter-order ceiling; Filter classes should only use static_assert(N >= 1)
Use unsigned int for template parameter order N; implicit conversion with consteig::Size (size_t) is acceptable; do not add casts unless compiler error requires it
Follow a[0] = 1 convention for discrete transfer functions (monic-denominator); Filter stores _a with a[0] included for uniformity but DF2T equations assume a[0] = 1
Implement stability check via CheckStab template parameter (default true) to reject Stability::Unstable systems while accepting Stable and MarginallyStable; check not yet implemented (issue#14)
Use DirectForm II Transposed (DF2T) implementation for runtime filter with both sample-by-sample real-time path and constexpr-capable batch path
Support Butterworth filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options
Support Elliptic filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options; use Nome-series Cauer design (see docs/elliptic.md)
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hppinclude/constfilt/butterworth.hppinclude/constfilt/elliptic.hpp
**/*.hpp
⚙️ CodeRabbit configuration file
**/*.hpp: This is a constexpr C++17 header-only filter design library. Focus on:
- Correctness of constexpr computations and static_assert usage
- Type consistency (float vs double tolerances must match the filter's numeric type)
- No standard library dependencies in core library code (freestanding constraint)
- Correct use of vendored consteig (matrix ops) and gcem (math functions) via the wrapper
- Headers named *_reference.hpp contain precomputed test data — do not suggest algorithmic changes to their values
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hppinclude/constfilt/butterworth.hppinclude/constfilt/elliptic.hpp
🔇 Additional comments (4)
include/constfilt/discretize.hpp (1)
106-119: LGTM!Also applies to: 181-250, 287-308, 685-834, 958-986
include/constfilt/analog_filter.hpp (1)
92-100: LGTM!Also applies to: 143-151
include/constfilt/butterworth.hpp (1)
37-39: LGTM!Also applies to: 95-132
include/constfilt/elliptic.hpp (1)
46-47: LGTM!Also applies to: 311-349, 493-543, 545-604
There was a problem hiding this comment.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
include/constfilt/discretize.hpp (1)
951-952: 🎯 Functional Correctness | 🟠 Major | 🏗️ Heavy liftAvoid routing factored ZOH back through
ss_to_tf.
ss_to_tf(...)callschar_poly(sys_d.A), andchar_polyrecovers the denominator viaconsteig::eigenvalues. That reintroduces the eigendecomposition/root-recovery step the factored path is meant to bypass for preserving analytic poles. Build the discrete denominator directly fromexp(Ts * factored_tf.poles[i]), then reusemarkov_numeratorforb.🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@include/constfilt/discretize.hpp` around lines 951 - 952, The factored ZOH path in `zoh_discretize_with_evals` is being sent back through `ss_to_tf`, which reintroduces denominator reconstruction via `char_poly` and `consteig::eigenvalues`. Update the code to build the discrete denominator directly from the analytic poles in `factored_tf` using `exp(Ts * factored_tf.poles[i])`, and keep the numerator computation on the existing `markov_numerator` path for `b` so the factored path stays factor-preserving.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Outside diff comments:
In `@include/constfilt/discretize.hpp`:
- Around line 951-952: The factored ZOH path in `zoh_discretize_with_evals` is
being sent back through `ss_to_tf`, which reintroduces denominator
reconstruction via `char_poly` and `consteig::eigenvalues`. Update the code to
build the discrete denominator directly from the analytic poles in `factored_tf`
using `exp(Ts * factored_tf.poles[i])`, and keep the numerator computation on
the existing `markov_numerator` path for `b` so the factored path stays
factor-preserving.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: 905d8f8a-f4f6-417b-87af-933eafe0a70e
📒 Files selected for processing (1)
include/constfilt/discretize.hpp
📜 Review details
⏰ Context from checks skipped due to timeout. (8)
- GitHub Check: gcc
- GitHub Check: cross.arm-clang
- GitHub Check: cross.arm-gcc
- GitHub Check: clang
- GitHub Check: cross.arm-gcc
- GitHub Check: gcc
- GitHub Check: cross.arm-clang
- GitHub Check: clang
🧰 Additional context used
📓 Path-based instructions (5)
**/*.{hpp,h,cpp}
📄 CodeRabbit inference engine (AGENTS.md)
**/*.{hpp,h,cpp}: Use C++17 only; do not use C++20 features (no std::is_constant_evaluated, no concepts)
Enforce clang-format code formatting via CI; run make container.make.format before committing
Files:
include/constfilt/discretize.hpp
include/constfilt/**
📄 CodeRabbit inference engine (AGENTS.md)
Maintain header-only library structure; do not add .cpp translation units to the library itself
Files:
include/constfilt/discretize.hpp
**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
Implement constexpr everywhere possible; coefficient computation must work at compile time; Filter::operator() batch path must remain constexpr
Files:
include/constfilt/discretize.hpp
include/constfilt/**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
include/constfilt/**/*.{hpp,h}: Do not introduce hard-coded filter-order ceiling; Filter classes should only use static_assert(N >= 1)
Use unsigned int for template parameter order N; implicit conversion with consteig::Size (size_t) is acceptable; do not add casts unless compiler error requires it
Follow a[0] = 1 convention for discrete transfer functions (monic-denominator); Filter stores _a with a[0] included for uniformity but DF2T equations assume a[0] = 1
Implement stability check via CheckStab template parameter (default true) to reject Stability::Unstable systems while accepting Stable and MarginallyStable; check not yet implemented (issue#14)
Use DirectForm II Transposed (DF2T) implementation for runtime filter with both sample-by-sample real-time path and constexpr-capable batch path
Support Butterworth filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options
Support Elliptic filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options; use Nome-series Cauer design (see docs/elliptic.md)
Files:
include/constfilt/discretize.hpp
**/*.hpp
⚙️ CodeRabbit configuration file
**/*.hpp: This is a constexpr C++17 header-only filter design library. Focus on:
- Correctness of constexpr computations and static_assert usage
- Type consistency (float vs double tolerances must match the filter's numeric type)
- No standard library dependencies in core library code (freestanding constraint)
- Correct use of vendored consteig (matrix ops) and gcem (math functions) via the wrapper
- Headers named *_reference.hpp contain precomputed test data — do not suggest algorithmic changes to their values
Files:
include/constfilt/discretize.hpp
🔇 Additional comments (1)
include/constfilt/discretize.hpp (1)
122-193: LGTM!Also applies to: 225-232, 278-290, 679-816, 956-960
| // analytically from the known poles, bypassing the | ||
| // polynomial->eigenvalue roundtrip that degrades accuracy | ||
| // at high orders. MatchedZ maps the supplied poles/zeros | ||
| // directly via z=exp(p*Ts), skipping companion-matrix |
There was a problem hiding this comment.
This is erroneously perspective. The comment is true but reading this you have no context so remove it.
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
include/constfilt/discretize.hpp (1)
516-550: 🎯 Functional Correctness | 🟠 Major | ⚡ Quick winMatch gain on the imaginary axis.
Matched-Z frequency matching should evaluate
H_c(jω)andH_d(e^{jωTs}); the current code evaluates along the real axis (s = ω,z = e^{ωTs}). Highpass filters with an origin zero will move to a nonzerow_c, making the gain scaling wrong. As per coding guidelines: “Support Butterworth filters with lowpass and highpass modes…” and “Support Elliptic filters with lowpass and highpass modes…”Proposed fix
- const T dr = w_c - poles[i].real; - const T di = poles[i].imag; + const T dr = poles[i].real; + const T di = w_c - poles[i].imag; @@ - const T dr = w_c - zeros[i].real; - const T di = zeros[i].imag; + const T dr = zeros[i].real; + const T di = w_c - zeros[i].imag; @@ - const Complex w_c_cx{w_c, static_cast<T>(0)}; + const Complex w_c_cx{static_cast<T>(0), w_c}; const Complex w_d_cx = - consteig::exp(w_c_cx * Complex{Ts, static_cast<T>(0)}); + consteig::exp(Complex{static_cast<T>(0), w_c * Ts});🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@include/constfilt/discretize.hpp` around lines 516 - 550, The gain-matching logic in the matched-Z path is using the real axis instead of the imaginary axis, so update the frequency evaluation in the Step 4/5 block to compute H_c(jω) and H_d(e^{jωTs}) rather than H_c(ω) and H_d(e^{ωTs}). Keep the collision-avoidance search for w_c, but ensure the final gain scaling in discretize.hpp uses the imaginary-axis complex frequencies via the existing poles/zeros handling and the discrete gain computation around w_c_cx and w_d_cx.Source: Coding guidelines
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Inline comments:
In `@include/constfilt/discretize.hpp`:
- Around line 459-464: Remove the stray standalone template declaration before
matched_z_assemble, and keep only the actual function template so the
declaration parses correctly. In matched_z_assemble, fix the gain-matching logic
to evaluate the continuous and discrete transfer functions at the complex
frequency jωc and z = e^{jωcTs} rather than using the real-valued wc/exp(wcTs);
update the gain scaling path in that function so the matched-Z normalization is
computed from the correct complex response.
---
Outside diff comments:
In `@include/constfilt/discretize.hpp`:
- Around line 516-550: The gain-matching logic in the matched-Z path is using
the real axis instead of the imaginary axis, so update the frequency evaluation
in the Step 4/5 block to compute H_c(jω) and H_d(e^{jωTs}) rather than H_c(ω)
and H_d(e^{ωTs}). Keep the collision-avoidance search for w_c, but ensure the
final gain scaling in discretize.hpp uses the imaginary-axis complex frequencies
via the existing poles/zeros handling and the discrete gain computation around
w_c_cx and w_d_cx.
🪄 Autofix (Beta)
Fix all unresolved CodeRabbit comments on this PR:
- Push a commit to this branch (recommended)
- Create a new PR with the fixes
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: eb0fabc8-cd97-42f4-b706-3053143e9dee
📒 Files selected for processing (1)
include/constfilt/discretize.hpp
📜 Review details
⏰ Context from checks skipped due to timeout. (8)
- GitHub Check: cross.arm-gcc
- GitHub Check: clang
- GitHub Check: gcc
- GitHub Check: cross.arm-clang
- GitHub Check: cross.arm-gcc
- GitHub Check: gcc
- GitHub Check: clang
- GitHub Check: cross.arm-clang
🧰 Additional context used
📓 Path-based instructions (5)
**/*.{hpp,h,cpp}
📄 CodeRabbit inference engine (AGENTS.md)
**/*.{hpp,h,cpp}: Use C++17 only; do not use C++20 features (no std::is_constant_evaluated, no concepts)
Enforce clang-format code formatting via CI; run make container.make.format before committing
Files:
include/constfilt/discretize.hpp
include/constfilt/**
📄 CodeRabbit inference engine (AGENTS.md)
Maintain header-only library structure; do not add .cpp translation units to the library itself
Files:
include/constfilt/discretize.hpp
**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
Implement constexpr everywhere possible; coefficient computation must work at compile time; Filter::operator() batch path must remain constexpr
Files:
include/constfilt/discretize.hpp
include/constfilt/**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
include/constfilt/**/*.{hpp,h}: Do not introduce hard-coded filter-order ceiling; Filter classes should only use static_assert(N >= 1)
Use unsigned int for template parameter order N; implicit conversion with consteig::Size (size_t) is acceptable; do not add casts unless compiler error requires it
Follow a[0] = 1 convention for discrete transfer functions (monic-denominator); Filter stores _a with a[0] included for uniformity but DF2T equations assume a[0] = 1
Implement stability check via CheckStab template parameter (default true) to reject Stability::Unstable systems while accepting Stable and MarginallyStable; check not yet implemented (issue#14)
Use DirectForm II Transposed (DF2T) implementation for runtime filter with both sample-by-sample real-time path and constexpr-capable batch path
Support Butterworth filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options
Support Elliptic filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options; use Nome-series Cauer design (see docs/elliptic.md)
Files:
include/constfilt/discretize.hpp
**/*.hpp
⚙️ CodeRabbit configuration file
**/*.hpp: This is a constexpr C++17 header-only filter design library. Focus on:
- Correctness of constexpr computations and static_assert usage
- Type consistency (float vs double tolerances must match the filter's numeric type)
- No standard library dependencies in core library code (freestanding constraint)
- Correct use of vendored consteig (matrix ops) and gcem (math functions) via the wrapper
- Headers named *_reference.hpp contain precomputed test data — do not suggest algorithmic changes to their values
Files:
include/constfilt/discretize.hpp
🔇 Additional comments (1)
include/constfilt/discretize.hpp (1)
106-118: LGTM!Also applies to: 122-194, 196-233, 271-291, 460-503, 547-704, 829-857
| template <typename T, consteig::Size N> | ||
| constexpr TransferFunction<T, N + 1u, N + 1u> matched_z_discretize_tf( | ||
| const T (&b_c)[N + 1u], const T (&a_c)[N + 1u], T Ts, MatchedZ /*tag*/) | ||
| // Shared kernel: given analog poles, finite zeros, zero count, continuous gain, | ||
| // and sample period, maps to z-domain, pads missing zeros at z=-1, matches | ||
| // gain, and assembles the discrete TF. | ||
| template <typename T, consteig::Size N> | ||
| constexpr TransferFunction<T, N + 1u, N + 1u> matched_z_assemble( |
There was a problem hiding this comment.
🎯 Functional Correctness | 🔴 Critical | ⚡ Quick win
🧩 Analysis chain
🏁 Script executed:
#!/bin/bash
# Verify there is only one template head before matched_z_assemble.
sed -n '453,466p' include/constfilt/discretize.hppRepository: MitchellThompkins/constfilt
Length of output: 948
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== matched_z_assemble definition =="
sed -n '459,560p' include/constfilt/discretize.hpp
echo
echo "== matched_z_assemble call sites =="
rg -n "matched_z_assemble\\(" include/constfilt/discretize.hpp
echo
echo "== surrounding matched-Z helpers =="
sed -n '560,760p' include/constfilt/discretize.hppRepository: MitchellThompkins/constfilt
Length of output: 11003
🏁 Script executed:
#!/bin/bash
set -euo pipefail
echo "== matched-z gain-matching block =="
sed -n '500,545p' include/constfilt/discretize.hpp
echo
echo "== docs / references mentioning matched-Z =="
rg -n "matched-z|matched Z|c2d__|test frequency|gain matching" -S include docs .Repository: MitchellThompkins/constfilt
Length of output: 3446
Remove the stray template head and fix matched-Z gain evaluation.
- Line 459 leaves a standalone
template <typename T, consteig::Size N>before the actual declaration, which is a compile error. - The gain-matching block evaluates at real
w_c/exp(w_c T_s)instead ofjω_c/e^{jω_c T_s}, so the matched-Z scaling is off.
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
In `@include/constfilt/discretize.hpp` around lines 459 - 464, Remove the stray
standalone template declaration before matched_z_assemble, and keep only the
actual function template so the declaration parses correctly. In
matched_z_assemble, fix the gain-matching logic to evaluate the continuous and
discrete transfer functions at the complex frequency jωc and z = e^{jωcTs}
rather than using the real-valued wc/exp(wcTs); update the gain scaling path in
that function so the matched-Z normalization is computed from the correct
complex response.
There was a problem hiding this comment.
On the second point:
The gain-matching block evaluates at real w_c/exp(w_c T_s) instead of jω_c/e^{jω_c T_s}, so the matched-Z scaling is off.
I took the implementation from Octave.
## Copyright (C) 2009-2016 Lukas F. Reichlin
##
## 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 TF model into its discrete-time equivalent.
## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
## Created: October 2009
## Version: 0.3
function sys = __c2d__ (sys, tsam, method = "zoh", w0 = 0)
if (strncmpi (method, "i", 1)) # "impulse invariant"
[~,~,~,D] = ssdata (sys);
if (any (D(:)))
error ("c2d: impuls invariant discrete-time models only supported for systems without direct feedthrough\n");
endif
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
[z_c, p_c, k_c] = zpkdata (sys, "vector");
p_d = exp (p_c * tsam);
z_d = exp (z_c * tsam);
if (any (! isfinite (p_d)) || any (! isfinite (z_d)))
error ("tf: 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 a certain frequency (w_c, w_d)
## to continuous-time gain k_c. the dc gain is taken (w_c=0, w_d=1) unless
## there are continuous-time poles/zeros near 0. then w_c=1/tsam is taken.
w_c = 0; # dc gain
tol = sqrt (eps); # poles/zeros below tol are assumed to be zero
while (any (abs ([p_c; z_c] - w_c) < tol))
w_c += 0.1 / tsam;
endwhile
w_d = exp (w_c * tsam);
k_d = real (k_c * prod (w_c - z_c) / prod (w_c - p_c) * prod (w_d - p_d) / prod (w_d - z_d));
tmp = zpk (z_d, p_d, k_d, tsam);
sys.num = tmp.num;
sys.den = tmp.den;
else
[p, m] = size (sys);
for i = 1 : p
for j = 1 : m
idx = substruct ("()", {i, j});
tmp = subsref (sys, idx);
tmp = c2d (ss (tmp), tsam, method, w0);
[num, den] = tfdata (tmp, "tfpoly");
sys.num(i, j) = num;
sys.den(i, j) = den;
endfor
endfor
endif
sys.tfvar = "z";
endfunction
I think my implementation aligns with theirs. That's not saying Octave is right, but I can't find a reference that says that the gain-matching block at w_c/exp(w_c T_s) is wrong (though I agree it's weird). #51 documents this.
There was a problem hiding this comment.
@MitchellThompkins Good catch — thanks for pulling the Octave source.
You're right: the implementation in include/constfilt/discretize.hpp is aligned with Octave's matched discretization path, so I don't have enough basis to keep pushing on the jω_c / e^{jω_c T_s} point here. I'll drop that part of the review.
That leaves only the separate stray-template compile issue as the actionable item from my original comment. Tracking the theory/reference question in #51 is the right move.
(◕ᵕ◕)ᐢ
✏️ Learnings added
Learnt from: MitchellThompkins
Repo: MitchellThompkins/constfilt PR: 47
File: include/constfilt/discretize.hpp:459-464
Timestamp: 2026-06-27T06:04:41.991Z
Learning: In `include/constfilt/discretize.hpp`, the matched-Z gain-matching logic in `matched_z_assemble(...)` intentionally follows Octave control package `tf/__c2d__.m`, using real `w_c` and `exp(w_c * Ts)` rather than `jω_c` / `e^{jω_c Ts}`; any theoretical reconsideration of that choice is tracked in GitHub issue `#51`.
Note: Learnings are effective only in the context of similar code segments. To apply general coding guidelines, please configure review instructions. You can manage existing learnings in the app.
|
Note Docstrings generation - SKIPPED |
Docstrings generation was requested by @MitchellThompkins. * #47 (comment) The following files were modified: * `include/constfilt/analog_filter.hpp` * `include/constfilt/butterworth.hpp` * `include/constfilt/discretize.hpp` * `include/constfilt/elliptic.hpp`
|
@coderabbitai resume |
✅ Action performedReviews resumed. |
There was a problem hiding this comment.
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
include/constfilt/discretize.hpp (1)
116-116: 🩺 Stability & Availability | 🟠 Major | ⚡ Quick winBound
nzbefore indexing factored zeros.Line 116 exposes
nzthrough the publicFactoredTF, and Line 701 forwards it directly. If a direct caller suppliesnz > N, thezeros[]loops read past the array andpad = N - num_degcan underflow. Bound the count once before all indexing.Proposed fix
constexpr TransferFunction<T, N + 1u, N + 1u> matched_z_assemble( const consteig::Complex<T> (&poles)[N], const consteig::Complex<T> (&zeros)[N], consteig::Size nz, T k_c, T Ts) { using Complex = consteig::Complex<T>; + const consteig::Size nz_used = (nz < N) ? nz : N; @@ - for (consteig::Size k = 0; k < nz; ++k) + for (consteig::Size k = 0; k < nz_used; ++k) @@ - const consteig::Size n_extra = (nz + 1u < N) ? (N - nz - 1u) : 0u; + const consteig::Size n_extra = + (nz_used + 1u < N) ? (N - nz_used - 1u) : 0u; @@ - const consteig::Size cur_deg = nz + e; + const consteig::Size cur_deg = nz_used + e; @@ - const consteig::Size num_deg = nz + n_extra; + const consteig::Size num_deg = nz_used + n_extra; @@ - for (consteig::Size i = 0; i < nz && !collision; ++i) + for (consteig::Size i = 0; i < nz_used && !collision; ++i) @@ - for (consteig::Size i = 0; i < nz; ++i) + for (consteig::Size i = 0; i < nz_used; ++i) @@ - for (consteig::Size i = 0; i < nz; ++i) + for (consteig::Size i = 0; i < nz_used; ++i)Also applies to: 485-509, 526-548, 566-594, 701-702
🤖 Prompt for AI Agents
Verify each finding against current code. Fix only still-valid issues, skip the rest with a brief reason, keep changes minimal, and validate. In `@include/constfilt/discretize.hpp` at line 116, The FactoredTF zero-count field nz is used directly in multiple zeros-processing paths, so a caller can pass a value larger than the available zeros and trigger out-of-bounds reads or underflow. Add a single clamp for nz to the valid range before any indexing or pad calculations, and ensure the bounded value is the one forwarded from FactoredTF into the zero-handling loops in the related discretize paths and the caller at FactoredTF::operator()/forwarding site.
🤖 Prompt for all review comments with AI agents
Verify each finding against current code. Fix only still-valid issues, skip the
rest with a brief reason, keep changes minimal, and validate.
Outside diff comments:
In `@include/constfilt/discretize.hpp`:
- Line 116: The FactoredTF zero-count field nz is used directly in multiple
zeros-processing paths, so a caller can pass a value larger than the available
zeros and trigger out-of-bounds reads or underflow. Add a single clamp for nz to
the valid range before any indexing or pad calculations, and ensure the bounded
value is the one forwarded from FactoredTF into the zero-handling loops in the
related discretize paths and the caller at FactoredTF::operator()/forwarding
site.
ℹ️ Review info
⚙️ Run configuration
Configuration used: Path: .coderabbit.yaml
Review profile: ASSERTIVE
Plan: Pro
Run ID: c56c2540-ba9a-41b6-8f73-23c9f966c643
📒 Files selected for processing (2)
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
📜 Review details
🧰 Additional context used
📓 Path-based instructions (5)
**/*.{hpp,h,cpp}
📄 CodeRabbit inference engine (AGENTS.md)
**/*.{hpp,h,cpp}: Use C++17 only; do not use C++20 features (no std::is_constant_evaluated, no concepts)
Enforce clang-format code formatting via CI; run make container.make.format before committing
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
include/constfilt/**
📄 CodeRabbit inference engine (AGENTS.md)
Maintain header-only library structure; do not add .cpp translation units to the library itself
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
Implement constexpr everywhere possible; coefficient computation must work at compile time; Filter::operator() batch path must remain constexpr
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
include/constfilt/**/*.{hpp,h}
📄 CodeRabbit inference engine (AGENTS.md)
include/constfilt/**/*.{hpp,h}: Do not introduce hard-coded filter-order ceiling; Filter classes should only use static_assert(N >= 1)
Use unsigned int for template parameter order N; implicit conversion with consteig::Size (size_t) is acceptable; do not add casts unless compiler error requires it
Follow a[0] = 1 convention for discrete transfer functions (monic-denominator); Filter stores _a with a[0] included for uniformity but DF2T equations assume a[0] = 1
Implement stability check via CheckStab template parameter (default true) to reject Stability::Unstable systems while accepting Stable and MarginallyStable; check not yet implemented (issue#14)
Use DirectForm II Transposed (DF2T) implementation for runtime filter with both sample-by-sample real-time path and constexpr-capable batch path
Support Butterworth filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options
Support Elliptic filters with lowpass and highpass modes, arbitrary order (N >= 1), and ZOH or Matched-Z discretization options; use Nome-series Cauer design (see docs/elliptic.md)
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
**/*.hpp
⚙️ CodeRabbit configuration file
**/*.hpp: This is a constexpr C++17 header-only filter design library. Focus on:
- Correctness of constexpr computations and static_assert usage
- Type consistency (float vs double tolerances must match the filter's numeric type)
- No standard library dependencies in core library code (freestanding constraint)
- Correct use of vendored consteig (matrix ops) and gcem (math functions) via the wrapper
- Headers named *_reference.hpp contain precomputed test data — do not suggest algorithmic changes to their values
Files:
include/constfilt/analog_filter.hppinclude/constfilt/discretize.hpp
🔇 Additional comments (2)
include/constfilt/discretize.hpp (1)
120-234: LGTM!Also applies to: 271-292, 455-484, 511-525, 535-541, 600-690, 695-700, 703-704, 828-856
include/constfilt/analog_filter.hpp (1)
41-46: LGTM!Also applies to: 89-97, 140-148
1 participant
Resolves #33.
Description
The motivation behind this PR was to avoid doing repeat work that I suspected led to deterioration in numerical accuracy (that part actually wasn't completely accurate and is documented in the Accuracy results below).
AnalogFiltertook a polynomial, so the butterworth and elliptic filter designs analytically determined their poles but then turned that into a polynomial. However the pole/zero formulas in Butterworth and Elliptic design algorithms are closed-form expressions for the eigenvalues of the system. I designed the system so I have the eigenvalues. I turned it into a polynomial b/c it felt standard but didn't really think too hard about "oh I need to recover those".Thus the original design was such that
AnalogFiltertook an arbitrary continuous-time transfer function in the form of polynomial coefficients (like a filter you design in matlab), converted to state-space (CCF), then let each discretization method derive whatever it needs from that. For Tustin that is pure matrix algebra (still on CCF matrix), for Matched-Z it needs eigenvalues (roots of the polynomial), for ZOH it needs the eigenvalues and eigenvectors to compute the matrix exponential.If the filter design is already doing the heavy lifting of finding the analytical poles, there's no reason to then collapse that info into a state-space form only to re-extract it. It was clean and mathematically consistent but incurred unnecessary extra work and introduced numerical instability. The solution was to add new constructors that took the analytically derived poles/zeros (again, the eigenvalues) directly.
The other functions still exist because they have no reason not to. It's still perfectly valid to define an
AnalogFilterin terms of a polynomial, and if you do that, you have no choice but to perform eigendecomposition (at least for ZOH and Matched-Z).More on this below, but regarding numerical accuracy I thought I was going to remove eigendecomposition for all discretization methods and see a huge improvement but I didn't really. I started with trying #25 but saw basically no change. The only impactful change was really the removal of the
consteig::eigenvector()call in the ZOH method when designing poles analytically.Compile Time Results
Compile times improved. This follows naturally because I'm just not doing iterative matrix operations that the QR algorithm requires to get the eigenvalues.
Before
After
Accuracy Results
There actually wasn't a huge shift here like I had hoped. ZOH at high orders definitely got better (though the elliptic filter still seems like it may have an accuracy issue which is tracked by #54). But that's about it. In fact, this didn't largely improve at all until I removed the eigenvector call in the ZOH discretization method. Even if you have the eigenvalues analytically, you need to solve the eigenvectors. The discrete
B_dis computed asA_c^(-1) * (e^(A_c*Ts) - I) * B_c, which requirese^(A_c*Ts). That's the full matrix exponential so eigenvectors are unavoidable for ZOH regardless of whether you have analytic eigenvalues.consteig::eigenvectorsuses inverse iteration on the CCF A matrix to find the eigenvectors numerically and the CCF A matrix has the polynomial coefficients in its last row. At high orders those coefficients are very large numbers, making the matrix ill-conditioned. Inverse iteration on an ill-conditioned matrix produces inaccurate eigenvectors, which drives error in the ultimate derived coefficients.There is a quick fix to this that I used a Vandermonde matrix to find the eigenvectors. This works only because the A matrix is already in CCF. The CCF has a special structure which you can verify algebraically that multiplying the CCF A matrix by
[1, lambda, lambda^2, ..., lambda^(N-1)]gives backlambda * [1, lambda, lambda^2, ..., lambda^(N-1)]. When we use the constructed eigenvectorsVin conjunction withA_Cto findB_d, ourVis only valid specifically because of how we definedA_c.The definition of an eigenvector is
A·v = lambda·vwhich means applying the matrix to the vector just scales it by the eigenvalue. For CCF, if you takev = [1, lambda, lambda^2, ..., lambda^(N-1)]and multiply by the CCF A matrix, you getlambda * [1, lambda, lambda^2, ..., lambda^(N-1)]. That's exactlyA·v = lambda·vwhich is exactly the definition of the eigenvector. So for CCF specifically, the Vandermonde vectors are eigenvectors by construction. For a general matrix this wouldn't hold.Looking at it in detail, here is the CCF A matrix:
The last row contains the characteristic polynomial coefficients. Because the poles lambda_k are roots of that polynomial, multiplying A by the Vandermonde vector
v_k = [1, lambda_k, lambda_k^2, ..., lambda_k^(N-1)]gives backlambda_k * v_k, satisfyingA*v_k = lambda_k*v_kwhich again is the definition of an eigenvector.For each pole lambda_k, the corresponding eigenvector column is therefore:
So the matrix of eigenvectors V becomes:
Then the same matrix exponential calculation applies as it did before I made these changes:
Before
After
Summary