diff --git a/Project.toml b/Project.toml index 378659a..5324fa1 100644 --- a/Project.toml +++ b/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"] diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index b4d9f37..2336ced 100644 --- a/docs/src/tutorials/noise.md +++ b/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. \ No newline at end of file + +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. diff --git a/src/ModelingToolkitSampledData.jl b/src/ModelingToolkitSampledData.jl index 58ce736..c578b54 100644 --- a/src/ModelingToolkitSampledData.jl +++ b/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 diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl index 8d89fc6..f6683ed 100644 --- a/src/discrete_blocks.jl +++ b/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 \ No newline at end of file +# 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) diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index 262a2e2..d81a568 100644 --- a/test/test_discrete_blocks.jl +++ b/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 + +