diff --git a/.travis.yml b/.travis.yml
index 2f95bd7..b231de1 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -3,13 +3,30 @@ language: julia
 os:
   - linux
   - osx
+  - windows
+
+arch:
+  - x64
+  - x86
+  - arm64
 
 julia:
   - 0.5
   - 0.6
+  - 0.7
+  - 1.1
   - nightly
 
 matrix:
+  exclude:
+    - os: osx
+      arch: x86
+    - os: osx
+      arch: arm64
+    - os: windows
+      arch: arm64
+    - julia: nightly
+      arch: arm64
   fast_finish: true
   allow_failures:
     - julia: nightly
@@ -23,11 +40,11 @@ notifications:
     on_failure: change
     on_start: never
 
-script:
-  - julia -e 'Pkg.clone(pwd()); Pkg.build("LTISystems")'
-  - julia --check-bounds=yes --inline=no -e 'Pkg.test("LTISystems",
-    coverage=true)'
+#script:
+#  - julia -e 'Pkg.clone(pwd()); Pkg.build("LTISystems")'
+#  - julia --check-bounds=yes --inline=no -e 'Pkg.test("LTISystems",
+#    coverage=true)'
 
 after_success:
-  - julia -e 'cd(Pkg.dir("LTISystems")); Pkg.add("Coverage"); using Coverage;
+  - julia -e 'using Pkg; Pkg.dir("LTISystems"); Pkg.add("Coverage"); using Coverage;
     Coveralls.submit(process_folder()); Codecov.submit(process_folder())'
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..44ae996
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,20 @@
+name = "LTISystems"
+uuid = "67f030dc-aa52-5f83-ac42-e4024659685c"
+authors = ["Arda Aytekin <aytekin@kth.se>, Niklas Everitt <neveritt@kth.se>"]
+
+[deps]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+PolynomialMatrices = "68ae48f1-ae9d-5916-b516-bab924cb4d98"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+RationalFunctions = "a89fc88f-fc1e-5208-8949-7a3f8ddd21cd"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+
+[compat]
+Compat = "≥ 0.18.0"
+OrdinaryDiffEq = "≥ 1.0.0"
+PolynomialMatrices = "≥ 0.1.1"
+RationalFunctions = "≥ 0.1.1"
+julia = "≥ 0.7.0"
diff --git a/src/LTISystems.jl b/src/LTISystems.jl
index 10364b2..d787e99 100644
--- a/src/LTISystems.jl
+++ b/src/LTISystems.jl
@@ -1,4 +1,4 @@
-__precompile__(true)
+#__precompile__(true)
 
 module LTISystems
 
@@ -17,6 +17,9 @@ using PolynomialMatrices
 # Polynomials-related things
 using Polynomials
 
+# Linear Algebra related things
+using LinearAlgebra
+
 # RationalFunctions-related things
 using RationalFunctions
 import RationalFunctions: poles
@@ -31,16 +34,16 @@ import Base: convert, promote_rule
 import Base: one, zero
 
 # Import num, den for getting numerator and denominator polynomials
-import Base: num, den
+import Base: numerator, denominator
 
 # Import inv and zeros
 import Base: inv, zeros
 
 # Indexing
-import Base: getindex, setindex!, endof
+import Base: getindex, setindex!, firstindex, lastindex
 
 # Import iteration interface functions
-import Base: start, next, done, eltype, length, size
+import Base: iterate, IteratorSize, IteratorEltype, eltype, length, size
 
 # Import printing functions
 import Base: summary
@@ -49,7 +52,7 @@ import Base: summary
 import Base: step
 
 # Import mathematical operations for overloading
-import Base: +, .+, -, .-, *, .*, /, ./, ==, !=, isapprox, transpose
+import Base: +, -, *, /, ==, !=, isapprox, transpose
 
 # Export only the useful functions
 export
diff --git a/src/conversions/mfd2ss.jl b/src/conversions/mfd2ss.jl
index 937457a..2d89e4c 100644
--- a/src/conversions/mfd2ss.jl
+++ b/src/conversions/mfd2ss.jl
@@ -4,7 +4,7 @@
 # realization (Kailath, 1980, Section 6.4.1).
 # NOTE: Base this function on the SLICOT routine TC04AD.
 # NOTE: Is it necessary to make a SISO version of this function?
-function _mfd2ss{S,M1,M2}(mfd::MatrixFractionDescription{Val{:mimo},S,Val{:rfd},M1,M2})
+function _mfd2ss(mfd::MatrixFractionDescription{Val{:mimo},S,Val{:rfd},M1,M2}) where {S,M1,M2}
   Den, Num  = colred(mfd.D, mfd.N)
   kden, Phc = high_col_deg_matrix(Den)
   knum      = col_degree(Num)
