From 882e89cb089e7dbfab76a7a4e4baed4ac22131f5 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Wed, 4 Sep 2024 07:43:28 +0200
Subject: [PATCH] add moving-average filter

---
 docs/src/tutorials/noise.md       |  2 +-
 src/ModelingToolkitSampledData.jl |  3 ++-
 src/discrete_blocks.jl            | 30 ++++++++++++++++++++++++++++++
 test/test_discrete_blocks.jl      | 26 ++++++++++++++++++++++++++
 4 files changed, 59 insertions(+), 2 deletions(-)

diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md
index ed26b5a..f2733b5 100644
--- a/docs/src/tutorials/noise.md
+++ b/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.
diff --git a/src/ModelingToolkitSampledData.jl b/src/ModelingToolkitSampledData.jl
index c578b54..7edf2ed 100644
--- a/src/ModelingToolkitSampledData.jl
+++ b/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
diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl
index f6683ed..f546fe6 100644
--- a/src/discrete_blocks.jl
+++ b/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
diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl
index c01a9c0..cf78bb1 100644
--- a/test/test_discrete_blocks.jl
+++ b/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