Quantization (#10)

* add stable noise sources and tutorial

* add quantization

* add stats

* add quantization tutorial

* add quantization mode midtread

Author: baggepinnen

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
+
+