Merge pull request #16 from JuliaComputing/expfilt

add exponential filter component

Author: baggepinnen

docs/src/examples/dc_motor_pi.md

@@ -104,7 +104,7 @@ Plots.plot(p1, p2, layout = (2, 1))
 
 ## Discrete-time controller
 
-Until now, we have modeled both the physical part of the system, the DC motor, and the control systems, in continuous time. In practice, it is common to implement control systems on a computer operating at a fixed sample rate, i.e, in discrete time. A system containing both continuous-time parts and discrete-time parts is often referred to as a "sampled-data system", and the ModelingToolkit standard library contains several components to model such systems.
+Until now, we have modeled both the physical part of the system, the DC motor, and the control systems, in continuous time. In practice, it is common to implement control systems on a computer operating at a fixed sample rate, i.e, in discrete time. A system containing both continuous-time parts and discrete-time parts is often referred to as a "sampled-data system", and the ModelingToolkit standard library together with ModelingToolkitSampledData.jl contain several components to model such systems.
 
 Below, we re-model the system, this time with a discrete-time controller: [`DiscretePIDStandard`](@ref). To interface between the continuous and discrete parts of the model, we make use of a [`Sampler`](@ref) and [`ZeroOrderHold`](@ref) blocks. Apart from the three aforementioned components, the model is the same as before.
 

docs/src/tutorials/noise.md

@@ -1,11 +1,12 @@
 # Measurement noise and corruption
-Measurement noise is practically always present in signals originating from real-world sensors. In a sampled-data system, analyzing the influence of measurement noise using simulation is relatively straight forward. Below, we add Gaussian white noise to the speed sensor signal in the DC motor example. The noise is added using the [`NormalNoise`](@ref) block.
+Measurement noise is practically always present in signals originating from real-world sensors. In a sampled-data system, analyzing the influence of measurement noise using simulation is relatively straight forward. Below, we add Gaussian white noise to the speed sensor signal in the DC-motor example from [DC Motor with PI-controller](@ref). The noise is added using the [`NormalNoise`](@ref) block.
 
 This block has two modes of operation
 1. If `additive = false` (default), the block has the connector `output` only, and this output is the noise signal.
 2. If `additive = true`, the block has the connectors `input` and `output`, and the output is the sum of the input and the noise signal, i.e., the noise is _added_ to the input signal. This mode makes it convenient to add noise to a signal in a sampled-data system.
 
 ## Example: Noise
+This example is a continuation of the DC-motor example from [DC Motor with PI-controller](@ref). We add Gaussian white noise with ``σ^2 = 0.1^2`` to the speed sensor signal.
 ```@example NOISE
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t
@@ -84,9 +85,9 @@ plot(figy, figu, plot_title = "DC Motor with Discrete-time Speed Controller")
 ```
 
 ## Noise filtering
-No discrete-time filter components are available yet. You may, e.g.
-- Add exponential filtering using `xf(k) ~ (1-α)xf(k-1) + α*x(k)`, where `α` is the filter coefficient and `x` is the signal to be filtered.
-- Add moving average filtering using `xf(k) ~ 1/N sum(i->x(k-i), i=0:N-1)`, where `N` is the number of samples to average over.
+You may, e.g.
+- Use [`ExponentialFilter`](@ref) to add exponential filtering using `y(k) ~ (1-a)y(k-1) + a*u(k)`, where `a` is the filter coefficient and `u` is the signal to be filtered.
+- Add moving average filtering using `y(k) ~ 1/N sum(i->u(k-i), i=0:N-1)`, where `N` is the number of samples to average over.
 
 ## Colored noise
 Colored noise can be achieved by filtering white noise through a filter with the desired spectrum. No components are available for this yet.

src/ModelingToolkitSampledData.jl

@@ -7,7 +7,8 @@ export get_clock
 export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler,
        ClockChanger,
        DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
-       DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization
+       DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization,
+       ExponentialFilter
 include("discrete_blocks.jl")
 
 end

src/discrete_blocks.jl

@@ -888,6 +888,7 @@ A quantization block that quantizes the input signal to a specified number of bi
 - `y_min`: Lower limit of output
 - `bits`: Number of bits of quantization
 - `quantized`: If quantization effects shall be computed. If false, the output is equal to the input, which may be useful for, e.g., linearization.
+- `midrise`: (structural) If true (default), the quantizer is a midrise quantizer, otherwise it is a midtread quantizer. See [Docs: Quantization](https://juliacomputing.github.io/ModelingToolkitSampledData.jl/dev/tutorials/noise/#Quantization) for more details.
 
 # Connectors:
 - `input`
@@ -940,3 +941,33 @@ function quantize(u, bits, y_min, y_max, midrise)
 end
 
 @register_symbolic quantize(u::Real, bits::Real, y_min::Real, y_max::Real, midrise::Bool)
+
+"""
+    ExponentialFilter(a = 0.1)
+
+Exponential filtering with input-output relation ``y(z) ~ (1 - a) y(z-1) + a u(z-1)``
+
+# Parameters:
+- `a`: Filter parameter `[0, 1]`, a small value implies stronger filtering. 
+
+# Variables:
+- `u`: Input signal
+- `y`: Output signal
+
+# Connectors:
+- `input::RealInput`: Input signal
+- `output::RealOutput`: Output signal
+"""
+@mtkmodel ExponentialFilter begin
+    @extend u, y = siso = SISO()
+    @structural_parameters begin
+        z = ShiftIndex()
+    end
+    @parameters begin
+        a = 0.1, [description = "Filter parameter"]
+    end
+    @equations begin
+        y(z) ~ (1 - a) * y(z-1) + a * u(z)
+    end
+end
+

test/test_discrete_blocks.jl

@@ -462,3 +462,29 @@ end
 end
 
 
+@testset "ExponentialFilter" begin
+    @info "Testing ExponentialFilter"
+    z = ShiftIndex(Clock(0.1))
+    @mtkmodel ExponentialFilterModel begin
+        @components begin
+            input = Step(start_time=1, smooth=false)
+            filter = ExponentialFilter(; a = 0.1, z)
+        end
+        @variables begin
+            x(t) = 0 # Dummy variable to workaround JSCompiler bug
+        end
+        @equations begin
+            connect(input.output, filter.input)
+            D(x) ~ 0
+        end
+    end
+    @named m = ExponentialFilterModel()
+    m = complete(m)
+    ssys = structural_simplify(IRSystem(m))
+    prob = ODEProblem(ssys, [m.filter.y(z-1) => 0], (0.0, 10.0))
+    sol = solve(prob, Tsit5(), dtmax=0.1)
+    @test sol(10, idxs=m.filter.y) ≈ 1 atol=0.001
+    @test sol(0.999, idxs=m.filter.y) == 0
+    @test sol(1.1, idxs=m.filter.y) > 0
+
+end