diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index d2b115a..eb293ce 100644 --- a/docs/src/tutorials/noise.md +++ b/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,38 @@ 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 an ideal [`Sampler`](@ref), a [`NormalNoise](@ref) and a [`Quantization`](@ref) block to simulate the undesirable but practically occurring effects of 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.1, noisy = true, quantized=true, 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)) +sol = solve(prob, Tsit5()) +plot(sol, idxs=m.input.output.u) +plot!(sol, idxs=m.sampling.y, label="AD converted output") +``` + +Both quantization and noise addition are optional and turned off by default. In the example above, we turn them on with keywords `noisy = true` and `quantized = true`. The noise is Gaussian white noise with standard deviation `sigma`, and the quantization is a 3-bit midrise quantizer (8 output levels) with limits `y_min` and `y_max`. Limits have to be provided when quantization is used. The `dt` parameter is the sampling time, if left unspecified, it will be inferred from context. + +Things to notice in the plot: +- The sampled signal is saturated at the quantization limits ±1. +- The noise is added to the signal before quantization, which means that the sampled signal has ``2^\text{bits}`` distinct output levels only. +- 0 is not a possible output value. In situations where 0 is an important value (such as in the presence of integration of a quantized value that is expected to be close to 0), the mid-tread quantizer should be used instead by passing `midrise = false`. + diff --git a/src/ModelingToolkitSampledData.jl b/src/ModelingToolkitSampledData.jl index 6e6039c..4a6b4c0 100644 --- a/src/ModelingToolkitSampledData.jl +++ b/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 diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl index 197bba7..bf5172d 100644 --- a/src/discrete_blocks.jl +++ b/src/discrete_blocks.jl @@ -910,11 +910,9 @@ This block is not differentiable, the derivative is zero everywhere exect for at @parameters begin y_max = 1, [description = "Upper limit of output"] y_min = -1, [description = "Lower limit of output"] - bits::Int = 8, [description = "Number of bits of quantization"] + bits = 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(; bits, 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 \ No newline at end of file diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index c6c0dbd..c3f4446 100644 --- a/test/test_discrete_blocks.jl +++ b/test/test_discrete_blocks.jl @@ -514,4 +514,31 @@ end @test sol(0.999, idxs=m.filter.y) == 0 @test sol(1.1, idxs=m.filter.y) > 0 +end + +@testset "sampling with AD effects" begin + @info "Testing sampling with AD effects" + z = ShiftIndex() + @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.1, noisy = true, quantized=true, 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)) + sol = solve(prob, Tsit5()) + + @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 \ No newline at end of file