From ab0f3d0c77360990cda79b1e021a14915abd635b Mon Sep 17 00:00:00 2001
From: VaiTon <eyadlorenzo@gmail.com>
Date: Thu, 11 Jan 2024 17:07:35 +0100
Subject: [PATCH] Allow bode and bodeplot to use a range of frequencies and a
freq to calculate the freq. vector
---
lib/ControlSystemsBase/src/freqresp.jl | 28 +-
lib/ControlSystemsBase/src/plotting.jl | 5 +
lib/ControlSystemsBase/src/utilities.jl | 52 ++-
lib/ControlSystemsBase/test/test_freqresp.jl | 367 ++++++++++---------
4 files changed, 255 insertions(+), 197 deletions(-)
diff --git a/lib/ControlSystemsBase/src/freqresp.jl b/lib/ControlSystemsBase/src/freqresp.jl
index ce8148e26..a578afd59 100644
--- a/lib/ControlSystemsBase/src/freqresp.jl
+++ b/lib/ControlSystemsBase/src/freqresp.jl
@@ -296,13 +296,27 @@ function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle
end
"""
- mag, phase, w = bode(sys[, w]; unwrap=true)
+ mag, phase, w = bode(sys; unwrap=true)
+ mag, phase, w = bode(sys, w; unwrap=true)
+ mag, phase, w = bode(sys, (w1, w2); freq=100, unwrap=true)
-Compute the magnitude and phase parts of the frequency response of system `sys`
-at frequencies `w`. See also [`bodeplot`](@ref)
+Compute the magnitude and phase parts of the frequency response of system `sys`.
-`mag` and `phase` has size `(ny, nu, length(w))`.
-If `unwrap` is true (default), the function `unwrap!` will be applied to the phase angles. This procedure is costly and can be avoided if the unwrapping is not required.
+If the frequency vector is not provided, it is automatically generated based on
+the poles and zeros of the system.
+
+If a frequency vector is provided, the magnitude and phase are computed at the
+frequencies in the vector.
+
+If a frequency range is provided, the magnitude and phase are computed at `freq`
+logarithmically spaced points per decade between `w1` and `w2`.
+
+If `unwrap` is true (default), the function `unwrap!` will be applied to the phase angles.
+This procedure is costly and can be avoided if the unwrapping is not required.
+
+`mag` and `phase` have size `(ny, nu, length(w))`.
+
+To plot a Bode diagram, see [`bodeplot`](@ref).
For higher performance, see the function [`bodemag!`](@ref) that computes the magnitude only.
"""
@@ -313,8 +327,10 @@ For higher performance, see the function [`bodemag!`](@ref) that computes the ma
@. angles = rad2deg(angles)
return abs.(resp), angles, w
end
-@autovec (1, 2) bode(sys::LTISystem) = bode(sys, _default_freq_vector(sys, Val{:bode}()))
+@autovec (1, 2) bode(sys::LTISystem, w::Tuple{<:Real,<:Real}, freq=100; unwrap=true) =
+ bode(sys, _tuple_to_vec(w, freq), unwrap=unwrap)
+@autovec (1, 2) bode(sys::LTISystem) = bode(sys, _default_freq_vector(sys, Val{:bode}()))
# Performance difference between bode and bodemag for tf. Note how expensive the phase unwrapping is.
# using ControlSystemsBase
# G = tf(ssrand(2,2,5))
diff --git a/lib/ControlSystemsBase/src/plotting.jl b/lib/ControlSystemsBase/src/plotting.jl
index 667567266..5deacf65d 100644
--- a/lib/ControlSystemsBase/src/plotting.jl
+++ b/lib/ControlSystemsBase/src/plotting.jl
@@ -231,6 +231,11 @@ _processfreqplot(plottype, system::LTISystem, args...) =
_processfreqplot(plottype, [system], args...)
# Catch when system is not vector, with and without frequency input
+_processfreqplot(plottype, systems::AbstractVector{<:LTISystem},
+ t::Tuple{Real,Real}, freq::Real) =
+ _processfreqplot(plottype, systems, _tuple_to_vec(t, freq))
+
+
# Catch correct form
function _processfreqplot(plottype, systems::AbstractVector{<:LTISystem},
w = _default_freq_vector(systems, plottype))
diff --git a/lib/ControlSystemsBase/src/utilities.jl b/lib/ControlSystemsBase/src/utilities.jl
index b9166a344..238fe42ed 100644
--- a/lib/ControlSystemsBase/src/utilities.jl
+++ b/lib/ControlSystemsBase/src/utilities.jl
@@ -1,15 +1,15 @@
# Are only really needed for cases when we accept general LTISystem
# we should either use them consistently, with a good definition, or remove them
basetype(sys::LTISystem) = basetype(typeof(sys))
-basetype(::Type{T}) where T <: LTISystem = T.name.wrapper
-numeric_type(::Type{SisoRational{T}}) where T = T
-numeric_type(::Type{<:SisoZpk{T}}) where T = T
+basetype(::Type{T}) where {T<:LTISystem} = T.name.wrapper
+numeric_type(::Type{SisoRational{T}}) where {T} = T
+numeric_type(::Type{<:SisoZpk{T}}) where {T} = T
numeric_type(sys::SisoTf) = numeric_type(typeof(sys))
numeric_type(::Type{TransferFunction{TE,S}}) where {TE,S} = numeric_type(S)
numeric_type(::Type{<:StateSpace{TE,T}}) where {TE,T} = T
numeric_type(::Type{<:HeteroStateSpace{TE,AT}}) where {TE,AT} = eltype(AT)
-numeric_type(::Type{<:LFTSystem{<:Any, T}}) where {T} = T
+numeric_type(::Type{<:LFTSystem{<:Any,T}}) where {T} = T
numeric_type(sys::AbstractStateSpace) = eltype(sys.A)
numeric_type(sys::LTISystem) = numeric_type(typeof(sys))
@@ -18,12 +18,12 @@ to_matrix(T, A::AbstractVector) = Matrix{T}(reshape(A, length(A), 1))
to_matrix(T, A::AbstractMatrix) = convert(Matrix{T}, A) # Fallback
to_matrix(T, A::Number) = fill(T(A), 1, 1)
# Handle Adjoint Matrices
-to_matrix(T, A::Adjoint{R, MT}) where {R<:Number, MT<:AbstractMatrix} = to_matrix(T, MT(A))
+to_matrix(T, A::Adjoint{R,MT}) where {R<:Number,MT<:AbstractMatrix} = to_matrix(T, MT(A))
to_abstract_matrix(A::AbstractMatrix) = A
function to_abstract_matrix(A::AbstractArray)
try
- return convert(AbstractMatrix,A)
+ return convert(AbstractMatrix, A)
catch
@warn "Could not convert $(typeof(A)) to `AbstractMatrix`. A HeteroStateSpace must consist of AbstractMatrix."
rethrow()
@@ -35,7 +35,7 @@ to_abstract_matrix(A::Number) = fill(A, 1, 1)
# Do no sorting of eigenvalues
eigvalsnosort(args...; kwargs...) = eigvals(args...; sortby=nothing, kwargs...)
-eigvalsnosort(A::StaticArraysCore.StaticMatrix; kwargs...) = eigvalsnosort(Matrix(A); kwargs...)
+eigvalsnosort(A::StaticArraysCore.StaticMatrix; kwargs...) = eigvalsnosort(Matrix(A); kwargs...)
roots(args...; kwargs...) = Polynomials.roots(args...; sortby=nothing, kwargs...)
issemiposdef(A) = ishermitian(A) && minimum(real.(eigvals(A))) >= 0
@@ -44,7 +44,7 @@ issemiposdef(A::UniformScaling) = real(A.λ) >= 0
""" f = printpolyfun(var)
`fun` Prints polynomial in descending order, with variable `var`
"""
-printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printpoly(io, Polynomial(p.coeffs, var), mimetype, descending_powers=true, mulsymbol="")
+printpolyfun(var) = (io, p, mimetype=MIME"text/plain"()) -> Polynomials.printpoly(io, Polynomial(p.coeffs, var), mimetype, descending_powers=true, mulsymbol="")
# NOTE: Tolerances for checking real-ness removed, shouldn't happen from LAPACK?
# TODO: This doesn't play too well with dual numbers..
@@ -53,14 +53,14 @@ printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printp
# This function rely on that the every complex roots is followed by its exact conjugate,
# and that the first complex root in each pair has positive imaginary part. This format is always
# returned by LAPACK routines for eigenvalues.
-function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
+function roots2real_poly_factors(roots::Vector{cT}) where {cT<:Number}
T = real(cT)
poly_factors = Vector{Polynomial{T,:x}}()
- for k=1:length(roots)
+ for k = 1:length(roots)
r = roots[k]
if isreal(r)
- push!(poly_factors,Polynomial{T,:x}([-real(r),1]))
+ push!(poly_factors, Polynomial{T,:x}([-real(r), 1]))
else
if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
continue
@@ -70,7 +70,7 @@ function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
throw(ArgumentError("Found pole without matching conjugate."))
end
- push!(poly_factors,Polynomial{T,:x}([real(r)^2+imag(r)^2, -2*real(r), 1]))
+ push!(poly_factors, Polynomial{T,:x}([real(r)^2 + imag(r)^2, -2 * real(r), 1]))
# k += 1 # Skip one iteration in the loop
end
end
@@ -78,7 +78,7 @@ function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
return poly_factors
end
# This function should handle both Complex as well as symbolic types
-function roots2poly_factors(roots::Vector{T}) where T <: Number
+function roots2poly_factors(roots::Vector{T}) where {T<:Number}
return Polynomial{T,:x}[Polynomial{T,:x}([-r, 1]) for r in roots]
end
@@ -96,7 +96,7 @@ function _fix_conjugate_pairs!(v::AbstractVector{<:Complex})
# Do nothing
else
if isapprox(v[k], conj(v[k+1]), rtol=1e-15)
- z = (v[k] + conj(v[k+1]))/2
+ z = (v[k] + conj(v[k+1])) / 2
v[k] = z
v[k+1] = conj(z)
k += 1
@@ -114,14 +114,14 @@ poly2vec(p::Polynomial) = p.coeffs[1:end]
function unwrap!(M::Array, dim=1)
- alldims(i) = ntuple(n->n==dim ? i : (1:size(M,n)), ndims(M))
+ alldims(i) = ntuple(n -> n == dim ? i : (1:size(M, n)), ndims(M))
π2 = eltype(M)(2π)
for i = 2:size(M, dim)
#This is a copy of slicedim from the JuliaLang but enables us to write to it
#The code (with dim=1) is equivalent to
# d = M[i,:,:,...,:] - M[i-1,:,...,:]
# M[i,:,:,...,:] -= floor((d+π) / (2π)) * 2π
- d = M[alldims(i)...] - M[alldims(i-1)...]
+ d = M[alldims(i)...] - M[alldims(i - 1)...]
M[alldims(i)...] -= @. floor((d + π) / π2) * π2
end
return M
@@ -196,7 +196,7 @@ macro autovec(indices, f)
end
end
result = $(esc(fname))($(argnames...);
- $(map(x->Expr(:(=), esc(x), esc(x)), kwargnames)...))
+ $(map(x -> Expr(:(=), esc(x), esc(x)), kwargnames)...))
return $return_exp
end
end
@@ -205,3 +205,21 @@ end
function extractvarname(a)
typeof(a) == Symbol ? a : extractvarname(a.args[1])
end
+
+"""
+ _tuple_to_range(t, freq)
+
+Helper function to convert a tuple of frequencies to a range of frequencies.
+The vector is created by taking `freq` points logarithmically spaced every decade.
+
+# Example
+
+ _tuple_to_range((1e-1, 1e5), 100)
+
+will create a range from 1e-1 to 1e5 with 100 points logarithmically spaced every decade.
+"""
+function _tuple_to_vec(t::Tuple{Real,Real}, freq::L) where {L<:Real}
+ w = log10.(t)
+ n_points = round(Int, freq * (w[2] - w[1]))
+ exp10.(LinRange(w[1], w[2], n_points))
+end
\ No newline at end of file
diff --git a/lib/ControlSystemsBase/test/test_freqresp.jl b/lib/ControlSystemsBase/test/test_freqresp.jl
index 89a385e75..37075b59f 100644
--- a/lib/ControlSystemsBase/test/test_freqresp.jl
+++ b/lib/ControlSystemsBase/test/test_freqresp.jl
@@ -1,178 +1,197 @@
@testset "test_freqresp" begin
-## EVALFR ##
-H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])]
-G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
- [0 3 0 0; -2 0 0 1], [0 0; 1 0])
-
-
-
-@test evalfr(H, -6) == [0.0 -0.45; 5.0 Inf]
-@test evalfr(H, -5) == [0.0 -0.5; Inf 2.0]
-@test evalfr(H, -1) == [0.0 -0.3; 0.0 0.4]
-@test evalfr(H, 0) ≈ [0.0 0.0; 0.2 1/3]
-@inferred evalfr(H, 0)
-
-@test (@test_logs (:warn, "Got exception SingularException(4), returning Inf") evalfr(G, -6)) == [Inf Inf; Inf Inf]
-@test (@test_logs (:warn, "Got exception SingularException(1), returning Inf") evalfr(G, -5)) == [Inf Inf; Inf Inf]
-@test evalfr(G, -1) ≈ [0.0 -0.3; 0.0 0.4]
-@test evalfr(G, 0) ≈ [0.0 0.0; 0.2 1/3]
-
-## Constant system
-w = exp10.(range(-1, stop=1, length=50))
-
-sys1 = ss(1)
-sys1s = HeteroStateSpace(sparse([1]))
-G1 = tf(1)
-H1 = zpk(1)
-resp1 = ones(ComplexF64, 1, 1, length(w))
-
-@test evalfr(sys1, im*w[1]) == fill(resp1[1], 1, 1)
-@test evalfr(G1, im*w[1]) == fill(resp1[1], 1, 1)
-@test evalfr(H1, im*w[1]) == fill(resp1[1], 1, 1)
-
-@test freqresp(sys1, w) == resp1
-@test freqresp(sys1s, w) == resp1
-@test freqresp(G1, w) == resp1
-@test freqresp(H1, w) == resp1
-@inferred freqresp(sys1, w)
-@inferred freqresp(G1, w)
-@inferred freqresp(H1, w)
-
-for G in [2, 2I, 2I(2), randn(2,2)]
- @test dcgain(G) == G
- !(G isa UniformScaling) && @test freqresp(G, [1,2]) == cat(G,G, dims=3) # I does not have size
- @test freqresp(G, 1) == G
- @test evalfr(G, 1im) == G
-end
-
-## First order system
-
-sys2 = ss(-1, 1, 1, 1)
-sys2s = HeteroStateSpace(sparse([-1;;]),sparse([1;;]),sparse([1;;]),sparse([1;;])) # test that freqresp works for matrix types that don't support hessenberg
-G2 = tf([1, 2], [1,1])
-H2 = zpk([-2], [-1.0], 1.0)
-resp2 = reshape((im*w .+ 2)./(im*w .+ 1), 1, 1, length(w))
-
-@test evalfr(sys2, im*w[1]) ≈ fill(resp2[1], 1, 1)
-@test evalfr(G2, im*w[1]) == fill(resp2[1], 1, 1)
-@test evalfr(H2, im*w[1]) == fill(resp2[1], 1, 1)
-
-@test freqresp(sys2, w) ≈ resp2 rtol=1e-15
-@test freqresp(sys2s, w) ≈ resp2 rtol=1e-15
-@test freqresp(G2, w) == resp2
-@test freqresp(H2, w) == resp2
-
-@inferred freqresp(sys2, w)
-@inferred freqresp(sys2s, w)
-@inferred freqresp(G2, w)
-@inferred freqresp(H2, w)
-@test (@allocated freqresp(sys2, w)) < 1.2*56672 # allow 20% increase due to compiler variations
-@test (@allocated freqresp(G2, w)) < 1.2*976 # allow 20% increase due to compiler variations
-
-## Complex-coefficient system
-sys3 = ss(-1+im, 1, (1-im), (1-im))
-G3 = (1-im)*tf([1,2-im], [1,1-im])
-H3 = zpk([-2+im], [-1+im], (1-im))
-resp3 = reshape((1 .- im)*(2 .+ im*(w .- 1))./(1 .+ im*(w .- 1)), 1, 1, length(w))
-
-@test evalfr(sys3, im*w[1]) ≈ fill(resp3[1], 1, 1) rtol=1e-16
-@test evalfr(G3, im*w[1]) == fill(resp3[1], 1, 1)
-@test evalfr(H3, im*w[1]) == fill(resp3[1], 1, 1)
-
-@test freqresp(sys3, w) ≈ resp3 rtol=1e-15
-@test freqresp(G3, w) ≈ resp3 rtol=1e-16
-@test freqresp(H3, w) == resp3
-
-
-
-
-## Shortcut notation for evalfr ##
-F = tf([1],[1,0.5],-1)
-omega = 2
-z = 0.5(1+im)
-# This test is not correct if we dont have a sample time
-# @test_throws F(omega,true)[1] == 1/(exp(-im*2)+0.5)
-@test F(z,false)[1] == 1/(z+0.5)
-@test_throws ErrorException F(z,true)
-
-F = [tf([1],[1,0.5],2.0) 3*tf([1],[1,0.5],2.0)]
-omegas = [1,2]
-z = 0.5(1+im)
-@test F(omegas[1],true) ≈ [1 3].*1/(exp(im*2)+0.5) atol=1e-14
-@test F(omegas[2],true) == [1 3].*1/(exp(2*im*2)+0.5)
-@test F(omegas,true) ≈ [k/(exp(omega*im*2)+0.5) for o=1:1, k=[1,3], omega=omegas] atol=1e-14
-@test F(omegas,false) ≈ [k/(omega+0.5) for o=1:1, k=[1,3], omega=omegas] atol=1e-14
-@test F(z,false)[1] == 1/(z+0.5)
-@test_throws ErrorException F(z,true)
-
-## Test bode, nyquist and sigma
-sys = [tf([1,-1], [1,1,1]) 0; 0 tf([1],[1,1])]
-f(s) = [(s-1)./(s.^2+s+1) 0; 0 1 ./(1+s)]
-ws = exp10.(range(-2, stop=2, length=50))
-resp = Array{ComplexF64}(undef,2,2,50)
-for (i,w) in enumerate(ws)
- resp[:,:,i] = f(im*w)
-end
-
-# A little test helper
-Base.isapprox(t1::Tuple, t2::Tuple) = t1[1] ≈ t2[1] && t1[2] ≈ t2[2]
-
-@test bode(sys, ws)[1:2] ≈ (abs.(resp), rad2deg.(angle.(resp)))
-workspace = BodemagWorkspace(sys, ws)
-@test bode(sys, ws)[1] == bodemag!(workspace, sys, ws)
-@test nyquist(sys, ws)[1:2] ≈ (real(resp), imag(resp))
-sigs = Array{Float64}(undef, 2, 50)
-for i in eachindex(ws)
- sigs[:, i] = svdvals(resp[:,:,i])
-end
-@test sigma(sys, ws)[1] ≈ sigs
-
-
-## Test sigma when either ny = 1 or nu = 1
-sys = ssrand(1,2,1)
-resp = freqresp(sys, ws)
-sigs = Array{Float64}(undef, 1, 50)
-for i in eachindex(ws)
- sigs[:, i] = svdvals(resp[:,:,i])
-end
-@test sigma(sys, ws)[1] ≈ sigs
-
-sys = ssrand(2,1,1)
-resp = freqresp(sys, ws)
-sigs = Array{Float64}(undef, 1, 50)
-for i in eachindex(ws)
- sigs[:, i] = svdvals(resp[:,:,i])
-end
-@test sigma(sys, ws)[1] ≈ sigs
-
-
-# test unwrap
-P = tf(1, [1, 1, 0]) * tf(1, [1, 1])
-w = exp10.(LinRange(-2, 4, 200))
-mag, ph = bode(P, w)
-@test ph[end] ≈ -(180 + 90) rtol = 1e-2
-
-# test explicit size and type constructor
-ws2 = BodemagWorkspace{Float32}(2,2,length(ws))
-@test size(ws2.mag) == (2,2,length(ws))
-@test size(ws2.R) == (2,2,length(ws))
-@test ws2.mag isa Array{Float32, 3}
-@test ws2.R isa Array{Complex{Float32}, 3}
-
-#Test default freq vector contains at least half a decade more than all features
-p, z = 100, 1/1000
-sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
-mag, phase, ws2 = bode(sys2)
-@test maximum(ws2) >= 2max(p,z)
-@test minimum(ws2) <= 0.5min(p,z)
-
-#Test default freq vector not too small
-p, z = 1, 1.2
-sys2 = [tf([1],[1/p,1]) tf([1/z, 2/z, 1],[1])]
-mag, mag, ws2 = bode(sys2)
-@test maximum(ws2) >= 5max(p,z)
-@test minimum(ws2) <= 0.2min(p,z)
-@test length(ws2) > 100
+ ## EVALFR ##
+ @testset "evalfr" begin
+ H = [tf(0) tf([3, 0], [1, 1, 10]); tf([1, 1], [1, 5]) tf([2], [1, 6])]
+ G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
+ [0 3 0 0; -2 0 0 1], [0 0; 1 0])
+
+ @test evalfr(H, -6) == [0.0 -0.45; 5.0 Inf]
+ @test evalfr(H, -5) == [0.0 -0.5; Inf 2.0]
+ @test evalfr(H, -1) == [0.0 -0.3; 0.0 0.4]
+ @test evalfr(H, 0) ≈ [0.0 0.0; 0.2 1/3]
+ @inferred evalfr(H, 0)
+
+ @test (@test_logs (:warn, "Got exception SingularException(4), returning Inf") evalfr(G, -6)) == [Inf Inf; Inf Inf]
+ @test (@test_logs (:warn, "Got exception SingularException(1), returning Inf") evalfr(G, -5)) == [Inf Inf; Inf Inf]
+ @test evalfr(G, -1) ≈ [0.0 -0.3; 0.0 0.4]
+ @test evalfr(G, 0) ≈ [0.0 0.0; 0.2 1/3]
+ end
+
+ w = exp10.(range(-1, stop=1, length=50))
+
+ @testset "Constant system" begin
+
+ sys1 = ss(1)
+ sys1s = HeteroStateSpace(sparse([1]))
+ G1 = tf(1)
+ H1 = zpk(1)
+ resp1 = ones(ComplexF64, 1, 1, length(w))
+
+ @test evalfr(sys1, im * w[1]) == fill(resp1[1], 1, 1)
+ @test evalfr(G1, im * w[1]) == fill(resp1[1], 1, 1)
+ @test evalfr(H1, im * w[1]) == fill(resp1[1], 1, 1)
+
+ @test freqresp(sys1, w) == resp1
+ @test freqresp(sys1s, w) == resp1
+ @test freqresp(G1, w) == resp1
+ @test freqresp(H1, w) == resp1
+ @inferred freqresp(sys1, w)
+ @inferred freqresp(G1, w)
+ @inferred freqresp(H1, w)
+
+ for G in [2, 2I, 2I(2), randn(2, 2)]
+ @test dcgain(G) == G
+ !(G isa UniformScaling) && @test freqresp(G, [1, 2]) == cat(G, G, dims=3) # I does not have size
+ @test freqresp(G, 1) == G
+ @test evalfr(G, 1im) == G
+ end
+
+ end
+
+ @testset "First order system" begin
+
+ sys2 = ss(-1, 1, 1, 1)
+ sys2s = HeteroStateSpace(sparse([-1;;]), sparse([1;;]), sparse([1;;]), sparse([1;;])) # test that freqresp works for matrix types that don't support hessenberg
+ G2 = tf([1, 2], [1, 1])
+ H2 = zpk([-2], [-1.0], 1.0)
+ resp2 = reshape((im * w .+ 2) ./ (im * w .+ 1), 1, 1, length(w))
+
+ @test evalfr(sys2, im * w[1]) ≈ fill(resp2[1], 1, 1)
+ @test evalfr(G2, im * w[1]) == fill(resp2[1], 1, 1)
+ @test evalfr(H2, im * w[1]) == fill(resp2[1], 1, 1)
+
+ @test freqresp(sys2, w) ≈ resp2 rtol = 1e-15
+ @test freqresp(sys2s, w) ≈ resp2 rtol = 1e-15
+ @test freqresp(G2, w) == resp2
+ @test freqresp(H2, w) == resp2
+
+ @inferred freqresp(sys2, w)
+ @inferred freqresp(sys2s, w)
+ @inferred freqresp(G2, w)
+ @inferred freqresp(H2, w)
+ @test (@allocated freqresp(sys2, w)) < 1.2 * 56672 # allow 20% increase due to compiler variations
+ @test (@allocated freqresp(G2, w)) < 1.2 * 976 # allow 20% increase due to compiler variations
+
+ end
+
+ ## Complex-coefficient system
+ sys3 = ss(-1 + im, 1, (1 - im), (1 - im))
+ G3 = (1 - im) * tf([1, 2 - im], [1, 1 - im])
+ H3 = zpk([-2 + im], [-1 + im], (1 - im))
+ resp3 = reshape((1 .- im) * (2 .+ im * (w .- 1)) ./ (1 .+ im * (w .- 1)), 1, 1, length(w))
+
+ @test evalfr(sys3, im * w[1]) ≈ fill(resp3[1], 1, 1) rtol = 1e-16
+ @test evalfr(G3, im * w[1]) == fill(resp3[1], 1, 1)
+ @test evalfr(H3, im * w[1]) == fill(resp3[1], 1, 1)
+
+ @test freqresp(sys3, w) ≈ resp3 rtol = 1e-15
+ @test freqresp(G3, w) ≈ resp3 rtol = 1e-16
+ @test freqresp(H3, w) == resp3
+
+
+ ## Shortcut notation for evalfr ##
+ F = tf([1], [1, 0.5], -1)
+ omega = 2
+ z = 0.5(1 + im)
+ # This test is not correct if we dont have a sample time
+ # @test_throws F(omega,true)[1] == 1/(exp(-im*2)+0.5)
+ @test F(z, false)[1] == 1 / (z + 0.5)
+ @test_throws ErrorException F(z, true)
+
+ F = [tf([1], [1, 0.5], 2.0) 3 * tf([1], [1, 0.5], 2.0)]
+ omegas = [1, 2]
+ z = 0.5(1 + im)
+ @test F(omegas[1], true) ≈ [1 3] .* 1 / (exp(im * 2) + 0.5) atol = 1e-14
+ @test F(omegas[2], true) == [1 3] .* 1 / (exp(2 * im * 2) + 0.5)
+ @test F(omegas, true) ≈ [k / (exp(omega * im * 2) + 0.5) for o = 1:1, k = [1, 3], omega = omegas] atol = 1e-14
+ @test F(omegas, false) ≈ [k / (omega + 0.5) for o = 1:1, k = [1, 3], omega = omegas] atol = 1e-14
+ @test F(z, false)[1] == 1 / (z + 0.5)
+ @test_throws ErrorException F(z, true)
+
+ @testset "Test bode, nyquist and sigma" begin
+ # A little test helper
+ Base.isapprox(t1::Tuple, t2::Tuple) = t1[1] ≈ t2[1] && t1[2] ≈ t2[2]
+
+ sys = [tf([1, -1], [1, 1, 1]) 0; 0 tf([1], [1, 1])]
+ f(s) = [(s-1)./(s .^ 2+s+1) 0; 0 1 ./ (1+s)]
+
+ freq = 20
+ start_dec = 10^-2
+ stop_dec = 10^2
+
+ npoints = round(Int, (log10(stop_dec) - log10(start_dec)) * freq)
+ ws = exp10.(LinRange(log10(start_dec), log10(stop_dec), npoints))
+
+ resp = Array{ComplexF64}(undef, 2, 2, npoints)
+ for (i, w) in enumerate(ws)
+ resp[:, :, i] = f(im * w)
+ end
+
+ @testset "bode" begin
+ @test bode(sys, (start_dec, stop_dec), freq) ≈ (abs.(resp), rad2deg.(angle.(resp)), ws)
+ @test bode(sys, ws)[1:2] ≈ (abs.(resp), rad2deg.(angle.(resp)))
+
+ workspace = BodemagWorkspace(sys, ws)
+ @test bode(sys, ws)[1] == bodemag!(workspace, sys, ws)
+ end
+
+ @test nyquist(sys, ws)[1:2] ≈ (real(resp), imag(resp))
+
+ sigs = Array{Float64}(undef, 2, npoints)
+ for i in eachindex(ws)
+ sigs[:, i] = svdvals(resp[:, :, i])
+ end
+ @test sigma(sys, ws)[1] ≈ sigs
+
+
+ ## Test sigma when either ny = 1 or nu = 1
+ sys = ssrand(1, 2, 1)
+ resp = freqresp(sys, ws)
+ sigs = Array{Float64}(undef, 1, npoints)
+ for i in eachindex(ws)
+ sigs[:, i] = svdvals(resp[:, :, i])
+ end
+ @test sigma(sys, ws)[1] ≈ sigs
+
+ sys = ssrand(2, 1, 1)
+ resp = freqresp(sys, ws)
+ sigs = Array{Float64}(undef, 1, npoints)
+ for i in eachindex(ws)
+ sigs[:, i] = svdvals(resp[:, :, i])
+ end
+ @test sigma(sys, ws)[1] ≈ sigs
+
+ end
+
+ # test unwrap
+ @testset "unwrap" begin
+ P = tf(1, [1, 1, 0]) * tf(1, [1, 1])
+ w = exp10.(LinRange(-2, 4, 200))
+ mag, ph = bode(P, w)
+ @test ph[end] ≈ -(180 + 90) rtol = 1e-2
+ end
+
+ # test explicit size and type constructor
+ ws2 = BodemagWorkspace{Float32}(2, 2, length(ws))
+ @test size(ws2.mag) == (2, 2, length(ws))
+ @test size(ws2.R) == (2, 2, length(ws))
+ @test ws2.mag isa Array{Float32,3}
+ @test ws2.R isa Array{Complex{Float32},3}
+
+ #Test default freq vector contains at least half a decade more than all features
+ p, z = 100, 1 / 1000
+ sys2 = [tf([1], [1 / p, 1]) tf([1 / z, 2 / z, 1], [1])]
+ mag, phase, ws2 = bode(sys2)
+ @test maximum(ws2) >= 2max(p, z)
+ @test minimum(ws2) <= 0.5min(p, z)
+
+ #Test default freq vector not too small
+ p, z = 1, 1.2
+ sys2 = [tf([1], [1 / p, 1]) tf([1 / z, 2 / z, 1], [1])]
+ mag, mag, ws2 = bode(sys2)
+ @test maximum(ws2) >= 5max(p, z)
+ @test minimum(ws2) <= 0.2min(p, z)
+ @test length(ws2) > 100
end