updates to julia1

REQUIRE

@@ -5,3 +5,4 @@ Polynomials
 ProgressMeter
 JLD
 Documenter
+DSP

src/BallAndBeam.jl

+using Pkg
 installed_packages = Pkg.installed()
 if "LabProcesses" ∉ keys(installed_packages)
 	Pkg.clone("https://gitlab.control.lth.se/processes/LabProcesses.jl.git")
@@ -19,7 +20,7 @@ bopl, bopl!, nypl, nypl!, plot, fbdesign, ffdesign, opendoc
 export SysFilter, run_control_2DOF # For documentation
 
 
-using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter
+using LabProcesses, Plots, Polynomials, ControlSystems, ProgressMeter, DSP, LinearAlgebra, Statistics
 
 include("arx.jl")
 
@@ -30,7 +31,7 @@ Perform fra-experiemnt For a single frequency `ω`. Called from inside `fra`
 """
 function run_experiment(P::Beam, ω, duration, settling_time, amplitude)
 	h       = sampletime(P)
-	data    = zeros(0:h:duration)
+	data    = zeros(length(0:h:duration))
 	LabProcesses.initialize(P)
 	for t = 0:h:settling_time
 		@periodically h begin
@@ -51,7 +52,7 @@ end
 
 function run_experiment(P::BeamSimulator, ω, duration, settling_time, amplitude)
 	h       = sampletime(P)
-	data    = zeros(0:h:duration)
+	data    = zeros(length(0:h:duration))
 	LabProcesses.initialize(P)
 	for t = 0:h:settling_time
 		u = amplitude*sin(ω*t)
@@ -139,10 +140,16 @@ end
 
 @recipe function Nypl(p::Nypl)
 	G = p.args[1]
+	@series begin
+		linestyle := :dash
+		color := :black
+		α = LinRange(0,2π,100)
+		cos.(α),sin.(α)
+	end
 	xlabel --> "\$Re(G)\$"
 	ylabel --> "\$Im(G)\$"
 	title  --> "Nyquist plot"
-	real.(G[:,2]), imag.(G[:,2])
+	@series real.(G[:,2]), imag.(G[:,2])
 end
 
 
@@ -175,7 +182,7 @@ function fbdesign(G::AbstractMatrix, polevect, zerovect, gain)
 	sysFB = ss(zpk(zerovect,polevect,gain*prod(abs.(polevect))/pzv))
 	C = Number[ω freqresp(sysFB, ω)]
 	L = Number[ω G[:,2].*C[:,2]]
-	T = Number[ω L[:,2]./(1+L[:,2])]
+	T = Number[ω L[:,2]./(1 .+L[:,2])]
 	sysFB,L,T,C
 end
 
@@ -201,8 +208,8 @@ poles2poly(poles) = reduce((r,l) -> conv(r,[1, l]), [1], poles) #|> reverse |> P
 Build and launch the documentation website.
 """
 function opendoc()
-	include(Pkg.dir("BallAndBeam", "docs","make.jl"))
-	docpath = Pkg.dir("BallAndBeam", "docs","build","index.html")
+	include(joinpath(dirname(pathof("BallAndBeam")), "docs","make.jl"))
+	docpath = joinpath(dirname(pathof("BallAndBeam")), "docs","build","index.html")
 	run(`xdg-open $docpath`)
 end
 

src/FRTN35_lab1.jl

-using BallAndBeam, ControlSystems, JLD, LabProcesses, Plots
+using BallAndBeam, ControlSystems, LabProcesses, Plots, JLD
 # @load "workspace.jld" # Run this command to restore a saved workspace
 
 P    = LabProcesses.Beam(bias = 0.)  # Replace for BeamSimulator() to run simulations
@@ -11,15 +11,15 @@ nbr_of_periods = 3
 # and run three experiments. You may modify the freqency vectors
 # any way you want and add/remove experiments as needed.
 
