## Experiments in Certified Digital Audio Processing

### Emilio J. Gallego Arias (joint work with Pierre Jouvelot)

#### MINES ParisTech, PSL Research University, France

##### 11 Avril 2016 - Specfun - Palaiseau

Real Time Signal Processing
$\cap$
Programming Language Theory
$\cap$
Theorem Proving

FEEVER ANR Project

## Audio Objects

###### Recursive 2nd order Filter

Image Credits: Julius Orion Smith III

###### Waveguide Resonator

Image Credits: Julius Orion Smith III

#### Questions?

• Implementation ease?
• Composability?
• Mathematical Semantics?
• FP Semantics?
• PL Semantics!
• Domain practice.
• Frequency response.
• State Space.
• Optimization.

Julius Orion Smith III audio-focused book's series:

• Mathematics of the Discrete Fourier Transform (DFT)
• [Partially handled + missing inf. sum.]
• Introduction to Digital Filters
• [Focus of this talk]
• Physical Audio Signal Processing
• [Minor bits in this talk + missing ODE]
• Spectral Audio Signal Processing
• [Future work]
##### Main Points of the Series:
• Free mix of mathematics and a little computation.
• Linear algebra is pervasive.
• Proofs of uneven difficulty, not constructive-friendly.

# Where to Look?

## DSP & Theorem Proving

### Isabelle/HOL[Akbarpour, Tahar et al.]

Focus on fixed systems, numerical issues.

### Lucid, Velus[Boulmé, Auger, ...]

Safety Property, some work unpublished?

### Kahn Networks[Paulin-Mohring]

Domain Theory point of view, setoid-rewriting based.

### VeriDrone Project[Ricketts, Machela, Lerner, ...]

Work in progress, from LTL to Lyapunov!

### FFT in Coq [Capretta]

Main successor of our current world view.

## DSP & Programming Languages

### Arrows/FRP [Elliot, Hughes, Hudak, Nilsson, ...]

Too abstract, not convenient for sample-based DSP view.

### Synchronous Languages [Berry, Caspi, Halbwachs, Pouzet, ...]

