A Coq Formalization of Digital Filters Diane Gallois-Wong, Sylvie - - PowerPoint PPT Presentation

A Coq Formalization of Digital Filters

Diane Gallois-Wong, Sylvie Boldo and Thibault Hilaire

Universit´ e Paris-Sud, LRI (Orsay), Inria Saclay

Calculemus - August 15, 2018

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 1 / 11

Signal Processing and Digital Filters

Signal: audio, video, various physical measurements Applications: communication, robotics, automotive, aeronautics, etc.

ime

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 2 / 11

Signal Processing and Digital Filters

Signal: audio, video, various physical measurements Applications: communication, robotics, automotive, aeronautics, etc. Digital signal:

Time

discrete time k ∈ Z − → value u(k) ∈ R

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 2 / 11

Signal Processing and Digital Filters

Signal: audio, video, various physical measurements Applications: communication, robotics, automotive, aeronautics, etc. Digital signal:

Time

discrete time k ∈ Z − → value u(k) ∈ R Digital filter H : H u(k)

input signal

y(k)

• utput signal

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 2 / 11

LTI Digital Filters

Digital filter H : H u(k)

input signal

y(k)

• utput signal

y(k) depends on u(k) but also on the past Example: y(k) = u(k) − 3y(k − 1)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 3 / 11

LTI Digital Filters

Digital filter H : H u(k)

input signal

y(k)

• utput signal

y(k) depends on u(k) but also on the past Example: y(k) = u(k) − 3y(k − 1) We are interested in Linear Time-Invariant (LTI) filters:

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 3 / 11

LTI Digital Filters

Digital filter H : H u(k)

input signal

y(k)

• utput signal

y(k) depends on u(k) but also on the past Example: y(k) = u(k) − 3y(k − 1) We are interested in Linear Time-Invariant (LTI) filters: valid operations are addition and scaling

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 3 / 11

LTI Digital Filters

Digital filter H : H u(k)

input signal

y(k)

• utput signal

y(k) depends on u(k) but also on the past Example: y(k) = u(k) − 3y(k − 1) We are interested in Linear Time-Invariant (LTI) filters: valid operations are addition and scaling if the input is delayed, then the output is delayed as well

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 3 / 11

Theory and Practice: Finite Precision

Theory: mathematical definition ∀k ∈ Z. y(k) = u(k) − 3y(k − 1) infinite precision (real numbers R) ↓ Practice: implementation foreach k do y(k) ← u(k) − 3 ∗ y(k − 1) end finite precision

(floating- or fixed-point numbers)

֒ → rounding errors

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 4 / 11

Theory and Practice: Finite Precision

Theory: mathematical definition ∀k ∈ Z. y(k) = u(k) − 3y(k − 1) infinite precision (real numbers R) ↓ Practice: implementation foreach k do y(k) ← u(k) − 3 ∗ y(k − 1) end finite precision

(floating- or fixed-point numbers)

֒ → rounding errors which propagate and may accumulate: y ∗(k) ← u(k) ⊖ 3 ⊗ y ∗(k − 1)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 4 / 11

Theory and Practice: Finite Precision

Theory: mathematical definition ∀k ∈ Z. y(k) = u(k) − 3y(k − 1) infinite precision (real numbers R) ↓ Practice: implementation foreach k do y(k) ← u(k) − 3 ∗ y(k − 1) end finite precision

(floating- or fixed-point numbers)

֒ → rounding errors which propagate and may accumulate: y ∗(k) ← u(k) ⊖ 3 ⊗ y ∗(k − 1) y ∗(k) ← ֓ y ∗(k − 1) ← ֓ y ∗(k − 2) ← ֓ ...

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 4 / 11

Rounding Errors in Digital Filters

Digital filters in embedded systems use fixed-point numbers: consume less energy, less expensive than floating-point numbers. Optimisation: trying to use as few bits as possible ֒ → difficult, especially as errors may accumulate unexpectedly

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 5 / 11

Rounding Errors in Digital Filters

