Adventures in the not so complex space.

Un petit exercise en la theorie de DSP

Emilio Jesús Gallego Arias, Pierre Jouvelot

Mines PARISTECH

Motivations

Mechanized Models for Sound Processing:

Program verification:

Education

Forget about numerical issues for now.

The Discrete Fourier Transform

Assume finite length signals $x \in \mathbb{C}^N$.

Then the DFT is:

$$X(k) = \sum_{n=0}^{N-1} x(n) * e^{-i 2 \pi n k / N} \quad k = 0,1,2,\ldots,N-1$$

and the inverse DFT is:

$$x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) * e^{i 2 \pi n k / N} \quad k = 0,1,2,\ldots,N-1$$ We ignore sampling interval, etc... as standard.

The DFT in Coq

$$ X(k) = \cssId{dft-sum}{\class{slide}{\sum_{n=0}^{N-1}}} \cssId{dft-signal}{\class{slide}{x(n)}} * \cssId{dft-unity-root}{\class{slide}{e^{-i 2 \pi n k / N}}} $$

Bigops: mature library.

No problem found, all the proofs straigforward.

matrix.v + ssrnum.v

No problem either, a couple of tweaks here and there.

Complex exponentiation?

Slight problem. Q: Do we need to resort to Coq (classical) reals? A: Look at what we need for the required theorems.

Complex Numbers in Coq

What we could find:

Our choice: algC

A Cool Theorem

C_prim_root_exists:
 <- closed_field_poly_normal [poly.v]
    {poly is factorizable}
 <- root_prod_XsubC          [poly.v]
    {factors are roots}
 <- has_prim_root            [cyclic.v]
    {be n' (>=n) diff n.-unity roots, one is primitive}
 <- separable_prod_XsubC     [separable.v]
    {separable factored form = uniqueness of factors}
 <- separable_Xn_sub_1       [cyclotomic.v]
    {X^n - 1 is separable}

We have a primitive unity root of type algC !

Cyclic Exponentiation


Definition expz_def (z : F) (i : 'I_n) := z ^+ i.
Definition expz := nosimpl expz_def.

Notation "x @+ n" := (expz x n)    : ring_scope.
Notation "x @- n" := (expz x (-n)) : ring_scope.

Hypothesis zP : N.-primitive_root z.

Lemma expzD m n : z @+ (m + n) = z @+ m * z @+ n.
Lemma expzM m n : z @+ (m * n) = z @+ m @+ n.
Lemma expzV n : (z @+ n)^-1 = z @- n.

Primitive Roots


Notation n := N.+2.
Hypothesis zP : n.-primitive_root z.

Lemma sum_expr_z_zero : \sum_(0 <= j < n) z ^+ j = 0.
Lemma norm_primX k    : `|z @+ k| = 1.
Lemma prim_conjC      : z^* = z^-1.
Lemma expz_conjC      : (z @+ k)^* = z @- k.
Lemma z_ortho: \sum_k (z @+ (i * k)) * (z @+ (j * k))^* = n * (i == j).

Basic Properties


Definition dft x k := \sum_n (x n 0) * 'ω ^+ (k*n).

Lemma dft_scale a x k : a * dft x k = dft (a *: x) k.
Lemma dft_sum x y k : dft (x + y) k = dft x k + dft y k.
Lemma dft_shifts x m k : dft (shifts x m) k = 'ω ^+ (k * m) * dft x k.
Lemma dft_flips x (k : 'I_N) : dft (flips x) k = dft x (-k).

(Circular) Convolution

$$(x \circledast y) (k) = \sum_{m=0}^{N-1} x(m) y(k-m)$$ $$\mathsf{DFT}_k~(x \circledast y) = X_k Y_k$$ 

Definition convs x y := \col_k \sum_m x m 0 * y (k-m) 0.

Lemma convsC : commutative convs.
Lemma dft_convs x y k : dft (convs x y) k = dft x k * dft y k.
[5 lines of proof...]

More Theorems

The next interesting theorem is Parseval's

$$\sum_{n=0}^{N-1} |x(n)|^2 = \frac{1}{N} \sum_{n=0}^{N-1} |X(n)|^2 $$
But things being to feel a little ad-hoc.

Geometric Signal Theory

An inner product over a vector space is an operation satisfying:
$$ \newcommand{\ip}[2]{\langle{#1},{#2}\rangle} \begin{array}{l} \ip{x}{y} = \overline{\ip{y}{x}} \\ \ip{ax}{y} = a\ip{x}{y} \\ \ip{x}{x} \ge 0\\ \ip{x}{x} = 0 => x = 0 \end{array} $$
A norm can be defined:
$$ \|x\| = \sqrt{\ip{x}{x}} $$

Satisfying the usual identities, with $algC^N$ we have a pre-Hilbert space.

Geometric Signal Theory

Now the DFT becomes

$$ 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.

The DFT in matrix form:

Definition W := \matrix_(i < n, j < n) 'ω ^+ (i*j).
Lemma dftE x k : (W *m x) k 0 = dft x k.
Lemma dftDr x y : W *m (x + y) = W *m x + W *m y.
Proof. exact: linearD. Qed.

In an Ideal World...

We'd use a theory of unitary matrices to derive the rest of the results.

Equivalent conditions (taken from WPD).

In an Ideal World...

We are not so far from it.

Problem to generalize: Hierarchy of classes.

Going Further...

Algebraic Signal Processing [Püschel et al. 2006-2015]

Filters from a ring $A$:

The operation of filtering defines an algebra.

More interesting examples in the paper.

The End

← →

/