@@ -82,7 +82,7 @@ function _mfd2ss{S,M1,M2}(mfd::MatrixFractionDescription{Val{:mimo},S,Val{:rfd},
   return A, B, C, D
 end
 
-function _mfd2ss{S,M1,M2}(mfd::MatrixFractionDescription{Val{:mimo},S,Val{:lfd},M1,M2})
+function _mfd2ss(mfd::MatrixFractionDescription{Val{:mimo},S,Val{:lfd},M1,M2}) where {S,M1,M2}
   Den, Num  = rowred(mfd.D,mfd.N)
   kden, Phr = high_row_deg_matrix(Den)
   knum      = row_degree(Num)
diff --git a/src/conversions/ss2mfd.jl b/src/conversions/ss2mfd.jl
index 8f3bf00..58b7594 100644
--- a/src/conversions/ss2mfd.jl
+++ b/src/conversions/ss2mfd.jl
@@ -1,12 +1,12 @@
 # Generates a right MatrixFractionDescription whose transfer function coincides with that
 # of an observable state-space model.
-function _ss2rfd{T,S}(sys::StateSpace{Val{T},Val{S}}, var::Symbol)
+function _ss2rfd(sys::StateSpace{Val{T},Val{S}}, var::Symbol) where {T,S}
   n     = numstates(sys)
   nu    = numinputs(sys)
   ny    = numoutputs(sys)
 
   # Construct system matrix (Rosenbrock, 1970)
-  Dₗ    = PolyMatrix(vcat(sys.A, -eye(sys.A)), (n,n), Val{_pmvaltype(sys)})
+  Dₗ    = PolyMatrix(vcat(sys.A, -I(sys.A)), (n,n), Val{_pmvaltype(sys)})
   Nₗ    = PolyMatrix(sys.B, (n,nu), Val{_pmvaltype(sys)})
 
   # Write (sI-A)^{-1} B as a left MatrixFractionDescription and convert to right MatrixFractionDescription
@@ -33,13 +33,13 @@ rfd(sys::StateSpace{Val{:mimo},Val{:cont}})  =
 rfd(sys::StateSpace{Val{:mimo},Val{:disc}}) =
   rfd(_ss2rfd(sys,:z)...,samplingtime(sys))
 
-function _ss2lfd{T,S}(sys::StateSpace{Val{T},Val{S}}, var::Symbol)
+function _ss2lfd(sys::StateSpace{Val{T},Val{S}}, var::Symbol) where {T,S}
   n     = numstates(sys)
   nu    = numinputs(sys)
   ny    = numoutputs(sys)
 
   # Construct system matrix (Rosenbrock, 1970)
-  Dᵣ    = PolyMatrix(vcat(-sys.A, eye(sys.A)), (n,n), Val{_pmvaltype(sys)})
+  Dᵣ    = PolyMatrix(vcat(-sys.A, I(sys.A)), (n,n), Val{_pmvaltype(sys)})
   Nᵣ    = PolyMatrix(sys.C, (ny,n), Val{_pmvaltype(sys)})
   # Write C(sI-A)^{-1} as a right MatrixFractionDescription and convert to left MatrixFractionDescription
   R,U   = rtriang(vcat(-Dᵣ, Nᵣ))
@@ -64,5 +64,5 @@ lfd(sys::StateSpace{Val{:mimo},Val{:cont}})  =
 lfd(sys::StateSpace{Val{:mimo},Val{:disc}}) =
   lfd(_ss2lfd(sys,:z)...,samplingtime(sys))
 
-_pmvaltype{T}(s::LtiSystem{Val{T},Val{:cont}}) = :s
-_pmvaltype{T}(s::LtiSystem{Val{T},Val{:disc}}) = :z
+_pmvaltype(s::LtiSystem{Val{T},Val{:cont}}) where {T} = :s
+_pmvaltype(s::LtiSystem{Val{T},Val{:disc}}) where {T} = :z
diff --git a/src/conversions/tf2mfd.jl b/src/conversions/tf2mfd.jl
index 57026e6..25a0a4f 100644
--- a/src/conversions/tf2mfd.jl
+++ b/src/conversions/tf2mfd.jl
@@ -8,47 +8,47 @@
 
 # s = R(s)*inv(D(s)) where D(s) is diagonal with the product of all denominators
 # of column i in D[i,i]
-# R[i,j] is num(s[i,j]) times all denominators in den(s[:,j]) except den(s[i,j])
-function _tf2rfd(s::LTISystems.TransferFunction)
+# R[i,j] is num(s[i,j]) times all denominators in denominator(s[:,j]) except denominator(s[i,j])
+function _tf2rfd(s::TransferFunction)
   mat = s.mat
   n,m = size(mat)
-  dp  = fill(zero(den(mat[1])), m)
-  R   = map(num, mat)
-  map(num,s.mat)
+  dp  = fill(zero(denominator(mat[1])), m)
+  R   = map(numerator, mat)
+  map(numerator,s.mat)
   for i in 1:m
-    dp[i]   = Base.reduce(*, one(den(s.mat[1,i])), map(den,s.mat[:,i]))
+    dp[i] = Base.reduce(*, map(denominator,s.mat[:,i]), init=one(denominator(s.mat[1,i])))
     for j in 1:n
       for k in 1:n
         k == j && continue
-        R[j,i] *= den(s.mat[k,i])
+        R[j,i] *= denominator(s.mat[k,i])
       end
     end
   end
   R = PolyMatrix(R)
-  D = PolyMatrix(diagm(dp))
+  D = PolyMatrix(Matrix{eltype(dp)}(Diagonal(dp)))
   R, D
 end
 
 # s = inv(D(s))*R(s) where D(s) is diagonal with the product of all denominators
 # of column i in D[i,i]
-# R[i,j] is num(s[i,j]) times all denominators in den(s[j,:]) except den(s[i,j])
-function _tf2lfd(s::LTISystems.TransferFunction)
+# R[i,j] is num(s[i,j]) times all denominators in denominator(s[j,:]) except denominator(s[i,j])
+function _tf2lfd(s::TransferFunction)
   mat = s.mat
   n,m = size(mat)
-  dp  = fill(zero(den(mat[1])), n)
-  R   = map(num, mat)
-  map(num,s.mat)
+  dp  = fill(zero(denominator(mat[1])), n)
+  R   = map(numerator, mat)
+  map(numerator,s.mat)
   for i in 1:n
-    dp[i] = Base.reduce(*, one(den(s.mat[i,1])), map(den,s.mat[i,:]))
+    dp[i] = Base.reduce(*, map(denominator,s.mat[i,:]), init=one(denominator(s.mat[i,1])))
     for j in 1:m
       for k in 1:m
         k == j && continue
-        R[i,j] *= den(s.mat[i,k])
+        R[i,j] *= denominator(s.mat[i,k])
       end
     end
   end
   R = PolyMatrix(R)
-  D = PolyMatrix(diagm(dp))
+  D = PolyMatrix(Matrix{eltype(dp)}(Diagonal(dp)))
   R, D
 end
 
diff --git a/src/conversions/tf2ss.jl b/src/conversions/tf2ss.jl
index 856b2a8..4adb1b1 100644
--- a/src/conversions/tf2ss.jl
+++ b/src/conversions/tf2ss.jl
@@ -8,8 +8,8 @@ ss(s::TransferFunction{Val{:siso},Val{:disc}}) = minreal(ss(_tf2ss(s)..., s.Ts))
 # TODO num(s) instead of num(s.mat[1])
 # controllable canonical form
 function _tf2ss(s::TransferFunction{Val{:siso}})
-  nump  = RationalFunctions.num(s.mat[1])
-  denp  = RationalFunctions.den(s.mat[1])
+  nump  = RationalFunctions.numerator(s.mat[1])
+  denp  = RationalFunctions.denominator(s.mat[1])
   nump /= denp[end]
   denp /= denp[end]
   m     = degree(nump)
diff --git a/src/methods/feedback.jl b/src/methods/feedback.jl
index 3edea1d..73da664 100644
--- a/src/methods/feedback.jl
+++ b/src/methods/feedback.jl
@@ -18,18 +18,18 @@ Continuous time rational transfer function model
 x^2 + 4x + 7
 ```
 """
-function feedback{T1<:LtiSystem, T2<:LtiSystem}(s1::T1, s2::T2)
+function feedback(s1::T1, s2::T2) where {T1<:LtiSystem, T2<:LtiSystem}
   /(s1, parallel(one(s1), series(s1,s2)))
 end
 
-feedback{T1<:LtiSystem, T2<:Real}(s1::T1, n::T2) =
+feedback(s1::T1, n::T2) where {T1<:LtiSystem, T2<:Real} =
   /(s1, parallel(one(s1), series(s1, n)))
-feedback{T1<:LtiSystem, T2<:Real}(n::T2, s1::T1) =
+feedback(n::T2, s1::T1) where {T1<:LtiSystem, T2<:Real} =
   /(n, parallel(one(s1), series(n, s1)))
 
-feedback{T1<:LtiSystem, T2<:Real}(s1::T1, n::Matrix{T2}) =
+feedback(s1::T1, n::Matrix{T2}) where {T1<:LtiSystem, T2<:Real} =
   /(s1, parallel(one(s1), series(s1, n)))
-feedback{T1<:LtiSystem, T2<:Real}(n::Matrix{T2}, s1::T1) =
+feedback(n::Matrix{T2}, s1::T1) where {T1<:LtiSystem, T2<:Real} =
   /(n, parallel(one(s1), series(n, s1)))
 
 # # TODO: Implement the specific cases for speedup purposes. Take into account the
diff --git a/src/methods/freqresp.jl b/src/methods/freqresp.jl
index 94684c2..32bad45 100644
--- a/src/methods/freqresp.jl
+++ b/src/methods/freqresp.jl
@@ -28,13 +28,13 @@ systems.
 
 **See also:** `bode`, `nyquist`, `samplingtime`.
 """
-freqresp{S}(sys::LtiSystem{Val{S},Val{:cont}}, x::Real)                     =
+freqresp(sys::LtiSystem{Val{S},Val{:cont}}, x::Real) where {S}                     =
   sys(1im*x)
-freqresp{S,M<:Real}(sys::LtiSystem{Val{S},Val{:cont}}, X::AbstractArray{M}) =
+freqresp(sys::LtiSystem{Val{S},Val{:cont}}, X::AbstractArray{M}) where {S,M<:Real} =
   sys(1im*X)
-freqresp{S}(sys::LtiSystem{Val{S},Val{:disc}}, x::Real)                     =
+freqresp(sys::LtiSystem{Val{S},Val{:disc}}, x::Real) where {S}                     =
   sys(exp(1im*samplingtime(sys)*x))
-freqresp{S,M<:Real}(sys::LtiSystem{Val{S},Val{:disc}}, X::AbstractArray{M}) =
+freqresp(sys::LtiSystem{Val{S},Val{:disc}}, X::AbstractArray{M}) where {S,M<:Real} =
   sys(exp(1im*samplingtime(sys)*X))
 
 # Implements algorithm found in:
@@ -45,7 +45,7 @@ function _eval(sys::StateSpace, x::Number)
   convert(AbstractMatrix{Complex128}, sys.D + sys.C*(R\sys.B))
 end
 
-function _eval{M<:Number}(sys::StateSpace, X::AbstractArray{M})
+function _eval(sys::StateSpace, X::AbstractArray{M}) where {M<:Number}
   resp = Array{Complex128}(size(sys)..., size(X)...)
   for idx in CartesianRange(indices(X))
     resp[:, :, idx] = sys.D + sys.C*((X[idx]*I - sys.A)\sys.B)
diff --git a/src/methods/minreal.jl b/src/methods/minreal.jl
index 27d90dc..e848331 100644
--- a/src/methods/minreal.jl
+++ b/src/methods/minreal.jl
@@ -60,8 +60,8 @@ julia> Dm
  0.0
 ```
 """
-function minreal{M1<:AbstractMatrix,M2<:AbstractMatrix,M3<:AbstractMatrix,
-  M4<:AbstractMatrix}(A::M1, B::M2, C::M3, D::M4, tol::Float64 = zero(Float64))
+function minreal(A::M1, B::M2, C::M3, D::M4, tol::Float64 = zero(Float64)) where {
+  M1<:AbstractMatrix,M2<:AbstractMatrix,M3<:AbstractMatrix,M4<:AbstractMatrix}
   @assert size(A,1) == size(A,2) "minreal: A must be square"
   @assert size(A,1) == size(B,1) "minreal: A and B must have the same row size"
   @assert size(A,1) == size(C,2) "minreal: A and C must have the same column size"
@@ -72,7 +72,7 @@ function minreal{M1<:AbstractMatrix,M2<:AbstractMatrix,M3<:AbstractMatrix,
   Tc, c = rosenbrock(A, B, tol)
   ATemp = (Tc'*A*Tc)[n-c+1:n,n-c+1:n]
   CTemp = (C*Tc)[:,n-c+1:n]
-  To, o = rosenbrock(ATemp.', CTemp.', tol)
+  To, o = rosenbrock(transpose(ATemp), transpose(CTemp), tol)
 
   T     = view(Tc, :, n-c+1:n)*To
 
diff --git a/src/methods/parallel.jl b/src/methods/parallel.jl
index 39dcfa8..9cffe20 100644
--- a/src/methods/parallel.jl
+++ b/src/methods/parallel.jl
@@ -18,10 +18,10 @@ Continuous time rational transfer function model
 x + 2
 ```
 """
-parallel{T1<:LtiSystem, T2<:LtiSystem}(s1::T1, s2::T2)     = +(s1, s2)
+parallel(s1::T1, s2::T2) where {T1<:LtiSystem, T2<:LtiSystem}     = +(s1, s2)
 
-parallel{T1<:LtiSystem, T2<:Real}(s1::T1, n::T2)           = +(s1, n)
-parallel{T1<:LtiSystem, T2<:Real}(n::T2, s1::T1)           = +(n, s1)
+parallel(s1::T1, n::T2) where {T1<:LtiSystem, T2<:Real}           = +(s1, n)
+parallel(n::T2, s1::T1) where {T1<:LtiSystem, T2<:Real}           = +(n, s1)
 
-parallel{T1<:LtiSystem, M2<:AbstractMatrix}(s1::T1, n::M2) = +(s1, n)
-parallel{T1<:LtiSystem, M2<:AbstractMatrix}(n::M2, s1::T1) = +(n, s1)
+parallel(s1::T1, n::M2) where {T1<:LtiSystem, M2<:AbstractMatrix} = +(s1, n)
+parallel(n::M2, s1::T1) where {T1<:LtiSystem, M2<:AbstractMatrix} = +(n, s1)
diff --git a/src/methods/pzmap.jl b/src/methods/pzmap.jl
index 1892808..7f6846c 100644
--- a/src/methods/pzmap.jl
+++ b/src/methods/pzmap.jl
@@ -1,10 +1,10 @@
-immutable PZData{T}
-  p::Vector{Complex128}
-  z::Vector{Complex128}
+struct PZData{T}
+  p::Vector{Complex{Float64}}
+  z::Vector{Complex{Float64}}
 
-  function (::Type{PZData{Val{U}}}){T<:Number,S<:Number,U}(p::AbstractVector{T},
-    z::AbstractVector{S})
-    new{Val{U}}(convert(Vector{Complex128}, p), convert(Vector{Complex128}, z))
+  function (::Type{PZData{Val{U}}})(p::AbstractVector{T},
+    z::AbstractVector{S}) where {T<:Number,S<:Number,U}
+    new{Val{U}}(convert(Vector{Complex{Float64}}, p), convert(Vector{Complex{Float64}}, z))
   end
 end
 
@@ -37,7 +37,7 @@ However, plotting recipes are defined for `pzdata`, which allows for
 
 **See also:** `poles`, `zeros`, `tzeros`.
 """
-pzmap{T,S}(sys::LtiSystem{Val{T},Val{S}}) = PZData{Val{S}}(_pzmap(sys)...)
+pzmap(sys::LtiSystem{Val{T},Val{S}}) where {T,S} = PZData{Val{S}}(_pzmap(sys)...)
 
 @recipe function f{T}(pzdata::PZData{Val{T}}; stability = true)
   @assert isa(stability, Bool) "plot(pzdata): `stability` must be a `Bool`"
diff --git a/src/methods/reduce.jl b/src/methods/reduce.jl
index 899c50c..60820fb 100644
--- a/src/methods/reduce.jl
+++ b/src/methods/reduce.jl
@@ -1,13 +1,13 @@
-function reduce{M1<:AbstractMatrix, M2<:AbstractMatrix, M3<:AbstractMatrix,
-  M4<:AbstractMatrix}(A::M1, B::M2, C::M3, D::M4, tol::Float64 = zero(Float64))
+function reduce(A::M1, B::M2, C::M3, D::M4, tol::Float64 = zero(Float64)) where {M1<:AbstractMatrix, M2<:AbstractMatrix, M3<:AbstractMatrix,
+  M4<:AbstractMatrix}
   @assert size(A,1) == size(A,2) "reduce: A must be square"
   @assert size(A,1) == size(B,1) "reduce: A and B must have the same row size"
   @assert size(A,1) == size(C,2) "reduce: A and C must have the same column size"
   @assert size(B,2) == size(D,2) "reduce: B and D must have the same column size"
   @assert size(C,1) == size(D,1) "reduce: C and D must have the same row size"
 
-  n1, n2  = size(B,1,2)
-  n3, n4  = size(C,1,2)
+  n1, n2  = size(B)
+  n3, n4  = size(C)
 
   ν₁, ν₂  = n1, n1
   σ₁, σ₂  = n3, n3
@@ -27,8 +27,8 @@ function reduce{M1<:AbstractMatrix, M2<:AbstractMatrix, M3<:AbstractMatrix,
   tol     = max(tol, 10*max(n3,n2+n4)*norm(temp1,1)*eps(Float64) + eps(Float64))
 
   while τ ≠ 0 && ρ ≠ 0 && ν₂ ≠ 0
-    n3, n4  = size(Ctemp,1,2)
-    svdobj  = svdfact(Dtemp, thin = false)
+    n3, n4  = size(Ctemp)
+    svdobj  = svd(Dtemp, full = true)
     σ₂      = sum(svdobj.S .>= tol)
     τ       = n3 - σ₂
 
@@ -38,14 +38,14 @@ function reduce{M1<:AbstractMatrix, M2<:AbstractMatrix, M3<:AbstractMatrix,
     D̄       = temp2[1:σ₂, n4+1:end]
 
     if τ ≠ 0
-      svdobj  = svdfact(view(temp2, σ₂+1:n3, 1:n4), thin = false)
+      svdobj  = svd(view(temp2, σ₂+1:n3, 1:n4), full = true)
       ρ       = sum(svdobj.S .>= tol)
       ν₂      = n4 - ρ
 
       if ρ ≠ 0 && ν₂ ≠ 0
         μ     = ρ + σ₂
         δ     += ρ
-        V     = flipdim(svdobj.Vt', 2)
+        V     = reverse(svdobj.Vt', dims=2)
         temp3 = [V'*Atemp*V ; C̄*V   ]
         temp4 = [V'*Btemp   ; D̄ ]
 
diff --git a/src/methods/rosenbrock.jl b/src/methods/rosenbrock.jl
index ea382f4..092a510 100644
--- a/src/methods/rosenbrock.jl
+++ b/src/methods/rosenbrock.jl
@@ -59,14 +59,13 @@ julia> Cc = C*Tc; Dc = D;
 -  [2]: H.H. Rosenbrock, State-Space and Multivariable Theory. London: Thomas
         Nelson and Sons, 1970.
 """
-function rosenbrock{M1<:AbstractMatrix,M2<:AbstractMatrix}(A::M1, B::M2,
-  tol::Float64 = zero(Float64))
-  n1, n2  = size(A,1,2)
-  n3, n4  = size(B,1,2)
+function rosenbrock(A::M1, B::M2, tol::Float64 = zero(Float64)) where {M1<:AbstractMatrix,M2<:AbstractMatrix}
+  n1, n2  = size(A)
+  n3, n4  = size(B)
   @assert n1 == n2 "rosenbrock: A must be square"
   @assert n1 == n3 "rosenbrock: A and B must have same row sizes"
   c::Int              = 0
-  T::Matrix{Float64}  = eye(n1)
+  T::Matrix{Float64}  = I(n1)
   τ₁, τ₂  = (n1,1)
   ρ₁, ρ₂  = (0,1)
   Atemp   = A
@@ -74,8 +73,8 @@ function rosenbrock{M1<:AbstractMatrix,M2<:AbstractMatrix}(A::M1, B::M2,
   tol     = max(tol, 10*max(n3,n4)*norm(B,1)*eps(Float64) + eps(Float64))
   while ρ₂ != 0 && τ₂ != 0
     n       = size(Btemp,1)
-    svdobj  = svdfact(Btemp, thin = false)
-    U       = flipdim(svdobj.U, 2)
+    svdobj  = svd(Btemp, full = true)
+    U       = reverse(svdobj.U, dims=2)
     ρ₂      = sum(svdobj.S .>= tol)
     τ₂      = n - ρ₂
 
@@ -84,7 +83,7 @@ function rosenbrock{M1<:AbstractMatrix,M2<:AbstractMatrix}(A::M1, B::M2,
       Btemp   = view(Atemp, 1:τ₂, τ₂+1:n)
       Atemp   = view(Atemp, 1:τ₂, 1:τ₂)
       T[:]    = T*vcat(hcat(U, zeros(size(U, 1),c)),
-                        hcat(zeros(c, size(U, 2)), eye(c)))
+                        hcat(zeros(c, size(U, 2)), I(c)))
 
       c       = c + ρ₂
       ρ₁, τ₁  = (ρ₂, τ₂)
diff --git a/src/methods/series.jl b/src/methods/series.jl
index 5a555d0..68c177d 100644
--- a/src/methods/series.jl
+++ b/src/methods/series.jl
@@ -18,10 +18,10 @@ Continuous time rational transfer function model
 x^2 + 4x + 4
 ```
 """
-series{T1<:LtiSystem, T2<:LtiSystem}(s1::T1, s2::T2)     = *(s1,s2)
+series(s1::T1, s2::T2) where {T1<:LtiSystem, T2<:LtiSystem}     = *(s1,s2)
 
-series{T1<:LtiSystem, T2<:Real}(s1::T1, n::T2)           = *(s1, n)
-series{T1<:LtiSystem, T2<:Real}(n::T2, s1::T1)           = *(n ,s1)
+series(s1::T1, n::T2) where {T1<:LtiSystem, T2<:Real}           = *(s1, n)
+series(n::T2, s1::T1) where {T1<:LtiSystem, T2<:Real}           = *(n ,s1)
 
-series{T1<:LtiSystem, M2<:AbstractMatrix}(s1::T1, n::M2) = *(s1, n)
-series{T1<:LtiSystem, M2<:AbstractMatrix}(n::M2, s1::T1) = *(n ,s1)
+series(s1::T1, n::M2) where {T1<:LtiSystem, M2<:AbstractMatrix} = *(s1, n)
+series(n::M2, s1::T1) where {T1<:LtiSystem, M2<:AbstractMatrix} = *(n ,s1)
diff --git a/src/methods/simulate.jl b/src/methods/simulate.jl
index aa9afb3..96c1f69 100644
--- a/src/methods/simulate.jl
+++ b/src/methods/simulate.jl
@@ -1,7 +1,7 @@
 # TODO: For now, keep SimType{T} in `dense` format. Later on, think of relaxing
 #       this thing, *i.e.*, SimType{T<:Real,V1<:AbstractVector{T},etc.}.
 
-type SimType{T<:Real} <: DEDataVector{T}
+mutable struct SimType{T<:Real} <: DEDataVector{T}
   x::Vector{T}
   y::Vector{T}
   u::Vector{T}
@@ -12,22 +12,22 @@ function (::Type{SimType})(x::Vector, y::Vector, u::Vector)
   SimType(convert(Vector{T}, x), convert(Vector{T}, y), convert(Vector{T}, u))
 end
 
-function simulate{T}(sys::LtiSystem{Val{T},Val{:cont}}, tspan;
-  input = (t,x)->zeros(numinputs(sys)), alg::AbstractODEAlgorithm = Tsit5(),
-  initial::AbstractVector = zeros(numstates(sys)), tstops = Float64[], kwargs...)
+function simulate(sys::LtiSystem{Val{T},Val{:cont}}, tspan;
+  input = (t,x)->zeros(numinputs(sys)), alg::DiffEqBase.AbstractODEAlgorithm = Tsit5(),
+  initial::AbstractVector = zeros(numstates(sys)), tstops = Float64[], kwargs...) where {T}
 
-  f     = (t,x,dx)->sys(t,x,dx,input)
-  simvar= SimType(initial, zeros(numoutputs(sys)), zeros(numinputs(sys)))
-  f(tspan[1], simvar, zeros(initial))
+  f      = (dx,x,p,t)->sys(t,x,dx,input)
+  simvar = SimType(initial, zeros(numoutputs(sys)), zeros(numinputs(sys)))
+  f(zero(initial), simvar, 0, tspan[1])
 
-  tstops= [tstops..., discontinuities(input, tspan)...]
-  prob  = ODEProblem(f, simvar, tspan)
-  sln   = OrdinaryDiffEq.solve(prob, alg; tstops = tstops, kwargs...)
+  tstops = [tstops..., discontinuities(input, tspan)...]
+  prob   = ODEProblem(f, simvar, tspan)
+  sln    = OrdinaryDiffEq.solve(prob, alg; tstops = tstops, kwargs...)
 
-  tvals = sln.t
-  xvals = Matrix{eltype(sln[1].x)}(length(tvals), numstates(sys))
-  yvals = Matrix{eltype(sln[1].y)}(length(tvals), numoutputs(sys))
-  uvals = Matrix{eltype(sln[1].u)}(length(tvals), numinputs(sys))
+  tvals  = sln.t
+  xvals  = Matrix{eltype(sln[1].x)}(undef, length(tvals), numstates(sys))
+  yvals  = Matrix{eltype(sln[1].y)}(undef, length(tvals), numoutputs(sys))
+  uvals  = Matrix{eltype(sln[1].u)}(undef, length(tvals), numinputs(sys))
   for idx in 1:length(tvals)
     xvals[idx,:] = sln[idx]
     yvals[idx,:] = sln[idx].y
@@ -37,14 +37,14 @@ function simulate{T}(sys::LtiSystem{Val{T},Val{:cont}}, tspan;
   TimeResponse(tvals,xvals,yvals,uvals)
 end
 
-function simulate{T}(sys::LtiSystem{Val{T},Val{:disc}}, tspan;
+function simulate(sys::LtiSystem{Val{T},Val{:disc}}, tspan;
   input = (t,x)->zeros(numinputs(sys)),
-  initial::AbstractVector = zeros(numstates(sys)), kwargs...)
+  initial::AbstractVector = zeros(numstates(sys)), kwargs...) where {T}
 
-  f      = (t,x,dx)->sys(t,x,dx,input)
+  f      = (dx,x,p,t)->sys(t,x,dx,input)
   simvar = SimType(initial, zeros(numoutputs(sys)), zeros(numinputs(sys)))
-  temp   = zeros(simvar.x)
-  f(0.0, simvar, temp)  # ensure first step is correct
+  temp   = zero(simvar.x)
+  f(temp, simvar, 0, 0.0)  # ensure first step is correct
   simvar.x = temp       # alternative is adding extra time step, see below
 
   dt    = samplingtime(sys)
@@ -53,9 +53,9 @@ function simulate{T}(sys::LtiSystem{Val{T},Val{:disc}}, tspan;
     kwargs..., dt = dt)
 
   tvals = sln.t
-  xvals = Matrix{eltype(sln[1].x)}(length(tvals), numstates(sys))
-  yvals = Matrix{eltype(sln[1].y)}(length(tvals), numoutputs(sys))
-  uvals = Matrix{eltype(sln[1].u)}(length(tvals), numinputs(sys))
+  xvals = Matrix{eltype(sln[1].x)}(undef, length(tvals), numstates(sys))
+  yvals = Matrix{eltype(sln[1].y)}(undef, length(tvals), numoutputs(sys))
+  uvals = Matrix{eltype(sln[1].u)}(undef, length(tvals), numinputs(sys))
   for idx in 1:length(tvals)
     xvals[idx,:] = sln[idx]
     yvals[idx,:] = sln[idx].y
diff --git a/src/methods/unwrap.jl b/src/methods/unwrap.jl
index ea2dc53..958b1f5 100644
--- a/src/methods/unwrap.jl
+++ b/src/methods/unwrap.jl
@@ -11,8 +11,8 @@ place and returns the modified version, whereas `unwrap` does not modify `Θ₀`
 
 **See also:* `freqresp`, `bode`, `nyquist`.
 """
-function unwrap!{T<:AbstractFloat}(Θ::AbstractArray{T}, dim::Integer = ndims(Θ);
-  tol::Real = π)
+function unwrap!(Θ::AbstractArray{T}, dim::Integer = ndims(Θ);
+  tol::Real = π) where {T<:AbstractFloat}
   size(Θ, dim) < 2 && return Θ
   settol  = max(float(tol), π)
   setrng  = 2settol
@@ -26,5 +26,5 @@ function unwrap!{T<:AbstractFloat}(Θ::AbstractArray{T}, dim::Integer = ndims(Θ
   return Θ
 end
 
-unwrap{T<:Real}(Θ::AbstractArray{T}, dim::Integer = ndims(Θ);
-  tol::Real = π) = unwrap!(map(float,Θ), dim; tol = tol)
+unwrap(Θ::AbstractArray{T}, dim::Integer = ndims(Θ);
+  tol::Real = π) where {T<:Real} = unwrap!(map(float,Θ), dim; tol = tol)
diff --git a/src/types/response/boderesponse.jl b/src/types/response/boderesponse.jl
index fe00225..e42cad6 100644
--- a/src/types/response/boderesponse.jl
+++ b/src/types/response/boderesponse.jl
@@ -1,10 +1,10 @@
-immutable BodeResponse{T<:Real} <: SystemResponse
+struct BodeResponse{T<:Real} <: SystemResponse
   freqs::Vector{T}  # rad/sec
   mag::Array{T,3}   # abs, i.e., |G|
   phase::Array{T,3} # radians
 
-  function (::Type{BodeResponse}){T1<:Real,T2<:Real,T3<:Real}(
-    freqs::AbstractVector{T1}, mag::AbstractArray{T2,3}, phase::AbstractArray{T3,3})
+  function (::Type{BodeResponse})(freqs::AbstractVector{T1}, mag::AbstractArray{T2,3},
+    phase::AbstractArray{T3,3}) where {T1<:Real,T2<:Real,T3<:Real}
     if size(mag) ≠ size(phase)
       warn("BodeResponse(freqs, mag, phase): `mag` and `phase` must have same dimensions")
       throw(DomainError())
@@ -22,9 +22,9 @@ immutable BodeResponse{T<:Real} <: SystemResponse
   end
 end
 
-_bode{T<:Real}(sys::LtiSystem{Val{:siso}}, ω::AbstractVector{T})  =
+_bode(sys::LtiSystem{Val{:siso}}, ω::AbstractVector{T}) where {T<:Real} =
   reshape(freqresp(sys, ω), size(sys)..., length(ω))
-_bode{T<:Real}(sys::LtiSystem{Val{:mimo}}, ω::AbstractVector{T})  =
+_bode(sys::LtiSystem{Val{:mimo}}, ω::AbstractVector{T}) where {T<:Real} =
   freqresp(sys, ω)
 
 """
@@ -58,13 +58,13 @@ when `using Plots`.
 
 **See also:** `freqresp`, `nyquist`.
 """
-function bode{T}(sys::LtiSystem, ω::AbstractVector{T})
+function bode(sys::LtiSystem, ω::AbstractVector{T}) where {T}
   # TODO: should we check for ω ≥ 0 ?
   fr = _bode(sys, ω)
   BodeResponse(ω, abs(fr), unwrap!(angle(fr)))
 end
 
-bode{T}(sys::LtiSystem{Val{T},Val{:disc}}) =
+bode(sys::LtiSystem{Val{T},Val{:disc}}) where {T} =
   bode(sys, logspace(-6, log10(π/samplingtime(sys)), 1000))
 
 @recipe function f(br::BodeResponse; iopairs = Tuple{Int,Int}[], freqs = :rads,
diff --git a/src/types/response/nyquistresponse.jl b/src/types/response/nyquistresponse.jl
index 7565b90..ed01131 100644
--- a/src/types/response/nyquistresponse.jl
+++ b/src/types/response/nyquistresponse.jl
@@ -1,10 +1,10 @@
-immutable NyquistResponse{T<:Real} <: SystemResponse
+struct NyquistResponse{T<:Real} <: SystemResponse
   freqs::Vector{T}  # rad/sec
   real::Array{T,3}   # -
   imag::Array{T,3} # -
 
-  function (::Type{NyquistResponse}){T1<:Real,T2<:Real,T3<:Real}(
-    freqs::AbstractVector{T1}, real::AbstractArray{T2,3}, imag::AbstractArray{T3,3})
+  function (::Type{NyquistResponse})(freqs::AbstractVector{T1}, real::AbstractArray{T2,3},
+    imag::AbstractArray{T3,3}) where {T1<:Real,T2<:Real,T3<:Real}
     if size(real) ≠ size(imag)
       warn("NyquistResponse(freqs, real, imag): `real` and `imag` must have same dimensions")
       throw(DomainError())
@@ -19,9 +19,9 @@ immutable NyquistResponse{T<:Real} <: SystemResponse
   end
 end
 
-_nyquist{T<:Real}(sys::LtiSystem{Val{:siso}}, ω::AbstractVector{T})  =
+_nyquist(sys::LtiSystem{Val{:siso}}, ω::AbstractVector{T}) where {T<:Real} =
   reshape(freqresp(sys, ω), size(sys)..., length(ω))
-_nyquist{T<:Real}(sys::LtiSystem{Val{:mimo}}, ω::AbstractVector{T})  =
+_nyquist(sys::LtiSystem{Val{:mimo}}, ω::AbstractVector{T}) where {T<:Real} =
   freqresp(sys, ω)
 
 """
@@ -51,13 +51,13 @@ when `using Plots`.
 
 **See also:** `freqresp`, `bode`.
 """
-function nyquist{T}(sys::LtiSystem, ω::AbstractVector{T})
+function nyquist(sys::LtiSystem, ω::AbstractVector{T}) where {T}
   # TODO: should we check for ω ≥ 0 ?
   fr = _nyquist(sys, ω)
   NyquistResponse(ω, real(fr), imag(fr))
 end
 
-nyquist{T}(sys::LtiSystem{Val{T},Val{:disc}}) =
+nyquist(sys::LtiSystem{Val{T},Val{:disc}}) where {T} =
   nyquist(sys, logspace(-6, log10(π/samplingtime(sys)), 1000))
 
 @recipe function f(nr::NyquistResponse; iopairs = Tuple{Int,Int}[], freqs = :rads,
@@ -151,25 +151,25 @@ nyquist{T}(sys::LtiSystem{Val{T},Val{:disc}}) =
 end
 
 # args specified in degrees
-function _iso_phases{T1<:Real,T2<:Real}(args::AbstractVector{T1},
-  w::AbstractVector{T2}=linspace(0, 2π, 200))
+function _iso_phases(args::AbstractVector{T1},
+  w::AbstractVector{T2}=linspace(0, 2π, 200)) where {T1<:Real,T2<:Real}
 
   N = tan(args*π/180)
   radius = sqrt(1+N.^2)./2N
   xc = -1/2
-  yc = 1./2N
+  yc = 1/2N
 
   c = cos(w)
   s = sin(w)
 
-  reals = xc*ones(c*radius.') + c*radius.'
-  imags = ones(s)*yc.' + s*radius.'
+  reals = xc*ones(c*transpose(radius)) + c*transpose(radius)
+  imags = ones(s)*transpose(yc) + s*transpose(radius)
   reals, imags
 end
 
 # modules specified in db
-function _iso_modules{T1<:Real,T2<:Real}(modules::AbstractVector{T1},
-  w::AbstractVector{T2}=linspace(0, 2π, 200))
+function _iso_modules(modules::AbstractVector{T1},
+  w::AbstractVector{T2}=linspace(0, 2π, 200)) where {T1<:Real,T2<:Real}
 
   M = exp10(modules/20)
   abs2(M)
@@ -181,7 +181,7 @@ function _iso_modules{T1<:Real,T2<:Real}(modules::AbstractVector{T1},
   c = cos(w)
   s = sin(w)
 
-  reals = ones(c)*xc.' + c*radius.'
-  imags = yc*ones(s*radius.') + s*radius.'
+  reals = ones(c)*transpose(xc) + c*transpose(radius)
+  imags = yc*ones(s*transpose(radius)) + s*transpose(radius)
   reals, imags
 end
diff --git a/src/types/response/timeresponse.jl b/src/types/response/timeresponse.jl
index bb77c82..b930f5c 100644
--- a/src/types/response/timeresponse.jl
+++ b/src/types/response/timeresponse.jl
@@ -1,12 +1,12 @@
-type TimeResponse{T,S} <: SystemResponse
+mutable struct TimeResponse{T,S} <: SystemResponse
   t::Vector{T}
   x::Matrix{T}
   y::Matrix{T}
   u::Matrix{T}
   tsloc::Int
 
-  function (::Type{TimeResponse}){S}(t::Vector, x::Matrix, y::Matrix, u::Matrix,
-    ::Type{Val{S}} = Val{:cont})
+  function (::Type{TimeResponse})(t::Vector, x::Matrix, y::Matrix, u::Matrix,
+    ::Type{Val{S}} = Val{:cont}) where {S}
     (length(t) == size(x, 1) == size(y, 1) == size(u, 1)) ||
       DimensionMismatch("TimeResponse(t,x,y,u): dimension mistmatch")
 
diff --git a/src/types/signals/abstractsignal.jl b/src/types/signals/abstractsignal.jl
index a108e84..b9626cc 100644
--- a/src/types/signals/abstractsignal.jl
+++ b/src/types/signals/abstractsignal.jl
@@ -1,9 +1,9 @@
 @compat abstract type AbstractSignal{T<:Real,N} end
 
-discontinuities{T1<:Real,T2<:Real}(u, tspan::Tuple{T1,T2}) = promote_type(T1,T2)[]
+discontinuities(u, tspan::Tuple{T1,T2}) where {T1<:Real,T2<:Real} = promote_type(T1,T2)[]
 
-_eltype{T,N}(s::AbstractSignal{T,N}) = T
-_length{T,N}(s::AbstractSignal{T,N}) = N
+_eltype(s::AbstractSignal{T,N}) where {T,N} = T
+_length(s::AbstractSignal{T,N}) where {T,N} = N
 
-+{T1,T2,N}(s1::AbstractSignal{T1,N}, s2::AbstractSignal{T2,N}) = SumOfSignals(s1,  s2)
--{T1,T2,N}(s1::AbstractSignal{T1,N}, s2::AbstractSignal{T2,N}) = SumOfSignals(s1, -s2)
++(s1::AbstractSignal{T1,N}, s2::AbstractSignal{T2,N}) where {T1,T2,N} = SumOfSignals(s1,  s2)
+-(s1::AbstractSignal{T1,N}, s2::AbstractSignal{T2,N}) where {T1,T2,N} = SumOfSignals(s1, -s2)
diff --git a/src/types/signals/prbs.jl b/src/types/signals/prbs.jl
index 81dbac2..ba70347 100644
--- a/src/types/signals/prbs.jl
+++ b/src/types/signals/prbs.jl
@@ -1,11 +1,10 @@
-immutable PRBS{T,N,S1,S2} <: AbstractSignal{T,N}
+struct PRBS{T,N,S1,S2} <: AbstractSignal{T,N}
   mag::Vector{T}
   dc::Vector{T}
   internal::Vector{Tuple{S1,S1,S1,S1}}
 
-  function (::Type{PRBS}){T1<:Real,T2<:Real,T3<:Integer,N}(mag::Vector{T1},
-    dc::Vector{T2} = zeros(mag), seed::Vector{T3} = ones(Int, dc),
-    ::Type{Val{N}} = Val{7})
+  function (::Type{PRBS})(mag::Vector{T1}, dc::Vector{T2} = zeros(mag), seed::Vector{T3} = ones(Int, dc),
+    ::Type{Val{N}} = Val{7}) where {T1<:Real,T2<:Real,T3<:Integer,N}
     if !(length(mag) == length(dc) == length(seed))
       warn("PRBS(mag, dc, seed, Val{N}): `mag`, `dc` and `seed` must all have the same lengths")
       throw(DomainError())
@@ -26,8 +25,8 @@ immutable PRBS{T,N,S1,S2} <: AbstractSignal{T,N}
       internal)
   end
 
-  function (::Type{PRBS}){T1<:Real,T2<:Real,S1<:Integer,N}(mag::Vector{T1},
-    dc::Vector{T2}, internal::Vector{Tuple{S1,S1,S1,S1}}, ::Type{Val{N}})
+  function (::Type{PRBS})(mag::Vector{T1}, dc::Vector{T2},
+    internal::Vector{Tuple{S1,S1,S1,S1}}, ::Type{Val{N}}) where {T1<:Real,T2<:Real,S1<:Integer,N}
     if !(length(mag) == length(dc) == length(internal))
       warn("PRBS(mag, dc, internal, Val{N}): `mag`, `dc` and `internal` must all have the same lengths")
       throw(DomainError())
@@ -50,11 +49,11 @@ immutable PRBS{T,N,S1,S2} <: AbstractSignal{T,N}
     new{T,length(internal),S1,Val{N}}(convert(Vector{T}, mag), convert(Vector{T}, dc), internal)
   end
 
-  (::Type{PRBS}){N}(mag::Real, dc::Real = zero(mag), seed::Integer = 1,
-    t::Type{Val{N}} = Val{7}) = PRBS([mag], [dc], [seed], t)
+  (::Type{PRBS})(mag::Real, dc::Real = zero(mag), seed::Integer = 1,
+    t::Type{Val{N}} = Val{7}) where {N} = PRBS([mag], [dc], [seed], t)
   (::Type{PRBS})() = PRBS(1.)
 
-  function (p::PRBS{T,N,S1}){T,N,S1}()
+  function (p::PRBS{T,N,S1})() where {T,N,S1}
     temp = zeros(S1, N)
     for idx in 1:length(p.internal)
       a, b, val, _    = p.internal[idx]
@@ -66,7 +65,7 @@ immutable PRBS{T,N,S1,S2} <: AbstractSignal{T,N}
     temp
   end
 
-  function (p::PRBS{T,N,S1,Val{S2}}){T,N,S1,S2}(t::Real, x = nothing)
+  function (p::PRBS{T,N,S1,Val{S2}})(t::Real, x = nothing) where {T,N,S1,S2}
     temp  = zeros(Bool, N)
     for idx in 1:length(temp)
       a, b, val, k = p.internal[idx]
@@ -106,20 +105,20 @@ _coeffs(::Type{Val{31}})  = (30 , 27)
 # discontinuities{T1<:Real,T2<:Real}(p::PRBS, tspan::Tuple{T1,T2}) =
 #   promote_type(T1, T2)[]
 
-convert{T,N,S1,S2}(::Type{AbstractSignal{T,N}}, p::PRBS{T,N,S1,Val{S2}})        = p
-convert{T1,T2,N,S1,S2}(::Type{AbstractSignal{T1,N}}, p::PRBS{T2,N,S1,Val{S2}})  =
+convert(::Type{AbstractSignal{T,N}}, p::PRBS{T,N,S1,Val{S2}}) where {T,N,S1,S2} = p
+convert(::Type{AbstractSignal{T1,N}}, p::PRBS{T2,N,S1,Val{S2}}) where {T1,T2,N,S1,S2} =
   PRBS(convert(Vector{T1}, p.mag), convert(Vector{T1}, p.dc), internal, Val{S2})
 
--{T,N,S1,S2}(p::PRBS{T,N,S1,Val{S2}}) = PRBS(-p.mag, -p.dc, p.internal, Val{S2})
+-(p::PRBS{T,N,S1,Val{S2}}) where {T,N,S1,S2} = PRBS(-p.mag, -p.dc, p.internal, Val{S2})
 
 # Relation between `Real`s
 
-+{T1<:Real,T2,N,S1,S2}(x::Union{T1,AbstractVector{T1}}, p::PRBS{T2,N,S1,Val{S2}}) =
++(x::Union{T1,AbstractVector{T1}}, p::PRBS{T2,N,S1,Val{S2}}) where {T1<:Real,T2,N,S1,S2} =
   PRBS(p.mag, p.dc + x, p.internal, Val{S2})
-+{T<:Real}(p::PRBS, x::Union{T,AbstractVector{T}})  = +(x,  p)
--{T<:Real}(x::Union{T,AbstractVector{T}}, p::PRBS)  = +(x, -p)
--{T<:Real}(p::PRBS, x::Union{T,AbstractVector{T}})  = +(p, -x)
++(p::PRBS, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  p)
+-(x::Union{T,AbstractVector{T}}, p::PRBS) where {T<:Real} = +(x, -p)
+-(p::PRBS, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(p, -x)
 
-*{T,N,S1,S2}(x::Real, p::PRBS{T,N,S1,Val{S2}}) = PRBS(x*p.mag, x*p.dc, p.internal, Val{S2})
+*(x::Real, p::PRBS{T,N,S1,Val{S2}}) where {T,N,S1,S2} = PRBS(x*p.mag, x*p.dc, p.internal, Val{S2})
 *(p::PRBS, x::Real) = *(x,   p)
 /(p::PRBS, x::Real) = *(p, 1/x)
diff --git a/src/types/signals/ramp.jl b/src/types/signals/ramp.jl
index fe907a1..20bcac2 100644
--- a/src/types/signals/ramp.jl
+++ b/src/types/signals/ramp.jl
@@ -1,10 +1,10 @@
-immutable Ramp{T,N} <: AbstractSignal{T,N}
+struct Ramp{T,N} <: AbstractSignal{T,N}
   time::Vector{T}
   rate::Vector{T}
   dc::Vector{T}
 
-  function (::Type{Ramp}){T1<:Real,T2<:Real,T3<:Real}(rt::Vector{T1}, rr::Vector{T2},
-    dc::Vector{T3} = zeros(rr))
+  function (::Type{Ramp})(rt::Vector{T1}, rr::Vector{T2},
+    dc::Vector{T3} = zeros(rr)) where {T1<:Real,T2<:Real,T3<:Real}
     if length(rt) ≠ length(rr) || length(rt) ≠ length(dc)
       warn("Ramp(rt, rr, dc): `rt`, `rr` and `dc` must all have the same lengths")
       throw(DomainError())
@@ -16,7 +16,7 @@ immutable Ramp{T,N} <: AbstractSignal{T,N}
   (::Type{Ramp})(rt::Real, rr::Real, dc::Real = zero(rr)) = Ramp([rt], [rr], [dc])
   (::Type{Ramp})()                                        = Ramp([0.], [1.], [0.])
 
-  function (r::Ramp{T,N}){T,N}(t::Real, x = nothing)
+  function (r::Ramp{T,N})(t::Real, x = nothing) where {T,N}
     res = one(T) + one(t)*one(T)
     R   = typeof(res)
 
@@ -24,11 +24,11 @@ immutable Ramp{T,N} <: AbstractSignal{T,N}
   end
 end
 
-discontinuities{T1<:Real,T2<:Real}(r::Ramp, tspan::Tuple{T1,T2}) =
+discontinuities(r::Ramp, tspan::Tuple{T1,T2}) where {T1<:Real,T2<:Real} =
   [rt for rt in r.time if tspan[1] ≤ rt ≤ tspan[2]]
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, r::Ramp{T,N})       = r
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, r::Ramp{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, r::Ramp{T,N}) where {T,N} = r
+convert(::Type{AbstractSignal{T1,N}}, r::Ramp{T2,N}) where {T1,T2,N} =
   Ramp(convert(Vector{T1}, r.time), convert(Vector{T1}, r.rate),
     convert(Vector{T1}, r.dc))
 
@@ -36,10 +36,10 @@ convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, r::Ramp{T2,N}) =
 
 # Relation between `Real`s
 
-+{T<:Real}(x::Union{T,AbstractVector{T}}, r::Ramp) = Ramp(r.time, r.rate, r.dc + x)
-+{T<:Real}(r::Ramp, x::Union{T,AbstractVector{T}}) = +(x,  r)
--{T<:Real}(x::Union{T,AbstractVector{T}}, r::Ramp) = +(x, -r)
--{T<:Real}(r::Ramp, x::Union{T,AbstractVector{T}}) = +(r, -x)
++(x::Union{T,AbstractVector{T}}, r::Ramp) where {T<:Real} = Ramp(r.time, r.rate, r.dc + x)
++(r::Ramp, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  r)
+-(x::Union{T,AbstractVector{T}}, r::Ramp) where {T<:Real} = +(x, -r)
+-(r::Ramp, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(r, -x)
 
 *(x::Real, r::Ramp) = Ramp(r.time, x*r.rate, x*r.dc)
 *(r::Ramp, x::Real) = *(x,   r)
diff --git a/src/types/signals/sinusoid.jl b/src/types/signals/sinusoid.jl
index e4657c7..1b5c34d 100644
--- a/src/types/signals/sinusoid.jl
+++ b/src/types/signals/sinusoid.jl
@@ -1,13 +1,13 @@
-immutable Sinusoid{T,N} <: AbstractSignal{T,N}
+struct Sinusoid{T,N} <: AbstractSignal{T,N}
   magnitude::Vector{Vector{T}}
   dc::Vector{Vector{T}}
   shift::Vector{Vector{T}}
   frequency::Vector{Vector{T}}
 
-  function (::Type{Sinusoid}){T1<:Real,T2<:Real,T3<:Real,T4<:Real}(mag::Vector{Vector{T1}},
+  function (::Type{Sinusoid})(mag::Vector{Vector{T1}},
     dc::Vector{Vector{T2}}    = [zeros(elem) for elem in mag],
     shift::Vector{Vector{T3}} = [zeros(elem) for elem in mag],
-    freq::Vector{Vector{T4}}  = [ones(elem) for elem in mag])
+    freq::Vector{Vector{T4}}  = [ones(elem) for elem in mag]) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}
     if !(length(mag) == length(dc) == length(shift) == length(freq))
       warn("Sinusoid(mag, dc, shift, freq): `mag`, `dc`, `shift` and `freq` must all have the same lengths")
       throw(DomainError())
@@ -26,10 +26,9 @@ immutable Sinusoid{T,N} <: AbstractSignal{T,N}
           convert(Vector{Vector{T}}, shift), convert(Vector{Vector{T}}, freq))
   end
 
-  (::Type{Sinusoid}){T1<:Real,T2<:Real,T3<:Real,T4<:Real}(mag::Vector{T1},
-    dc::Vector{T2} = zeros(mag), shift::Vector{T3} = zeros(mag),
-    freq::Vector{T4} = ones(mag)) = Sinusoid([[elem] for elem in mag],
-    [[elem] for elem in dc], [[elem] for elem in shift], [[elem] for elem in freq])
+  (::Type{Sinusoid})(mag::Vector{T1}, dc::Vector{T2} = zeros(mag),
+    shift::Vector{T3} = zeros(mag), freq::Vector{T4} = ones(mag)) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real} =
+    Sinusoid([[elem] for elem in mag], [[elem] for elem in dc], [[elem] for elem in shift], [[elem] for elem in freq])
 
   (::Type{Sinusoid})(mag::Real, dc::Real = zero(mag), shift::Real = zero(mag),
     freq::Real = one(mag)) = Sinusoid([[mag]], [[dc]], [[shift]], [[freq]])
@@ -40,13 +39,13 @@ immutable Sinusoid{T,N} <: AbstractSignal{T,N}
   end
 end
 
-function discontinuities{T1<:Real,T2<:Real}(s::Sinusoid, tspan::Tuple{T1,T2})
+function discontinuities(s::Sinusoid, tspan::Tuple{T1,T2}) where {T1<:Real,T2<:Real}
   resp = s(tspan[1])
   return (resp ≈ zeros(resp)) ? T1[] : [tspan[1]]
 end
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, s::Sinusoid{T,N})       = s
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Sinusoid{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, s::Sinusoid{T,N}) where {T,N}       = s
+convert(::Type{AbstractSignal{T1,N}}, s::Sinusoid{T2,N}) where {T1,T2,N} =
   Sinusoid(convert(Vector{Vector{T1}}, s.magnitude),
            convert(Vector{Vector{T1}}, s.dc),
            convert(Vector{Vector{T1}}, s.shift),
@@ -56,11 +55,11 @@ convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Sinusoid{T2,N}) =
 
 # Relation between `Real`s
 
-+{T<:Real}(x::Union{T,AbstractVector{T}}, s::Sinusoid) =
++(x::Union{T,AbstractVector{T}}, s::Sinusoid) where {T<:Real} =
   Sinusoid(s.magnitude, s.dc + x, s.shift, s.frequency)
-+{T<:Real}(s::Sinusoid, x::Union{T,AbstractVector{T}}) = +(x,  s)
--{T<:Real}(x::Union{T,AbstractVector{T}}, s::Sinusoid) = +(x, -s)
--{T<:Real}(s::Sinusoid, x::Union{T,AbstractVector{T}}) = +(s, -x)
++(s::Sinusoid, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  s)
+-(x::Union{T,AbstractVector{T}}, s::Sinusoid) where {T<:Real} = +(x, -s)
+-(s::Sinusoid, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(s, -x)
 
 *(x::Real, s::Sinusoid) = Sinusoid(x*s.magnitude, x*s.dc, s.shift, s.frequency)
 *(s::Sinusoid, x::Real) = *(x,   s)
diff --git a/src/types/signals/square.jl b/src/types/signals/square.jl
index 7a4b1e5..edfa9a3 100644
--- a/src/types/signals/square.jl
+++ b/src/types/signals/square.jl
@@ -1,13 +1,13 @@
-immutable Square{T,N} <: AbstractSignal{T,N}
+struct Square{T,N} <: AbstractSignal{T,N}
   magnitude::Vector{T}
   width::Vector{T}
   dc::Vector{T}
   shift::Vector{T}
   period::Vector{T}
 
-  function (::Type{Square}){T1,T2,T3,T4,T5}(mag::Vector{T1},
+  function (::Type{Square})(mag::Vector{T1},
     width::Vector{T2} = ones(mag), dc::Vector{T3} = zeros(mag),
-    shift::Vector{T4} = zeros(mag), per::Vector{T5} = width)
+    shift::Vector{T4} = zeros(mag), per::Vector{T5} = width) where {T1,T2,T3,T4,T5}
     if !(length(mag) == length(width) == length(dc) == length(shift) == length(per))
       warn("Square(mag, width, dc, shift, per): all inputs must have the same lengths")
       throw(DomainError())
@@ -30,7 +30,7 @@ immutable Square{T,N} <: AbstractSignal{T,N}
     [shift], [per])
   (::Type{Square})() = Square(1.)
 
-  function (s::Square{T}){T}(t::Real, x = nothing)
+  function (s::Square{T})(t::Real, x = nothing) where {T}
     res   = divrem.(t - s.shift, s.period)
     tnorm = map(y->y[2], res)
     [(t < zero(t)) ? dc : (t < 0.5w) ? dc + m : (t < w) ? dc - m : dc for
@@ -38,7 +38,7 @@ immutable Square{T,N} <: AbstractSignal{T,N}
   end
 end
 
-function discontinuities{T1<:Real,T2<:Real,T3<:Real}(s::Square{T1}, tspan::Tuple{T2,T3})
+function discontinuities(s::Square{T1}, tspan::Tuple{T2,T3}) where {T1<:Real,T2<:Real,T3<:Real}
   tstops = Set{T1}()
 
   for i1 in 1:length(s.period)
@@ -62,8 +62,8 @@ function discontinuities{T1<:Real,T2<:Real,T3<:Real}(s::Square{T1}, tspan::Tuple
   [tstop for tstop in tstops if tspan[1] ≤ tstop ≤ tspan[2]]
 end
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, s::Square{T,N})       = s
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Square{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, s::Square{T,N}) where {T,N} = s
+convert(::Type{AbstractSignal{T1,N}}, s::Square{T2,N}) where {T1,T2,N} =
   Square(convert(Vector{T1}, s.magnitude), convert(Vector{T1}, s.width),
     convert(Vector{T1}, s.dc), convert(Vector{T1}, s.shift),
     convert(Vector{T1}, s.period))
@@ -72,11 +72,11 @@ convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Square{T2,N}) =
 
 # Relation between `Real`s
 
-+{T<:Real}(x::Union{T,AbstractVector{T}}, s::Square)  =
++(x::Union{T,AbstractVector{T}}, s::Square) where {T<:Real} =
   Square(s.magnitude, s.width, s.dc + x, s.shift, s.period)
-+{T<:Real}(s::Square, x::Union{T,AbstractVector{T}})  = +(x,  r)
--{T<:Real}(x::Union{T,AbstractVector{T}}, s::Square)  = +(x, -r)
--{T<:Real}(s::Square, x::Union{T,AbstractVector{T}})  = +(r, -x)
++(s::Square, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  r)
+-(x::Union{T,AbstractVector{T}}, s::Square) where {T<:Real} = +(x, -r)
+-(s::Square, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(r, -x)
 
 *(x::Real, s::Square) = Square(x*s.magnitude, s.width, x*s.dc, s.shift, s.period)
 *(s::Square, x::Real) = *(x,   r)
diff --git a/src/types/signals/step.jl b/src/types/signals/step.jl
index 0b190dd..94e0cd9 100644
--- a/src/types/signals/step.jl
+++ b/src/types/signals/step.jl
@@ -1,10 +1,10 @@
-immutable Step{T,N} <: AbstractSignal{T,N}
+struct Step{T,N} <: AbstractSignal{T,N}
   time::Vector{T}
   step::Vector{T}
   dc::Vector{T}
 
-  function (::Type{Step}){T1<:Real,T2<:Real,T3<:Real}(st::Vector{T1}, sv::Vector{T2},
-    dc::Vector{T3} = zeros(sv))
+  function (::Type{Step})(st::Vector{T1}, sv::Vector{T2},
+    dc::Vector{T3} = zeros(sv)) where {T1<:Real,T2<:Real,T3<:Real}
     if length(st) ≠ length(sv) || length(st) ≠ length(dc)
       warn("Step(st, sv, dc): `st`, `sv` and `dc` must all have the same lengths")
       throw(DomainError())
@@ -21,11 +21,11 @@ immutable Step{T,N} <: AbstractSignal{T,N}
   end
 end
 
-discontinuities{T1<:Real,T2<:Real}(s::Step, tspan::Tuple{T1,T2}) =
+discontinuities(s::Step, tspan::Tuple{T1,T2}) where {T1<:Real,T2<:Real} =
   [st for st in s.time if tspan[1] ≤ st ≤ tspan[2]]
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, s::Step{T,N})       = s
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Step{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, s::Step{T,N}) where {T,N}       = s
+convert(::Type{AbstractSignal{T1,N}}, s::Step{T2,N}) where {T1,T2,N} =
   Step(convert(Vector{T1}, s.time), convert(Vector{T1}, s.step),
     convert(Vector{T1}, s.dc))
 
@@ -33,10 +33,10 @@ convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, s::Step{T2,N}) =
 
 # Relation between `Real`s
 
-+{T<:Real}(x::Union{T,AbstractVector{T}}, s::Step) = Step(s.time, s.step, s.dc + x)
-+{T<:Real}(s::Step, x::Union{T,AbstractVector{T}}) = +(x,  s)
--{T<:Real}(x::Union{T,AbstractVector{T}}, s::Step) = +(x, -s)
--{T<:Real}(s::Step, x::Union{T,AbstractVector{T}}) = +(s, -x)
++(x::Union{T,AbstractVector{T}}, s::Step) where {T<:Real} = Step(s.time, s.step, s.dc + x)
++(s::Step, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  s)
+-(x::Union{T,AbstractVector{T}}, s::Step) where {T<:Real} = +(x, -s)
+-(s::Step, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(s, -x)
 
 *(x::Real, s::Step) = Step(s.time, x*s.step, x*s.dc)
 *(s::Step, x::Real) = *(x,   s)
diff --git a/src/types/signals/sumofsignals.jl b/src/types/signals/sumofsignals.jl
index b122e12..4fc5550 100644
--- a/src/types/signals/sumofsignals.jl
+++ b/src/types/signals/sumofsignals.jl
@@ -1,7 +1,7 @@
-immutable SumOfSignals{T,N} <: AbstractSignal{T,N}
+struct SumOfSignals{T,N} <: AbstractSignal{T,N}
   signals::Vector{AbstractSignal{T,N}}
 
-  function (::Type{SumOfSignals}){T1<:AbstractSignal}(signals::Vector{T1})
+  function (::Type{SumOfSignals})(signals::Vector{T1}) where {T1<:AbstractSignal}
     if isempty(signals)
       warn("SumOfSignals(signals): `signals` is empty")
       throw(DomainError())
@@ -23,7 +23,7 @@ immutable SumOfSignals{T,N} <: AbstractSignal{T,N}
     SumOfSignals(signals)
   end
 
-  function (sos::SumOfSignals{T,N}){T,N}(t::Real, x = nothing)
+  function (sos::SumOfSignals{T,N})(t::Real, x = nothing) where {T,N}
     temp = zeros(T,N)
     for signal in sos.signals
       temp += signal(t, x)
@@ -32,8 +32,8 @@ immutable SumOfSignals{T,N} <: AbstractSignal{T,N}
   end
 end
 
-function discontinuities{T1,T2<:Real,T3<:Real}(sos::SumOfSignals{T1},
-  tspan::Tuple{T2,T3})
+function discontinuities(sos::SumOfSignals{T1},
+  tspan::Tuple{T2,T3}) where {T1,T2<:Real,T3<:Real}
   tstops = Set{T1}()
   for signal in sos.signals
     push!(tstops, discontinuities(signal, tspan)...)
@@ -41,21 +41,21 @@ function discontinuities{T1,T2<:Real,T3<:Real}(sos::SumOfSignals{T1},
   [tstop for tstop in tstops]
 end
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, sos::SumOfSignals{T,N})       = sos
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, sos::SumOfSignals{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, sos::SumOfSignals{T,N}) where {T,N}       = sos
+convert(::Type{AbstractSignal{T1,N}}, sos::SumOfSignals{T2,N}) where {T1,T2,N} =
   SumOfSignals(convert(Vector{AbstractSignal{T1,N}}, sos.signals))
 
 -(sos::SumOfSignals) = SumOfSignals(-signals)
 
 # Relation between `Real`s
 
-function +{T<:Real}(x::Union{T,AbstractVector{T}}, sos::SumOfSignals)
+function +(x::Union{T,AbstractVector{T}}, sos::SumOfSignals) where {T<:Real}
   temp = [signal + x for signal in sos.signals]
   SumOfSignals(temp)
 end
-+{T<:Real}(sos::SumOfSignals, x::Union{T,AbstractVector{T}}) = +(x,  sos)
--{T<:Real}(x::Union{T,AbstractVector{T}}, sos::SumOfSignals) = +(x, -sos)
--{T<:Real}(sos::SumOfSignals, x::Union{T,AbstractVector{T}}) = +(sos, -x)
++(sos::SumOfSignals, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(x,  sos)
+-(x::Union{T,AbstractVector{T}}, sos::SumOfSignals) where {T<:Real} = +(x, -sos)
+-(sos::SumOfSignals, x::Union{T,AbstractVector{T}}) where {T<:Real} = +(sos, -x)
 
 function *(x::Real, sos::SumOfSignals)
   temp = [x*signal for signal in sos.signals]
@@ -66,15 +66,15 @@ end
 
 # Relationship between `AbstractSignal`s
 
-+{T1,T2,N}(sos::SumOfSignals{T1,N}, s::AbstractSignal{T2,N}) =
++(sos::SumOfSignals{T1,N}, s::AbstractSignal{T2,N}) where {T1,T2,N} =
   SumOfSignals(push!(copy(sos.signals), s))
-+{T1,T2,N}(s::AbstractSignal{T2,N}, sos::SumOfSignals{T1,N}) = +(sos,  s)
--{T1,T2,N}(s::AbstractSignal{T2,N}, sos::SumOfSignals{T1,N}) = +(s, -sos)
--{T1,T2,N}(sos::SumOfSignals{T1,N}, s::AbstractSignal{T2,N}) = +(sos, -s)
++(s::AbstractSignal{T2,N}, sos::SumOfSignals{T1,N}) where {T1,T2,N} = +(sos,  s)
+-(s::AbstractSignal{T2,N}, sos::SumOfSignals{T1,N}) where {T1,T2,N} = +(s, -sos)
+-(sos::SumOfSignals{T1,N}, s::AbstractSignal{T2,N}) where {T1,T2,N} = +(sos, -s)
 
 # Relationship among `SumOfSignals`
 
-+{T1,T2,N}(sos1::SumOfSignals{T1,N}, sos2::SumOfSignals{T2,N}) =
++(sos1::SumOfSignals{T1,N}, sos2::SumOfSignals{T2,N}) where {T1,T2,N} =
   SumOfSignals(AbstractSignal{promote_type(T1,T2),N}[sos1.signals..., sos2.signals...])
--{T1,T2,N}(sos1::SumOfSignals{T1,N}, sos2::SumOfSignals{T2,N}) =
+-(sos1::SumOfSignals{T1,N}, sos2::SumOfSignals{T2,N}) where {T1,T2,N} =
   SumOfSignals(AbstractSignal{promote_type(T1,T2),N}[sos1.signals..., -sos2.signals...])
diff --git a/src/types/signals/triangle.jl b/src/types/signals/triangle.jl
index 09e6c6a..4b0ec4f 100644
--- a/src/types/signals/triangle.jl
+++ b/src/types/signals/triangle.jl
@@ -1,13 +1,13 @@
-immutable Triangle{T,N} <: AbstractSignal{T,N}
+struct Triangle{T,N} <: AbstractSignal{T,N}
   magnitude::Vector{T}
   width::Vector{T}
   dc::Vector{T}
   shift::Vector{T}
   period::Vector{T}
 
-  function (::Type{Triangle}){T1<:Real,T2<:Real,T3<:Real,T4<:Real,T5<:Real}(mag::Vector{T1},
+  function (::Type{Triangle})(mag::Vector{T1},
     width::Vector{T2} = ones(mag), dc::Vector{T3} = zeros(mag),
-    shift::Vector{T4} = zeros(mag), per::Vector{T5} = width)
+    shift::Vector{T4} = zeros(mag), per::Vector{T5} = width) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real,T5<:Real}
     if !(length(mag) == length(width) == length(dc) == length(shift) == length(per))
       warn("Triangle(mag, width, dc, shift, per): all inputs must have the same lengths")
       throw(DomainError())
@@ -30,7 +30,7 @@ immutable Triangle{T,N} <: AbstractSignal{T,N}
     [shift], [per])
   (::Type{Triangle})() = Triangle(1.)
 
-  function (tr::Triangle{T}){T}(t::Real, x = nothing)
+  function (tr::Triangle{T})(t::Real, x = nothing) where {T}
     res   = divrem.(t - tr.shift, tr.period)
     tnorm = map(y->y[2], res)
     [(t < zero(t)) ? dc : (t < 0.25w) ? dc + 4m/w*t : (t < 0.75w) ? dc + m - 4m/w*(t-0.25w) :
@@ -39,7 +39,7 @@ immutable Triangle{T,N} <: AbstractSignal{T,N}
   end
 end
 
-function discontinuities{T1<:Real,T2<:Real,T3<:Real}(t::Triangle{T1}, tspan::Tuple{T2,T3})
+function discontinuities(t::Triangle{T1}, tspan::Tuple{T2,T3}) where {T1<:Real,T2<:Real,T3<:Real}
   tstops = Set{T1}()
 
   for i1 in 1:length(t.period)
@@ -64,8 +64,8 @@ function discontinuities{T1<:Real,T2<:Real,T3<:Real}(t::Triangle{T1}, tspan::Tup
   [tstop for tstop in tstops if tspan[1] ≤ tstop ≤ tspan[2]]
 end
 
-convert{T,N}(::Type{AbstractSignal{T,N}}, t::Triangle{T,N})       = t
-convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, t::Triangle{T2,N}) =
+convert(::Type{AbstractSignal{T,N}}, t::Triangle{T,N}) where {T,N}       = t
+convert(::Type{AbstractSignal{T1,N}}, t::Triangle{T2,N}) where {T1,T2,N} =
   Triangle(convert(Vector{T1}, t.magnitude), convert(Vector{T1}, t.width),
     convert(Vector{T1}, t.dc), convert(Vector{T1}, t.shift),
     convert(Vector{T1}, t.period))
@@ -74,11 +74,11 @@ convert{T1,T2,N}(::Type{AbstractSignal{T1,N}}, t::Triangle{T2,N}) =
 
 # Relation between `Real`s
 
-+{T<:Real}(x::Union{T,AbstractVector{T}}, r::Triangle)  =
++(x::Union{T,AbstractVector{T}}, r::Triangle) where {T<:Real}  =
   Triangle(t.magnitude, t.width, t.dc + x, t.shift, t.period)
-+{T<:Real}(t::Triangle, x::Union{T,AbstractVector{T}})  = +(x,  t)
--{T<:Real}(x::Union{T,AbstractVector{T}}, t::Triangle)  = +(x, -t)
--{T<:Real}(t::Triangle, x::Union{T,AbstractVector{T}})  = +(t, -x)
++(t::Triangle, x::Union{T,AbstractVector{T}}) where {T<:Real}  = +(x,  t)
+-(x::Union{T,AbstractVector{T}}, t::Triangle) where {T<:Real}  = +(x, -t)
+-(t::Triangle, x::Union{T,AbstractVector{T}}) where {T<:Real}  = +(t, -x)
 
 *(x::Real, t::Triangle) = Triangle(x*t.magnitude, t.width, x*t.dc, t.shift, t.period)
 *(t::Triangle, x::Real) = *(x,   t)
diff --git a/src/types/system/ltisystem.jl b/src/types/system/ltisystem.jl
index cba1086..a75681d 100644
--- a/src/types/system/ltisystem.jl
+++ b/src/types/system/ltisystem.jl
@@ -16,7 +16,7 @@ function mimo(s::LtiSystem{Val{:siso}})
 end
 
 # Sampling-time
-iscontinuous{T}(::LtiSystem{Val{T},Val{:cont}}) = true
-iscontinuous{T}(::LtiSystem{Val{T},Val{:disc}}) = false
-isdiscrete(s::LtiSystem)                        = !iscontinuous(s)
-samplingtime{T}(::LtiSystem{Val{T},Val{:disc}}) = zero(Float64)
+iscontinuous(::LtiSystem{Val{T},Val{:cont}}) where {T} = true
+iscontinuous(::LtiSystem{Val{T},Val{:disc}}) where {T} = false
+isdiscrete(s::LtiSystem)                               = !iscontinuous(s)
+samplingtime(::LtiSystem{Val{T},Val{:disc}}) where {T} = zero(Float64)
diff --git a/src/types/system/mfd.jl b/src/types/system/mfd.jl
index cea1126..ced59a1 100644
--- a/src/types/system/mfd.jl
+++ b/src/types/system/mfd.jl
@@ -12,7 +12,7 @@
 #
 # NOTE: Should SISO MatrixFractionDescriptions be based on 1x1 polynomial matrices (similarly to `TransferFunction`)?
 # NOTE: Should the constructors verify that D is nonsingular and that the MatrixFractionDescription is proper?
-immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
+struct MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   N::M1
   D::M2
   nu::Int
@@ -20,7 +20,7 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   Ts::Float64
 
   # Continuous-time, single-input-single-output MatrixFractionDescription model
-  function (::Type{MatrixFractionDescription}){L,M1<:Polynomials.Poly,M2<:Polynomials.Poly}(N::M1, D::M2, ::Type{Val{L}})
+  function (::Type{MatrixFractionDescription})(N::M1, D::M2, ::Type{Val{L}}) where {L,M1<:Polynomials.Poly,M2<:Polynomials.Poly}
     mfdcheck(N,D)
     nN,nD = length(N),length(D)
     _N = PolyMatrix(reshape(coeffs(N),1,1,nN), (1,1,nN), Val{:s})
@@ -29,7 +29,7 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   end
 
   # Discrete-time, single-input-single-output MatrixFractionDescription model
-  function (::Type{MatrixFractionDescription}){L,M1<:Polynomials.Poly,M2<:Polynomials.Poly}(N::M1, D::M2, Ts::Real, ::Type{Val{L}})
+  function (::Type{MatrixFractionDescription})(N::M1, D::M2, Ts::Real, ::Type{Val{L}}) where {L,M1<:Polynomials.Poly,M2<:Polynomials.Poly}
     mfdcheck(N,D,Ts)
     nN,nD = length(N),length(D)
     _N = PolyMatrix(reshape(coeffs(N),1,1,nN), (1,1,nN), Val{:s})
@@ -38,8 +38,8 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   end
 
   # Continuous-time, multi-input-multi-output MatrixFractionDescription model
-  function (::Type{MatrixFractionDescription}){L,M1<:PolynomialMatrices.PolyMatrix,
-    M2<:PolynomialMatrices.PolyMatrix}(N::M1, D::M2, ::Type{Val{L}})
+  function (::Type{MatrixFractionDescription})(N::M1, D::M2, ::Type{Val{L}}) where {L,M1<:PolynomialMatrices.PolyMatrix,
+    M2<:PolynomialMatrices.PolyMatrix}
     ny, nu = mfdcheck(N, D, Val{L})
     _N = PolyMatrix(coeffs(N), size(N), Val{:s})
     _D = PolyMatrix(coeffs(D), size(D), Val{:s})
@@ -49,8 +49,8 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   end
 
   # Discrete-time, multi-input-multi-output MatrixFractionDescription model
-  function (::Type{MatrixFractionDescription}){L,M1<:PolynomialMatrices.PolyMatrix,
-    M2<:PolynomialMatrices.PolyMatrix}(N::M1, D::M2, Ts::Real, ::Type{Val{L}})
+  function (::Type{MatrixFractionDescription})(N::M1, D::M2, Ts::Real, ::Type{Val{L}}) where {L,M1<:PolynomialMatrices.PolyMatrix,
+    M2<:PolynomialMatrices.PolyMatrix}
     ny, nu = mfdcheck(N, D, Val{L}, Ts)
     _N = PolyMatrix(coeffs(N), size(N), Val{:z})
     _D = PolyMatrix(coeffs(D), size(D), Val{:z})
@@ -62,12 +62,11 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   # Function calls
   # TODO needs to be properly set up according to e.g. Kailath Linear systems 6.4
   # to properly handle all column/row reduced cases
-  function (sys::MatrixFractionDescription{T,S,Val{:lfd}}){T,S}(t::Real, x::DiffEqBase.DEDataArray, dx::AbstractVector, u)
+  function (sys::MatrixFractionDescription{T,S,Val{:lfd}})(t::Real, x::DiffEqBase.DEDataArray, dx::AbstractVector, u) where {T,S}
     ucalc = u(t, x.y)
     ucalc = isa(ucalc, Real) ? [ucalc] : ucalc
     if !isa(ucalc, AbstractVector) || length(ucalc) ≠ sys.nu || !(eltype(ucalc) <: Real)
-      warn("sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values")
-      throw(DomainError())
+      throw(DomainError((u, t), "sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values"))
     end
 
     ny  = numoutputs(sys)
@@ -78,25 +77,26 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
     Dhci = inv(Dhc)
     D = Dhci*D
 
-    xidx = (cumsum(vcat(0,degs))+1)[1:end-1]
-    dxidx = cumsum(degs,1)[:]
+    xidx = broadcast(+, 1, cumsum(vcat(0,degs), dims=1))[1:end-1]
+    dxidx = cumsum(degs, dims=1)[:]
 
-    dx[:] = zeros(dx)
+    dx[:] = zero(dx)
     # "A matrix"
     for (k,v) in coeffs(D)
-      for row in find(degs .> k)
+      for row in findall(degs[:] .> k)
         dx[dxidx[row]-k] += -dot(v[row,:], x.x[xidx])
       end
     end
 
     for row in 1:ny
-      nrowx = degs[row]-1
-      dx[xidx[row]-1+(1:nrowx)] += x.x[xidx[row]+(1:nrowx)]
+      for krowx in 1:(degs[row]-1)
+        dx[xidx[row]-1+krowx] += x.x[xidx[row]+krowx]
+      end
     end
 
     # "B matrix"
     for (k,v) in coeffs(N)
-      for row in find(degs .> k)
+      for row in findall(degs[:] .> k)
         dx[dxidx[row]-k] += dot(v[row,:], ucalc)
       end
     end
@@ -105,7 +105,7 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
     if degree(N) == degree(D)
       tmp = coeffs(N)[degree(N)]*ucalc
       for (k,v) in coeffs(D)
-        for row in find(degs .> k)
+        for row in findall(degs[:] .> k)
           dx[dxidx[row]-k] += -dot(v[row,:], tmp[xidx])
         end
       end
@@ -118,12 +118,11 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
     x.u[:]  = ucalc
   end
 
-  function (sys::MatrixFractionDescription{T,S,Val{:rfd}}){T,S}(t::Real, x::DiffEqBase.DEDataArray, dx::AbstractVector, u)
+  function (sys::MatrixFractionDescription{T,S,Val{:rfd}})(t::Real, x::DiffEqBase.DEDataArray, dx::AbstractVector, u) where {T,S}
     ucalc = u(t, x.y)
     ucalc = isa(ucalc, Real) ? [ucalc] : ucalc
     if !isa(ucalc, AbstractVector) || length(ucalc) ≠ sys.nu || !(eltype(ucalc) <: Real)
-      warn("sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values")
-      throw(DomainError())
+      throw(DomainError((u, t), "sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values"))
     end
 
     ny  = numoutputs(sys)
@@ -135,21 +134,22 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
     Dhci = inv(Dhc)
     D = Dhci*D
 
-    xidx = cumsum(degs,2)[:]
-    dxidx = (cumsum(hcat(0,degs),2)+1)[1:end-1]
+    xidx = cumsum(degs, dims=2)[:]
+    dxidx = broadcast(+, 1, cumsum(hcat(0,degs), dims=2))[1:end-1]
 
-    x.y = zeros(x.y)
-    dx[:] = zeros(dx)
+    x.y = zero(x.y)
+    dx[:] = zero(dx)
 
     # "A matrix"
     for (k,v) in coeffs(D)
-      for col in find(degs .> k)
+      for col in findall(degs[:] .> k)
         dx[dxidx] += -v[:,col]*x.x[xidx[col]-k]
       end
     end
     for col in 1:nu
-      ncolx = degs[col]-1
-      dx[dxidx[col]+(1:ncolx)] += x.x[dxidx[col]-1+(1:ncolx)]
+      for kcolx in 1:(degs[col]-1)
+        dx[dxidx[col]+kcolx] += x.x[dxidx[col]-1+kcolx]
+      end
     end
 
     # "B matrix"
@@ -157,7 +157,7 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
 
     # "C matrix"
     for (k,v) in coeffs(N)
-      for col in find(degs .> k)
+      for col in findall(degs[:] .> k)
         x.y[:] += v[:,col]*x.x[xidx[col]-k]
       end
     end
@@ -168,7 +168,7 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
       b0 = coeffs(N)[degree(N)]
       tmp = ucalc
       for (k,v) in coeffs(D)
-        for col in find(degs .> k)
+        for col in findall(degs[:] .> k)
           tmp += -v[:,col]*x.x[xidx[col]-k]
         end
       end
@@ -177,13 +177,13 @@ immutable MatrixFractionDescription{T,S,L,M1,M2}  <: LtiSystem{T,S}
   end
 end
 
-function mfdcheck{T<:Real,S<:Real}(N::Poly{T}, D::Poly{S}, Ts::Real = zero(Float64))
+function mfdcheck(N::Poly{T}, D::Poly{S}, Ts::Real = zero(Float64)) where {T<:Real,S<:Real}
   @assert Ts ≥ zero(Ts) && !isinf(Ts) "MatrixFractionDescription: Ts must be non-negative real number"
 end
 
 isproper(s::MatrixFractionDescription{Val{:siso}})  = degree(s.D) ≥ degree(s.N)
 
-function isproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
+function isproper(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}}) where {S}
   if is_col_proper(s.D)
     return all(col_degree(s.D) .>= col_degree(s.N))
   else
@@ -191,7 +191,7 @@ function isproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
   end
 end
 
-function isproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}})
+function isproper(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}}) where {S}
   if is_row_proper(s.D)
     return all(row_degree(s.D) .>= row_degree(s.N))
   else
@@ -201,7 +201,7 @@ end
 
 isstrictlyproper(s::MatrixFractionDescription{Val{:siso}}) = degree(s.D) > degree(s.N)
 
-function isstrictlyproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
+function isstrictlyproper(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}}) where {S}
   if is_col_proper(s.D)
     return all(col_degree(s.D) .>= col_degree(s.N))
   else
@@ -209,7 +209,7 @@ function isstrictlyproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{
   end
 end
 
-function isstrictlyproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}})
+function isstrictlyproper(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}}) where {S}
   if is_row_proper(s.D)
     return all(row_degree(s.D) .>= row_degree(s.N))
   else
@@ -218,22 +218,19 @@ function isstrictlyproper{S}(s::MatrixFractionDescription{Val{:mimo},Val{S},Val{
 end
 
 # Enforce rational transfer function type invariance
-function mfdcheck{M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}(
-  N::M1, D::M2, ::Type{Val{:lfd}})
+function mfdcheck(N::M1, D::M2, ::Type{Val{:lfd}}) where {M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}
   @assert size(N,1) == size(D,1) "MatrixFractionDescription: size(N,1) ≠ size(D,1)"
   @assert size(D,1) == size(D,2) "MatrixFractionDescription: size(D,1) ≠ size(D,2)"
 
   return size(N,1), size(N,2)
 end
-function mfdcheck{M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}(
-  N::M1, D::M2, ::Type{Val{:rfd}})
+function mfdcheck(N::M1, D::M2, ::Type{Val{:rfd}}) where {M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}
   @assert size(N,2) == size(D,2) "MatrixFractionDescription: size(N,2) ≠ size(D,2)"
   @assert size(D,1) == size(D,2) "MatrixFractionDescription: size(D,1) ≠ size(D,2)"
 
   return size(N,1), size(N,2)
 end
-function mfdcheck{T,M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}(
-  N::M1, D::M2, ::Type{Val{T}}, Ts::Real)
+function mfdcheck(N::M1, D::M2, ::Type{Val{T}}, Ts::Real) where {T,M1<:PolynomialMatrices.PolyMatrix,M2<:PolynomialMatrices.PolyMatrix}
   @assert Ts ≥ zero(Ts) && !isinf(Ts) "MatrixFractionDescription: Ts must be non-negative real number"
   return mfdcheck(N, D, Val{T})
 end
@@ -241,9 +238,9 @@ end
 # Outer constructors
 lfd(N::Poly, D::Poly)           = MatrixFractionDescription(N, D, Val{:lfd})
 lfd(N::Poly, D::Poly, Ts::Real) = MatrixFractionDescription(N, D, convert(Float64, Ts), Val{:lfd})
-lfd{M1<:PolyMatrix,M2<:PolyMatrix}(N::M1, D::M2) =
+lfd(N::M1, D::M2) where {M1<:PolyMatrix,M2<:PolyMatrix} =
   MatrixFractionDescription(N, D, Val{:lfd})
-lfd{M1<:PolyMatrix,M2<:PolyMatrix}(N::M1, D::M2, Ts::Real) =
+lfd(N::M1, D::M2, Ts::Real) where {M1<:PolyMatrix,M2<:PolyMatrix} =
   MatrixFractionDescription(N, D, convert(Float64, Ts), Val{:lfd})
 
 rfd(N::Poly, D::Poly) = MatrixFractionDescription(N, D, Val{:rfd})
@@ -252,22 +249,21 @@ rfd(N::PolyMatrix, D::PolyMatrix) = MatrixFractionDescription(N, D, Val{:rfd})
 rfd(N::PolyMatrix, D::PolyMatrix, Ts::Real) = MatrixFractionDescription(N, D, convert(Float64, Ts), Val{:rfd})
 
 # Vector constructors
-lfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2}) =
+lfd(N::AbstractVector{T1}, D::AbstractVector{T2}) where {T1<:Real, T2<:Real} =
   lfd(Poly(reverse(N), :s), Poly(reverse(D), :s))
-lfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real) =
+lfd(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real) where {T1<:Real, T2<:Real} =
   lfd(Poly(reverse(N), :z), Poly(reverse(D), :z), Ts)
-rfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2}) =
+rfd(N::AbstractVector{T1}, D::AbstractVector{T2}) where {T1<:Real, T2<:Real} =
   rfd(Poly(reverse(N), :s), Poly(reverse(D), :s))
-rfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real) =
+rfd(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real) where {T1<:Real, T2<:Real} =
   rfd(Poly(reverse(N), :z), Poly(reverse(D), :z), Ts)
 
-function lfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2},
-  Ts::Real, var::Symbol)
+function lfd(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real, var::Symbol) where {T1<:Real, T2<:Real}
   vars    = [:z̄,:q̄,:qinv,:zinv]
   @assert var ∈ vars string("tf: var ∉ ", vars)
 
-  nlast       = findlast(N)
-  dlast       = findlast(D)
+  nlast       = findlast(!iszero, N)
+  dlast       = findlast(!iszero, D)
   order       = max(nlast, dlast)
   N_          = zeros(T1, order)
   N_[1:nlast] = N[1:nlast]
@@ -277,13 +273,12 @@ function lfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2},
   lfd(Poly(reverse(N_), :z), Poly(reverse(D_), :z), Ts)
 end
 
-function rfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2},
-  Ts::Real, var::Symbol)
+function rfd(N::AbstractVector{T1}, D::AbstractVector{T2}, Ts::Real, var::Symbol) where {T1<:Real, T2<:Real}
   vars    = [:z̄,:q̄,:qinv,:zinv]
   @assert var ∈ vars string("tf: var ∉ ", vars)
 
-  nlast       = findlast(N)
-  dlast       = findlast(D)
+  nlast       = findlast(!iszero, N)
+  dlast       = findlast(!iszero, D)
   order       = max(nlast, dlast)
   N_          = zeros(T1, order)
   N_[1:nlast] = N[1:nlast]
@@ -294,16 +289,16 @@ function rfd{T1<:Real, T2<:Real}(N::AbstractVector{T1}, D::AbstractVector{T2},
 end
 
 # Interfaces
-samplingtime(s::MatrixFractionDescription)              = s.Ts
-islfd{T,S,L}(s::MatrixFractionDescription{T,S,Val{L}})  = false
-islfd{T,S}(s::MatrixFractionDescription{T,S,Val{:lfd}}) = true
-isrfd{T,S,L}(s::MatrixFractionDescription{T,S,Val{L}})  = !islfd(s)
-num(s::MatrixFractionDescription)                       = s.N
-den(s::MatrixFractionDescription)                       = s.D
+samplingtime(s::MatrixFractionDescription)                      = s.Ts
+islfd(s::MatrixFractionDescription{T,S,Val{L}}) where {T,S,L}   = false
+islfd(s::MatrixFractionDescription{T,S,Val{:lfd}}) where {T,S}  = true
+isrfd(s::MatrixFractionDescription{T,S,Val{L}}) where {T,S,L}   = !islfd(s)
+num(s::MatrixFractionDescription)                               = s.N
+den(s::MatrixFractionDescription)                               = s.D
 
 # Think carefully about how to implement numstates
-numstates{T,S}(s::MatrixFractionDescription{T,S,Val{:lfd}}) = sum(row_degree(s.D))
-numstates{T,S}(s::MatrixFractionDescription{T,S,Val{:rfd}}) = sum(col_degree(s.D))
+numstates(s::MatrixFractionDescription{T,S,Val{:lfd}}) where {T,S} = sum(row_degree(s.D))
+numstates(s::MatrixFractionDescription{T,S,Val{:rfd}}) where {T,S} = sum(col_degree(s.D))
 # Currently, we only allow for proper systems
 numinputs(s::MatrixFractionDescription)               = s.nu
 numoutputs(s::MatrixFractionDescription)              = s.ny
@@ -314,27 +309,25 @@ size(s::MatrixFractionDescription, d)                 = size(s.N, d)
 length(s::MatrixFractionDescription{Val{:mimo}})      = length(s.N)
 
 # conversion between 1×1 mimo and siso
-function siso{L}(s::MatrixFractionDescription{Val{:mimo},Val{:cont},Val{L}})
+function siso(s::MatrixFractionDescription{Val{:mimo},Val{:cont},Val{L}}) where {L}
   if size(s) != (1,1)
-    warn("siso(s): system is not 1×1")
-    throw(DomainError())
+    throw(DomainError(s, "siso(s): system is not 1×1"))
   end
   MatrixFractionDescription(s.N[1], s.D[1], Val{L})
 end
 
-function siso{L}(s::MatrixFractionDescription{Val{:mimo},Val{:disc},Val{L}})
+function siso(s::MatrixFractionDescription{Val{:mimo},Val{:disc},Val{L}}) where {L}
   if size(s) != (1,1)
-    warn("siso(s): system is not 1×1")
-    throw(DomainError())
+    throw(DomainError(s, "siso(s): system is not 1×1"))
   end
   MatrixFractionDescription(s.N[1], s.D[1], s.Ts, Val{L})
 end
 
-function mimo{L}(s::MatrixFractionDescription{Val{:siso},Val{:cont},Val{L}})
+function mimo(s::MatrixFractionDescription{Val{:siso},Val{:cont},Val{L}}) where {L}
   MatrixFractionDescription(PolyMatrix(s.N, (1,1), Val{:s}) , PolyMatrix(s.D, (1,1), Val{:s}), Val{L})
 end
 
-function mimo{L}(s::MatrixFractionDescription{Val{:siso},Val{:disc},Val{L}})
+function mimo(s::MatrixFractionDescription{Val{:siso},Val{:disc},Val{L}}) where {L}
   MatrixFractionDescription(PolyMatrix(s.N, (1,1), Val{:z}) , PolyMatrix(s.D, (1,1), Val{:z}), s.Ts, Val{L})
 end
 
@@ -357,12 +350,10 @@ end
 #   rfd(ss(s)[I...])
 # end
 #
-# endof(s::MatrixFractionDescription{Val{:mimo}}) = endof(s.N)
-#
 # Conversion and promotion
-promote_rule{T<:Real,S,L}(::Type{T}, ::Type{MatrixFractionDescription{Val{:siso},S,L}}) =
+promote_rule(::Type{T}, ::Type{MatrixFractionDescription{Val{:siso},S,L}}) where {T<:Real,S,L} =
   MatrixFractionDescription{Val{:siso},S,L}
-promote_rule{T<:AbstractMatrix,S,L}(::Type{T}, ::Type{MatrixFractionDescription{Val{:mimo},S,L}}) =
+promote_rule(::Type{T}, ::Type{MatrixFractionDescription{Val{:mimo},S,L}}) where {T<:AbstractMatrix,S,L} =
   MatrixFractionDescription{Val{:mimo},S,L}
 
 convert(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd}}}, g::Real)            =
@@ -370,47 +361,47 @@ convert(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd}}}, g::R
 convert(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd}}}, g::Real)            =
   lfd(Poly(g,:z), Poly(one(g),:z), zero(Float64))
 convert(::Type{MatrixFractionDescription{Val{:mimo},Val{:cont},Val{:lfd}}}, g::AbstractMatrix)  =
-  lfd(PolyMatrix(g, Val{:s}), PolyMatrix(eye(eltype(g), size(g,1), size(g,1))))
+  lfd(PolyMatrix(g, Val{:s}), PolyMatrix(I(eltype(g), size(g,1), size(g,1))))
 convert(::Type{MatrixFractionDescription{Val{:mimo},Val{:disc},Val{:lfd}}}, g::AbstractMatrix)  =
-  lfd(PolyMatrix(g, Val{:z}), PolyMatrix(eye(eltype(g), size(g,1), size(g,1))), zero(Float64))
+  lfd(PolyMatrix(g, Val{:z}), PolyMatrix(I(eltype(g), size(g,1), size(g,1))), zero(Float64))
 
 convert(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd}}}, g::Real)            =
   rfd(Poly(g,:s), Poly(one(g),:s))
 convert(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd}}}, g::Real)            =
   rfd(Poly(g,:z), Poly(one(g),:z), zero(Float64))
 convert(::Type{MatrixFractionDescription{Val{:mimo},Val{:cont},Val{:rfd}}}, g::AbstractMatrix)  =
-  rfd(PolyMatrix(g, Val{:s}), PolyMatrix(eye(eltype(g), size(g,2), size(g,2))))
+  rfd(PolyMatrix(g, Val{:s}), PolyMatrix(I(eltype(g), size(g,2), size(g,2))))
 convert(::Type{MatrixFractionDescription{Val{:mimo},Val{:disc},Val{:rfd}}}, g::AbstractMatrix)  =
-  rfd(PolyMatrix(g, Val{:z}), PolyMatrix(eye(eltype(g), size(g,2), size(g,2))), zero(Float64))
+  rfd(PolyMatrix(g, Val{:z}), PolyMatrix(I(eltype(g), size(g,2), size(g,2))), zero(Float64))
 
 # conversions between lfd and rfd
-convert{S,T,L}(::Type{MatrixFractionDescription{Val{S},Val{T},Val{:lfd}}},
-  s::MatrixFractionDescription{Val{S},Val{T},Val{L}}) = lfd(s)
-convert{S,T,L}(::Type{MatrixFractionDescription{Val{S},Val{T},Val{:rfd}}},
-  s::MatrixFractionDescription{Val{S},Val{T},Val{L}}) = rfd(s)
+convert(::Type{MatrixFractionDescription{Val{S},Val{T},Val{:lfd}}},
+  s::MatrixFractionDescription{Val{S},Val{T},Val{L}}) where {S,T,L} = lfd(s)
+convert(::Type{MatrixFractionDescription{Val{S},Val{T},Val{:rfd}}},
+  s::MatrixFractionDescription{Val{S},Val{T},Val{L}}) where {S,T,L} = rfd(s)
 
 # Multiplicative and additive identities (meaningful only for SISO)
-one{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd},M1,M2}})  =
+one(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd},M1,M2}}) where {M1,M2} =
   lfd(one(Int8))
-one{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd},M1,M2}})  =
+one(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd},M1,M2}}) where {M1,M2} =
   lfd(one(Int8), zero(Float64))
-zero{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd},M1,M2}}) =
+zero(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd},M1,M2}}) where {M1,M2} =
   lfd(zero(Int8))
-zero{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd},M1,M2}}) =
+zero(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd},M1,M2}}) where {M1,M2} =
   lfd(zero(Int8), zero(Float64))
 
-one{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd},M1,M2}})  =
+one(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd},M1,M2}}) where {M1,M2} =
   rfd(one(Int8))
-one{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd},M1,M2}})  =
+one(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd},M1,M2}}) where {M1,M2} =
   rfd(one(Int8), zero(Float64))
-zero{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd},M1,M2}}) =
+zero(::Type{MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd},M1,M2}}) where {M1,M2} =
   rfd(zero(Int8))
-zero{M1,M2}(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd},M1,M2}}) =
+zero(::Type{MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd},M1,M2}}) where {M1,M2} =
   tf(zero(Int8), zero(Float64))
 
-one{T,S,L,M1,M2}(s::MatrixFractionDescription{Val{T},Val{S},Val{L},M1,M2})   =
+one(s::MatrixFractionDescription{Val{T},Val{S},Val{L},M1,M2}) where {T,S,L,M1,M2} =
   one(typeof(s))
-zero{T,S,L,M1,M2}(s::MatrixFractionDescription{Val{T},Val{S},Val{L},M1,M2})  =
+zero(s::MatrixFractionDescription{Val{T},Val{S},Val{L},M1,M2}) where {T,S,L,M1,M2} =
   zero(typeof(s))
 
 # conversions between lfd and rfd
@@ -418,15 +409,15 @@ lfd(s::MatrixFractionDescription{Val{:siso},Val{:cont},Val{:rfd}}) = lfd(s.N,s.D
 lfd(s::MatrixFractionDescription{Val{:siso},Val{:disc},Val{:rfd}}) = lfd(s.N,s.D,s.Ts)
 lfd(s::MatrixFractionDescription{Val{:mimo},Val{:cont},Val{:rfd}}) = lfd(_rfd2lfd(s)...)
 lfd(s::MatrixFractionDescription{Val{:mimo},Val{:disc},Val{:rfd}}) = lfd(_rfd2lfd(s)...,s.Ts)
-lfd{T,S}(s::MatrixFractionDescription{T,S,Val{:lfd}})              = s
+lfd(s::MatrixFractionDescription{T,S,Val{:lfd}}) where {T,S}             = s
 
 rfd(s::MatrixFractionDescription{Val{:siso},Val{:cont},Val{:lfd}}) = rfd(s.N,s.D)
 rfd(s::MatrixFractionDescription{Val{:siso},Val{:disc},Val{:lfd}}) = rfd(s.N,s.D,s.Ts)
 rfd(s::MatrixFractionDescription{Val{:mimo},Val{:cont},Val{:lfd}}) = rfd(_lfd2rfd(s)...)
 rfd(s::MatrixFractionDescription{Val{:mimo},Val{:disc},Val{:lfd}}) = rfd(_lfd2rfd(s)...,s.Ts)
-rfd{T,S}(s::MatrixFractionDescription{T,S,Val{:rfd}})              = s
+rfd(s::MatrixFractionDescription{T,S,Val{:rfd}}) where {T,S}       = s
 
-function _lfd2rfd{T,S}(s::MatrixFractionDescription{T,S,Val{:lfd}})
+function _lfd2rfd(s::MatrixFractionDescription{T,S,Val{:lfd}}) where {T,S}
   n,m = size(s)
   p   = hcat(-s.N, s.D)
   _,U = ltriang(p)
@@ -436,7 +427,7 @@ function _lfd2rfd{T,S}(s::MatrixFractionDescription{T,S,Val{:lfd}})
   return N,D
 end
 
-function _rfd2lfd{T,S}(s::MatrixFractionDescription{T,S,Val{:rfd}})
+function _rfd2lfd(s::MatrixFractionDescription{T,S,Val{:rfd}}) where {T,S}
   n,m = size(s)
   p   = vcat(-s.D, s.N)
   _,U = rtriang(p)
@@ -446,13 +437,13 @@ function _rfd2lfd{T,S}(s::MatrixFractionDescription{T,S,Val{:rfd}})
   return N,D
 end
 
-function inv{M<:MatrixFractionDescription}(s::M)
+function inv(s::MatrixFractionDescription)
   _mfdinvcheck(s)
   _inv(s)
 end
 
-_inv{T,L}(s::MatrixFractionDescription{Val{T},Val{:cont},Val{L}}) = MatrixFractionDescription(copy(s.D), copy(s.N), Val{L})
-_inv{T,L}(s::MatrixFractionDescription{Val{T},Val{:disc},Val{L}}) = MatrixFractionDescription(copy(s.D), copy(s.N), s.Ts, Val{L})
+_inv(s::MatrixFractionDescription{Val{T},Val{:cont},Val{L}}) where {T,L} = MatrixFractionDescription(copy(s.D), copy(s.N), Val{L})
+_inv(s::MatrixFractionDescription{Val{T},Val{:disc},Val{L}}) where {T,L} = MatrixFractionDescription(copy(s.D), copy(s.N), s.Ts, Val{L})
 
 function _mfdinvcheck(s::MatrixFractionDescription)
   if s.ny ≠ s.nu
@@ -467,14 +458,14 @@ function _mfdinvcheck(s::MatrixFractionDescription)
 end
 
 # Negative of a transfer-function model
--{T}(s::MatrixFractionDescription{Val{T},Val{:cont},Val{:lfd}}) = lfd(-s.N, copy(s.D))
--{T}(s::MatrixFractionDescription{Val{T},Val{:disc},Val{:lfd}}) = lfd(-s.N, copy(s.D), s.Ts)
--{T}(s::MatrixFractionDescription{Val{T},Val{:cont},Val{:rfd}}) = rfd(-s.N, copy(s.D))
--{T}(s::MatrixFractionDescription{Val{T},Val{:disc},Val{:rfd}}) = rfd(-s.N, copy(s.D), s.Ts)
+-(s::MatrixFractionDescription{Val{T},Val{:cont},Val{:lfd}}) where {T} = lfd(-s.N, copy(s.D))
+-(s::MatrixFractionDescription{Val{T},Val{:disc},Val{:lfd}}) where {T} = lfd(-s.N, copy(s.D), s.Ts)
+-(s::MatrixFractionDescription{Val{T},Val{:cont},Val{:rfd}}) where {T} = rfd(-s.N, copy(s.D))
+-(s::MatrixFractionDescription{Val{T},Val{:disc},Val{:rfd}}) where {T} = rfd(-s.N, copy(s.D), s.Ts)
 
 # Addition (parallel)
-function _mfdparallelcheck{T1,T2,S,L}(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L}},
-  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L}})
+function _mfdparallelcheck(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L}},
+  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L}}) where {T1,T2,S,L}
   if s₁.Ts ≉ s₂.Ts && s₁.Ts ≠ zero(s₁.Ts) && s₂.Ts ≠ zero(s₂.Ts)
     warn("parallel(s₁,s₂): Sampling time mismatch")
     throw(DomainError())
@@ -487,8 +478,8 @@ function _mfdparallelcheck{T1,T2,S,L}(s₁::MatrixFractionDescription{Val{T1},Va
 end
 
 # siso version
-function _mfdparallel{S,L}(s₁::MatrixFractionDescription{Val{:siso},Val{S},Val{:L}},
-  s₂::MatrixFractionDescription{Val{:siso},Val{S},Val{L}})
+function _mfdparallel(s₁::MatrixFractionDescription{Val{:siso},Val{S},Val{:L}},
+  s₂::MatrixFractionDescription{Val{:siso},Val{S},Val{L}}) where {S,L}
   R   = gcd(s₁.D, s₂.D)
   D₁  = div(s₁.D, R)
   D   = D₁*s₂.D        # only include common part R once
@@ -498,8 +489,8 @@ function _mfdparallel{S,L}(s₁::MatrixFractionDescription{Val{:siso},Val{S},Val
 end
 
 # mimo lfd version
-function _mfdparallel{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}})
+function _mfdparallel(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}}) where {S}
   R, V₁, V₂ = gcrd(s₁.D, s₂.D)
   detV₁, adjV₁ = inv(V₁)
   detV₂, adjV₂ = inv(V₂)
@@ -509,8 +500,8 @@ function _mfdparallel{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:
 end
 
 # mimo rfd version
-function _mfdparallel{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
+function _mfdparallel(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}}) where {S}
   L, V₁, V₂ = gcld(s₁.D, s₂.D)
   detV₁, adjV₁ = inv(V₁)
   detV₂, adjV₂ = inv(V₂)
@@ -520,57 +511,47 @@ function _mfdparallel{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:
 end
 
 # mimo mixed lfd/rfd version
-function _mfdparallel{S,L1,L2}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{L1}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{L2}})
+function _mfdparallel(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{L1}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{L2}}) where {S,L1,L2}
   _mfdparallel(lfd(s₁), lfd(s₂))
 end
 
 # siso and mimo of dimensions 1×1
-function _mfdparallel{T1,T2,S,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L2}})
+function _mfdparallel(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L2}}) where {T1,T2,S,L1,L2}
   _mfdparallel(mimo(s₁), mimo(s₂))
 end
 
-function +{T1,T2,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{:cont},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{:cont},Val{L2}})
+function +(s₁::MatrixFractionDescription{Val{T1},Val{:cont},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{:cont},Val{L2}}) where {T1,T2,L1,L2}
   _mfdparallelcheck(s₁, s₂)
   N, D, _ = _mfdparallel(s₁, s₂)
   MatrixFractionDescription(N, D, Val{L1})
 end
 
-function +{T1,T2,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{:disc},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{:disc},Val{L2}})
+function +(s₁::MatrixFractionDescription{Val{T1},Val{:disc},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{:disc},Val{L2}}) where {T1,T2,L1,L2}
   _mfdparallelcheck(s₁, s₂)
   N, D, Ts = _mfdparallel(s₁, s₂)
   MatrixFractionDescription(N, D, Ts, Val{L1})
 end
 
-.+(s₁::MatrixFractionDescription{Val{:siso}}, s₂::MatrixFractionDescription{Val{:siso}}) = +(s₁, s₂)
-
-+{T,S}(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) =
++(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) where {T,S} =
   +(s, convert(typeof(s), g))
-+{T,S}(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) =
++(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) where {T,S} =
   +(convert(typeof(s), g), s)
 
-.+(s::MatrixFractionDescription{Val{:siso}}, g::Real) = +(s, g)
-.+(g::Real, s::MatrixFractionDescription{Val{:siso}}) = +(g, s)
-
 # Subtraction
 -(s₁::MatrixFractionDescription, s₂::MatrixFractionDescription) = +(s₁, -s₂)
 
-.-(s₁::MatrixFractionDescription{Val{:siso}}, s₂::MatrixFractionDescription{Val{:siso}}) = -(s₁, s₂)
-
--{T,S}(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) =
+-(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) where {T,S} =
   -(s, convert(typeof(s), g))
--{T,S}(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) =
+-(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) where {T,S} =
   -(convert(typeof(s), g), s)
 
-.-(s::MatrixFractionDescription{Val{:siso}}, g::Real)    = -(s, g)
-.-(g::Real, s::MatrixFractionDescription{Val{:siso}})    = -(g, s)
-
 # Multiplication
-function _mfdseriescheck{T1,T2,S}(s₁::MatrixFractionDescription{Val{T1},Val{S}},
-  s₂::MatrixFractionDescription{Val{T2},Val{S}})
+function _mfdseriescheck(s₁::MatrixFractionDescription{Val{T1},Val{S}},
+  s₂::MatrixFractionDescription{Val{T2},Val{S}}) where {T1,T2,S}
   # Remark: s₁*s₂ implies u -> s₂ -> s₁ -> y
 
   if s₁.Ts ≉ s₂.Ts && s₁.Ts ≠ zero(s₁.Ts) && s₂.Ts ≠ zero(s₂.Ts)
@@ -585,8 +566,8 @@ function _mfdseriescheck{T1,T2,S}(s₁::MatrixFractionDescription{Val{T1},Val{S}
 end
 
 # siso version
-function _mfdseries{S,L1,L2}(s₁::MatrixFractionDescription{Val{:siso},Val{S},Val{L1}},
-  s₂::MatrixFractionDescription{Val{:siso},Val{S},Val{L2}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{:siso},Val{S},Val{L1}},
+  s₂::MatrixFractionDescription{Val{:siso},Val{S},Val{L2}}) where {S,L1,L2}
   R₁  = gcd(s₁.D, s₂.N)
   R₂  = gcd(s₂.D, s₁.N)
   D   = div(s₁.D, R₁)*div(s₂.D, R₂)
@@ -595,8 +576,8 @@ function _mfdseries{S,L1,L2}(s₁::MatrixFractionDescription{Val{:siso},Val{S},V
 end
 
 # mimo lfd version
-function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{lfd}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{lfd}}) where {S}
   # D₁^-1 N₁ D₂^-1 N₂
   sᵢ  = lfd(rfd(s₁.N, s₂.D))
   # D₁^-1 Dᵢ^-1 Nᵢ N₂
@@ -608,8 +589,8 @@ function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lf
 end
 
 # mimo rfd version
-function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}}) where {S}
   # N₁ D₁^-1 N₂ D₂^-1
   sᵢ  = rfd(lfd(s₂.N, s₁.D))
   # N₁ Nᵢ Dᵢ^-1 D₂^-1
@@ -621,8 +602,8 @@ function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rf
 end
 
 # mimo mixed lfd/rfd versions
-function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}}) where {S}
   # N₁ D₁^-1 D₂^-1 N₂
   Dᵢ  = s₂.D*s₁.D
   sᵢ  = lfd(rfd(s₂.N, Dᵢ))
@@ -633,8 +614,8 @@ function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lf
   V₁, V₂, max(s₁.Ts, s₂.Ts)
 end
 
-function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
-  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rfd}},
+  s₂::MatrixFractionDescription{Val{:mimo},Val{S},Val{:lfd}}) where {S}
   # D₁^-1 N₁ N₂ D₂^-1
   Nᵢ  = s₁.N*s₂.N
   sᵢ  = rfd(lfd(Nᵢ, s₁.D))
@@ -646,46 +627,41 @@ function _mfdseries{S}(s₁::MatrixFractionDescription{Val{:mimo},Val{S},Val{:rf
 end
 
 # siso and mimo of dimensions 1×1
-function _mfdseries{T1,T2,S,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L2}})
+function _mfdseries(s₁::MatrixFractionDescription{Val{T1},Val{S},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{S},Val{L2}}) where {T1,T2,S,L1,L2}
   _mfdseries(mimo(s₁), mimo(s₂))
 end
 
-function *{T1,T2,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{:cont},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{:cont},Val{L2}})
+function *(s₁::MatrixFractionDescription{Val{T1},Val{:cont},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{:cont},Val{L2}}) where {T1,T2,L1,L2}
   _mfdseriescheck(s₁, s₂)
   N, D, _ = _mfdseries(s₁, s₂)
   MatrixFractionDescription(N, D, Val{L1})
 end
 
-function *{T1,T2,L1,L2}(s₁::MatrixFractionDescription{Val{T1},Val{:disc},Val{L1}},
-  s₂::MatrixFractionDescription{Val{T2},Val{:disc},Val{L2}})
+function *(s₁::MatrixFractionDescription{Val{T1},Val{:disc},Val{L1}},
+  s₂::MatrixFractionDescription{Val{T2},Val{:disc},Val{L2}}) where {T1,T2,L1,L2}
   _mfdseriescheck(s₁, s₂)
   N, D, Ts = _mfdseries(s₁, s₂)
   MatrixFractionDescription(N, D, Ts, Val{L1})
 end
 
-.*(s₁::MatrixFractionDescription{Val{:siso}}, s₂::MatrixFractionDescription{Val{:siso}}) = *(s₁, s₂)
-
-*{T,S}(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) =
+*(s::MatrixFractionDescription{Val{T},Val{S}}, g::Union{Real,AbstractMatrix}) where {T,S} =
   *(s, convert(typeof(s), g))
-*{T,S}(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) =
+*(g::Union{Real,AbstractMatrix}, s::MatrixFractionDescription{Val{T},Val{S}}) where {T,S} =
   *(convert(typeof(s), g), s)
 
-.*(s::MatrixFractionDescription{Val{:siso}}, g::Real)    = *(s, g)
-.*(g::Real, s::MatrixFractionDescription{Val{:siso}})    = *(g, s)
-
 ## Comparison
-=={T,S,L}(s₁::MatrixFractionDescription{T,S,L}, s₂::MatrixFractionDescription{T,S,L}) =
+==(s₁::MatrixFractionDescription{T,S,L}, s₂::MatrixFractionDescription{T,S,L}) where {T,S,L} =
   (s₁.N == s₂.N) && (s₁.D == s₂.D) && (s₁.Ts == s₂.Ts)
-=={T1,S1,L1,T2,S2,L2}(s₁::MatrixFractionDescription{T1,S1,L1}, s₂::MatrixFractionDescription{T2,S2,L2}) = false
+==(s₁::MatrixFractionDescription{T1,S1,L1}, s₂::MatrixFractionDescription{T2,S2,L2}) where {T1,S1,L1,T2,S2,L2} = false
 
 hash(s::MatrixFractionDescription, h::UInt)     = hash(s.D, hash(S.N, hash(S.Ts, h)))
 isequal(s₁::MatrixFractionDescription, s₂::MatrixFractionDescription) = (hash(s₁) == hash(s₂))
 
-function isapprox{T,S,L1,L2,M1,M2,M3,M4}(s₁::MatrixFractionDescription{T,S,L1,M1,M2}, s₂::MatrixFractionDescription{T,S,L2,M3,M4};
+function isapprox(s₁::MatrixFractionDescription{T,S,L1,M1,M2}, s₂::MatrixFractionDescription{T,S,L2,M3,M4};
   rtol::Real=Base.rtoldefault(promote_type(eltype(mattype(s₁.N)),eltype(mattype(s₁.D)),eltype(mattype(s₂.N)),eltype(mattype(s₂.D)))),
-  atol::Real=0, norm::Function=vecnorm)
+  atol::Real=0, norm::Function=vecnorm) where {T,S,L1,L2,M1,M2,M3,M4}
   isapprox(s₁.Ts, s₂.Ts) || return false # quick exit
   lfd1 = lfd(s₁)
   lfd2 = lfd(s₂)
diff --git a/src/types/system/printing.jl b/src/types/system/printing.jl
index 5a60c96..0d89894 100644
--- a/src/types/system/printing.jl
+++ b/src/types/system/printing.jl
@@ -1,28 +1,28 @@
-_dimensions(s::LtiSystem{Val{:siso}})     = ""
-_dimensions(s::LtiSystem{Val{:mimo}})     = "$(numoutputs(s))×$(numinputs(s))"
-_time{T}(::LtiSystem{Val{T},Val{:cont}}) = "continuous"
-_time{T}(::LtiSystem{Val{T},Val{:disc}}) = "discrete"
-_timevar{T}(::LtiSystem{Val{T},Val{:cont}}) = "t"
-_timevar{T}(::LtiSystem{Val{T},Val{:disc}}) = "k"
-_printsamplingtime{T}(s::LtiSystem{Val{T},Val{:cont}}) = ""
-_printsamplingtime{T}(s::LtiSystem{Val{T},Val{:disc}}) = "+$(samplingtime(s))"
+_dimensions(s::LtiSystem{Val{:siso}})              = ""
+_dimensions(s::LtiSystem{Val{:mimo}})              = "$(numoutputs(s))×$(numinputs(s))"
+_time(::LtiSystem{Val{T},Val{:cont}}) where {T}    = "continuous"
+_time(::LtiSystem{Val{T},Val{:disc}}) where {T}    = "discrete"
+_timevar(::LtiSystem{Val{T},Val{:cont}}) where {T} = "t"
+_timevar(::LtiSystem{Val{T},Val{:disc}}) where {T} = "k"
+_printsamplingtime(s::LtiSystem{Val{T},Val{:cont}}) where {T} = ""
+_printsamplingtime(s::LtiSystem{Val{T},Val{:disc}}) where {T} = "+$(samplingtime(s))"
 
 import Base: summary
 
 # StateSpace printing
-_printxupdate{T}(s::StateSpace{Val{T},Val{:cont}})  = "ẋ(t)"
-_printxupdate{T}(s::StateSpace{Val{T},Val{:disc}})  = "x(k+1)"
+_printxupdate(s::StateSpace{Val{T},Val{:cont}}) where {T} = "ẋ(t)"
+_printxupdate(s::StateSpace{Val{T},Val{:disc}}) where {T} = "x(k+1)"
 
-summary{T}(s::StateSpace{Val{:siso},Val{T}})  = "StateSpace"
-summary{T}(s::StateSpace{Val{:mimo},Val{T}})  = "$(_dimensions(s)) StateSpace"
+summary(s::StateSpace{Val{:siso},Val{T}}) where {T} = "StateSpace"
+summary(s::StateSpace{Val{:mimo},Val{T}}) where {T} = "$(_dimensions(s)) StateSpace"
 
 # Compact representations
-function _compact{T,S}(stream, ::MIME"text/plain", s::StateSpace{Val{T},Val{S}})
+function _compact(stream, ::MIME"text/plain", s::StateSpace{Val{T},Val{S}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "G($(var))")
 end
 
-function _compact{T,S}(stream, ::MIME"text/latex", s::StateSpace{Val{T},Val{S}})
+function _compact(stream, ::MIME"text/latex", s::StateSpace{Val{T},Val{S}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "\$")
   print(stream, "G($(var))")
@@ -32,7 +32,7 @@ end
 # TODO: Think about text/html
 
 # Full representations
-function _full{T,S}(stream, m::MIME"text/plain", s::StateSpace{Val{T},Val{S}})
+function _full(stream, m::MIME"text/plain", s::StateSpace{Val{T},Val{S}}) where {T,S}
   tvar = _timevar(s)
   println(stream, summary(s))
   println(stream, "$(_printxupdate(s)) = Ax($tvar) + Bu($tvar,x)")
@@ -40,7 +40,7 @@ function _full{T,S}(stream, m::MIME"text/plain", s::StateSpace{Val{T},Val{S}})
   println(stream, "with $(numstates(s)) states in $(_time(s)) time.")
 end
 
-function _full{T,S}(stream, m::MIME"text/latex", s::StateSpace{Val{T},Val{S}})
+function _full(stream, m::MIME"text/latex", s::StateSpace{Val{T},Val{S}}) where {T,S}
   tvar = _timevar(s)
   println(stream, summary(s))
   print(stream, "\\begin{align*}")
@@ -61,16 +61,16 @@ end
   get(stream, :compact, false) ? _compact(stream, mime, s) : _full(stream, mime, s)
 
 ## TransferFunction printing
-summary{T}(s::TransferFunction{Val{:siso},Val{T}})  = "TransferFunction"
-summary{T}(s::TransferFunction{Val{:mimo},Val{T}})  = "$(_dimensions(s)) TransferFunction"
+summary(s::TransferFunction{Val{:siso},Val{T}}) where {T} = "TransferFunction"
+summary(s::TransferFunction{Val{:mimo},Val{T}}) where {T} = "$(_dimensions(s)) TransferFunction"
 
 # Compact representations
-function _compact{T,S}(stream, ::MIME"text/plain", s::TransferFunction{Val{T},Val{S}})
+function _compact(stream, ::MIME"text/plain", s::TransferFunction{Val{T},Val{S}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "n($(var))/d($(var))")
 end
 
-function _compact{T,S}(stream, ::MIME"text/latex", s::TransferFunction{Val{T},Val{S}})
+function _compact(stream, ::MIME"text/latex", s::TransferFunction{Val{T},Val{S}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "\$")
   print(stream, "\\tfrac{\\mathrm{n}($(var))}{\\mathrm{d}($(var))}")
@@ -80,7 +80,7 @@ end
 # TODO: Think about text/html
 
 # Full representations
-function _full{T,S}(stream, m::MIME"text/plain", s::TransferFunction{Val{T},Val{S}})
+function _full(stream, m::MIME"text/plain", s::TransferFunction{Val{T},Val{S}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
@@ -88,7 +88,7 @@ function _full{T,S}(stream, m::MIME"text/plain", s::TransferFunction{Val{T},Val{
   println(stream, "in $(_time(s)) time.")
 end
 
-function _full{T,S}(stream, m::MIME"text/latex", s::TransferFunction{Val{T},Val{S}})
+function _full(stream, m::MIME"text/latex", s::TransferFunction{Val{T},Val{S}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
@@ -110,32 +110,32 @@ end
 
 ## MatrixFractionDescription printing
 
-_mfdtype{S,T,L}(::MatrixFractionDescription{Val{S},Val{T},Val{L}}) = ifelse(L == :lfd, "Left", "Right")
+_mfdtype(::MatrixFractionDescription{Val{S},Val{T},Val{L}}) where {S,T,L} = ifelse(L == :lfd, "Left", "Right")
 
-summary{T,L}(s::MatrixFractionDescription{Val{:siso},Val{T},Val{L}}) =
+summary(s::MatrixFractionDescription{Val{:siso},Val{T},Val{L}}) where {T,L} =
   "$(_mfdtype(s)) MatrixFractionDescription"
-summary{T,L}(s::MatrixFractionDescription{Val{:mimo},Val{T},Val{L}}) =
+summary(s::MatrixFractionDescription{Val{:mimo},Val{T},Val{L}}) where {T,L} =
   "$(_dimensions(s)) $(_mfdtype(s)) MatrixFractionDescription"
 
 # Compact representations
-function _compact{T,S}(stream, ::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}})
+function _compact(stream, ::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "d($(var))", "\\", "n($(var))")
 end
 
-function _compact{T,S}(stream, ::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}})
+function _compact(stream, ::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "n($(var))/d($(var))")
 end
 
-function _compact{T,S}(stream, ::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}})
+function _compact(stream, ::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "\$")
   print(stream, "d($(var))", "\\", "n($(var))")
   print(stream, "\$")
 end
 
-function _compact{T,S}(stream, ::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}})
+function _compact(stream, ::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}}) where {T,S}
   var = ifelse(S == :cont, "s", "z")
   print(stream, "\$")
   print(stream, "n($(var))/d($(var))")
@@ -145,7 +145,7 @@ end
 # TODO: Think about text/html
 
 # Full representations
-function _full{T,S}(stream, m::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}})
+function _full(stream, m::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
@@ -153,7 +153,7 @@ function _full{T,S}(stream, m::MIME"text/plain", s::MatrixFractionDescription{Va
   println(stream, "in $(_time(s)) time.")
 end
 
-function _full{T,S}(stream, m::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}})
+function _full(stream, m::MIME"text/plain", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
@@ -161,7 +161,7 @@ function _full{T,S}(stream, m::MIME"text/plain", s::MatrixFractionDescription{Va
   println(stream, "in $(_time(s)) time.")
 end
 
-function _full{T,S}(stream, m::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}})
+function _full(stream, m::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:lfd}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
@@ -171,7 +171,7 @@ function _full{T,S}(stream, m::MIME"text/latex", s::MatrixFractionDescription{Va
   println(stream, "in $(_time(s)) time.")
 end
 
-function _full{T,S}(stream, m::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}})
+function _full(stream, m::MIME"text/latex", s::MatrixFractionDescription{Val{T},Val{S},Val{:rfd}}) where {T,S}
   var  = ifelse(S == :cont, "s", "z")
   tvar = _timevar(s)
   println(stream, summary(s))
diff --git a/src/types/system/statespace.jl b/src/types/system/statespace.jl
index ce41de0..42ec399 100644
--- a/src/types/system/statespace.jl
+++ b/src/types/system/statespace.jl
@@ -1,4 +1,4 @@
-immutable StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
+struct StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
   A::M1
   B::M2
   C::M3
@@ -10,16 +10,16 @@ immutable StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
 
   # Constructors
   # Continuous-time, single-input-single-output state-space model
-  function (::Type{StateSpace}){M1<:AbstractMatrix,M2<:AbstractMatrix,
-    M3<:AbstractMatrix}(A::M1, B::M2, C::M3, D::Real)
+  function (::Type{StateSpace})(A::M1, B::M2, C::M3, D::Real) where {M1<:AbstractMatrix,M2<:AbstractMatrix,
+    M3<:AbstractMatrix}
     d = fill(D,1,1)
     nx, nu, ny = _sscheck(A, B, C, d)
     new{Val{:siso},Val{:cont},M1,M2,M3,typeof(d)}(A, B, C, d, nx, nu, ny, zero(Float64))
   end
 
   # Discrete-time, single-input-single-output state-space model
-  function (::Type{StateSpace}){M1<:AbstractMatrix,M2<:AbstractMatrix,
-    M3<:AbstractMatrix}(A::M1, B::M2, C::M3, D::Real, Ts::Real)
+  function (::Type{StateSpace})(A::M1, B::M2, C::M3, D::Real, Ts::Real) where {M1<:AbstractMatrix,M2<:AbstractMatrix,
+    M3<:AbstractMatrix}
     d = fill(D,1,1)
     nx, nu, ny = _sscheck(A, B, C, d, Ts)
     new{Val{:siso},Val{:disc},M1,M2,M3,typeof(d)}(A, B, C, d, nx, nu, ny,
@@ -27,15 +27,15 @@ immutable StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
   end
 
   # Continuous-time, multi-input-multi-output state-space model
-  function (::Type{StateSpace}){M1<:AbstractMatrix,M2<:AbstractMatrix,
-    M3<:AbstractMatrix,M4<:AbstractMatrix}(A::M1, B::M2, C::M3, D::M4)
+  function (::Type{StateSpace})(A::M1, B::M2, C::M3, D::M4) where {M1<:AbstractMatrix,M2<:AbstractMatrix,
+    M3<:AbstractMatrix,M4<:AbstractMatrix}
     nx, nu, ny = _sscheck(A, B, C, D)
     new{Val{:mimo},Val{:cont},M1,M2,M3,M4}(A, B, C, D, nx, nu, ny, zero(Float64))
   end
 
   # Discrete-time, multi-input-multi-output state-space model
-  function (::Type{StateSpace}){M1<:AbstractMatrix,M2<:AbstractMatrix,
-    M3<:AbstractMatrix,M4<:AbstractMatrix}(A::M1, B::M2, C::M3, D::M4, Ts::Real)
+  function (::Type{StateSpace})(A::M1, B::M2, C::M3, D::M4, Ts::Real) where {M1<:AbstractMatrix,M2<:AbstractMatrix,
+    M3<:AbstractMatrix,M4<:AbstractMatrix}
     nx, nu, ny = _sscheck(A, B, C, D, Ts)
     new{Val{:mimo},Val{:disc},M1,M2,M3,M4}(A, B, C, D, nx, nu, ny,
       convert(Float64, Ts))
@@ -58,8 +58,7 @@ immutable StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
     ucalc = u(t, x.x)
     ucalc = isa(ucalc, Real) ? [ucalc] : ucalc
     if !isa(ucalc, AbstractVector) || length(ucalc) ≠ sys.nu || !(eltype(ucalc) <: Real)
-      warn("sys(t,x,dx,u): u(t,x) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values")
-      throw(DomainError())
+      throw(DomainError(u(t,x), "sys(t,x,dx,u): u(t,x) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values"))
     end
     x.u   = ucalc
     x.y   = C*x.x + D*x.u
@@ -69,16 +68,15 @@ immutable StateSpace{T,S,M1,M2,M3,M4} <: LtiSystem{T,S}
   ## Frequency response
   (sys::StateSpace{Val{:siso}})(x::Number)                                    =
     _eval(sys, x)[1]
-  (sys::StateSpace{Val{:siso}}){M<:Number}(X::AbstractArray{M})               =
+  (sys::StateSpace{Val{:siso}})(X::AbstractArray{M}) where {M<:Number}        =
     reshape(_eval(sys, X), size(X))
   (sys::StateSpace{Val{:mimo}})(x::Number)                                    =
     _eval(sys, x)
-  (sys::StateSpace{Val{:mimo}}){M<:Number}(X::AbstractArray{M})               =
+  (sys::StateSpace{Val{:mimo}})(X::AbstractArray{M}) where {M<:Number}        =
     _eval(sys, X)
-  function (sys::StateSpace){T<:Real}(; ω::Union{T, AbstractArray{T}} = Float64[])
+  function (sys::StateSpace)(; ω::Union{T, AbstractArray{T}} = Float64[]) where {T<:Real}
     if isempty(ω)
-      warn("sys(): Provide an argument for the function call. Refer to `?freqresp`.")
-      throw(DomainError())
+      throw(DomainError(w, "sys(): Provide an argument for the function call. Refer to `?freqresp`."))
     end
     freqresp(sys, ω)
   end
@@ -93,48 +91,39 @@ function _sscheck(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix,
   nd, md  = size(D)
 
   if Ts < zero(Ts) || isinf(Ts)
-    warn("StateSpace: Ts must be non-negative real number")
-    throw(DomainError())
+    throw(DomainError(Ts, "StateSpace: Ts must be non-negative real number"))
   end
 
   if na ≠ ma || !(eltype(A) <: Real)
-    warn("StateSpace: A must be a square matrix of real numbers")
-    throw(DomainError())
+    throw(DomainError(A, "StateSpace: A must be a square matrix of real numbers"))
   end
 
   if !(eltype(B) <: Real)
-    warn("StateSpace: B must be a matrix of real numbers")
-    throw(DomainError())
+    throw(DomainError(B, "StateSpace: B must be a matrix of real numbers"))
   end
 
   if !(eltype(C) <: Real)
-    warn("StateSpace: C must be a matrix of real numbers")
-    throw(DomainError())
+    throw(DomainError(C, "StateSpace: C must be a matrix of real numbers"))
   end
 
   if !(eltype(D) <: Real)
-    warn("StateSpace: D must be a matrix of real numbers")
-    throw(DomainError())
+    throw(DomainError(D, "StateSpace: D must be a matrix of real numbers"))
   end
 
   if na ≠ nb
-    warn("StateSpace: A and B must have the same number of rows")
-    throw(DomainError())
+    throw(DomainError((A, B), "StateSpace: A and B must have the same number of rows"))
   end
 
   if ma ≠ mc
-    warn("StateSpace: A and C must have the same number of columns")
-    throw(DomainError())
+    throw(DomainError((A, C), "StateSpace: A and C must have the same number of columns"))
   end
 
   if nc ≠ nd || nc < 1
-    warn("StateSpace: C and D must have the same number (≥1) of rows")
-    throw(DomainError())
+    throw(DomainError((C, D), "StateSpace: C and D must have the same number (≥1) of rows"))
   end
 
   if mb ≠ md || mb < 1
-    warn("StateSpace: B and D must have the same number (≥1) of columns")
-    throw(DomainError())
+    throw(DomainError((B, D), "StateSpace: B and D must have the same number (≥1) of columns"))
   end
 
   return na, mb, nc
@@ -148,8 +137,8 @@ ss(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::Real = zero(Float
 ss(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::Real, Ts::Real)    =
   StateSpace(A, B, C, D, Ts)
 
-ss{T<:Real}(D::T)           = StateSpace(zeros(T,0,0), zeros(T,0,1), zeros(T,1,0), D)
-ss{T<:Real}(D::T, Ts::Real) = StateSpace(zeros(T,0,0), zeros(T,0,1), zeros(T,1,0), D, Ts)
+ss(D::T) where {T<:Real}           = StateSpace(zeros(T,0,0), zeros(T,0,1), zeros(T,1,0), D)
+ss(D::T, Ts::Real) where {T<:Real} = StateSpace(zeros(T,0,0), zeros(T,0,1), zeros(T,1,0), D, Ts)
 
 # MIMO
 ss(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix)    =
@@ -172,9 +161,9 @@ function ss(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix,
   StateSpace(A, B, C, d, Ts)
 end
 
-ss{T<:Real}(D::AbstractMatrix{T})           = StateSpace(zeros(T,0,0),
+ss(D::AbstractMatrix{T}) where {T<:Real}           = StateSpace(zeros(T,0,0),
   zeros(T,0,size(D,2)), zeros(T,size(D,1),0), D)
-ss{T<:Real}(D::AbstractMatrix{T}, Ts::Real) = StateSpace(zeros(T,0,0),
+ss(D::AbstractMatrix{T}, Ts::Real) where {T<:Real} = StateSpace(zeros(T,0,0),
   zeros(T,0,size(D,2)), zeros(T,size(D,1),0), D, Ts)
 
 # # Catch-all for convenience when dealing with scalars, vectors, etc.
@@ -204,7 +193,7 @@ start(s::StateSpace{Val{:mimo}})        = start(s.D)
 next(s::StateSpace{Val{:mimo}}, state)  = (s[state], state+1)
 done(s::StateSpace{Val{:mimo}}, state)  = done(s.D, state)
 
-eltype{S,M1}(::Type{StateSpace{Val{:mimo},Val{S},M1}}) =
+eltype(::Type{StateSpace{Val{:mimo},Val{S},M1}}) where {S,M1} =
   StateSpace{Val{:siso},Val{S},M1}
 
 length(s::StateSpace{Val{:mimo}}) = length(s.D)
@@ -258,19 +247,20 @@ getindex(s::StateSpace{Val{:mimo}}, rows, ::Colon)    = s[rows, 1:s.nu]
 getindex(s::StateSpace{Val{:mimo}}, ::Colon, cols)    = s[1:s.ny, cols]
 getindex(s::StateSpace{Val{:mimo}}, ::Colon)          = s[1:end]
 getindex(s::StateSpace{Val{:mimo}}, ::Colon, ::Colon) = s[1:s.ny,1:s.nu]
-endof(s::StateSpace{Val{:mimo}})                      = endof(s.D)
+firstindex(s::StateSpace{Val{:mimo}})                 = firstindex(s.D)
+lastindex(s::StateSpace{Val{:mimo}})                  = lastindex(s.D)
 
 # Multiplicative and additive identities (meaningful only for SISO)
-one{M1,M2,M3,M4}(::Type{StateSpace{Val{:siso},Val{:cont},M1,M2,M3,M4}})   =
+one(::Type{StateSpace{Val{:siso},Val{:cont},M1,M2,M3,M4}}) where {M1,M2,M3,M4}  =
   StateSpace(zeros(eltype(M1),0,0), zeros(eltype(M2),0,1), zeros(eltype(M3),1,0),
   one(eltype(M4)))
-one{M1,M2,M3,M4}(::Type{StateSpace{Val{:siso},Val{:disc},M1,M2,M3,M4}})   =
+one(::Type{StateSpace{Val{:siso},Val{:disc},M1,M2,M3,M4}}) where {M1,M2,M3,M4}  =
   StateSpace(zeros(eltype(M1),0,0), zeros(eltype(M2),0,1), zeros(eltype(M3),1,0),
   one(eltype(M4)), zero(Float64))
-zero{M1,M2,M3,M4}(::Type{StateSpace{Val{:siso},Val{:cont},M1,M2,M3,M4}})  =
+zero(::Type{StateSpace{Val{:siso},Val{:cont},M1,M2,M3,M4}}) where {M1,M2,M3,M4} =
   StateSpace(zeros(eltype(M1),0,0), zeros(eltype(M2),0,1), zeros(eltype(M3),1,0),
   zero(eltype(M4)))
-zero{M1,M2,M3,M4}(::Type{StateSpace{Val{:siso},Val{:disc},M1,M2,M3,M4}})  =
+zero(::Type{StateSpace{Val{:siso},Val{:disc},M1,M2,M3,M4}}) where {M1,M2,M3,M4} =
   StateSpace(zeros(eltype(M1),0,0), zeros(eltype(M2),0,1), zeros(eltype(M3),1,0),
   zero(eltype(M4)), zero(Float64))
 
@@ -324,21 +314,21 @@ function zeros(s::StateSpace)
   if nr == 0
     return Complex{Float64}[]
   end
-  Arc, Brc, Crc, Drc, mrc, nrc, prc = reduce(Ar.', Cr.', Br.', Dr.')
+  Arc, Brc, Crc, Drc, mrc, nrc, prc = reduce(transpose(Ar), transpose(Cr), transpose(Br), transpose(Dr))
   if nrc == 0
     return Complex{Float64}[]
   end
 
-  svdobj  = svdfact([Crc Drc], thin = false)
-  W       = flipdim(svdobj.Vt', 2)
+  svdobj  = svd([Crc Drc], full = true)
+  W       = reverse(adjoint(svdobj.Vt), dims=2)
   Af      = [Arc Brc]*W[:, 1:nrc]
 
   if mrc == 0
-    zerovalues = eigfact(Af).values
+    zerovalues = eigen(Af).values
     return zerovalues
   else
     Bf    = W[1:nrc,1:nrc]
-    zerovalues = eigfact(Af, Bf).values
+    zerovalues = eigen(Af, Bf).values
     return zerovalues
   end
 end
@@ -349,7 +339,7 @@ tzeros(s::StateSpace) = zeros(minreal(s))
 # Poles of a state-space model
 function poles(s::StateSpace)
   Aₘ, _, _, _ = minreal(s.A, s.B, s.C, s.D)
-  return eigfact(Aₘ).values
+  return eigen!(Aₘ).values
 end
 
 # Negative of a state-space model
@@ -359,8 +349,8 @@ end
 -(s::StateSpace{Val{:mimo},Val{:disc}}) = StateSpace(s.A, s.B, -s.C, -s.D, s.Ts)
 
 # Addition
-function _ssparallel{T1,T2,S}(s1::StateSpace{Val{T1},Val{S}},
-  s2::StateSpace{Val{T2},Val{S}})
+function _ssparallel(s1::StateSpace{Val{T1},Val{S}},
+  s2::StateSpace{Val{T2},Val{S}}) where {T1,T2,S}
   if s1.Ts ≉ s2.Ts && s1.Ts ≠ zero(s1.Ts) && s2.Ts ≠ zero(s2.Ts)
     warn("parallel(s1,s2): Sampling time mismatch")
     throw(DomainError())
@@ -393,46 +383,36 @@ function +(s1::StateSpace{Val{:siso},Val{:disc}},
   StateSpace(a, b, c, d[1], Ts)
 end
 
-function +{T1,T2}(s1::StateSpace{Val{T1},Val{:cont}},
-  s2::StateSpace{Val{T2},Val{:cont}})
+function +(s1::StateSpace{Val{T1},Val{:cont}},
+  s2::StateSpace{Val{T2},Val{:cont}}) where {T1,T2}
   a, b, c, d, _ = _ssparallel(s1, s2)
   StateSpace(a, b, c, d)
 end
 
-function +{T1,T2}(s1::StateSpace{Val{T1},Val{:disc}},
-  s2::StateSpace{Val{T2},Val{:disc}})
+function +(s1::StateSpace{Val{T1},Val{:disc}},
+  s2::StateSpace{Val{T2},Val{:disc}}) where {T1,T2}
   a, b, c, d, Ts = _ssparallel(s1, s2)
   StateSpace(a, b, c, d, Ts)
 end
 
-.+(s1::StateSpace{Val{:siso}}, s2::StateSpace{Val{:siso}}) = +(s1, s2)
-
-+{T}(s::StateSpace{Val{T},Val{:disc}}, g::Union{Real,AbstractMatrix}) =
++(s::StateSpace{Val{T},Val{:disc}}, g::Union{Real,AbstractMatrix}) where {T} =
   +(s, ss(g, samplingtime(s)))
-+{T}(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:disc}}) =
++(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:disc}}) where {T} =
   +(ss(g, samplingtime(s)), s)
-+{T}(s::StateSpace{Val{T},Val{:cont}}, g::Union{Real,AbstractMatrix}) =
++(s::StateSpace{Val{T},Val{:cont}}, g::Union{Real,AbstractMatrix}) where {T} =
   +(s, ss(g))
-+{T}(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:cont}}) =
++(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:cont}}) where {T} =
   +(ss(g), s)
 
-.+(s::StateSpace{Val{:siso}}, g::Real)    = +(s, g)
-.+(g::Real, s::StateSpace{Val{:siso}})    = +(g, s)
-
 # Subtraction
 -(s1::StateSpace, s2::StateSpace) = +(s1, -s2)
 
-.-(s1::StateSpace{Val{:siso}}, s2::StateSpace{Val{:siso}}) = -(s1, s2)
-
 -(s::StateSpace, g::Union{Real,AbstractMatrix}) = +(s, -g)
 -(g::Union{Real,AbstractMatrix}, s::StateSpace) = +(g, -s)
 
-.-(s::StateSpace{Val{:siso}}, g::Real)    = -(s, g)
-.-(g::Real, s::StateSpace{Val{:siso}})    = -(g, s)
-
 # Multiplication
-function _ssseries{T1,T2,S}(s1::StateSpace{Val{T1},Val{S}},
-  s2::StateSpace{Val{T2},Val{S}})
+function _ssseries(s1::StateSpace{Val{T1},Val{S}},
+  s2::StateSpace{Val{T2},Val{S}}) where {T1,T2,S}
   # Remark: s1*s2 implies u -> s2 -> s1 -> y
 
   if s1.Ts ≉ s2.Ts && s1.Ts ≠ zero(s1.Ts) && s2.Ts ≠ zero(s2.Ts)
@@ -468,41 +448,31 @@ function *(s1::StateSpace{Val{:siso},Val{:disc}},
   StateSpace(a, b, c, d[1], Ts)
 end
 
-function *{T1,T2}(s1::StateSpace{Val{T1},Val{:cont}},
-  s2::StateSpace{Val{T2},Val{:cont}})
+function *(s1::StateSpace{Val{T1},Val{:cont}},
+  s2::StateSpace{Val{T2},Val{:cont}}) where {T1,T2}
   a, b, c, d, _ = _ssseries(s1, s2)
   StateSpace(a, b, c, d)
 end
 
-function *{T1,T2}(s1::StateSpace{Val{T1},Val{:disc}},
-  s2::StateSpace{Val{T2},Val{:disc}})
+function *(s1::StateSpace{Val{T1},Val{:disc}},
+  s2::StateSpace{Val{T2},Val{:disc}}) where {T1,T2}
   a, b, c, d, Ts = _ssseries(s1, s2)
   StateSpace(a, b, c, d, Ts)
 end
 
-.*(s1::StateSpace{Val{:siso}}, s2::StateSpace{Val{:siso}}) = *(s1, s2)
-
-*{T}(s::StateSpace{Val{T},Val{:disc}}, g::Union{Real,AbstractMatrix}) =
+*(s::StateSpace{Val{T},Val{:disc}}, g::Union{Real,AbstractMatrix}) where {T} =
   *(s, ss(g, samplingtime(s)))
-*{T}(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:disc}}) =
+*(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:disc}}) where {T} =
   *(ss(g, samplingtime(s)), s)
-*{T}(s::StateSpace{Val{T},Val{:cont}}, g::Union{Real,AbstractMatrix}) =
+*(s::StateSpace{Val{T},Val{:cont}}, g::Union{Real,AbstractMatrix}) where {T} =
   *(s, ss(g))
-*{T}(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:cont}}) =
+*(g::Union{Real,AbstractMatrix}, s::StateSpace{Val{T},Val{:cont}}) where {T} =
   *(ss(g), s)
 
-.*(s::StateSpace{Val{:siso}}, g::Real)    = *(s, g)
-.*(g::Real, s::StateSpace{Val{:siso}})    = *(g, s)
-
 # Division
 /(s1::StateSpace, s2::StateSpace)         = *(s1, inv(s2))
 
-./(s1::StateSpace{Val{:siso}}, s2::StateSpace{Val{:siso}}) = /(s1, s2)
-
 /(s::StateSpace, g::Union{Real,AbstractMatrix}) =
   *(s, inv(g))
 /(g::Union{Real,AbstractMatrix}, s::StateSpace) =
   *(g, inv(s))
-
-./(s::StateSpace{Val{:siso}}, g::Real)    = /(s, g)
-./(g::Real, s::StateSpace{Val{:siso}})    = /(g, s)
diff --git a/src/types/system/transferfunction.jl b/src/types/system/transferfunction.jl
index 9e64ab2..3591927 100644
--- a/src/types/system/transferfunction.jl
+++ b/src/types/system/transferfunction.jl
@@ -1,46 +1,45 @@
-immutable TransferFunction{T,S,U,V} <: LtiSystem{T,S}
+struct TransferFunction{T,S,U,V} <: LtiSystem{T,S}
   mat::V
   nu::Int
   ny::Int
   Ts::Float64
 
   # Continuous-time, single-input-single-output rational transfer function model
-  function (::Type{TransferFunction}){S,U<:Real,V<:Real}(r::RationalFunction{
-    Val{:s},Val{S},U,V})
+  function (::Type{TransferFunction})(r::RationalFunction{Val{:s},Val{S},U,V}) where {S,U<:Real,V<:Real}
     mat     = fill(r, 1, 1)
     ny, nu  = _tfcheck(mat)
     new{Val{:siso},Val{:cont},Val{S},typeof(mat)}(mat, nu, ny, zero(Float64))
   end
 
   # Discrete-time, single-input-single-output rational transfer function model
-  function (::Type{TransferFunction}){S,U<:Real,V<:Real}(r::RationalFunction{
-    Val{:z},Val{S},U,V}, Ts::Real)
+  function (::Type{TransferFunction})(r::RationalFunction{Val{:z},Val{S},U,V},
+    Ts::Real) where {S,U<:Real,V<:Real}
     mat     = fill(r, 1, 1)
     ny, nu  = _tfcheck(mat, Ts)
     new{Val{:siso},Val{:disc},Val{S},typeof(mat)}(mat, nu, ny, convert(Float64, Ts))
   end
 
   # Continuous-time, multi-input-multi-output rational transfer function model
-  function (::Type{TransferFunction}){S,U<:Real,V<:Real}(
-    mat::AbstractMatrix{RationalFunction{Val{:s},Val{S},U,V}})
+  function (::Type{TransferFunction})(
+    mat::AbstractMatrix{RationalFunction{Val{:s},Val{S},U,V}}) where {S,U<:Real,V<:Real}
     ny, nu  = _tfcheck(mat)
     new{Val{:mimo},Val{:cont},Val{S},typeof(mat)}(mat, nu, ny, zero(Float64))
   end
-  function (::Type{TransferFunction}){T<:Real,S}(mat::AbstractMatrix{T},
-    t::Type{Val{S}} = Val{:notc})
+  function (::Type{TransferFunction})(mat::AbstractMatrix{T},
+    t::Type{Val{S}} = Val{:notc}) where {T<:Real,S}
     m = map(x->RationalFunction(x, :s, t), mat)
     ny, nu  = _tfcheck(m)
     new{Val{:mimo},Val{:cont},t,typeof(m)}(m, nu, ny, zero(Float64))
   end
 
   # Discrete-time, multi-input-multi-output rational transfer function model
-  function (::Type{TransferFunction}){S,U<:Real,V<:Real}(
-    mat::AbstractMatrix{RationalFunction{Val{:z},Val{S},U,V}}, Ts::Real)
+  function (::Type{TransferFunction})(
+    mat::AbstractMatrix{RationalFunction{Val{:z},Val{S},U,V}}, Ts::Real) where {S,U<:Real,V<:Real}
     ny, nu  = _tfcheck(mat, Ts)
     new{Val{:mimo},Val{:disc},Val{S},typeof(mat)}(mat, nu, ny, convert(Float64, Ts))
   end
-  function (::Type{TransferFunction}){T<:Real,S}(mat::AbstractMatrix{T},
-    Ts::Real, t::Type{Val{S}} = Val{:notc})
+  function (::Type{TransferFunction})(mat::AbstractMatrix{T},
+    Ts::Real, t::Type{Val{S}} = Val{:notc}) where {T<:Real,S}
     m = map(x->RationalFunction(x, :z, t), mat)
     ny, nu  = _tfcheck(m)
     new{Val{:mimo},Val{:disc},t,typeof(m)}(m, nu, ny, convert(Float64, Ts))
@@ -51,8 +50,7 @@ immutable TransferFunction{T,S,U,V} <: LtiSystem{T,S}
     ucalc = u(t, x.y)
     ucalc = isa(ucalc, Real) ? [ucalc] : ucalc
     if !isa(ucalc, AbstractVector) || length(ucalc) ≠ sys.nu || !(eltype(ucalc) <: Real)
-      warn("sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values")
-      throw(DomainError())
+      throw(DomainError((t,x,dx,u), "sys(t,x,dx,u): u(t,y) has to be an `AbstractVector` of length $(sys.nu), containing `Real` values"))
     end
 
     sidx = 0 # number of states visited so far
@@ -71,8 +69,8 @@ immutable TransferFunction{T,S,U,V} <: LtiSystem{T,S}
   function (sys::TransferFunction)(x::DiffEqBase.DEDataArray, dx::AbstractVector, u, i, j, idx)
     # @assert degree(nump) ≤ degree(denp)
     r     = sys.mat[i,j]
-    nump  = num(r)
-    denp  = den(r)
+    nump  = numerator(r)
+    denp  = denominator(r)
     nump /= denp[end]
     denp /= denp[end]
 
@@ -82,15 +80,15 @@ immutable TransferFunction{T,S,U,V} <: LtiSystem{T,S}
     nb    = db == da ? db-1          : db
 
     if da > 0
-      dx[idx+(da:-1:1)]     = -denp[0:end-1]*x.x[idx+1]
-      dx[idx+(da:-1:da-nb)] += nump[0:nb]*u[j]
+      dx[(idx+da):-1:(idx+1)]      = -denp[0:end-1]*x.x[idx+1]
+      dx[(idx+da):-1:(idx+da-nb)] +=  nump[0:nb]*u[j]
       for k = 2:da
         dx[idx+k-1] += x.x[idx+k]
       end
       out += x.x[idx+1]
       # take care of direct term
       if db == da
-        dx[idx+(da:-1:1)] += -denp[0:end-1]*nump[end]*u[j]
+        dx[(idx+da):-1:(idx+1)] += -denp[0:end-1]*nump[end]*u[j]
       end
     end
     out, da
@@ -99,16 +97,15 @@ immutable TransferFunction{T,S,U,V} <: LtiSystem{T,S}
   # Evaluate system models at given complex numbers
   (sys::TransferFunction{Val{:siso}})(x::Number)                                    =
     (f = sys.mat[1]; convert(Complex128, f(x)))
-  (sys::TransferFunction{Val{:siso}}){M<:Number}(X::AbstractArray{M})               =
+  (sys::TransferFunction{Val{:siso}})(X::AbstractArray{M}) where {M<:Number}        =
     (f = sys.mat[1]; Complex128[f(x) for x in X])
   (sys::TransferFunction{Val{:mimo}})(x::Number)                                    =
     Complex128[f(x) for f in sys.mat]
-  (sys::TransferFunction{Val{:mimo}}){M<:Number}(X::AbstractArray{M})               =
-    Complex128[f(x) for f in sys.mat, x in X]
-  function (sys::TransferFunction){T<:Real}(; ω::Union{T, AbstractArray{T}} = Float64[])
+  (sys::TransferFunction{Val{:mimo}})(X::AbstractArray{M}) where {M<:Number}        =
+     Complex128[f(x) for f in sys.mat, x in X]
+  function (sys::TransferFunction)(; ω::Union{T, AbstractArray{T}} = Float64[]) where {T<:Real}
     if isempty(ω)
-      warn("sys(): Provide an argument for the function call. Refer to `?freqresp`.")
-      throw(DomainError())
+      throw(DomainError("sys(): Provide an argument for the function call. Refer to `?freqresp`."))
     end
     freqresp(sys, ω)
   end
@@ -116,32 +113,27 @@ end
 
 # Warn the user in other type constructions
 function (::Type{TransferFunction})(r::RationalFunction)
-  warn("tf(r): `r` can only be a real-coefficient `RationalFunction` of variable `:s`")
-  throw(DomainError())
+  throw(DomainError("tf(r): `r` can only be a real-coefficient `RationalFunction` of variable `:s`"))
 end
 
 function (::Type{TransferFunction})(r::RationalFunction, Ts::Real)
-  warn("tf(r, Ts): `r` can only be a real-coefficient `RationalFunction` of variable `:z`")
-  throw(DomainError())
+  throw(DomainError("tf(r, Ts): `r` can only be a real-coefficient `RationalFunction` of variable `:z`"))
 end
 
 function (::Type{TransferFunction})(m::AbstractMatrix)
-  warn("tf(m): `m` can only be an `AbstractMatrix` of real-coefficient `RationalFunction` objects of variable `:s`")
-  throw(DomainError())
+  throw(DomainError("tf(m): `m` can only be an `AbstractMatrix` of real-coefficient `RationalFunction` objects of variable `:s`"))
 end
 
 function (::Type{TransferFunction})(m::AbstractMatrix, Ts::Real)
-  warn("tf(m, Ts): `m` can only be an `AbstractMatrix` of real-coefficient `RationalFunction` objects of variable `:z`")
-  throw(DomainError())
+  throw(DomainError("tf(m, Ts): `m` can only be an `AbstractMatrix` of real-coefficient `RationalFunction` objects of variable `:z`"))
 end
 
 # Enforce rational transfer function type invariance
-function _tfcheck{T,S,U<:Real,V<:Real}(mat::AbstractMatrix{RationalFunction{Val{T},
-  Val{S},U,V}}, Ts::Real = zero(Float64))
+function _tfcheck(mat::AbstractMatrix{RationalFunction{Val{T},
+  Val{S},U,V}}, Ts::Real = zero(Float64)) where {T,S,U<:Real,V<:Real}
   # Check sampling time
   if Ts < zero(Ts) || isinf(Ts)
-    warn("tf(m, Ts): `Ts` must be a non-negative real number")
-    throw(DomainError())
+    throw(DomainError((m, Ts), "tf(m, Ts): `Ts` must be a non-negative real number"))
   end
 
   # Check input-output dimensions
@@ -198,8 +190,8 @@ function tf(num::AbstractVector, den::AbstractVector, Ts::Real, var::Symbol)
     throw(DomainError())
   end
 
-  numlast         = findlast(num)
-  denlast         = findlast(den)
+  numlast         = findlast(!iszero, num)
+  denlast         = findlast(!iszero, den)
   order           = max(numlast, denlast)
   num_            = zeros(eltype(num), order)
   num_[1:numlast] = num[1:numlast]
@@ -212,19 +204,19 @@ function tf(num::AbstractVector, den::AbstractVector, Ts::Real, var::Symbol)
 end
 
 # Continuous-time, multi-input-multi-output rational transfer function model
-tf{S,U,V}(mat::AbstractMatrix{RationalFunction{Val{:s},Val{S},U,V}})            =
+tf(mat::AbstractMatrix{RationalFunction{Val{:s},Val{S},U,V}}) where {S,U,V} =
   TransferFunction(mat)
 
 # Discrete-time, multi-input-multi-output rational transfer function model
-tf{S,U,V}(mat::AbstractMatrix{RationalFunction{Val{:z},Val{S},U,V}}, Ts::Real)  =
+tf(mat::AbstractMatrix{RationalFunction{Val{:z},Val{S},U,V}}, Ts::Real) where {S,U,V} =
   TransferFunction(mat, Ts)
 
 # Continuous-time, multi-input-multi-output rational transfer function model
-tf{T<:Real,S}(mat::AbstractMatrix{T}, t::Type{Val{S}} = Val{:notc})             =
+tf(mat::AbstractMatrix{T}, t::Type{Val{S}} = Val{:notc}) where {T<:Real,S} =
   TransferFunction(mat, t)
 
 # Discrete-time, multi-input-multi-output rational transfer function model
-tf{T<:Real,S}(mat::AbstractMatrix{T}, Ts::Real, t::Type{Val{S}} = Val{:notc})   =
+tf(mat::AbstractMatrix{T}, Ts::Real, t::Type{Val{S}} = Val{:notc}) where {T<:Real,S} =
   TransferFunction(mat, Ts, t)
 
 # Interfaces
@@ -233,17 +225,16 @@ numstates(s::TransferFunction)    = sum(x->degree(x)[2], s.mat)
 numinputs(s::TransferFunction)    = s.nu
 numoutputs(s::TransferFunction)   = s.ny
 
-num(s::TransferFunction{Val{:siso}})  = num(s.mat[1])
-num(s::TransferFunction{Val{:mimo}})  = map(num, s.mat)
-den(s::TransferFunction{Val{:siso}})  = den(s.mat[1])
-den(s::TransferFunction{Val{:mimo}})  = map(den, s.mat)
+numerator(s::TransferFunction{Val{:siso}})    = numerator(s.mat[1])
+numerator(s::TransferFunction{Val{:mimo}})    = map(numerator, s.mat)
+denominator(s::TransferFunction{Val{:siso}})  = denominator(s.mat[1])
+denominator(s::TransferFunction{Val{:mimo}})  = map(denominator, s.mat)
 
 # Iteration interface
-start(s::TransferFunction{Val{:mimo}})        = start(s.mat)
-next(s::TransferFunction{Val{:mimo}}, state)  = (s[state], state+1)
-done(s::TransferFunction{Val{:mimo}}, state)  = done(s.mat, state)
+iterate(s::TransferFunction{Val{:mimo}})        = iterate(s.mat)
+iterate(s::TransferFunction{Val{:mimo}}, state) = iterate(smat, state)
 
-eltype{S,U,V}(::Type{TransferFunction{Val{:mimo},Val{S},Val{U},V}}) =
+eltype(::Type{TransferFunction{Val{:mimo},Val{S},Val{U},V}}) where {S,U,V} =
   TransferFunction{Val{:siso},Val{S},Val{U},Matrix{eltype(V)}}
 
 length(s::TransferFunction{Val{:mimo}}) = length(s.mat)
@@ -297,35 +288,36 @@ getindex(s::TransferFunction{Val{:mimo}}, rows, ::Colon)    = s[rows, 1:s.nu]
 getindex(s::TransferFunction{Val{:mimo}}, ::Colon, cols)    = s[1:s.ny, cols]
 getindex(s::TransferFunction{Val{:mimo}}, ::Colon)          = s[1:end]
 getindex(s::TransferFunction{Val{:mimo}}, ::Colon, ::Colon) = s[1:s.ny,1:s.nu]
-endof(s::TransferFunction{Val{:mimo}})                      = endof(s.mat)
+firstindex(s::TransferFunction{Val{:mimo}})                 = firstindex(s.mat)
+lastindex(s::TransferFunction{Val{:mimo}})                  = lastindex(s.mat)
 
 # Conversion and promotion
-promote_rule{T1<:Real,T2,S,U,V}(::Type{T1}, ::Type{TransferFunction{Val{T2},Val{S},Val{U},V}}) =
+promote_rule(::Type{T1}, ::Type{TransferFunction{Val{T2},Val{S},Val{U},V}}) where {T1<:Real,T2,S,U,V} =
   TransferFunction{Val{T2},Val{S},Val{U}}
-promote_rule{T<:AbstractMatrix,S,U,V}(::Type{T}, ::Type{TransferFunction{Val{:mimo},Val{S},Val{U},V}}) =
+promote_rule(::Type{T}, ::Type{TransferFunction{Val{:mimo},Val{S},Val{U},V}}) where {T<:AbstractMatrix,S,U,V} =
   TransferFunction{Val{:mimo},Val{S},Val{U}}
 
-convert{U,V}(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}}, g::Real) =
+convert(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}}, g::Real) where {U,V} =
   tf(RationalFunction(g, :s, Val{U}))
-convert{U,V}(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}}, g::Real) =
+convert(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}}, g::Real) where {U,V} =
   tf(RationalFunction(g, :z, Val{U}), zero(Float64))
-convert{U,V}(::Type{TransferFunction{Val{:mimo},Val{:cont},Val{U},V}}, g::Real) =
+convert(::Type{TransferFunction{Val{:mimo},Val{:cont},Val{U},V}}, g::Real) where {U,V} =
   tf(fill(g,1,1), Val{U})
-convert{U,V}(::Type{TransferFunction{Val{:mimo},Val{:disc},Val{U},V}}, g::Real) =
+convert(::Type{TransferFunction{Val{:mimo},Val{:disc},Val{U},V}}, g::Real) where {U,V} =
   tf(fill(g,1,1), zero(Float64), Val{U})
-convert{U,V}(::Type{TransferFunction{Val{:mimo},Val{:cont},Val{U},V}}, g::AbstractMatrix) =
+convert(::Type{TransferFunction{Val{:mimo},Val{:cont},Val{U},V}}, g::AbstractMatrix) where {U,V} =
   tf(g)
-convert{U,V}(::Type{TransferFunction{Val{:mimo},Val{:disc},Val{U},V}}, g::AbstractMatrix) =
+convert(::Type{TransferFunction{Val{:mimo},Val{:disc},Val{U},V}}, g::AbstractMatrix) where {U,V} =
   tf(g, zero(Float64))
 
 # Multiplicative and additive identities (meaningful only for SISO)
-one{U,V}(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}})  =
+one(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}}) where {U,V}   =
   tf(one(eltype(V)))
-one{U,V}(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}})  =
+one(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}}) where {U,V}   =
   tf(one(eltype(V)), zero(Float64))
-zero{U,V}(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}}) =
+zero(::Type{TransferFunction{Val{:siso},Val{:cont},Val{U},V}}) where {U,V}  =
   tf(zero(eltype(V)))
-zero{U,V}(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}}) =
+zero(::Type{TransferFunction{Val{:siso},Val{:disc},Val{U},V}}) where {U,V}  =
   tf(zero(eltype(V)), zero(Float64))
 
 one(s::TransferFunction{Val{:siso},Val{:cont}})   = one(typeof(s))
@@ -395,8 +387,8 @@ end
 -(s::TransferFunction{Val{:mimo},Val{:disc}}) = TransferFunction(-s.mat, s.Ts)
 
 # Addition
-function _tfparallel{T1,T2,S,U}(s1::TransferFunction{Val{T1},Val{S},Val{U}},
-  s2::TransferFunction{Val{T2},Val{S},Val{U}})
+function _tfparallel(s1::TransferFunction{Val{T1},Val{S},Val{U}},
+  s2::TransferFunction{Val{T2},Val{S},Val{U}}) where {T1,T2,S,U}
   if s1.Ts ≉ s2.Ts && s1.Ts ≠ zero(s1.Ts) && s2.Ts ≠ zero(s2.Ts)
     warn("parallel(s1,s2): Sampling time mismatch")
     throw(DomainError())
@@ -410,56 +402,43 @@ function _tfparallel{T1,T2,S,U}(s1::TransferFunction{Val{T1},Val{S},Val{U}},
   return s1.mat + s2.mat, max(s1.Ts, s2.Ts)
 end
 
-function +{U}(s1::TransferFunction{Val{:siso},Val{:cont},Val{U}},
-  s2::TransferFunction{Val{:siso},Val{:cont},Val{U}})
+function +(s1::TransferFunction{Val{:siso},Val{:cont},Val{U}},
+  s2::TransferFunction{Val{:siso},Val{:cont},Val{U}}) where {U}
   mat, _ = _tfparallel(s1, s2)
   TransferFunction(mat[1])
 end
 
-function +{U}(s1::TransferFunction{Val{:siso},Val{:disc},Val{U}},
-  s2::TransferFunction{Val{:siso},Val{:disc},Val{U}})
+function +(s1::TransferFunction{Val{:siso},Val{:disc},Val{U}},
+  s2::TransferFunction{Val{:siso},Val{:disc},Val{U}}) where {U}
   mat, Ts = _tfparallel(s1, s2)
   TransferFunction(mat[1], Ts)
 end
 
-function +{T1,T2,U}(s1::TransferFunction{Val{T1},Val{:cont},Val{U}},
-  s2::TransferFunction{Val{T2},Val{:cont},Val{U}})
+function +(s1::TransferFunction{Val{T1},Val{:cont},Val{U}},
+  s2::TransferFunction{Val{T2},Val{:cont},Val{U}}) where {T1,T2,U}
   mat, _ = _tfparallel(s1, s2)
   TransferFunction(mat)
 end
 
-function +{T1,T2,U}(s1::TransferFunction{Val{T1},Val{:disc},Val{U}},
-  s2::TransferFunction{Val{T2},Val{:disc},Val{U}})
+function +(s1::TransferFunction{Val{T1},Val{:disc},Val{U}},
+  s2::TransferFunction{Val{T2},Val{:disc},Val{U}}) where {T1,T2,U}
   mat, Ts = _tfparallel(s1, s2)
   TransferFunction(mat, Ts)
 end
 
-.+(s1::TransferFunction{Val{:siso}}, s2::TransferFunction{Val{:siso}}) = +(s1, s2)
-
-+{T,S,U}(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix})  =
++(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix}) where {T,S,U} =
   +(s, convert(typeof(s), g))
-+{T,S,U}(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}})  =
++(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}}) where {T,S,U} =
   +(convert(typeof(s), g), s)
 
-.+(s::TransferFunction{Val{:siso}}, g::Real)    = +(s, g)
-.+(g::Real, s::TransferFunction{Val{:siso}})    = +(g, s)
-
-# Subtraction
--(s1::TransferFunction, s2::TransferFunction) = +(s1, -s2)
-
-.-(s1::TransferFunction{Val{:siso}}, s2::TransferFunction{Val{:siso}}) = -(s1, s2)
-
--{T,S,U}(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix})  =
+-(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix}) where {T,S,U} =
   -(s, convert(typeof(s), g))
--{T,S,U}(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}})  =
+-(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}}) where {T,S,U} =
   -(convert(typeof(s), g), s)
 
-.-(s::TransferFunction{Val{:siso}}, g::Real)    = -(s, g)
-.-(g::Real, s::TransferFunction{Val{:siso}})    = -(g, s)
-
 # Multiplication
-function _tfseries{T1,T2,S,U}(s1::TransferFunction{Val{T1},Val{S},Val{U}},
-  s2::TransferFunction{Val{T2},Val{S},Val{U}})
+function _tfseries(s1::TransferFunction{Val{T1},Val{S},Val{U}},
+  s2::TransferFunction{Val{T2},Val{S},Val{U}}) where {T1,T2,S,U}
   # Remark: s1*s2 implies u -> s2 -> s1 -> y
 
   if s1.Ts ≉ s2.Ts && s1.Ts ≠ zero(s1.Ts) && s2.Ts == zero(s2.Ts)
@@ -482,49 +461,39 @@ function _tfseries{T1,T2,S,U}(s1::TransferFunction{Val{T1},Val{S},Val{U}},
   return temp, max(s1.Ts, s2.Ts)
 end
 
-function *{U}(s1::TransferFunction{Val{:siso},Val{:cont},Val{U}},
-  s2::TransferFunction{Val{:siso},Val{:cont},Val{U}})
+function *(s1::TransferFunction{Val{:siso},Val{:cont},Val{U}},
+  s2::TransferFunction{Val{:siso},Val{:cont},Val{U}}) where {U}
   mat, _ = _tfseries(s1, s2)
   TransferFunction(mat[1])
 end
 
-function *{U}(s1::TransferFunction{Val{:siso},Val{:disc},Val{U}},
-  s2::TransferFunction{Val{:siso},Val{:disc},Val{U}})
+function *(s1::TransferFunction{Val{:siso},Val{:disc},Val{U}},
+  s2::TransferFunction{Val{:siso},Val{:disc},Val{U}}) where {U}
   mat, Ts = _tfseries(s1, s2)
   TransferFunction(mat[1], Ts)
 end
 
-function *{T1,T2,U}(s1::TransferFunction{Val{T1},Val{:cont},Val{U}},
-  s2::TransferFunction{Val{T2},Val{:cont},Val{U}})
+function *(s1::TransferFunction{Val{T1},Val{:cont},Val{U}},
+  s2::TransferFunction{Val{T2},Val{:cont},Val{U}}) where {T1,T2,U}
   mat, _ = _tfseries(s1, s2)
   TransferFunction(mat)
 end
 
-function *{T1,T2,U}(s1::TransferFunction{Val{T1},Val{:disc},Val{U}},
-  s2::TransferFunction{Val{T2},Val{:disc},Val{U}})
+function *(s1::TransferFunction{Val{T1},Val{:disc},Val{U}},
+  s2::TransferFunction{Val{T2},Val{:disc},Val{U}}) where {T1,T2,U}
   mat, Ts = _tfseries(s1, s2)
   TransferFunction(mat, Ts)
 end
 
-.*(s1::TransferFunction{Val{:siso}}, s2::TransferFunction{Val{:siso}}) = *(s1, s2)
-
-*{T,S,U}(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix})  =
+*(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix}) where {T,S,U} =
   *(s, convert(typeof(s), g))
-*{T,S,U}(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}})  =
+*(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}}) where {T,S,U} =
   *(convert(typeof(s), g), s)
 
-.*(s::TransferFunction{Val{:siso}}, g::Real)    = *(s, g)
-.*(g::Real, s::TransferFunction{Val{:siso}})    = *(g, s)
-
 # Division
 /(s1::TransferFunction, s2::TransferFunction)         = *(s1, inv(s2))
 
-./(s1::TransferFunction{Val{:siso}}, s2::TransferFunction{Val{:siso}}) = /(s1, s2)
-
-/{T,S,U}(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix})  =
+/(s::TransferFunction{Val{T},Val{S},Val{U}}, g::Union{Real,AbstractMatrix}) where {T,S,U} =
   /(s, convert(typeof(s), g))
-/{T,S,U}(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}})  =
+/(g::Union{Real,AbstractMatrix}, s::TransferFunction{Val{T},Val{S},Val{U}}) where {T,S,U} =
   /(convert(typeof(s), g), s)
-
-./(s::TransferFunction{Val{:siso}}, g::Real)    = /(s, g)
-./(g::Real, s::TransferFunction{Val{:siso}})    = /(g, s)
diff --git a/test/conversions/tf2ss.jl b/test/conversions/tf2ss.jl
index f9e2ac4..b3f0144 100644
--- a/test/conversions/tf2ss.jl
+++ b/test/conversions/tf2ss.jl
@@ -39,7 +39,7 @@ r34 = RationalFunction(2*Poly([34., 27., 5.], :s), poly([-1., -3., -5.], :s))
 Ptf = tf([r11 r12 r13 r14; r21 r22 r23 r24; r31 r32 r33 r34])
 
 Plfd = lfd(LTISystems._tf2lfd(Ptf)...)
-# ss(Plfd) # TODO find bug in _tf2lfd
+#ss(Plfd) # TODO find bug in _tf2lfd
 
 Prfd = rfd(LTISystems._tf2rfd(Ptf)...)
 ss(Prfd)
diff --git a/test/methods/simulation.jl b/test/methods/simulation.jl
index a625d6d..2279e71 100644
--- a/test/methods/simulation.jl
+++ b/test/methods/simulation.jl
@@ -1,4 +1,4 @@
-type TestSystem{M1,M2}
+struct TestSystem{M1,M2}
   v::Vector{LTISystems.LtiSystem}
   sol::M1
   input::M2
@@ -7,9 +7,8 @@ type TestSystem{M1,M2}
   end
 end
 
-Base.start(s::TestSystem)       = start(s.v)
-Base.next(s::TestSystem, state) = (s.v[state], state+1)
-Base.done(s::TestSystem, state) = done(s.v, state)
+Base.iterate(s::TestSystem)         = (s.v[1], 1)
+Base.iterate(s::TestSystem, i::Int) = i < lastindex(s.v) ? (s.v[i+1], i+1) : nothing
 
 T = 20.
 time = (0, T)
@@ -18,14 +17,14 @@ simtol = 1e-4
 function testsim(sim, sol, simtol)
   simsum = zero(sim.y[1])
   for (i,t) in enumerate(sim.t)
-    simsum += sum(abs2.(sol(t)-sim.y[i,:]))
+    simsum += sum(abs2.(broadcast(-, sol(t), sim.y[i,:])))
   end
   return simsum < simtol
 end
 
 # first continous system
 # s1 = 1/(s+2) step response 1/2 - 1/2*exp(-2t)
-c1v = Vector{LTISystems.LtiSystem}(0)
+c1v = Vector{LTISystems.LtiSystem}()
 push!(c1v, tf([1],[1, 2]))
 push!(c1v, lfd([1], [1, 2]))
 push!(c1v, rfd([1], [1, 2]))
@@ -117,7 +116,7 @@ c4rfd = rfd([1.,6.,9], [1.,3.,2.])
 A = [ 0.  1.;
      -2. -3.]
 
-B = [0. 1.].'
+B = transpose([0. 1.])
 
 C = [7. 3.]
 D = 1.
@@ -132,28 +131,28 @@ c4 = TestSystem(c4v, c4sol, c4inp)
 @testset "continuous simulation" begin
   @testset "first order SISO system" begin
     for sys in c1
-      sim = simulate(sys, time; input = c1.input, initial = zeros(numstates(sys)))
+      sim = simulate(sys, time; input = c1.input, initial = zeros(Float64, numstates(sys)))
       @test testsim(sim, c1.sol, simtol)
     end
   end
 
   @testset "second order MIMO system" begin
     for sys in c2
-      sim = simulate(sys, time; input = c2.input, initial = zeros(numstates(sys)))
+      sim = simulate(sys, time; input = c2.input, initial = zeros(Float64, numstates(sys)))
       @test testsim(sim, c2.sol, simtol)
     end
   end
 
   @testset "MIMO system with different row and col degrees" begin
     for sys in c3
-      sim = simulate(sys, time; input = c3.input, initial = zeros(numstates(sys)))
+      sim = simulate(sys, time; input = c3.input, initial = zeros(Float64, numstates(sys)))
       @test testsim(sim, c3.sol, simtol)
     end
   end
 
   @testset "second order SISO system with direct term" begin
     for sys in c4
-      sim = simulate(sys, time; input = c4.input, initial = zeros(numstates(sys)))
+      sim = simulate(sys, time; input = c4.input, initial = zeros(Float64, numstates(sys)))
       @test testsim(sim, c4.sol, simtol)
     end
   end
@@ -171,7 +170,7 @@ a = 0.2
 d1sol = t->fores(t,b,a)
 d1inp = Signals.Step([0.], [1.], [0.])
 
-d1v = Vector{LTISystems.LtiSystem}(0)
+d1v = Vector{LTISystems.LtiSystem}()
 push!(d1v, tf([b],[1, -a], 1))
 push!(d1v, lfd([b],[1, -a], 1))
 push!(d1v, rfd([b],[1, -a], 1))
diff --git a/test/runtests.jl b/test/runtests.jl
index 8379faa..23a9621 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,12 @@
 using Polynomials
 using RationalFunctions
 using PolynomialMatrices
+using LinearAlgebra
 
 using LTISystems
 
-using Base.Test
-using Base.Test.@testset
+using Test
+#using Test.@testset
 
 # conversions
 include("conversions/tf2ss.jl")
diff --git a/test/types/system/ltisystem.jl b/test/types/system/ltisystem.jl
index 6fb9e9a..d4d22ad 100644
--- a/test/types/system/ltisystem.jl
+++ b/test/types/system/ltisystem.jl
@@ -1,5 +1,5 @@
-info("Starting tests for `LtiSystem`...")
+print("Starting tests for `LtiSystem`...")
 
 @test 1 == 1
 
-info("Tests for `LtiSystem` finished.")
+print("Tests for `LtiSystem` finished.")
diff --git a/test/types/system/mfd.jl b/test/types/system/mfd.jl
index fb2b7f9..9619dc4 100644
--- a/test/types/system/mfd.jl
+++ b/test/types/system/mfd.jl
@@ -1,6 +1,4 @@
-## MatrixFractionDescriptions
-using Polynomials
-using PolynomialMatrices
+print("Starting tests for `MatrixFractionDescriptions`...")
 
 # Defining left MatrixFractionDescriptions
 s = Poly([0, 1],:s)
@@ -55,3 +53,5 @@ sys = ss(A,B,C,D)
 # ml.D
 #
 # @test isapprox(mr, ml; rtol=10)
+
+print("Tests for `MatrixFractionDescriptions` finished.")
diff --git a/test/types/system/statespace.jl b/test/types/system/statespace.jl
index 0bbf8ae..709f60d 100644
--- a/test/types/system/statespace.jl
+++ b/test/types/system/statespace.jl
@@ -1,5 +1,5 @@
-info("Starting tests for `StateSpace`...")
+print("Starting tests for `StateSpace`...")
 
 @test 1 == 1
 
-info("Tests for `StateSpace` finished.")
+print("Tests for `StateSpace` finished.")
diff --git a/test/types/system/transferfunction.jl b/test/types/system/transferfunction.jl
index 3dc3290..e5952dc 100644
--- a/test/types/system/transferfunction.jl
+++ b/test/types/system/transferfunction.jl
@@ -1,4 +1,4 @@
-info("Starting tests for `TransferFunction`...")
+print("Starting tests for `TransferFunction`...")
 
 # Rational transfer function constructions for SISO systems
 rc1   = RationalFunction([1], [2, 1], :s)     # 1/(s+2)
@@ -15,10 +15,10 @@ sysds4= tf(Poly(1, :z), [1, -0.5], 0.5)
 sysds5= tf(1, Poly([-0.5, 1], :z), 0.5)
 sysds6= tf([0, 1], [1, -0.5], 0.5, :z̄)
 
-@test num(syscs1) == num(syscs2) == num(syscs3) == num(syscs4) == num(syscs5)
-@test den(syscs1) == den(syscs2) == den(syscs3) == den(syscs4) == den(syscs5)
-@test num(sysds1) == num(sysds2) == num(sysds3) == num(sysds4) == num(sysds5) == num(sysds6)
-@test den(sysds1) == den(sysds2) == den(sysds3) == den(sysds4) == den(sysds5) == den(sysds6)
+@test numerator(syscs1) == numerator(syscs2) == numerator(syscs3) == numerator(syscs4) == numerator(syscs5)
+@test denominator(syscs1) == denominator(syscs2) == denominator(syscs3) == denominator(syscs4) == denominator(syscs5)
+@test numerator(sysds1) == numerator(sysds2) == numerator(sysds3) == numerator(sysds4) == numerator(sysds5) == numerator(sysds6)
+@test denominator(sysds1) == denominator(sysds2) == denominator(sysds3) == denominator(sysds4) == denominator(sysds5) == denominator(sysds6)
 
 # Rational transfer function constructions for MIMO systems
 rc2   = RationalFunction([1,1], [3,1], :s)
@@ -28,16 +28,16 @@ rd2   = RationalFunction([-0.1,1], [-0.3,1], :z)
 matd  = diagm([rd1, rd2])
 sysdm = tf(matd, 0.5)
 
-@test num(syscm) == [num(rc1) zero(num(rc1)); zero(num(rc1)) num(rc2)]
-@test den(syscm) == [den(rc1) one(num(rc1)); one(num(rc1)) den(rc2)]
-@test num(sysdm) == [num(rd1) zero(num(rd1)); zero(num(rd1)) num(rd2)]
-@test den(sysdm) == [den(rd1) one(num(rd1)); one(num(rd1)) den(rd2)]
+@test numerator(syscm) == [numerator(rc1) zero(numerator(rc1)); zero(numerator(rc1)) numerator(rc2)]
+@test denominator(syscm) == [denominator(rc1) one(numerator(rc1)); one(numerator(rc1)) denominator(rc2)]
+@test numerator(sysdm) == [numerator(rd1) zero(numerator(rd1)); zero(numerator(rd1)) numerator(rd2)]
+@test denominator(sysdm) == [denominator(rd1) one(numerator(rd1)); one(numerator(rd1)) denominator(rd2)]
 
-@test num(tf(diagm([1,2]))) == [Poly(1,:s) Poly(0,:s); Poly(0,:s) Poly(2,:s)]
-@test den(tf(diagm([1,2]))) == [Poly(1,:s) Poly(1,:s); Poly(1,:s) Poly(1,:s)]
+@test numerator(tf(diagm([1,2]))) == [Poly(1,:s) Poly(0,:s); Poly(0,:s) Poly(2,:s)]
+@test denominator(tf(diagm([1,2]))) == [Poly(1,:s) Poly(1,:s); Poly(1,:s) Poly(1,:s)]
 
-@test num(tf(diagm([1,2]), 0.5)) == [Poly(1,:z) Poly(0,:z); Poly(0,:z) Poly(2,:z)]
-@test den(tf(diagm([1,2]), 0.5)) == [Poly(1,:z) Poly(1,:z); Poly(1,:z) Poly(1,:z)]
+@test numerator(tf(diagm([1,2]), 0.5)) == [Poly(1,:z) Poly(0,:z); Poly(0,:z) Poly(2,:z)]
+@test denominator(tf(diagm([1,2]), 0.5)) == [Poly(1,:z) Poly(1,:z); Poly(1,:z) Poly(1,:z)]
 
 # Sampling time
 @test samplingtime(syscs1) == samplingtime(syscm) == zero(Float64)
@@ -48,4 +48,4 @@ sysdm = tf(matd, 0.5)
 @test numinputs(syscm)  == numinputs(sysdm)   == 2
 @test numoutputs(syscm) == numoutputs(sysdm)  == 2
 
-info("Tests for `TransferFunction` finished.")
+print("Tests for `TransferFunction` finished.")