Skip to content

Commit

Permalink
use map and comprehensions
Browse files Browse the repository at this point in the history
  • Loading branch information
baggepinnen committed Dec 11, 2018
1 parent 666e41b commit 08690a7
Showing 1 changed file with 10 additions and 39 deletions.
49 changes: 10 additions & 39 deletions src/freqresp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,22 +13,11 @@ function freqresp(sys::LTISystem, w_vec::AbstractVector{S}) where {S<:Real}
else
s_vec = im*w_vec
end
Tsys = numeric_type(sys)
T = promote_type(typeof(zero(Tsys)/norm(one(Tsys))), ComplexF32, S)
sys_fr = Array{T}(undef, length(w_vec), noutputs(sys), ninputs(sys))

if isa(sys, StateSpace)
sys = _preprocess_for_freqresp(sys)
end


for i=1:length(s_vec)
# TODO : This doesn't actually take advantage of Hessenberg structure
# for statespace version.
sys_fr[i, :, :] .= evalfr(sys, s_vec[i])
end

return sys_fr
ny,nu = noutputs(sys), ninputs(sys)
[evalfr(sys[i,j], s)[] for s in s_vec, i in 1:ny, j in 1:nu]
end

# Implements algorithm found in:
Expand Down Expand Up @@ -74,14 +63,7 @@ function evalfr(sys::StateSpace{T0}, s::Number) where {T0}
end

function evalfr(G::TransferFunction{<:SisoTf{T0}}, s::Number) where {T0}
T = promote_type(T0, typeof(one(T0)*one(typeof(s))/(one(T0)*one(typeof(s)))))
fr = Array{T}(undef, size(G))
for j = 1:ninputs(G)
for i = 1:noutputs(G)
fr[i, j] = evalfr(G.matrix[i, j], s)
end
end
return fr
map(m -> evalfr(m,s), G.matrix)
end

"""
Expand Down Expand Up @@ -109,11 +91,8 @@ function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle
@assert !iscontinuous(sys) "It makes no sense to call this function with continuous systems"
vals = sys.(z_or_omegas, map_to_unit_circle)# evalfr.(sys,exp.(evalpoints))
# Reshape from vector of evalfr matrizes, to (in,out,freq) Array
out = Array{eltype(eltype(vals)),3}(undef, length(z_or_omegas), size(sys)...)
for i in 1:length(vals)
out[i,:,:] .= vals[i]
end
return out
nu,ny = size(vals[1])
[v[i,j] for v in vals, i in 1:nu, j in 1:ny]
end

"""`mag, phase, w = bode(sys[, w])`
Expand Down Expand Up @@ -149,26 +128,18 @@ frequencies `w`
function sigma(sys::LTISystem, w::AbstractVector)
resp = freqresp(sys, w)
nw, ny, nu = size(resp)
sv = Array{numeric_type(sys)}(undef, nw, min(ny, nu))
for i=1:nw
sv[i, :] = svdvals(resp[i, :, :])
end
sv = dropdims(mapslices(svdvals, resp, dims=(2,3)),dims=3)
return sv, w
end
sigma(sys::LTISystem) = sigma(sys, _default_freq_vector(sys, :sigma))

function _default_freq_vector(systems::Vector{T}, plot::Symbol) where T<:LTISystem
min_pt_per_dec = 60
min_pt_total = 200
nsys = length(systems)
bounds = Array{numeric_type(systems[1])}(undef, 2, nsys)
for i=1:nsys
# TODO : For now we ignore the feature information. In the future,
# these can be used to improve the frequency vector near features.
bounds[:, i] = _bounds_and_features(systems[i], plot)[1]
end
w1 = minimum(bounds)
w2 = maximum(bounds)
bounds = map(sys -> _bounds_and_features(sys, plot)[1], systems)
w1 = minimum(minimum.(bounds))
w2 = maximum(maximum.(bounds))

nw = round(Int, max(min_pt_total, min_pt_per_dec*(w2 - w1)))
return exp10.(range(w1, stop=w2, length=nw))
end
Expand Down

0 comments on commit 08690a7

Please sign in to comment.