add moving-average filter

docs/src/tutorials/noise.md

@@ -86,7 +86,7 @@ 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.
+- Add moving average filtering using [`MovingAverageFilter`](@ref) according to `y(k) ~ 1/N * sum(u(k-i) for 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,
+       MovingAverageFilter
 include("discrete_blocks.jl")
 
 end

src/discrete_blocks.jl

@@ -940,3 +940,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)
+
+
+"""
+    MovingAverageFilter(N = 3)
+
+Exponential filtering with input-output relation ``y(z) ~ sum(u(z-i) for i in 0:N-1) / N``.
+
+Please note: this implementation of a moving average filter is not optimized for very large number of filter taps `N`.
+
+# Parameters:
+- `N`: (structural) Number of samples to average over
+
+# Variables:
+- `u`: Input signal
+- `y`: Output signal
+
+# Connectors:
+- `input::RealInput`: Input signal
+- `output::RealOutput`: Output signal
+"""
+@mtkmodel MovingAverageFilter begin
+    @extend u, y = siso = SISO()
+    @structural_parameters begin
+        z = ShiftIndex()
+        N = 3
+    end
+    @equations begin
+        y(z) ~ sum(u(z-i) for i in 0:N-1) / N
+    end
+end

test/test_discrete_blocks.jl

@@ -462,3 +462,29 @@ end
 end
 
 
+@testset "MovingAverageFilter" begin
+    @info "Testing MovingAverageFilter"
+    z = ShiftIndex(Clock(0.1))
+    @mtkmodel MovingAverageFilterModel begin
+        @components begin
+            input = Step(start_time=1, smooth=false)
+            filter = MovingAverageFilter(; N=3, 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 = MovingAverageFilterModel()
+    m = complete(m)
+    ssys = structural_simplify(IRSystem(m))
+    prob = ODEProblem(ssys, [m.filter.u(z-i) => 0 for i = 0:3], (0.0, 2.0))
+    sol = solve(prob, Tsit5(), dtmax=0.1)
+    # plot(sol, idxs=m.filter.y)
+    @test sol(1.5, idxs=m.filter.y) == 1
+    @test sol(0.999, idxs=m.filter.y) == 0
+    @test 0 < sol(1.1, idxs=m.filter.y) < 1
+end