Far from a "linear algebra" view, not particularly well suited for DSP, (state of the art: [Guatto]'s PhD).
Data-intensive vs control-intensive require quite different control techniques. [Berry, 2000]

### DataFlow Languages [Lee, Thies, ...]

IMHO most successful approach so far, weak on formalization?

### String Diagrams [Bonchi, Sobocinski, Zanasi]

Very interesting recent research, covers a limited subset of DSP.

### Other Approaches

VOBLA, Ziria, Halide, Darkroom, Julia, Spiral, FeldSpar
Varying degrees of formal semantics.

## A First Try: Mini-Faust & Coq

• Start from Faust [Orlarey 2002]: functional DSL for DSP based on arrow-like point-free combinators.
• Implement in Coq: step-indexing, math-comp.
• Works nicely for easy examples; supports self-composition, etc...

### Problems with the approach:

• Theorems in DSP not expressed in point-free combinators.
• PL methods awkward for DSP experts.
• Doesn't work so well for more complex examples.

## Example: Filter Stability

### In PL terms:

$$\fjud{x \in [a,b]}{\mathsf{*(1-c) : + \sim *(c)}}{x \in [a,b]}$$

### In DSP terms:

The impulse response decays to 0 as time goes to infinity.

High number of components make PL proof methods impractical.

## What now? Wagner

Simple λ-calculus with feedback, variables carry a time index.

### Types and Syntax:


## Wagner: Operational Semantics

We use a big-step, time-indexed style. Every program reduces to a value at a given time step $n$. Equivalent to synchronous dataflow.

$$\newcommand{\ev}[3]{{#1} \downarrow_{#2} {#3}} \begin{gather} \inferrule { } {\ev{x}{n}{x}} \qquad \inferrule {\ev{e_1}{n}{v_1} \quad … \quad \ev{e_k}{n}{v_k} \quad P(v_1, …, v_1) ↓ v } {\ev{P(e_1,…,e_k)}{n}{v}} \\[1em] \inferrule {\ev{eₐ}{n}{vₐ} \quad \ev{e_f[x/vₐ]}{n}{v}} {\ev{(λ~x.e_f)~eₐ}{n}{v}} \qquad \inferrule {\ev{eₐ}{n}{vₐ} \quad \ev{e_f[x/vₐ]}{n}{v}} {\ev{(Λ~x.e_f)~eₐ}{n}{v}} \\[1.5em] \inferrule { \ev{e_1}{n}{v_1} ~\ldots~ \ev{e_k}{n}{v_k} } { \ev{ \{e_1, \ldots, e_k \} }{n}{\{v_1, \ldots, v_k \} } } \qquad \inferrule { } {\ev{e_{\{k+1\}}}{0}{0_e}} \qquad \inferrule {\ev{e_{\{k\}}}{n}{v}} {\ev{e_{\{k+1\}}}{n+1}{v}} \\[2em] \inferrule {\ev{e[x/0_e]}{0}{v}} {\ev{\feed{x}{e}}{0}{v}} \qquad \inferrule {\ev{\feed{x}{e}}{n}{v_f} \quad \ev{e[x/v_f]}{n+1}{v}} {\ev{\feed{x}{e}}{n+1}{v}} \end{gather}$$

## Examples: DF-I Filter

$$\mathsf{df1} ≡ λ a b : R[k]. Λ x.~\feed{y}{ \{ b \cdot x + a \cdot y} \}$$

Note the overload of $x \cdot y$! The types of $a$ and $b$ will fix the op.

$$x \cdot y ≡ \sum_i x[i] ⬝ y[i]$$

## Examples: DF-II

$$\mathsf{df2} ≡ Λ x.~ \feedin{v}{ \{ x + a \cdot v \}}{ \{ b \cdot v \} }$$

Compare with df1:

$$\mathsf{df1} ≡ Λ x.~\feed{y}{ \{ b \cdot x + a \cdot y} \}$$

## Examples: WaveGuide OSC

$$\begin{array}{rcl} \text{feed}~x & = & C \cdot (G \cdot x + y) - y \\ y & = & C \cdot (G \cdot x + y) + G \cdot x \end{array}$$

or without sugar:

$$\feed{z}{ \{ (C \cdot (G \cdot z_1 + z_2) - z_2, C \cdot (G \cdot z_1 + z_2) + G \cdot z_1) \} }$$

## Typing and Equational Theory

Typing is almost standard: we restrict pre to variables and use coeffect-inspired annotations to keep track of history-access. We also subtype rated streams to arrays:

$$\begin{gather} \inferrule { } { \Gamma, x :_n T \vdash x_{\{n\}} : T} \qquad \inferrule { \Gamma, x :_n R@r \vdash x : R@r } { \Gamma, x :_n R@r \vdash x : R[rn] } \end{gather}$$

We can recover the general pre as:

$$\mathbf{pre}~e ≡ (Λ x.~x_{\{1\}})~e$$

Application has to distinguish between positive and negative types to correctly propagate environment annotations. For this talk we use a restricted functional type.

Note: Work in coeffects, focusing, and duality is relevant and moslty subsumes our rules.

## Now the fun starts

Lightweight embedding into Coq, use a mildy-typed syntax for a reduced subset (recall that: $R ⇒ R ⇒ R ≈ R ⊗ R → R$):

## Interpreting the Programs:

We want an easy embedding and "executable" embedding, our evaluation relation is not. Basically we are looking for something such as:

$$⟦ \_ ⟧ : ∀ t : ty, \mathsf{seq}~ℝ → \mathsf{expr}~t → ℕ → \mathsf{tyI}~t$$

However, what should the value of:

$$⟦ Γ ⊢ x : ℝ ⟧_n = ~??$$

Indeed, we need to equiq the environments with history, thus they become list of lists of values. A first try for the interpretation could be:

$$\begin{array}{lcl} \mathsf{tyI}(⊗~n) &= & \mathsf{cV[ℝ]_n} \\ \mathsf{tyI}(n_a ⇒ n_r) &= & \mathsf{seq~cV[ℝ]_{n_a} → cV[ℝ]_{n_r}} \\ \end{array}$$

## Step-Indexing

However it is not enough, as when interpreting λ we don't have enough history!

$$⟦ Γ ⊢ Λ x . e ⟧_n = \mathsf{fun}~x ⇒ ⟦ Γ, x ⊢ e ⟧_n$$

We must use a strong induction principle so the application rule can compute all previous values:

Many more refinements are possible!

## Case Study I: Linearity

The first application is to study the set of linear Wagner programs, as a first step towards the theory for linear systems and in particular, the Z-transform.

What should linearity mean for our programs? The proper definition involves something known as a logical relation:

$$\begin{array}{lcl} \mathsf{linR}_{⊗ n}(e) &≡& \forall Γ_1, Γ_2.~ ⟦Γ_1 ⊢ e⟧ - ⟦Γ_2 ⊢ e⟧ = ⟦Γ_1 - Γ_2 ⊢ e⟧ \\ \mathsf{linR}_{n1 ⇒ n2}(f) &≡& \forall Γ_1, Γ_2, r_1, r_2.\\ && \quad ⟦Γ_1 ⊢ f⟧ r_1 - ⟦Γ_2 ⊢ f⟧ r_2 = ⟦Γ_1 - Γ_2 ⊢ f⟧ (r_1 - r_2) \end{array}$$

The relation is simple but would also work at higher-types. If we restrict the syntax to linear programs, we every program satisfies the relation!

The proof proceeds first, by induction on the typing derivation, then by induction on the number of steps.

## Case Study I: Linearity

Compare: $$\begin{array}{l} λ x. Λ y . x \cdot y : R → S ⇒ T \\ Λ x. Λ y . x \cdot y : R ⇒ S ⇒ T \end{array}$$

Work in progress: Relate to the completeness theorem for signal dataflow! [Bonchi, Sobociński, Zanasi 2015/2016] [Baez]

## Case Study II: Frequency Domain

#### Quick Review

• The frequency response of a filter = output for the impulse response.
• A filter is stable = its Frequency Response is summable.
• Stability is implied by the existence of the DTFT.
• The DFT approximates the DTFT.
• The DTFT is the evaluation of the Z-transform for the unit circle.

#### Crucial Prerrequisites:

• Previous linearity theorem.
• A notion of DFT in Coq.

## Geometric Signal Theory

#### The DFT:

$$X(\omega_k) = \ip{x^T, \omega_k} ~~\text{where}~~ \omega_k = \omega^{k*0},\ldots,\omega^{k*(N-1)}$$

Most properties are an immediate consequency of matrix multiplication lemmas. Again life is easy here.

##### Primitive roots:

The constructive algebraic numbers in mathcomp provides us with a primitive root of the unity such that $ω^N = 1$.

## In an Ideal World...

#### ... we'd have everything we need in our library.

• We are not so far from it.
• Took the definitions from classfun.v and proved:

Energy theorems are trivial corollaries. Quite compact development for now.

## Transfer Functions

We want to relate programs to their frequency semantics; it is well known that:

$$H(e^{jωT}) = \text{DTFT}_{ωT}(h)$$

We'd like to relate programs to their transfer function:

## Z-interpretation of programs

We tried to build an effective procedure, however it is difficult.

## Z-interpretation of programs

We can thus prove:

We show then that evalution of the impulse response DFT corresponds to the polynomial plus a residual based on the poles.

Stability is implied if the magnitude of all the poles is less than 1. This also implies the existence of the DTFT.

## Open Topics

• Multirate processing: [quite a bit working]
• State Space analysis, control theory, MIMO systems
• Floating Point Interpretation.
• Schur-Cohn Stability Test
• Verification of transforms optimization. [Püschel]
• Computational Z-interpretation:
• Connection to more mainstream numerical approximation literature?

## Conclusions

• Implementation in Ocaml with a core type system. Full type system still not fully settled, interesting desing space.
• Same for Coq formalization, theorems are working but more experience is needed to choose some implementation details (for instance, choice of environment representation
• Next logical step: COLA/OLA.
• We are particularly excited to see what the relation with SPIRAL could be.
• Numerical issues are pervasive in this domain, but trying to get there. Library for floating point aritmetic in progress.
• Missing an integrated library for true real and complex numbers, recent effort by PY and Assia? Should we port Coquelicot?

## Context of the ANR project:

Practice of real time DSP still far from convenient.

##### Starting point, Faust:
• Faust [Orlarey 2002], a functional PL for audio programming.
• Abstracts away low-level complexity, efficient execution.
• A success, good number of users, high-interest topic. (CCRMA workshops, Kadenze course, etc...)
##### Faust's Future:
• Extend Faust to multirate processing.
• Reasoning about audio programs is not easy.

## A note on streams.

### The beggining of it all:

First exercise, port of Lucid paper to modern Coq/Ssreflect; // broken with coinductives! Dependent types not very scalable; decidability/extensionality of streams: bad.

Move to sequences!

"Inversion views", but better, wf_str is decidable now!

/