Allow bode and bodeplot to use a range of frequencies and a freq to calculate the freq. vector

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))

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))

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,15 +70,15 @@ 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
 
     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