add Sampler with AD effects

Author: baggepinnen

docs/src/tutorials/noise.md

@@ -117,6 +117,7 @@ z = ShiftIndex(Clock(0.1))
     @equations begin
         connect(input.output, quant.input)
         D(x) ~ 0 # Dummy equation
+        # D(x) ~ Hold(quant.y) # Dummy equation
     end
 end
 @named m = QuantizationModel()
@@ -147,3 +148,32 @@ plot(u, [y_mr y_mt], label=["Midrise" "Midtread"], xlabel="Input", ylabel="Outpu
 Note how the default mid-rise quantizer mode has a rise at the middle of the interval, while the mid-tread mode has a flat region (a tread) centered around the middle of the interval.
 
 The default option `midrise = true` includes both end points as possible output values, while `midrise = false` does not include the upper limit.
+
+
+## Sampling with AD effects
+The block [`SampleWithADEffects`](@ref) combines a [`Sampler`](@ref), a [`NormalNoise](@ref) and a [`Quantization`](@ref) block to simulate the effects of practical sampling, noise and quantization in an AD converter. The block has the connectors `input` and `output`, where the input is the continuous-time signal to be sampled, and the output is the quantized, noisy signal. Example
+
+```@example QUANT
+@mtkmodel PracticalSamplerModel begin
+    @components begin
+        input = Sine(amplitude=1.2, frequency=1, smooth=false)
+        sampling = SampleWithADEffects(; dt=0.03, bits=3, y_min = -1, y_max = 1, sigma = 0.05, noisy = false, quantized=false, midrise=true)
+    end
+    @variables begin
+        x(t) = 0 # Dummy variable to work around a bug for models without continuous-time state
+    end
+    @equations begin
+        connect(input.output, sampling.input)
+        D(x) ~ Hold(sampling.y) # Dummy equation
+    end
+end
+@named m = PracticalSamplerModel()
+m = complete(m)
+ssys = structural_simplify(IRSystem(m))
+# prob = ODEProblem(ssys, [m.sampling.noise.y(z-1) => 0], (0.0, 2.0))
+prob = ODEProblem(ssys, [], (0.0, 2.0))
+sol = solve(prob, Tsit5())
+plot(sol, idxs=m.input.output.u)
+plot!(sol, idxs=m.sampling.y, label="AD converted output")
+# plot!(sol, idxs=m.sampling.sampler.y, label="AD converted output")
+```

src/ModelingToolkitSampledData.jl

@@ -5,7 +5,7 @@ using StableRNGs
 
 export get_clock
 export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler,
-       ClockChanger,
+       ClockChanger, SampleWithADEffects,
        DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
        DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization,
        ExponentialFilter

src/discrete_blocks.jl

@@ -913,8 +913,6 @@ This block is not differentiable, the derivative is zero everywhere exect for at
         bits::Int = 8, [description = "Number of bits of quantization"]
         quantized::Bool = true, [description = "If quantization effects shall be computed."]
     end
-    begin
-    end
     @equations begin
         y(z) ~ ifelse(quantized == true, quantize(u(z), bits, y_min, y_max, midrise), u(z))
     end
@@ -1013,4 +1011,58 @@ Discrete-time On-Off controller with hysteresis. The controller switches between
             y(z) ~ k*(2*s(z) - 1)
         end
     end
+end
+
+"""
+    SampleWithADEffects(quantized = false, noisy = false)
+
+A sampler with additional effects that appear in practical systems, such as measurement noise and quantization.
+
+The operations occur in the order
+1. Sampling
+2. Noise addition
+3. Quantization
+
+# Structural parameters:
+- `quantized`: If true, the output is quantized. When this option is used, the output is quantized to the number of bits specified by the `bits` parameter. The quantization is midrise if `midrise = true`, otherwise it is midtread. The output is also limited to the range `[y_min, y_max]`.
+- `noisy`: If true, the output is corrupted by additive white Gaussian noise with standard deviation `sigma` (defaults to 0.1). If `noisy = false`, the noise block is a unit gain.
+- `dt`: Sample interval of the sampler. If not specified, the sample interval is inferred from the clock of the system.
+- `clock`: Clock signal of the system. If not specified, the sample interval is inferred from the clock of the system. If `clock` is specified, the parameter `dt` has no effect.
+
+# Parameters:
+- `y_min`: Lower limit of output, defaults to -1. Only used if `quantized = true`.
+- `y_max`: Upper limit of output, defaults to 1. Only used if `quantized = true`.
+- `bits`: Number of bits of quantization, defaults to 8 (256 output levels between `y_min` and `y_max`). Only used if `quantized = true`.
+- `sigma`: Standard deviation of the additive Gaussian noise, defaults to 0.1. Only used if `noisy = true`.
+"""
+@mtkmodel SampleWithADEffects begin
+    @extend input, output = siso = SISO()
+    @structural_parameters begin
+        dt = nothing
+        clock = (dt === nothing ? InferredDiscrete() : Clock(dt))
+        midrise = true
+        quantized = false
+        noisy = false
+    end
+    @parameters begin
+        y_min = -1, [description = "Lower limit of output"]
+        y_max = 1, [description = "Upper limit of output"]
+        bits = 8, [description = "Number of bits of quantization"]
+        sigma = 0.1, [description = "Standard deviation of the additive noise"]
+    end
+    @components begin
+        sampler = Sampler(; clock)
+        if noisy
+            noise = NormalNoise(; sigma, additive = true)
+        else
+            noise = Gain(; k = 1)
+        end
+        quantization = Quantization(; y_min, y_max, midrise, quantized)
+    end
+    @equations begin
+        connect(input, sampler.input)
+        connect(sampler.output, noise.input)
+        connect(noise.output, quantization.input)
+        connect(quantization.output, output)
+    end
 end