-w1_100 = logspace(log10(1),log10(300),8)
+w1_100 = exp10.(LinRange(log10(1),log10(300),8)) # exp10.(LinRange = logspace
 G1     = fra(P, w1_100, amplitude=1, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
 @save "workspace.jld"
 
-w1_200 = logspace(log10(50),log10(100),10)
+w1_200 = exp10.(LinRange(log10(50),log10(100),10))
 G2     = fra(P, w1_200, amplitude=2, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
 @save "workspace.jld"
 
-w1_300 = logspace(log10(100),log10(300),20)
+w1_300 = exp10.(LinRange(log10(100),log10(300),20))
 G3     = fra(P, w1_300, amplitude=2, nbr_of_periods=nbr_of_periods, settling_time=settling_time)
 @save "workspace.jld"
 
@@ -35,7 +35,7 @@ nypl(G123, m=:star)
 polevect = [-10]
 zerovect = []
 gain     = 1
-sysFBc,L,T,C = fbdesign(G, polevect, zerovect, gain)
+sysFBc,L,T,C = fbdesign(G123, polevect, zerovect, gain)
 
 polevect = [-10]
 zerovect = []
@@ -55,8 +55,8 @@ plot([y u r], lab = ["y" "u" "r"])
 # If you have time, estimate ARX model =====================================================
 prbs     = PRBSGenerator()
 function prbs_experiment(;amplitude = 1, duration = 10)
-    y = zeros(0:h:duration)
-    u = zeros(0:h:duration)
+    y = zeros(length(0:h:duration))
+    u = zeros(length(0:h:duration))
     LabProcesses.initialize(P)
     for (i,t) = enumerate(0:h:duration)
         @periodically h begin
@@ -71,12 +71,12 @@ function prbs_experiment(;amplitude = 1, duration = 10)
     y,u
 end
 y,u = prbs_experiment(duration = 10, amplitude = 1)
-plot([y u], lab=["y" "u"])
+plot([u y], lab=["u" "y"])
 
 na = 4 # Order of A polynomial
 nb = 2 # Order of B polynomial
 arxtf, Σ = arx(h, y, u, na, nb) # Estimate transfer function with ARX method
 
 bopl(G123, lab="Measured transfer function")
-bodeconfidence!(arxtf, Σ, ω = logspace(0,3,200), linecolor=:red) # The exclamation mark (!)
+bodeconfidence!(arxtf, Σ, ω = exp10.(LinRange(0,3,200)), linecolor=:red) # The exclamation mark (!)
 # adds to the previous plot. If there is no plot open, remove the !

src/arx.jl

@@ -8,10 +8,10 @@ fit(y,yh) = 100 * (1-rms(y-yh)./rms(y-mean(y)));
 aic(x,d) = log(sse(x)) + 2d/length(x)
 
 """
-    toeplitz{T}(c::AbstractArray{T},r::AbstractArray{T})
+    toeplitz(c::AbstractArray,r::AbstractArray)
 Returns a Toeplitz matrix where `c` is the first column and `r` is the first row.
 """
-function toeplitz{T}(c::AbstractVector{T},r::AbstractVector{T})
+function toeplitz(c::AbstractVector{T},r::AbstractVector{T}) where T
     nc = length(c)
     nr = length(r)
     A = zeros(T, nc, nr)
@@ -38,13 +38,13 @@ B     = [10,5] # B(z) coeffs
 u     = randn(100) # Simulate 100 time steps with Gaussian input
 y     = filt(B,A,u)
 yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
-x     = A\yr # Estimate model polynomials
+x     = A\\yr # Estimate model polynomials
 plot([yr A*x], lab=["Signal" "Prediction"])
 ```
 For nonlinear ARX-models, see [BasisFunctionExpansions.jl](https://github.com/baggepinnen/BasisFunctionExpansions.jl/)
 """
-function getARXregressor(y::AbstractVector{Float64},u::AbstractVecOrMat{Float64}, na, nb)
-    assert(length(nb) == size(u,2))
+function getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
+    @assert(length(nb) == size(u,2))
     m    = max(na+1,maximum(nb))
     n    = length(y) - m+1
     offs = m-na-1
@@ -86,7 +86,7 @@ function arx(h,y::AbstractVector{Float64}, u::AbstractVector{Float64}, na, nb;
     if λ == 0
         w = A\y_train
     else
-        w = (A'A + λ*eye(size(A,2)))\A'y_train
+        w = (A'A + λ*I)\A'y_train
     end
     a,b = params2poly(w,na,nb)
     model = tf(b,a,h)
@@ -116,27 +116,27 @@ function parameter_covariance(y_train, A, w, λ=0)
         ATAλ\ATA/ATAλ
     end
     iATA = (iATA+iATA')/2
-    Σ = σ²*iATA + sqrt(eps())*eye(iATA)
+    Σ = σ²*iATA + sqrt(eps())*I
 end
 
 """
-    bodeconfidence(arxtf::TransferFunction, Σ::Matrix; ω = logspace(0,3,200))
+    bodeconfidence(arxtf::TransferFunction, Σ::Matrix; ω = exp10.(LinRange(0,3,200)))
 Plot a bode diagram of a transfer function estimated with [`arx`](@ref) with confidence bounds on magnitude and phase.
 """
 bodeconfidence
 
 @userplot BodeConfidence
 
-@recipe function BodeConfidence(p::BodeConfidence; ω = logspace(-2,3,200))
+@recipe function BodeConfidence(p::BodeConfidence; ω = exp10.(LinRange(-2,3,200)))
     arxtfm = p.args[1]
     Σ      = p.args[2]
-    L      = chol(Hermitian(Σ))
+    L      = cholesky(Hermitian(Σ))
     am, bm = -reverse(denpoly(arxtfm)[1].a[1:end-1]), reverse(arxtfm.matrix[1].num.a)
     wm     = [am; bm]
     na,nb  = length(am), length(bm)
     mc     = 100
     res = map(1:mc) do _
-        w             = L'randn(size(L,1)) .+ wm
+        w             = L.L*randn(size(L,1)) .+ wm
         a,b           = params2poly(w,na,nb)
         arxtf         = tf(b,a,arxtfm.Ts)
         mag, phase, _ = bode(arxtf, ω)
@@ -144,8 +144,8 @@ bodeconfidence
     end
     magmc      = hcat(getindex.(res,1)...)
     phasemc    = hcat(getindex.(res,2)...)
-    mag        = mean(magmc,2)[:]
-    phase      = mean(phasemc,2)[:]
+    mag        = mean(magmc,dims=2)[:]
+    phase      = mean(phasemc,dims=2)[:]
     uppermag   = getpercentile(magmc,0.95)
     lowermag   = getpercentile(magmc,0.05)
     upperphase = getpercentile(phasemc,0.95)
@@ -176,7 +176,7 @@ bodeconfidence
 end
 
 function getpercentile(mag,p)
-    uppermag = mapslices(mag, 2) do magω
-        sort(magω)[round(Int,endof(magω)*p)]
+    uppermag = mapslices(mag, dims=2) do magω
+        sort(magω)[round(Int,length(magω)*p)]
     end
 end