.gitignore

 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+docs/build/
+docs/site/
+docs/.documenter

.gitlab-ci.yml

+# This file is a template, and might need editing before it works on your project.
+# An example .gitlab-ci.yml file to test (and optionally report the coverage
+# results of) your [Julia][1] packages. Please refer to the [documentation][2]
+# for more information about package development in Julia.
+
+
+# Below is the template to run your tests in Julia
+.test_template: &test_definition
+  # Uncomment below if you would like to run the tests on specific references
+  # only, such as the branches `master`, `development`, etc.
+  # only:
+  #   - master
+  #   - development
+  script:
+    # Let's run the tests. Substitute `coverage = false` below, if you do not
+    # want coverage results.
+    - julia -e 'Pkg.rm("LabProcesses");Pkg.clone(pwd()); Pkg.test("LabProcesses",coverage = false)'
+    #- julia -e 'Pkg.update();Pkg.test("LabProcesses", coverage = true)'
+    # Comment out below if you do not want coverage results.
+    #- julia -e 'Pkg.add("Coverage"); cd(Pkg.dir("LabProcesses"));
+    #  using Coverage; cl, tl = get_summary(process_folder());
+    #  println("(", cl/tl*100, "%) covered")'
+
+# Name a test and select an appropriate image.
+test:0.6.0:
+  image: julialang/julia:v0.6.0
+  <<: *test_definition
+
+
+# REMARK: Do not forget to enable the coverage feature for your project, if you
+# are using code coverage reporting above. This can be done by
+#
+# - Navigating to the `CI/CD Pipelines` settings of your project,
+# - Copying and pasting the default `Simplecov` regex example provided, i.e.,
+#   `\(\d+.\d+\%\) covered` in the `test coverage parsing` textfield.
+#
+# WARNING: This template is using the `julialang/julia` images from [Docker
+# Hub][3]. One can use custom Julia images and/or the official ones found
+# in the same place. However, care must be taken to correctly locate the binary
+# file (`/opt/julia/bin/julia` above), which is usually given on the image's
+# description page.
+#
+# [3]: http://hub.docker.com/

README.md

