Add Butterworth zeta constructor for uniform damping ratio

docs/butterworth.md

@@ -118,6 +118,86 @@ keeping the denominator monic.
 The numerator has $b[0] = 1$ and all other $b[k] = 0$, giving unity
 high-frequency gain: $H(\infty) = 1$.
 
+## Uniform damping ratio constructor
+
+The second constructor overload accepts an explicit damping ratio $\zeta$.
+
+```cpp
+Butterworth<T, N, Method, FilterType>(cutoff_hz, sample_rate_hz, zeta)
+```
+
+Every complex-conjugate pole pair is placed at the same $\zeta$, giving each
+quadratic factor the form $s^2 + 2\zeta\omega_c s + \omega_c^2$. This differs from
+the classical constructor, where each pair has a distinct damping ratio determined
+by its position on the unit circle.
+
+### Pole placement
+
+All $M = \lfloor N/2 \rfloor$ complex pairs share the same pole location:
+
+$$p = \omega_c(-\zeta \pm j\sqrt{1 - \zeta^2})$$
+
+The behavior depends on $\zeta$. For $0 < \zeta < 1$ the poles are complex
+conjugates (underdamped). For $\zeta = 1$ the two roots coincide at
+$-\omega_c$ (critically damped). For $\zeta > 1$, the equivalent real-pole form
+$p = \omega_c(-\zeta \pm \sqrt{\zeta^2 - 1})$ applies, so the poles are two
+distinct negative real values (overdamped). All three cases produce a stable
+filter with real polynomial coefficients, and `butterworth_poly_coeffs_zeta`
+handles them identically since it operates only on the real quadratic factor
+$s^2 + 2\zeta\omega_c s + \omega_c^2$.
+
+For odd $N$ the remaining real pole is placed at $-\omega_c$, identical to the
+classical case (a real pole has no damping ratio to control).
+
+### Denominator polynomial
+
+Because the poles are known analytically as quadratic factors, the denominator
+is built by real arithmetic directly, without complex types. Starting from the
+constant polynomial 1, each conjugate pair multiplies in a quadratic factor.
+
+```text
+for pair = 0 .. M-1:
+    cur = 2*pair
+    for j = cur+2 downto 2:
+        poly[j] += 2*zeta*poly[j-1] + poly[j-2]
+    poly[1] += 2*zeta*poly[0]
+if N is odd:
+    for j = N downto 1:
+        poly[j] += poly[j-1]
+```
+
+For example, with $\zeta = 0.5$:
+
+$$N=2: \quad s^2 + s + 1$$
+
+$$N=3: \quad s^3 + 2s^2 + 2s + 1$$
+
+$$N=4: \quad s^4 + 2s^3 + 3s^2 + 2s + 1$$
+
+Cutoff scaling and numerator construction follow the same rules as the classical
+constructor (Steps 3 and 4 above).
+
+### Highpass variant
+
+The highpass transform is identical to the classical case. The denominator
+coefficients are reversed and scaled by $\omega_c^k$, and the numerator is
+$b[0] = 1$ with all other entries zero.
+
+### ZOH order restriction
+
+The zeta constructor does not supply a `FactoredTF`, so ZOH discretization falls
+back to generic eigendecomposition via `matrix_exp`. With uniform $\zeta$ and
+$N \geq 4$, every complex pair lands at the same pole location, producing repeated
+eigenvalues in the controllable canonical form matrix. QR-based eigendecomposition
+is ill-conditioned for defective matrices of this kind. The eigenvector matrix is
+singular, so `matrix_exp` returns wrong results. MatchedZ and Tustin are
+unaffected because they do not use eigendecomposition.
+
+A `static_assert` rejects `ZOH` with `N >= 4` at compile time. See
+[issue 56](https://github.com/MitchellThompkins/constfilt/issues/56) for
+background. The classical constructor is not affected because its poles are
+always distinct by construction.
+
 ## Test reference
 
 `tests/butterworth_reference.hpp` is regenerated from Octave's stock `butter()`

include/constfilt/analog_filter.hpp

@@ -8,6 +8,14 @@
 namespace constfilt
 {
 
+struct LowPass
+{
+};
+
+struct HighPass
+{
+};
+
 // Discretize an analog (continuous-time, s-domain) transfer function into a
 // digital Filter.
 //

include/constfilt/butterworth.hpp

@@ -8,14 +8,6 @@
 namespace constfilt
 {
 
-struct LowPass
-{
-};
-
-struct HighPass
-{
-};
-
 // Template parameters:
 //   T          - floating-point scalar type
 //   N          - filter order (>= 1)
@@ -40,6 +32,25 @@ class Butterworth
     {
     }
 
+    // Construct with uniform damping ratio zeta across all complex pole pairs.
+    // Routes through generic eigendecomposition (no FactoredTF) to avoid the
+    // Vandermonde singularity that arises when N>=4 produces identical pole
+    // pairs.
+    constexpr Butterworth(T cutoff_hz, T sample_rate_hz, T zeta)
+        : AnalogFilter<T, N, BoundMethod>(
+              compute_continuous_tf_zeta(cutoff_hz, zeta), sample_rate_hz,
+              make_tustin_tag(cutoff_hz, BoundMethod{}))
+    {
+        // Uniform zeta produces repeated complex eigenvalue pairs for N>=4.
+        // consteig's QR iteration cannot split a defective eigenvalue block,
+        // so matrix_exp returns wrong results. MatchedZ and Tustin are fine.
+        static_assert(
+            !is_zoh_tag<Method>::value || N <= 3u,
+            "Butterworth zeta constructor: ZOH is unreliable for N>=4 due to "
+            "repeated eigenvalues. Use MatchedZ or TustinPW instead. "
+            "See https://github.com/MitchellThompkins/constfilt/issues/56");
+    }
+
   private:
     // Computes the continuous-time Butterworth transfer function.
     static constexpr TransferFunction<T, N + 1u, N + 1u> compute_continuous_tf(
@@ -173,6 +184,72 @@ class Butterworth
             result[i] = poly[N - i].real;
         }
     }
+
+    // Normalized Butterworth denominator coefficients for uniform damping
+    // ratio. Each complex pair contributes a quadratic factor (s^2 + 2*zeta*s +
+    // 1); an odd-order real pole contributes (s + 1). Pure real arithmetic, no
+    // complex types. Fills result in descending power order: [1, ..., 1].
+    static constexpr void butterworth_poly_coeffs_zeta(T (&result)[N + 1u],
+                                                       T zeta)
+    {
+        // poly[i] holds the coefficient of s^i (ascending order).
+        T poly[N + 1u]{};
+        poly[0] = static_cast<T>(1);
+        for (consteig::Size pair = 0u; pair < N / 2u; ++pair)
+        {
+            const consteig::Size cur = 2u * pair;
+            for (consteig::Size j = cur + 2u; j > 1u; --j)
+            {
+                poly[j] +=
+                    static_cast<T>(2) * zeta * poly[j - 1u] + poly[j - 2u];
+            }
+            poly[1] += static_cast<T>(2) * zeta * poly[0];
+        }
+        if (N % 2u == 1u)
+        {
+            for (consteig::Size j = N; j > 0u; --j)
+            {
+                poly[j] += poly[j - 1u];
+            }
+        }
+        for (consteig::Size i = 0u; i <= N; ++i)
+        {
+            result[i] = poly[N - i];
+        }
+    }
+
+    static constexpr void continuous_tf_zeta(T wc, T zeta, T (&b)[N + 1u],
+                                             T (&a)[N + 1u], LowPass)
+    {
+        b[N] = gcem::pow(wc, static_cast<int>(N));
+        T p[N + 1u]{};
+        butterworth_poly_coeffs_zeta(p, zeta);
+        for (consteig::Size k = 0; k <= N; ++k)
+        {
+            a[k] = p[k] * gcem::pow(wc, static_cast<int>(k));
+        }
+    }
+
+    static constexpr void continuous_tf_zeta(T wc, T zeta, T (&b)[N + 1u],
+                                             T (&a)[N + 1u], HighPass)
+    {
+        b[0] = static_cast<T>(1);
+        T p[N + 1u]{};
+        butterworth_poly_coeffs_zeta(p, zeta);
+        for (consteig::Size k = 0; k <= N; ++k)
+        {
+            a[k] = p[N - k] * gcem::pow(wc, static_cast<int>(k));
+        }
+    }
+
+    static constexpr TransferFunction<T, N + 1u, N + 1u>
+    compute_continuous_tf_zeta(T cutoff_hz, T zeta)
+    {
+        const T wc = static_cast<T>(2) * static_cast<T>(GCEM_PI) * cutoff_hz;
+        TransferFunction<T, N + 1u, N + 1u> tf{};
+        continuous_tf_zeta(wc, zeta, tf.b, tf.a, FilterType{});
+        return tf;
+    }
 };
 
 // Convenience aliases for first-order RC-equivalent filters.

include/constfilt/discretize.hpp

@@ -61,6 +61,16 @@ template <> struct is_tustinpw_tag<TustinPW>
     static constexpr bool value = true;
 };
 
