add moving-average filter

Author: baggepinnen

docs/src/tutorials/noise.md

@@ -87,7 +87,8 @@ plot(figy, figu, plot_title = "DC Motor with Discrete-time Speed Controller")
 ## Noise filtering
 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.
+No discrete-time filter components are available yet. You may, e.g.
+- 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

@@ -8,8 +8,9 @@ export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold,
        ClockChanger, SampleWithADEffects,
        DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
        DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization,
-       DiscreteSlewRateLimiter, ExponentialFilter
+       DiscreteSlewRateLimiter, ExponentialFilter, MovingAverageFilter
 export DiscreteOnOffController
+       
 include("discrete_blocks.jl")
 
 

src/discrete_blocks.jl

@@ -1039,6 +1039,34 @@ Y(z) = \\dfrac{a}{1 - (1 - a) z^{-1}} U(z)
     end
 end
 
+"""
+    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
+
 """
     DiscreteOnOffController(b = 0.1, bool = true)
 
@@ -1132,4 +1160,4 @@ The operations occur in the order
         connect(noise.output, quantization.input)
         connect(quantization.output, output)
     end
-end
+end

test/test_discrete_blocks.jl

@@ -512,7 +512,6 @@ end
     @test 0.89 <= sol(4, idxs=m.x) <= 1.11
 end
 
-
 @testset "ExponentialFilter" begin
     @info "Testing ExponentialFilter"
     z = ShiftIndex(Clock(0.1))
@@ -540,6 +539,34 @@ 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
+
 @testset "sampling with AD effects" begin
     @info "Testing sampling with AD effects"
     z = ShiftIndex()
@@ -565,4 +592,4 @@ end
     @test length(unique(sol[m.sampling.y])) == 8
     @test maximum(abs, sol(0:0.03:1, idxs=m.sampling.y) - sol(0:0.03:1, idxs=m.sampling.u)) < 0.3
 
-end
+end