diff --git a/src/weights.jl b/src/weights.jl index f625855a..e968c78d 100644 --- a/src/weights.jl +++ b/src/weights.jl @@ -6,7 +6,7 @@ Create a weighting function that goes from gain `low` at zero frequency, through # Arguments: - `low`: A number specifying the DC gain -- `mid`: A number specifying the frequency at which the gain is 1, or a tuple `(freq, gain)`. If `gain_mid` is not specified, the geometric mean of `high` and `low` is used. +- `mid`: A number specifying the frequency at which the gain is 1, or a tuple `(freq, gain)`. If `gain_mid` is not specified, it defaults to `1` whenever `low` and `high` straddle 1 (the usual case), and to the geometric mean `√(low*high)` otherwise. - `high`: A number specifying the gain at ∞ ```@example @@ -18,7 +18,14 @@ vline!([5], l=(:black, :dash), primary=false) ``` """ function makeweight(low, mid::Number, high) - makeweight(low, (mid, high < 1 ? √(high*low) : 1), high) + # Default `gain_mid` must lie strictly between `low` and `high` for the + # downstream formula to produce real coefficients. When `low` and `high` + # straddle 1 the historical default `gain_mid = 1` satisfies that and is + # preserved (back-compat). When both are on the same side of 1, fall back + # to the geometric mean `√(low*high)` so the previously-broken case where + # both `low > 1` and `high > 1` no longer produces complex/NaN poles. + gain_mid = (low <= 1 && high >= 1) ? one(promote_type(typeof(low), typeof(high))) : √(high*low) + makeweight(low, (mid, gain_mid), high) end diff --git a/test/test_weights.jl b/test/test_weights.jl index 865ae329..24061de0 100644 --- a/test/test_weights.jl +++ b/test/test_weights.jl @@ -12,6 +12,25 @@ w = gain_and_delay_uncertainty(1, 2, 1) w = makeweight(0.1, 1, 2) @test dcgain(w)[] ≈ 0.1 @test evalfr(w, 10000im)[] ≈ 2 atol=1e-3 +# `low` and `high` straddle 1, so the default gain_mid is 1 (preserved back-compat). +@test abs(evalfr(w, 1im)[]) ≈ 1 atol=1e-6 + +# Regression: previously, when both `low > 1` and `high > 1`, the default gain_mid was +# hard-coded to 1 which violates the formula's `low < mag < high` precondition and +# produced complex/NaN poles. Now falls back to √(low*high), which lies between them. +w_both_gt_1 = makeweight(2, 1, 3) +@test all(isreal, denvec(tf(w_both_gt_1))[1]) +@test all(isreal, numvec(tf(w_both_gt_1))[1]) +@test dcgain(w_both_gt_1)[] ≈ 2 +@test evalfr(w_both_gt_1, 10000im)[] ≈ 3 atol=1e-3 +@test abs(evalfr(w_both_gt_1, 1im)[]) ≈ √(2*3) atol=1e-6 + +# Symmetric same-side case: both below 1. Default gain_mid = √(low*high). +w_both_lt_1 = makeweight(0.5, 1, 0.1) +@test all(isreal, denvec(tf(w_both_lt_1))[1]) +@test dcgain(w_both_lt_1)[] ≈ 0.5 +@test evalfr(w_both_lt_1, 10000im)[] ≈ 0.1 atol=1e-3 +@test abs(evalfr(w_both_lt_1, 1im)[]) ≈ √(0.5*0.1) atol=1e-6 w = neglected_lag(1)