Digital filters in embedded systems use fixed-point numbers: consume less energy, less expensive than floating-point numbers. Optimisation: trying to use as few bits as possible ֒ → difficult, especially as errors may accumulate unexpectedly ֒ → more efficient algorithms have bigger risk of sizeable final error or overflow (exceeding the greatest representable value)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 5 / 11

Rounding Errors in Digital Filters

Digital filters in embedded systems use fixed-point numbers: consume less energy, less expensive than floating-point numbers. Optimisation: trying to use as few bits as possible ֒ → difficult, especially as errors may accumulate unexpectedly ֒ → more efficient algorithms have bigger risk of sizeable final error or overflow (exceeding the greatest representable value) How to ensure that rounding errors do not cause critical failures in robotics, automotive, aeronautics etc.?

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 5 / 11

Rounding Errors in Digital Filters

Digital filters in embedded systems use fixed-point numbers: consume less energy, less expensive than floating-point numbers. Optimisation: trying to use as few bits as possible ֒ → difficult, especially as errors may accumulate unexpectedly ֒ → more efficient algorithms have bigger risk of sizeable final error or overflow (exceeding the greatest representable value) How to ensure that rounding errors do not cause critical failures in robotics, automotive, aeronautics etc.? Error analysis with pen-and-paper proofs [Hilaire, Lopez 2013] etc. Formal methods [Akbarpour, Tahar 2007] [Siddique, Mahmoud, Tahar 2018] etc.

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 5 / 11

Contribution

Formalization in Coq: Definitions: signals and Linear Time-Invariant (LTI) filters Various realizations for filters and equivalences between them Theorem of the Error Filter to study propagation of errors Worst-Case Peak-Gain Theorem to bound the final error These are essential steps toward a fully proven rounding error analysis.

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 6 / 11

Defining Signals

Signal: function Z → R that takes the value 0 for all k < 0

Definition causal (x : Z → R) := (forall k : Z, (k < 0)%Z → x k = 0%R). Record signal := { signal_val :> Z → R ; signal_prop : causal signal_val }.

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 7 / 11

Defining Signals

Signal: function Z → R that takes the value 0 for all k < 0

Definition causal (x : Z → R) := (forall k : Z, (k < 0)%Z → x k = 0%R). Record signal := { signal_val :> Z → R ; signal_prop : causal signal_val }.

Why Z rather than nat?

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 7 / 11

Defining Signals

Signal: function Z → R that takes the value 0 for all k < 0

Definition causal (x : Z → R) := (forall k : Z, (k < 0)%Z → x k = 0%R). Record signal := { signal_val :> Z → R ; signal_prop : causal signal_val }.

Why Z rather than nat? + easier handling of initial conditions y(k) = u(k) − 3y(k − 1) + 5y(k − 3) in nat: separate cases for k ∈ {0, 1, 2}

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 7 / 11

Defining Signals

Signal: function Z → R that takes the value 0 for all k < 0

Definition causal (x : Z → R) := (forall k : Z, (k < 0)%Z → x k = 0%R). Record signal := { signal_val :> Z → R ; signal_prop : causal signal_val }.

Why Z rather than nat? + easier handling of initial conditions y(k) = u(k) − 3y(k − 1) + 5y(k − 3) in nat: separate cases for k ∈ {0, 1, 2} + better readability of theorems + more intuitive substraction (natural numbers in Coq: 3 − 5 = 0) − less library support (but often easy to adapt from nat to Z)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 7 / 11

One Realization to Represent them all

Filters can be defined under various forms, called realizations:

y(k) =

n

• i=0

biu(k − i) −

n

• i=1

aiy(k − i)        e(k) ← u(k) −

n

• i=1

aie(k − i) y(k) ←

n

• i=1

bie(k − i)

• x(k + 1)

= Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)

Coq definitions and proofs : Three realizations: Direct Forms I and II, and State-Space Corresponding filters are LTI Transformation functions between them that preserve the filter

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 8 / 11

One Realization to Represent them all

Filters can be defined under various forms, called realizations:

y(k) =

n

• i=0

biu(k − i) −

n

• i=1

aiy(k − i)        e(k) ← u(k) −

n

• i=1

aie(k − i) y(k) ←

n

• i=1

bie(k − i)

• x(k + 1)