+## News
+2018-12-07: Updated to julia v1.0, see commit eac09291 for last julia v0.6 version
+
+[![pipeline status](https://gitlab.control.lth.se/processes/LabProcesses.jl/badges/master/pipeline.svg)](https://gitlab.control.lth.se/processes/LabProcesses.jl/commits/master)
+[![coverage report](https://gitlab.control.lth.se/processes/LabProcesses.jl/badges/master/coverage.svg)](https://gitlab.control.lth.se/processes/LabProcesses.jl/commits/master)
+
 # LabProcesses
-This package contains an (programming- as well as connection-) interface to serve
-as a base for the implementation of lab-process software. The first exampel of
-an implementaiton of this interface is for the ball-and-beam process, which is
-used in Lab1 FRTN35: frequency response analysis of the beam. The lab is implemented
-in [BallAndBeam.jl](https://gitlab.control.lth.se/processes/BallAndBeam.jl), a
-package that makes use of `LabProcesses.jl` to handle the communication with
-the lab process and/or a simulated version thereof. This way, the code written
-for frequency response analysis of the beam can be run on another process
-implementing the same interface (or a simulated version) by changeing a single
-line of code :)
-
-## Installation
-1. Start julia
-2. Install LabProcesses.jl using command Pkg.clone("https://gitlab.control.lth.se/processes/LabProcesses.jl.git") Lots of packages will now be installed, this will take some time.
-
-## How to implement a new process
-1. Locate the file [interface.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface.jl)
-2. Copy all function definitions.
-3. Create a new file under `/interface_implementations` where you paste all the
-copied definitions and implement them. See [`/interface_implementations/ballandbeam.jl`](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface_implementations/ballandbeam.jl) for an example.
-4. Above all function implementations you must define the process type, e.g,
-```julia
-struct BallAndBeam <: PhysicalProcess
-    h::Float64
-end
-BallAndBeam() = BallAndBeam(0.01) # Constructor with default value of sample time
-```
-Make sure you inherit from `PhysicalProcess` or `SimulatedProcess` as appropriate.
-5. Documentation of all interface functions is available in the file [interface_documentation.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface_documentation.jl)
-
-## Control a process
-The interface `AbstractProcess` defines the functions `control(P, u)` and `measure(P)`.
-These functions can be used to implement your own control loops. A common loop
-with a feedback controller and a feedforward filter on the reference is implemented
-in the function [`run_control_2DOF`](@ref).
-
-The macro `@periodically` might come in handy if you want to implement your own loop.
-Consider the following example, in which the loop body will be run periodically
-with a sample time of `h` seconds.
-```julia
-for (i,t) = enumerate(0:h:duration)
-    @periodically h begin
-        y[i]       = measure(P)
-        r[i]       = reference(t)
-        u[i]       = control(i)
-        control(P, u[i])
-    end
-end
-```
-
-Often one finds the need to implement a stateful controller, i.e., a function
-that has a memory or state. To this end, the function [`sysfilter`](@ref) is
-provided. Usage is demonstrated below
-```julia
-stateC = init_sysfilter(C)
-function control(i)
-    e  = r[i]-y[i]
-    ui = sysfilter!(stateC, C, e)
-end
-```
-`C` must here be represented by a [`StateSpace`](@ref) type from [`ControlSystems.jl`](https://github.com/JuliaControl/ControlSystems.jl).
-`TransferFunction` types can easily be converted to a `StateSpace` by `Gss = ss(Gtf)`.
+Documentation available at [Documentation](http://processes.gitlab.control.lth.se/documentation/labprocesses/)
+
+
+# Automatiskt genererad dokumentation
+Det finns i skrivande stund två stycken mer eller mindre automatgenererade hemsidor med dokumentation, en sida till repot BallAndBeam.jl
+    http://processes.gitlab.control.lth.se/documentation/ballandbeam/
+
+och en sida till (detta) repot LabProcesses.jl
+    http://processes.gitlab.control.lth.se/documentation/labprocesses
+
+Mitt förslag är att alla som känner sig ansvariga för ett repo som är skapat med syfte att till slut överlämnas till framtida kollegor sätter upp liknande
+dokumentation för detta repo.
+
+## Generera dokumentation
+Jag har byggt dokumentationen med [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl). Detta verktyg accepterar en markdown-fil (docs/src/index.md)
+med lite text man skrivit samt macron som indikerar att man vill infoga docstrings från funktioner och typer i sitt juliapaket.
+Outputen är en html-sida med tillbehör som man ska se till att hosta på gitlab pages.
+
+Exempel på hur detta går till finns i repona LabProcesses.jl och BallAndBeam.jl. Man kan helt sonika kopiera docs/-mappen från ett av dessa repon och
+modifiera innehållet i alla filer så att det stämmer med ens nya repo.
+
+## Hosting
+För hosting finns repot
+https://gitlab.control.lth.se/processes/documentation/
+Det man behöver göra är att editera filen
+https://gitlab.control.lth.se/processes/documentation/blob/master/.gitlab-ci.yml
+och lägga till tre rader som
+
+- Bygger dokumentationen i repot man vill dokumentera
+- skapar en ny mapp med ett väl valt namn (myfoldername i exemplet nedan)
+- flyttar den byggda dokumentationen till den mappen
+- Det blir lätt att förstå hur man gör punkterna ovan när man kollar i filen .gitlab-ci.yml, bara kör mönstermatchning mot de repona som detta redan är gjort för.
+
+När .gitlab-ci.yml uppdateras i master triggas en pipline. Om denna lyckas kommer dokumentationen finnas under
+
+http://processes.gitlab.control.lth.se/documentation/myfoldername/

REQUIRE

-julia 0.6
+julia 0.7
 ControlSystems
+Parameters
+DSP

docs/make.jl

+using Documenter, LabProcesses
+# makedocs()
+# deploydocs(
+# deps   = Deps.pip("pygments", "mkdocs", "python-markdown-math", "mkdocs-cinder"),
+# repo   = "gitlab.control.lth.se/processes/LabProcesses.jl",
+# branch = "gh-pages",
+# julia  = "0.6",
+# osname = "linux"
+# )
+
+makedocs(
+    format = :html,
+    sitename = "LabProcesses",
+    pages = [
+        "index.md",
+    ]
+)

docs/mkdocs.yml

+site_name:        LabProcesses.jl
+repo_url:         https://gitlab.control.lth.se/labdev/software/tree/master/julia_ballandbeam/LabProcesses.jl
+site_description: Documentation for LabProcesses.jl
+site_author:      Fredri Bagge Carlson
+
+theme: cinder
+
+extra_css:
+  - assets/Documenter.css
+
+extra_javascript:
+  - https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML
+  - assets/mathjaxhelper.js
+
+markdown_extensions:
+  - extra
+  - tables
+  - fenced_code
+  - mdx_math:
+      enable_dollar_delimiter: True
+
+docs_dir: 'build'
+
+pages:
+  - Home: index.md

docs/src/index.md

+# LabProcesses
+
+```@contents
+Depth = 3
+```
+
+This package contains an (programming- as well as connection-) interface to serve
+as a base for the implementation of lab-process software. The first example of
+an implementaiton of this interface is for the ball-and-beam process, which is
+used in Lab1 FRTN35: frequency response analysis of the beam. The lab is implemented
+in [BallAndBeam.jl](https://gitlab.control.lth.se/processes/BallAndBeam.jl), a
+package that makes use of `LabProcesses.jl` to handle the communication with
+the lab process and/or a simulated version thereof. This way, the code written
+for frequency response analysis of the beam can be run on another process
+implementing the same interface (or a simulated version) by changeing a single
+line of code :)
+
+## Installation
+1. Start julia by typing `julia` in a terminal, make sure the printed info says it's
+`v0.6+` running. If not, visit [julialang.org](https://julialang.org/downloads/)
+to get the latest release.
+2. Install LabProcesses.jl using command `Pkg.clone("https://gitlab.control.lth.se/processes/LabProcesses.jl.git")` Lots of packages will now be installed, this will take some time. If this is your first time using Julia, you might have to run `julia> Pkg.init()` before you install any packages.
+
+# How to implement a new process
+### 1.
+Locate the file [interface.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface.jl). When the package is installed, you find its directory under `~/.julia/v0.6/LabProcesses/`, if not, run `julia> Pkg.dir("LabProcesses")` to locate the directory.
+(Alternatively, you can copy all definitions from [/interface_implementations/ballandbeam.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface_implementations/ballandbeam.jl) instead. Maybe it's easier to work from an existing implementaiton.)
+### 2.
+Copy all function definitions.
+### 3.
+Create a new file under `/interface_implementations` where you paste all the
+copied definitions and implement them. See [/interface_implementations/ballandbeam.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface_implementations/ballandbeam.jl) for an example.
+### 4.
+Above all function implementations you must define the process type, e.g,
+```julia
+struct BallAndBeam <: PhysicalProcess
+    h::Float64
+    bias::Float64
+end
+BallAndBeam() = BallAndBeam(0.01, 0.0) # Constructor with default value of sample time
+```
+Make sure you inherit from `PhysicalProcess` or `SimulatedProcess` as appropriate.
+This type must contains fields that hold information about everything that is
+relevant to a particular instance of the process. Different ballandbeam-process
+have different biases, hence this must be stored. A simulated process would have
+to keep track of its state etc. in order to implement the measure and control
+methods. See [Types in julia documentation](https://docs.julialang.org/en/stable/manual/types/#Composite-Types-1)
+for additional info regarding user defined types and (constructors)[https://docs.julialang.org/en/stable/manual/constructors/].
+### 5.
+Documentation of all interface functions is available in the file [interface_documentation.jl](https://gitlab.control.lth.se/processes/LabProcesses.jl/blob/master/src/interface_documentation.jl)
+
+# How to control a process
+The interface `AbstractProcess` defines the functions `control(P, u)` and `measure(P)`.
+These functions can be used to implement your own control loops. A common loop
+with a feedback controller and a feedforward filter on the reference is implemented
+in the function [`run_control_2DOF`](@ref), where the user can supply $G_1$ and $G_4$
+in the diagram below, with the process $P=G_2$.
+![block diagram](feedback4.png)
+
+The macro `@periodically` might come in handy if you want to implement your own loop.
+Consider the following example, in which the loop body will be run periodically
+with a sample time of `h` seconds.
+```julia
+for (i,t) = enumerate(0:h:duration)
+    @periodically h begin
+        y[i] = measure(P)
+        r[i] = reference(t)
+        u[i] = calc_control(y,r)
+        control(P, u[i])
+    end
+end
+```
+
+
+Often one finds the need to implement a stateful controller, i.e., a function
+that has a memory or state. To this end, the type [`SysFilter`](@ref) is
+provided. This type is used to implement control loops where a signal is
+filtered through a dynamical system, i.e., `U(z) = G1(z)E(z)`.
+Usage is demonstrated below, which is a simplified implementation of the block
+diagram above (transfer function- and signal names corresponds to the figure).
+First two `SysFilter` objects are created, these objects can now be used as
+functions of an input, and return the filtered output. The `SysFilter` type takes
+care of updating and remembering the state of the system when called.
+```julia
+G1f = SysFilter(G1)
+G4f = SysFilter(G4)
+function calc_control(y,r)
+    rf = G4f(r)
+    e  = rf-y
+    u  = G1f(e)
+end
+```
+`G1` and `G4` must here be represented by [`StateSpace`](http://juliacontrol.github.io/ControlSystems.jl/latest/lib/constructors/#ControlSystems.ss) types
+from [`ControlSystems.jl`](https://github.com/JuliaControl/ControlSystems.jl), e.g., `G1 = ss(A,B,C,D)`.
+`TransferFunction` types can easily be converted to a `StateSpace` by `Gss = ss(Gtf)`.
+Continuous time systems can be discretized using `Gd = c2d(Gc, h)[1]`. (The sample time of a process is available through `h = sampletime(P)`.)
+
+
+# How to implement a Simulated Process
+## Linear process
+This is very easy, just get a discrete time `StateSpace` model of your process
+(if you have a transfer function, `Gss = ss(Gtf)` will do the trick, if you have continuous time, `Gd = c2d(Gc,h)[1]` is your friend).
+
+You now have to implement the methods `control` and `measure` for your simulated type.
+The implementation for `BeamSimulator` is shown below
+```julia
+control(p::BeamSimulator, u) = p.Gf(u)
+measure(P) = vecdot(p.Gf.sys.C, p.Gf.state)
+```
+The `control` method accepts a control signal (`u`) and propagates the system state
+(`p.Gf.state`) forward using the statespace model (`p.Gf.sys`) of the beam. The object
+`Gf` (of type [`SysFilter`](@ref)) is familiar from the "Control" section above. What it does
+is essentially (simplified)
+```julia
+function Gf(input)
+    sys       = Gf.sys
+    Gf.state .= sys.A*Gf.state + sys.B*input
+    output    = sys.C*Gf.state + sys.D*input
+end
+```
+hence, it just performs one iteration of
+```math
+x' = Ax + Bu
+```
+```math
+y  = Cx + Du
+```
+
+The `measure` method performs the computation `y = Cx`, the reason for the call
+to `vecdot` is that `vecdot` produces a scalar output, whereas `C*x` produces a
+1-element `Matrix`. A scalar output is preferred in this case since the `Beam`
+is SISO.
+
+It should now be obvious which fields are required in the `BeamSimulator` type.
+It must know which sample time it has been discretized with, as well as its
+discrete-time system model. It must also remember the current state of the system.
+This is not needed in a physical process since it kind of remembers its own state.
+The system model and its state is conveniently covered by the type [`SysFilter`](@ref),
+which handles filtering of a signal through an LTI system.
+The full type specification for `BeamSimulator` is given below
+```julia
+struct BeamSimulator <: SimulatedProcess
+    h::Float64
+    Gf::SysFilter
+    BeamSimulator() = new(0.01, SysFilter(beam_system, 0.01))
+    BeamSimulator(h::Real) = new(Float64(h), SysFilter(beam_system, h))
+end
+```
+It contains three fields and two inner constructors. The constructors initializes
+the system filter by creating a [`SysFilter`](@ref).
+The variable `beam_system` is already defined outside the type specification.
+One of the constructors provides a default value for the sample time, in case
+the user is unsure about a reasonable value.
+
+## Non-linear process
+Your first option is to linearize the process and proceed like above.
+Other options include
+1. Make `control` perform forward Euler, i.e., `x[t+1] = x[t] + f(x[t],u[t])*h` for a general system model ``x' = f(x,u); y = g(x,u)`` and sample time ``h``.
+2. Integrate the system model using some fancy method like Runge-Kutta. See [DifferentialEquations.jl](http://docs.juliadiffeq.org/stable/types/discrete_types.html) for discrete-time solving of ODEs (don't be discouraged, this is almost as simple as forward Euler above).
+
+# Exported functions and types
+```@autodocs
+Modules = [LabProcesses]
+Private = false
+Pages   = ["LabProcesses.jl", "controllers.jl", "reference_generators.jl", "utilities.jl"]
+```
+
+# Process interface specification
+All processes must implement the following interface. See the existing
+implementations in the folder `interface_implementations` for guidance.
+
+```@autodocs
+Modules = [LabProcesses]
+Private = false
+Pages   = ["interface_documentation.jl"]
+```
+
+# Index
+
+```@index
+```

src/LabProcesses.jl

+# __precompile__()
+using Pkg
+installed_packages = Pkg.installed()
+if "LabConnections" ∉ keys(installed_packages)
+	Pkg.clone("https://gitlab.control.lth.se/cont-frb/LabConnections.jl")
+	# Pkg.checkout("LabConnection","v0.0.1")
+end
+
+
 module LabProcesses
 
-using ControlSystems
+using ControlSystems, LabConnections.Computer, Parameters, DSP, LinearAlgebra
+
+include("utilities.jl")
 
 include("interface.jl")
 include("interface_documentation.jl")
 include("interface_implementations/ballandbeam.jl")
+include("interface_implementations/eth_helicopter.jl")
 
-include("utilities.jl")
+include("reference_generators.jl")
 include("controllers.jl")
 
 end # module

src/controllers.jl

 export run_control_2DOF
 
 """
-	y,u,r = run_control(process, sysFB[, sysFF]; duration = 10, reference(t) = sign(sin(2π*t)))
-Perform control experiemnt where the feedback and feedforward controllers are given by
+	y,u,r = run_control_2DOF(process, sysFB[, sysFF]; duration = 10, reference(t) = sign(sin(2π*t)))
+Perform control experiemnt on process where the feedback and feedforward controllers are given by
 `sysFB` and `sysFF`, both of type `StateSpace`.
 
 `reference` is a reference generating function that accepts a scalar `t` (time in seconds) and outputs a scalar `r`, default is `reference(t) = sign(sin(2π*t))`.
 
 The outputs `y,u,r` are the beam angle, control signal and reference respectively.
+![block diagram](feedback4.png)
 """
 function run_control_2DOF(P::AbstractProcess,sysFB, sysFF=nothing; duration = 10, reference = t->sign(sin(2π*t)))
+	nu = num_inputs(P)
+	ny = num_outputs(P)
 	h  = sampletime(P)
-	y       = zeros(0:h:duration)
-	u       = zeros(0:h:duration)
-	r       = zeros(0:h:duration)
+	y  = zeros(ny, length(0:h:duration))
+	u  = zeros(nu, length(0:h:duration))
+	r  = zeros(ny, length(0:h:duration))
 
-	sysFB   = minreal(sysFB)
-	stateFB = init_sysfilter(sysFB)
+	Gfb = SysFilter(sysFB)
 	if sysFF != nothing
-		sysFF   = minreal(sysFF)
-		stateFF = init_sysfilter(sysFF)
+		Gff = SysFilter(sysFF)
 	end
 
-	function control(i)
-		rf = sysFF == nothing ? r[i] : sysfilter!(stateFF, sysFF, r[i])
-		e  = rf-y[i]
-		ui = sysfilter!(stateFB, sysFB, e)
+	function calc_control(i)
+		rf = sysFF == nothing ? r[:,i] : Gff(r[:,i])
+		e  = rf-y[:,i]
+		ui = Gfb(e)
+		ui .+ bias(P)
 	end
 
+	simulation = isa(P, SimulatedProcess)
 	initialize(P)
 	for (i,t) = enumerate(0:h:duration)
-		@periodically h begin
-			y[i]       = measure(P)
-			r[i]       = reference(t)
-			u[i]       = control(i) # y,r must be updated before u
-			control(P, u[i])
+		@periodically h simulation begin
+			y[:,i]    .= measure(P)
+			r[:,i]    .= reference(t)
+			u[:,i]    .= calc_control(i) # y,r must be updated before u
+			control(P, [clamp.(u[j,i], inputrange(P)[j]...) for j=1:nu])
 		end
 	end
 	finalize(P)
-	y,u,r
+	y',u',r'
 end

src/interface.jl

@@ -6,6 +6,7 @@ export  num_outputs,
         isstable,
         isasstable,
         sampletime,
+        bias,
         control,
         measure,
         initialize,
@@ -24,6 +25,7 @@ inputrange(p::AbstractProcess)  = error("Function not implemented for $(typeof(p
 isstable(p::AbstractProcess)    = error("Function not implemented for $(typeof(p))")
 isasstable(p::AbstractProcess)  = error("Function not implemented for $(typeof(p))")
 sampletime(p::AbstractProcess)  = error("Function not implemented for $(typeof(p))")
+bias(p::AbstractProcess)        = error("Function not implemented for $(typeof(p))")
 
 control(p::AbstractProcess, u)  = error("Function not implemented for $(typeof(p))")
 measure(p::AbstractProcess)     = error("Function not implemented for $(typeof(p))")

src/interface_documentation.jl

@@ -62,6 +62,14 @@ Return the sample time of the process in seconds.
 """
 sampletime
 
+"""
+    b = bias(P::AbstractProcess)
+Return an input bias for the process. This could be, i.e., the constant input u₀ around which
+a nonlinear system is linearized, or whatever other bias might exist on the input.
+`length(b) = num_inputs(P)`
+"""
+bias
+
 """
     control(P::AbstractProcess, u)
 Send a control signal to the process. `u` must have dimension equal to `num_inputs(P)`

src/interface_implementations/ballandbeam.jl

 # Interface implementation Ball And Beam ====================================================
 
-export BallAndBeam, BallAndBeamSimulator, BallAndBeamType
+# The ball and beam can be used in two modes, either just the Beam, in which case there is a
+# single output (measurement signal) or the BallAndBeam, in which case there are two.
+# There are a few union types defined for convenience, these are
+# AbstractBeam = Union{Beam, BeamSimulator}
+# AbstractBallAndBeam   = Union{BallAndBeam, BallAndBeamSimulator}
+# AbstractBeamOrBallAndBeam = All types
+# Although not Abstract per se, the names AbstractBeam etc. were chosen since this reflects
+# their usage in dispatch.
+
+export Beam, BeamSimulator, AbstractBeam, BallAndBeam, BallAndBeamSimulator, AbstractBeamOrBallAndBeam
+
+# @with_kw allows specification of default values for fields. If none is given, this value must be supplied by the user. replaces many constructors that would otherwise only supply default values.
+# Call constructor like Beam(bias=1.0) if you want a non-default value for bias
+"""
+    Beam(;kwargs...)
+Physical beam process
+
+#Arguments (fields)
+- `h::Float64 = 0.01`
+- `bias::Float64 = 0.0`
+- `stream::LabStream = ComediStream()`
+- `measure::AnalogInput10V = AnalogInput10V(0)`
+- `control::AnalogOutput10V = AnalogOutput10V(1)`
+"""
+struct Beam <: PhysicalProcess
+    h::Float64
+    bias::Float64
+    stream::LabStream
+    measure::AnalogInput10V
+    control::AnalogOutput10V
+end
+function Beam(;
+    h::Float64               = 0.01,
+    bias::Float64            = 0.,
+    stream::LabStream        = ComediStream(),
+    measure::AnalogInput10V  = AnalogInput10V(0),
+    control::AnalogOutput10V = AnalogOutput10V(1))
+    p = Beam(Float64(h),Float64(bias),stream,measure,control)
+    init_devices!(p.stream, p.measure, p.control)
+    p
+end
+
+include("define_beam_system.jl")
+const beam_system, nice_beam_controller = define_beam_system()
+# nice_beam_controller gives ϕₘ=56°, Aₘ=4, Mₛ = 1.6. Don't forget to discretize it before use
+struct BeamSimulator <: SimulatedProcess
+    h::Float64
+    s::SysFilter
+    BeamSimulator(;h::Real = 0.01, bias=0) = new(Float64(h), SysFilter(beam_system, h))
+end
 
 struct BallAndBeam <: PhysicalProcess
     h::Float64
+    bias::Float64
+    stream::LabStream
+    measure1::AnalogInput10V
+    measure2::AnalogInput10V
+    control::AnalogOutput10V
+end
+function BallAndBeam(;
+    h                        = 0.01,
+    bias                     = 0.,
+    stream                   = ComediStream(),
+    measure1::AnalogInput10V = AnalogInput10V(0),
+    measure2::AnalogInput10V = AnalogInput10V(1),
+    control::AnalogOutput10V = AnalogOutput10V(1))
+    p = BallAndBeam(h,bias,stream,measure1,measure2,control)
+    init_devices!(p.stream, p.measure1, p.measure2, p.control)
+    p
 end
-BallAndBeam() = BallAndBeam(0.01)
 
 struct BallAndBeamSimulator <: SimulatedProcess
     h::Float64
+    s::SysFilter
 end
-BallAndBeamSimulator() = BallAndBeamSimulator(0.01)
-
-const BallAndBeamType = Union{BallAndBeam, BallAndBeamSimulator}
-num_outputs(p::BallAndBeamType) = 2
-num_inputs(p::BallAndBeamType)  = 1
-outputrange(p::BallAndBeamType) = [(-10,10),(-1,1)] # Beam angle, Ball position
-inputrange(p::BallAndBeamType)  = [(-10,10)]
-isstable(p::BallAndBeamType)    = false
-isstable(p::BallAndBeamType)    = false
-sampletime(p::BallAndBeamType)  = p.h
-
-control(p::BallAndBeam, u)  = ccall((:comedi_write, comedipath),Int32,
-                                        (Int32,Int32,Int32,Int32),0,1,1,num2io(u[1]))
-measure(p::BallAndBeam)     = io2num(ccall((:comedi_read,comedipath), Int32,
-                                        (Int32,Int32,Int32), 0,0,0))
-
-control(p::BallAndBeamSimulator, u)  = error("Not yet implemented")
+
+const AbstractBeam              = Union{Beam, BeamSimulator}
+const AbstractBallAndBeam       = Union{BallAndBeam, BallAndBeamSimulator}
+const AbstractBeamOrBallAndBeam = Union{AbstractBeam, AbstractBallAndBeam}
+
+num_outputs(p::AbstractBeam)             = 1
+num_outputs(p::AbstractBallAndBeam)      = 2
+num_inputs(p::AbstractBeamOrBallAndBeam) = 1
+outputrange(p::AbstractBeam)             = [(-10,10)]
+outputrange(p::AbstractBallAndBeam)      = [(-10,10),(-1,1)] # Beam angle, Ball position
+inputrange(p::AbstractBeamOrBallAndBeam) = [(-10,10)]
+isstable(p::AbstractBeam)                = true
+isstable(p::AbstractBallAndBeam)         = false
+isasstable(p::AbstractBeamOrBallAndBeam) = false
+sampletime(p::AbstractBeamOrBallAndBeam) = p.h
+bias(p::AbstractBeamOrBallAndBeam)       = p.bias
+bias(p::BeamSimulator)                   = 0
+bias(p::BallAndBeamSimulator)            = 0
+
+function control(p::AbstractBeamOrBallAndBeam, u::AbstractArray)
+    length(u) == 1 || error("Process $(typeof(p)) only accepts one control signal, tried to send u=$u.")
+    control(p,u[1])
+end
+control(p::AbstractBeamOrBallAndBeam, u::Number) = send(p.control,u)
+control(p::BeamSimulator, u::Number)             = p.s(u)
+control(p::BallAndBeamSimulator, u::Number)      = error("Not yet implemented")
+
+measure(p::Beam)                         = read(p.measure)
+measure(p::BallAndBeam)                  = [read(p.measure1), read(p.measure2)]
+measure(p::BeamSimulator)                = dot(p.s.sys.C,p.s.state)
 measure(p::BallAndBeamSimulator)         = error("Not yet implemented")
 
-const comedipath = "../c/comedi_bridge.so"
-const conversion = 65535/20
-io2num(x) = x/conversion -10            # Converts from io to float
-num2io(x) = round(Int32,x*100 + 2050)   # Converts from regular number to io
 
-initialize(p::BallAndBeam)  = ccall((:comedi_start, comedipath),Int32,(Int32,), 0)
-finalize(p::BallAndBeam)    = (control(P,0);ccall((:comedi_stop, "../c/comedi_bridge.so"),Int32,(Int32,), 0))
+initialize(p::Beam)                    = nothing
+initialize(p::BallAndBeam)             = nothing
+finalize(p::AbstractBeamOrBallAndBeam) = foreach(close, p.stream.devices)
 initialize(p::BallAndBeamSimulator)    = nothing
 finalize(p::BallAndBeamSimulator)      = nothing
+initialize(p::BeamSimulator)           = p.s.state .*= 0
+finalize(p::BeamSimulator)             = nothing

src/interface_implementations/define_beam_system.jl

+using ControlSystems, DSP
+
+"""
+    beammodel, beamcontroller = define_beam_system(;doplot=false)
+
+Return a continuous time model of the beam (`::StateSpace`) and a feedback controller that
+gives ϕₘ = 56°, Aₘ = 4, Mₛ = 1.6.
+Don't forget to discretize both before use. Call this method with
+doplot = true to display a bunch of plots to convince yourself.
+
+The numbers in the code have been selected based on the fra performed and reported in
+`beamdetails.pdf`
+"""
+function define_beam_system(;doplot=false)
+    K  = 5    # gain of internal velocity control
+    I1 = 1    # intertia of motor mass
+    I2 = 1    # inertia of beam
+    C  = 5    # beam flexibility
+    D  = 0.1  # beam damping
+    T  = 0.01 # motor torque time constant
+    E  = 1    # motor viscous damping
+
+    # transfer function from input signal to motor position
+    G1 = 2*tf([I2, D, C]*K,conv([T, 1],conv([I1, (D+E), C],[I2, D, C]) - [0; 0; conv([D, C],[D, C])]) + K*[0, 0, I2, D, C, 0])
+
+    G2 = 2*tf([D, C]*K,conv([T, 1],conv([I1, (D+E), C],[I2, D, C]) - [0; 0; conv([D, C],[D, C])]) + K*[0, 0, I2, D, C, 0])
+    # pzmap(G1)
+    # bodeplot([G1,G2], logspace(-2,3,500))
+
+    # Above model is super inaccurate and I have no idea where it comes from, let's define
+    # a simple but more accurate model // Bagge
+
+    G1 = zpk([-145],[0, -60, -100], 4.5*60*100/145)
+    ωz = 145
+    ζz = 0.01
+    ωp = 150
+    ζp = 0.005
+    G1 *= tf([1,2ζz*ωz, ωz^2]*ωp^2, [1,2ζp*ωp, ωp^2]*ωz^2)
+
+
+    sysFBc = ss(zpk([], [-40, -40], 2*40*40)) # Controller suggesten in assistant manual, works fine =)
+    # gives ϕₘ = 56°, Aₘ = 4, Mₛ = 1.6. Don't forget to discretize it before use. Call this method with
+    # doplot = true to display a bunch of plots to convince yourself
+    if doplot
+        ControlSystems.setPlotScale("log10")
+        bodeplot(tf(G1), logspace(0,log10(200),300))
+        bodeplot([tf(G1), sysFBc*G1], logspace(0,log10(200),300))
+        marginplot(sysFBc*tf(G1))
+        gangoffourplot(tf(G1), tf(sysFBc))
+    end
+    ss(G1), sysFBc
+end

src/interface_implementations/eth_helicopter.jl

+# Interface implementation Ball And Beam ====================================================
+
+# There is a union type defined for convenience:
+# AbstractETHHelicopter = Union{ETHHelicopter, ETHHelicopterSimulator}
+# Although not Abstract per se, the names AbstractETHHelicopter etc. were chosen since this
+# reflects the usage in dispatch.
+
+export ETHHelicopter, ETHHelicopterSimulator, AbstractETHHelicopter
+
+# @with_kw allows specification of default values for fields. If none is given, this value must be supplied by the user. replaces many constructors that would otherwise only supply default values.
+# Call constructor like ETHHelicopter(bias=1.0) if you want a non-default value for bias
+"""
+    ETHHelicopter(;kwargs...)
+    Physical ETH helicopter process
+# Arguments (fields)
+- `h::Float64 = 0.05`
+- `bias::Float64 = 0.0`
+- `stream::LabStream = ComediStream()`
+- `measure1::AnalogInput10V = AnalogInput10V(0)`
+- `measure2::AnalogInput10V = AnalogInput10V(1)`
+- `control1::AnalogOutput10V = AnalogOutput10V(0)`
+- `control2::AnalogOutput10V = AnalogOutput10V(1)`
+"""
+@with_kw struct ETHHelicopter <: PhysicalProcess
+    h::Float64
+    bias::Float64
+    stream::LabStream
+    measure1::AnalogInput10V
+    measure2::AnalogInput10V
+    control1::AnalogOutput10V
+    control2::AnalogOutput10V
+end
+function ETHHelicopter(;
+    h                         = 0.05,
+    bias                      = 0.,
+    stream                    = ComediStream(),
+    measure1::AnalogInput10V  = AnalogInput10V(0),
+    measure2::AnalogInput10V  = AnalogInput10V(1),
+    control1::AnalogOutput10V = AnalogOutput10V(0),
+    control2::AnalogOutput10V = AnalogOutput10V(1))
+    p = ETHHelicopter(h,bias,stream,measure1,measure2,control1,control2)
+    init_devices!(p.stream, p.measure1, p.measure2, p.control1, p.control2)
+    p
+end
+
+
+struct ETHHelicopterSimulator <: SimulatedProcess
+    h::Float64
+    bias::Float64
+    state::Vector{Float64}
+end
+ETHHelicopterSimulator() = ETHHelicopterSimulator(0.01, zeros(4))
+
+const AbstractETHHelicopter = Union{ETHHelicopter, ETHHelicopterSimulator}
+num_outputs(p::AbstractETHHelicopter) = 2
+num_inputs(p::AbstractETHHelicopter)  = 2
+outputrange(p::AbstractETHHelicopter) = [(-10,10),(-10,10)]
+inputrange(p::AbstractETHHelicopter)  = [(-10,10),(-10,10)]
+isstable(p::AbstractETHHelicopter)    = false
+isasstable(p::AbstractETHHelicopter)  = false
+sampletime(p::AbstractETHHelicopter)  = p.h
+bias(p::AbstractETHHelicopter)        = p.bias
+
+
+function control(p::ETHHelicopter, u)
+    send(p.control1,u[1])
+    send(p.control2,u[2])
+end
+
+measure(p::ETHHelicopter) = [read(p.measure1), read(p.measure2)]
+control(p::ETHHelicopterSimulator, u)  = error("Not yet implemented")
+measure(p::ETHHelicopterSimulator)     = error("Not yet implemented")
+
+
+initialize(p::ETHHelicopter)          = nothing
+finalize(p::ETHHelicopter)            = foreach(close, p.stream.devices)
+initialize(p::ETHHelicopterSimulator) = nothing
+finalize(p::ETHHelicopterSimulator)   = nothing

src/reference_generators.jl

+export PRBSGenerator
+
+#PRBS
+"""
+    r = PRBSGenerator()
+Generates a pseudo-random binary sequence. Call like `random_input = r()`.
+"""
+mutable struct PRBSGenerator
+    state::Int
+end
+
+PRBSGenerator() = PRBSGenerator(Int(1))
+
+function (r::PRBSGenerator)(args...)
+    state = r.state
+    bit   = ((state >> 0) ⊻ (state >> 2) ⊻ (state >> 3) ⊻ (state >> 5) ) & 1
+    r.state = (state >> 1) | (bit << 15)
+    bit
+end

src/utilities.jl

-export @periodically, init_sysfilter, sysfilter!
+export @periodically, init_sysfilter, sysfilter!, SysFilter
 
 """
 	@periodically(h, body)
@@ -9,26 +9,73 @@ macro periodically(h, body)
 		local start_time = time()
 		$(esc(body))
 		local execution_time = time()-start_time
-		sleep(max(0,$h-execution_time))
+		Libc.systemsleep(max(0,$(esc(h))-execution_time))
 	end
 end
 
+"""
+	@periodically(h, simulation::Bool, body)
+Ensures that the body is run with an interval of `h >= 0.001` seconds.
+If `simulation == false`, no sleep is done
+"""
+macro periodically(h, simulation, body)
+	quote
+		local start_time = time()
+		$(esc(body))
+		local execution_time = time()-start_time
+		$(esc(simulation)) || Libc.systemsleep(max(0,$(esc(h))-execution_time))
+	end
+end
+
+
+"""
+	Csf = SysFilter(sys_discrete::StateSpace)
+	Csf = SysFilter(sys_continuous::StateSpace, sampletime)
+	Csf = SysFilter(sys::StateSpace, state::AbstractVector)
+Returns an object used for filtering signals through LTI systems.
+Create a SysFilter object that can be used to implement control loops and simulators
+with LTI systems, i.e., `U(z) = C(z)E(z)`. To filter a signal `u` through the filter,
+call like `y = Csf(u)`. Calculates the filtered output `y` in `y = Cx+Du, x'=Ax+Bu`
+"""
+struct SysFilter{T<:StateSpace}
+	sys::T
+	state::Vector{Float64}
+	function SysFilter(sys::StateSpace, state::AbstractVector)
+		@assert !ControlSystems.iscontinuous(sys) "Can not filter using continuous time model."
+		@assert length(state) == sys.nx "length(state) != sys.nx"
+		new{typeof(sys)}(sys, state)
+	end
+	function SysFilter(sys::StateSpace)
+		@assert !ControlSystems.iscontinuous(sys) "Can not filter using continuous time model. Supply sample time."
+		new{typeof(sys)}(sys, init_sysfilter(sys))
+	end
+	function SysFilter(sys::StateSpace, h::Real)
+		@assert ControlSystems.iscontinuous(sys) "Sample time supplied byt system model is already in discrete time."
+		sysd = c2d(sys, h)[1]
+		new{typeof(sysd)}(sysd, init_sysfilter(sysd))
+	end
+end
+(s::SysFilter)(input) = sysfilter!(s.state, s.sys, input)
+
 """
 	state = init_sysfilter(sys::StateSpace)
 Use together with [`sysfilter!`](@ref)
 """
 function init_sysfilter(sys::StateSpace)
- 	zeros(sys.nx,1)
+ 	zeros(sys.nx)
 end
 
 """
+	output = sysfilter!(s::SysFilter, input)
 	output = sysfilter!(state, sys::StateSpace, input)
 Returns the filtered output `y` in `y = Cx+Du, x'=Ax+Bu`
 
 This function is used to implement control loops where a signal is filtered through a
 dynamical system, i.e., `U(z) = C(z)E(z)`. Initialize `state` using [`init_sysfilter`](@ref).
 """
-function sysfilter!(state, sys::StateSpace, input)
-	state .= sys.A*state + sys.B*input
-	output = sys.C*state + sys.D*input
+function sysfilter!(state::AbstractVector, sys::StateSpace, input)
+	state .= vec(sys.A*state + sys.B*input)
+	output = vec(sys.C*state + sys.D*input)
 end
+
+sysfilter!(s::SysFilter, input) = sysfilter!(s.state, s.sys, input)

test/runtests.jl

-using LabProcesses
-using Base.Test
+using LabProcesses, ControlSystems, DSP
+using Test
 
-# write your own tests here
-@test 1 == 2
+# Reference generators
+r = PRBSGenerator(Int(4))
+seq = [r() for i = 1:10]
+@test all(seq .==  [1,0,1,0,0,0,0,0,0,0])
+foreach(r,1:10_000)
+
+function test_sysfilter()
+	N     = 10
+	u     = randn(N)
+	b     = [1, 1]
+	a     = [1, 0.1, 1]
+	sys   = ss(tf(b,a,1))
+	state = init_sysfilter(sys)
+	yf    = filt(b,a,u)
+	yff   = similar(yf)
+	for i in eachindex(u)
+		yff[i] = sysfilter!(state, sys, u[i])[1]
+	end
+	@test sum(abs,yf - yff) < √(eps())
+	sysfilt = SysFilter(sys)
+	for i in eachindex(u)
+		yff[i] = sysfilter!(sysfilt, u[i])[1]
+	end
+	@test sum(abs,yf - yff) < √(eps())
+	sysfilt = SysFilter(sys)
+	for i in eachindex(u)
+		yff[i] = sysfilt(u[i])[1]
+	end
+	@test sum(abs,yf - yff) < √(eps())
+end
+
+test_sysfilter()