From 2f167415c2f3a35289a312bc53f047bf536c4c18 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 4 Sep 2024 15:09:25 +0200 Subject: [PATCH 1/3] add Sampler with AD effects --- docs/src/tutorials/noise.md | 30 +++++++++++++++++ src/ModelingToolkitSampledData.jl | 2 +- src/discrete_blocks.jl | 56 +++++++++++++++++++++++++++++-- 3 files changed, 85 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index d2b115a..3199a2d 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,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") +``` \ No newline at end of file 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..748d4a6 100644 --- a/src/discrete_blocks.jl +++ b/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 \ No newline at end of file From c6dcd11c2c3c231d35c5902e4ab83d66c57feb91 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 5 Sep 2024 12:56:28 +0200 Subject: [PATCH 2/3] add test --- docs/src/tutorials/noise.md | 18 ++++++++++++------ src/discrete_blocks.jl | 4 ++-- test/test_discrete_blocks.jl | 27 +++++++++++++++++++++++++++ 3 files changed, 41 insertions(+), 8 deletions(-) diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index 3199a2d..eb293ce 100644 --- a/docs/src/tutorials/noise.md +++ b/docs/src/tutorials/noise.md @@ -151,13 +151,13 @@ The default option `midrise = true` includes both end points as possible output ## 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 +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.05, noisy = false, quantized=false, midrise=true) + 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 @@ -170,10 +170,16 @@ 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)) +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") -# plot!(sol, idxs=m.sampling.sampler.y, label="AD converted output") -``` \ No newline at end of file +``` + +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/discrete_blocks.jl b/src/discrete_blocks.jl index 748d4a6..bf5172d 100644 --- a/src/discrete_blocks.jl +++ b/src/discrete_blocks.jl @@ -910,7 +910,7 @@ 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 @equations begin @@ -1057,7 +1057,7 @@ The operations occur in the order else noise = Gain(; k = 1) end - quantization = Quantization(; y_min, y_max, midrise, quantized) + quantization = Quantization(; bits, y_min, y_max, midrise, quantized) end @equations begin connect(input, sampler.input) diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index c6c0dbd..2f92289 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" + + @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 From a8a5a659cc67f87ee3bbf578f05d5fdbc11247e5 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Fri, 6 Sep 2024 06:24:14 +0200 Subject: [PATCH 3/3] define z --- test/test_discrete_blocks.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index 2f92289..c3f4446 100644 --- a/test/test_discrete_blocks.jl +++ b/test/test_discrete_blocks.jl @@ -518,7 +518,7 @@ 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)