+template <typename M> struct is_zoh_tag
+{
+    static constexpr bool value = false;
+};
+
+template <> struct is_zoh_tag<ZOH>
+{
+    static constexpr bool value = true;
+};
+
 // Build the method tag from a cutoff frequency.
 // For TustinPWData<T>, fills in warp_omega = 2*pi*cutoff_hz.
 // For all other methods, returns a default-constructed tag (cutoff unused).

octave/generate_butterworth_tests.m

@@ -218,6 +218,166 @@ function emit_case(fid, sname, ord, fc, fs, b_d, a_d, y_step, y_imp, u_c, y_c)
               ord, fc, fs, b_tupw_hp, a_tupw_hp, y_step_tupw_hp, y_imp_tupw_hp, u_c, y_c_tupw_hp);
 end
 
+% =============================================================================
+% Uniform-zeta cases
+% =============================================================================
+%
+% Poles at wc*(-zeta +/- j*sqrt(1-zeta^2)) for each conjugate pair.
+% Avoids Octave's buttap (which uses classical Butterworth angles) so that the
+% Octave reference matches the C++ uniform-zeta constructor exactly.
+
+% There is no standard Octave signal package function for uniform-zeta
+% Butterworth design. buttap() produces classical Butterworth poles with
+% evenly spaced angles, which gives a different damping ratio per pair.
+% The poles must be constructed manually from zeta as below.
+function [b_s, a_s] = uniform_zeta_lp(N, wc, zeta)
+    omega = sqrt(1 - zeta^2);
+    poles = [];
+    for k = 1:floor(N/2)
+        poles = [poles; wc*(-zeta + 1j*omega); wc*(-zeta - 1j*omega)];
+    end
+    if mod(N,2) == 1; poles = [poles; -wc]; end
+    [b_s, a_s] = zp2tf([], poles, wc^N);
+end
+
+function [b_s, a_s] = uniform_zeta_hp(N, wc, zeta)
+    % HP poles via LP-to-HP transformation: wc^2 / p_lp.
+    % Since |p_lp| = wc and |p_lp|^2 = wc^2, this yields the complex conjugate
+    % of each LP pole, so re-derive directly.
+    omega = sqrt(1 - zeta^2);
+    poles_hp = [];
+    for k = 1:floor(N/2)
+        poles_hp = [poles_hp; wc*(-zeta - 1j*omega); wc*(-zeta + 1j*omega)];
+    end
+    if mod(N,2) == 1; poles_hp = [poles_hp; -wc]; end
+    [b_s, a_s] = zp2tf(zeros(N, 1), poles_hp, 1);
+end
+
+zeta_val = 0.5;
+zeta_tag = 'z50';
+
+zeta_lp_cases = {
+    {2, 100, 1000},
+    {3, 100, 1000},
+    {4, 100, 1000},
+    {5, 100, 1000},
+    {6, 100, 1000},
+    {7, 100, 1000},
+    {8, 100, 1000},
+};
+
+for ci = 1:length(zeta_lp_cases)
+    ord = zeta_lp_cases{ci}{1};
+    fc  = zeta_lp_cases{ci}{2};
+    fs  = zeta_lp_cases{ci}{3};
+    wc  = 2 * pi * fc;
+
+    [b_s, a_s] = uniform_zeta_lp(ord, wc, zeta_val);
+    sys_c = tf(b_s, a_s);
+    u = ones(1, STEP_LEN);
+    t_c = (0:CHIRP_LEN-1) / fs;
+    u_c = chirp(t_c, 0, t_c(end), fs/2);
+
+    % --- ZOH ---
+    sys_d = c2d(sys_c, 1.0/fs, 'zoh');
+    [b_d, a_d] = tfdata(sys_d, 'v');
+    while length(b_d) < length(a_d); b_d = [0, b_d]; end
+    y_step = filter(b_d, a_d, u);
+    y_imp  = filter(b_d, a_d, [1, zeros(1, STEP_LEN-1)]);
+    y_c    = filter(b_d, a_d, u_c);
+    emit_case(fid, sprintf('case_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_d, a_d, y_step, y_imp, u_c, y_c);
+
+    % --- Matched-Z ---
+    sys_mz = c2d(sys_c, 1.0/fs, 'matched');
+    [b_mz, a_mz] = tfdata(sys_mz, 'v');
+    while length(b_mz) < length(a_mz); b_mz = [0, b_mz]; end
+    y_step_mz = filter(b_mz, a_mz, u);
+    y_imp_mz  = filter(b_mz, a_mz, [1, zeros(1, STEP_LEN-1)]);
+    y_c_mz    = filter(b_mz, a_mz, u_c);
+    emit_case(fid, sprintf('case_mz_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_mz, a_mz, y_step_mz, y_imp_mz, u_c, y_c_mz);
+
+    % --- TustinNW ---
+    sys_tu = c2d(sys_c, 1.0/fs, 'tustin');
+    [b_tu, a_tu] = tfdata(sys_tu, 'v');
+    while length(b_tu) < length(a_tu); b_tu = [0, b_tu]; end
+    y_step_tu = filter(b_tu, a_tu, u);
+    y_imp_tu  = filter(b_tu, a_tu, [1, zeros(1, STEP_LEN-1)]);
+    y_c_tu    = filter(b_tu, a_tu, u_c);
+    emit_case(fid, sprintf('case_tu_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_tu, a_tu, y_step_tu, y_imp_tu, u_c, y_c_tu);
+
+    % --- TustinPW ---
+    sys_tupw = c2d(sys_c, 1.0/fs, 'prewarp', wc);
+    [b_tupw, a_tupw] = tfdata(sys_tupw, 'v');
+    while length(b_tupw) < length(a_tupw); b_tupw = [0, b_tupw]; end
+    y_step_tupw = filter(b_tupw, a_tupw, u);
+    y_imp_tupw  = filter(b_tupw, a_tupw, [1, zeros(1, STEP_LEN-1)]);
+    y_c_tupw    = filter(b_tupw, a_tupw, u_c);
+    emit_case(fid, sprintf('case_tupw_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_tupw, a_tupw, y_step_tupw, y_imp_tupw, u_c, y_c_tupw);
+end
+
+zeta_hp_cases = {
+    {2, 100, 1000},
+    {3, 100, 1000},
+    {4, 100, 1000},
+    {5, 100, 1000},
+    {6, 100, 1000},
+    {7, 100, 1000},
+    {8, 100, 1000},
+};
+
+for ci = 1:length(zeta_hp_cases)
+    ord = zeta_hp_cases{ci}{1};
+    fc  = zeta_hp_cases{ci}{2};
+    fs  = zeta_hp_cases{ci}{3};
+    wc  = 2 * pi * fc;
+
+    [b_s, a_s] = uniform_zeta_hp(ord, wc, zeta_val);
+    sys_c = tf(b_s, a_s);
+    u = ones(1, STEP_LEN);
+    t_c = (0:CHIRP_LEN-1) / fs;
+    u_c = chirp(t_c, 0, t_c(end), fs/2);
+
+    % --- ZOH ---
+    sys_d = c2d(sys_c, 1.0/fs, 'zoh');
+    [b_d, a_d] = tfdata(sys_d, 'v');
+    y_step = filter(b_d, a_d, u);
+    y_imp  = filter(b_d, a_d, [1, zeros(1, STEP_LEN-1)]);
+    y_c    = filter(b_d, a_d, u_c);
+    emit_case(fid, sprintf('case_hp_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_d, a_d, y_step, y_imp, u_c, y_c);
+
+    % --- Matched-Z ---
+    sys_mz = c2d(sys_c, 1.0/fs, 'matched');
+    [b_mz, a_mz] = tfdata(sys_mz, 'v');
+    y_step_mz = filter(b_mz, a_mz, u);
+    y_imp_mz  = filter(b_mz, a_mz, [1, zeros(1, STEP_LEN-1)]);
+    y_c_mz    = filter(b_mz, a_mz, u_c);
+    emit_case(fid, sprintf('case_mz_hp_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_mz, a_mz, y_step_mz, y_imp_mz, u_c, y_c_mz);
+
+    % --- TustinNW ---
+    sys_tu = c2d(sys_c, 1.0/fs, 'tustin');
+    [b_tu, a_tu] = tfdata(sys_tu, 'v');
+    y_step_tu = filter(b_tu, a_tu, u);
+    y_imp_tu  = filter(b_tu, a_tu, [1, zeros(1, STEP_LEN-1)]);
+    y_c_tu    = filter(b_tu, a_tu, u_c);
+    emit_case(fid, sprintf('case_tu_hp_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_tu, a_tu, y_step_tu, y_imp_tu, u_c, y_c_tu);
+
+    % --- TustinPW ---
+    sys_tupw = c2d(sys_c, 1.0/fs, 'prewarp', wc);
+    [b_tupw, a_tupw] = tfdata(sys_tupw, 'v');
+    y_step_tupw = filter(b_tupw, a_tupw, u);
+    y_imp_tupw  = filter(b_tupw, a_tupw, [1, zeros(1, STEP_LEN-1)]);
+    y_c_tupw    = filter(b_tupw, a_tupw, u_c);
+    emit_case(fid, sprintf('case_tupw_hp_zeta_%d_%dHz_%dHz_%s', ord, fc, fs, zeta_tag), ...
+              ord, fc, fs, b_tupw, a_tupw, y_step_tupw, y_imp_tupw, u_c, y_c_tupw);
+end
+
 fprintf(fid, '} // namespace bw_ref\n\n');
 fprintf(fid, '#endif // CONSTFILT_BUTTERWORTH_REFERENCE_HPP\n');
 fclose(fid);