= Ax(k) + Bu(k) y(k) = Cx(k) + Du(k)

Coq definitions and proofs : Three realizations: Direct Forms I and II, and State-Space Corresponding filters are LTI Transformation functions between them that preserve the filter ֒ → we can choose one of them to do the error analysis

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 8 / 11

Error Propagation: the Error Filter

Reminder

Digital filter H : H u(k)

input signal

y(k)

• utput signal

Example: y(k) = u(k) − 3y(k − 1) Implementation: y ∗(k) ← u(k) ⊖ 3 ⊗ y ∗(k − 1)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 9 / 11

Error Propagation: the Error Filter

Model: y(k) = u(k) − 3y(k − 1) Implementation: y ∗(k) = u(k) ⊖ 3 ⊗ y ∗(k − 1) (in reality, we use the previously chosen realization which is more complex, but the principle is the same)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 9 / 11

Error Propagation: the Error Filter

Model: y(k) = u(k) − 3y(k − 1) Implementation: y ∗(k) = u(k) ⊖ 3 ⊗ y ∗(k − 1) = u(k) − 3 × y ∗(k − 1) + ε(k)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 9 / 11

Error Propagation: the Error Filter

Model: y(k) = u(k) − 3y(k − 1) Implementation: y ∗(k) = u(k) ⊖ 3 ⊗ y ∗(k − 1) = u(k) − 3 × y ∗(k − 1) + ε(k)

Proven theorem

We can compute a filter Hε such that: Hε ε(k) ∆y(k) where ∆y(k) = y ∗(k) − y(k) is the final error after propagation

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 9 / 11

Error Propagation: the Error Filter

Model: y(k) = u(k) − 3y(k − 1) Implementation: y ∗(k) = u(k) ⊖ 3 ⊗ y ∗(k − 1) = u(k) − 3 × y ∗(k − 1) + ε(k)

Proven theorem

We can compute a filter Hε such that: Hε ε(k) ∆y(k) where ∆y(k) = y ∗(k) − y(k) is the final error after propagation Fixed-point principle: individual computations are correctly rounded ֒ → bound on ε(k)

???

− → bound on ∆y(k)

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 9 / 11

The Worst-Case Peak-Gain Theorem

Proven theorem

The Worst-Case Peak-Gain H = |D| + ∞

k=0

• CAkB
• f filter H

verifies: bound for the input

× H

• −

− − − − − − → bound for the output ∀k. |u(k)| ≤ M ∀k. |y(k)| ≤ H × M H The value of H is in R (it may be equal to +∞) and it is optimal

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 10 / 11

Conclusion and Future Work

Formalized1 in Coq: Equivalence between realizations, so we can choose one of them Theorem of the Error Filter to study propagation of errors Worst-Case Peak-Gain Theorem to bound the final error

1www.lri.fr/~gallois/code/coq-digital-filters-CICM18.tgz

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 11 / 11

Conclusion and Future Work

Formalized1 in Coq: Equivalence between realizations, so we can choose one of them Theorem of the Error Filter to study propagation of errors Worst-Case Peak-Gain Theorem to bound the final error To complete formal error analysis: Formalize two’s complement fixed-point arithmetic in Coq Prove that

H = |D| + ∞

k=0

• CAkB
• can be well approximated

Bound errors in correctly rounded individual computations

1www.lri.fr/~gallois/code/coq-digital-filters-CICM18.tgz

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 11 / 11

Conclusion and Future Work

Formalized1 in Coq: Equivalence between realizations, so we can choose one of them Theorem of the Error Filter to study propagation of errors Worst-Case Peak-Gain Theorem to bound the final error To complete formal error analysis: Formalize two’s complement fixed-point arithmetic in Coq Prove that

H = |D| + ∞

k=0

• CAkB
• can be well approximated

Bound errors in correctly rounded individual computations Longer-term goal: Automatically generate certified code together with formally proven error bound ...while reordering computations to minimize the rounding errors

1www.lri.fr/~gallois/code/coq-digital-filters-CICM18.tgz

Diane Gallois-Wong (U-PSud, LRI) A Coq Formalization of Digital Filters Calculemus - August 15, 2018 11 / 11