Quantization (#10)

* add stable noise sources and tutorial

* add quantization

* add stats

* add quantization tutorial

* add quantization mode midtread

Project.toml

@@ -14,18 +14,18 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [compat]
 DiffEqCallbacks = "~3.8"
-julia = "1.10"
-StableRNGs = "1"
 JuliaSimCompiler = "0.1.19"
 ModelingToolkit = "9"
 ModelingToolkitStandardLibrary = "2"
 OrdinaryDiffEq = "6.89"
 Random = "1"
-
+StableRNGs = "1"
+julia = "1.10"
 
 [extras]
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ControlSystemsBase", "Test"]
+test = ["ControlSystemsBase", "Test", "Statistics"]

docs/src/tutorials/noise.md

@@ -119,4 +119,50 @@ Internally, a random number generator from [StableRNGs.jl](https://github.com/Ju
 2. Multiple calls to the random number generator at the same time step all return the same number.
 
 ## Quantization
-Not yet available.
+
+A signal may be quantized to a fixed number of levels (e.g., 8-bit) using the [`Quantization`](@ref) block. This may be used to simulate, e.g., the quantization that occurs in a AD converter. Below, we have a simple example where a sine wave is quantized to 2 bits (4 levels), limited between -1 and 1:
+```@example QUANT
+using ModelingToolkit, ModelingToolkitSampledData, OrdinaryDiffEq, Plots
+using ModelingToolkit: t_nounits as t, D_nounits as D
+z = ShiftIndex(Clock(0.1))
+@mtkmodel QuantizationModel begin
+ @components begin
+ input = Sine(amplitude=1.5, frequency=1)
+ quant = Quantization(; z, bits=2, y_min = -1, y_max = 1)
+ 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, quant.input)
+ D(x) ~ 0 # Dummy equation
+ end
+end
+@named m = QuantizationModel()
+m = complete(m)
+ssys = structural_simplify(IRSystem(m))
+prob = ODEProblem(ssys, [], (0.0, 2.0))
+sol = solve(prob, Tsit5())
+plot(sol, idxs=m.input.output.u)
+plot!(sol, idxs=m.quant.y, label="Quantized output")
+```
+
+
+
+### Different quantization modes
+With the default option `midrise = true`, the output of the quantizer is always between `y_min` and `y_max` inclusive, and the number of distinct levels it can take is `2^bits`. The possible values are given by
+```@example
+bits = 2; y_min = -1; y_max = 1
+collect(range(y_min, stop=y_max, length=2^bits))
+```
+Notably, these possible levels _do not include 0_. If `midrise = false`, a mid-tread quantizer is used instead. The two options are visualized below:
+```@example QUANT
+y_min = -1; y_max = 1; bits = 2
+u = y_min:0.01:y_max
+y_mr = ModelingToolkitSampledData.quantize_midrise.(u, bits, y_min, y_max)
+y_mt = ModelingToolkitSampledData.quantize_midtread.(u, bits, y_min, y_max)
+plot(u, [y_mr y_mt], label=["Midrise" "Midtread"], xlabel="Input", ylabel="Output", framestyle=:zerolines, l=2, seriestype=:step)
+```
+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.

src/ModelingToolkitSampledData.jl

@@ -7,7 +7,7 @@ export get_clock
 export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler,
  ClockChanger,
  DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
- DiscreteTransferFunction, NormalNoise, UniformNoise
+ DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization
 include("discrete_blocks.jl")
 
 end

src/discrete_blocks.jl

@@ -876,4 +876,67 @@ See also [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK
 # push!(eqs, input.u ~ u)
 # push!(eqs, output.u ~ y)
 # compose(ODESystem(eqs, t, sts, pars; name = name), input, output)
-# end
+# end
+
+"""
+ Quantization
+
+A quantization block that quantizes the input signal to a specified number of bits.
+
+# Parameters:
+- `y_max`: Upper limit of output
+- `y_min`: Lower limit of output
+- `bits`: Number of bits of quantization
+- `quantized`: If quantization effects shall be computed. If false, the output is equal to the input, which may be useful for, e.g., linearization.
+
+# Connectors:
+- `input`
+- `output`
+
+# Variables
+- `y`: Output signal, equal to `output.u`
+- `u`: Input signal, equal to `input.u`
+"""
+@mtkmodel Quantization begin
+ @extend u, y = siso = SISO()
+ @structural_parameters begin
+ z = ShiftIndex()
+ midrise = true
+ end
+ @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"]
+ 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
+end
+
+function quantize_midrise(u, bits, y_min, y_max)
+ d = y_max - y_min
+ y1 = clamp(u, y_min, y_max)
+ y2 = (y1 - y_min) / d # between 0 and 1
+ Δ = 2^Int(bits)-1
+ y3 = round(y2 * Δ) / Δ # quantized between 0 and 1
+ y4 = y3*d + y_min
+ return y4
+end
+
+function quantize_midtread(u, bits, y_min, y_max)
+ Δ = (y_max - y_min) / (2^Int(bits)-1)
+ # clamp(Δ * floor(u / Δ + 0.5), y_min, y_max)
+ k = sign(u) * max(0, floor((abs(u) - Δ/2) / Δ + 1))
+ y0 = sign(k) * (Δ/2 + Δ*(abs(k)-1/2))
+ y1 = iszero(y0) ? zero(y0) : y0 # remove -0.0
+ return clamp(y1, y_min, y_max - Δ/2)
+end
+
+function quantize(u, bits, y_min, y_max, midrise)
+ midrise ? quantize_midrise(u, bits, y_min, y_max) : quantize_midtread(u, bits, y_min, y_max)
+end
+
+@register_symbolic quantize(u::Real, bits::Real, y_min::Real, y_max::Real, midrise::Bool)

test/test_discrete_blocks.jl

@@ -324,7 +324,33 @@ using Statistics
  @named m = NoiseModel()
  m = complete(m)
  ssys = structural_simplify(IRSystem(m))
- prob = ODEProblem(ssys, [], (0.0, 10.0))
+ prob = ODEProblem(ssys, [m.noise.y(k-1) => 0], (0.0, 10.0))
+ sol = solve(prob, Tsit5())
+ @test !all(iszero, sol.u)
+ tv = 0:k.clock.dt:sol.t[end]
+ @test std(sol(tv, idxs = m.plant.u)) ≈ 1 rtol=0.1
+ @test mean(sol(tv, idxs = m.plant.u)) ≈ 0 atol=0.08
+end
+
+@testset "UniformNoise" begin
+ k = ShiftIndex(Clock(0.01))
+
+ @mtkmodel NoiseModel begin
+ @components begin
+ noise = UniformNoise(z = k)
+ zoh = ZeroOrderHold(z = k)
+ plant = FirstOrder(T = 1e-4) # Included due to bug with only discrete-time systems
+ end
+ @equations begin
+ connect(noise.output, zoh.input)
+ connect(zoh.output, plant.input)
+ end
+ end
+
+ @named m = NoiseModel()
+ m = complete(m)
+ ssys = structural_simplify(IRSystem(m))
+ prob = ODEProblem(ssys, [m.noise.y(k-1) => 0], (0.0, 10.0))
  sol = solve(prob, Tsit5())
  @test !all(iszero, sol.u)
  tv = 0:k.clock.dt:sol.t[end]
@@ -385,3 +411,79 @@ end
 # @test reduce(vcat, sol((0:10) .+ 1e-2))[:]≈[zeros(2); 1; zeros(8)] atol=1e-2
 # end*
 
+
+@testset "quantization" begin
+ @info "Testing quantization"
+
+ function test_quant(y_min, y_max, bits)
+ u = y_min:(1/2^bits):y_max
+ y = ModelingToolkitSampledData.quantize_midrise.(u, bits, y_min, y_max)
+ uy = unique(y)
+ @test uy ≈ range(y_min, stop=y_max, length=2^bits)
+ end
+
+ test_quant(-1, 1, 2) # Symmetric
+ test_quant(-1, 2, 2) # Not symmetric
+ test_quant(-1, 1, 3) # Symmetric, uneven number of bits
+ test_quant(-5, -2, 2) # Only negative
+ test_quant(5, 12, 2) # Only positive
+ test_quant(-5, 12, 20) # Large number of bits
+
+
+
+ function test_quant2(y_min, y_max, bits)
+ u = y_min:(1/2^bits):y_max
+ y = ModelingToolkitSampledData.quantize_midtread.(u, bits, y_min, y_max)
+ uy = unique(y) 
+ # @test (2^bits - 2 <= length(uy) <= 2^bits) # This check is not reliable since there might be one ulp difference when the last clamp is applied
+ @test maximum(y) <= y_max
+ @test minimum(y) >= y_min
+ end
+
+ test_quant2(-1, 1, 2) # Symmetric
+ test_quant2(-1, 2, 2) # Not symmetric
+ test_quant2(-1, 1, 3) # Symmetric, uneven number of bits
+ test_quant2(-5, -2, 2) # Only negative
+ test_quant2(5, 12, 2) # Only positive
+ test_quant2(-5, 12, 20) # Large number of bits
+
+ z = ShiftIndex(Clock(0.1))
+ @mtkmodel QuantizationModel begin
+ @components begin
+ input = Sine(amplitude=1, frequency=1)
+ quant = Quantization(; z, bits=2)
+ end
+ @equations begin
+ connect(input.output, quant.input)
+ end
+ end
+ @named m = QuantizationModel()
+ m = complete(m)
+ ssys = structural_simplify(IRSystem(m))
+ prob = ODEProblem(ssys, [], (0.0, 10.0))
+ sol = solve(prob, Tsit5(), dtmax=0.01)
+ y = sol[m.quant.y]
+ uy = unique(y)
+ @test length(uy) == 4
+
+
+ @mtkmodel QuantizationModel2 begin
+ @components begin
+ input = Sine(amplitude=1, frequency=1)
+ quant = Quantization(; z, bits=2, midrise=false)
+ end
+ @equations begin
+ connect(input.output, quant.input)
+ end
+ end
+ @named m = QuantizationModel2()
+ m = complete(m)
+ ssys = structural_simplify(IRSystem(m))
+ prob = ODEProblem(ssys, [], (0.0, 10.0))
+ sol = solve(prob, Tsit5(), dtmax=0.01)
+ y = sol[m.quant.y]
+ uy = unique(y)
+ @test length(uy) == 4
+end
+
+