[PPT] - When computer arithmetic and signal processing meet Anastasia PowerPoint Presentation

SLIDE 1

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

When computer arithmetic and signal processing meet

Anastasia Volkova, Thibault Hilaire, Christoph Lauter

Metalbim workshop, March 12th, 2018

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

1/23

SLIDE 2

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Plan

1 Context and problematics 2 Determining the Most Significant Bit 3 Compute the quantized transfer function 4 WCPG

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

2/23

SLIDE 3

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Applications: reliable systems

Digital filters:

Algorithms that transform digital signals

Do not need guarantee in the majority of applications
A guarantee is necessary in other applications.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

3/23

SLIDE 4

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Applications: reliable systems

Digital filters:

Algorithms that transform digital signals

Do not need guarantee in the majority of applications
A guarantee is necessary in other applications.

We are interested in guarantees related to the implementation of numerical algorithms, especially in the embedded systems.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

3/23

SLIDE 5

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Automation

The implementation is done in several steps:

Specifications Constraints Software Numerical algorithm Code generation Mathematical algorithm

Numerous constraints:

performance
surface
accuracy
energy consumption
memory
etc.

We are interested in the automated process of reliable implementation.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

4/23

SLIDE 6

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Integer arithmetic:

y = Y

w 2w−1 20

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 7

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Integer arithmetic:

y = Y

w 2w−1 20

Fixed-Point arithmetic:

y = Y · 2ℓ where ℓ is an implicit factor

m + 1 −ℓ w −2m 20 2−1 2ℓ

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 8

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Integer arithmetic:

y = Y

w 2w−1 20

Fixed-Point arithmetic:

y = Y · 2ℓ where ℓ is an implicit factor

m + 1 −ℓ w −2m 20 2−1 2ℓ

Floating-Point arithmetic:

y = (−1)s · Y · 2e where e is an explicit factor

exponent mantissa s

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 9

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Integer arithmetic:

y = Y

w 2w−1 20

Fixed-Point arithmetic:

y = Y · 2ℓ where ℓ is an implicit factor

m + 1 −ℓ w −2m 20 2−1 2ℓ

Floating-Point arithmetic:

y = (−1)s · Y · 2e where e is an explicit factor

exponent mantissa s

Interval arithmetic:

[y, y] =

y ∈ R | y y y
Signal processing / computer arithmetic – Metalibm workshop
T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 10

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Integer arithmetic:

y = Y

w 2w−1 20

Fixed-Point arithmetic:

y = Y · 2ℓ where ℓ is an implicit factor

m + 1 −ℓ w −2m 20 2−1 2ℓ

Floating-Point arithmetic:

y = (−1)s · Y · 2e where e is an explicit factor

exponent mantissa s

Interval arithmetic:

[y, y] =

y ∈ R | y y y
Multiple-Precision arithmetic: the size
f the mantissa varies dynamically

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 11

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Arithmetics

Fixed-Point arithmetic:

y = Y · 2ℓ where ℓ is an implicit factor For the implementation

Floating-Point arithmetic:

y = (−1)s · Y · 2e where e is an explicit factor For the error analysis

Interval arithmetic:

[y, y] =

y ∈ R | y y y
For the error analysis
Multiple-Precision arithmetic: the size
f the mantissa varies dynamically

For the error analysis

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/23

SLIDE 12

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

State-space representation

State-space representation of an LTI filter H: H x(k + 1) = Ax(k) + bu(k) y(k) = cx(k) + du(k)

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

6/23

SLIDE 13

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Range of variables

Worst-Case Peak Gain

Input u(k) Stable filter Output y(k)

Temps Amplitude

H

Amplification/Attenuation

Amplitude Temps

∀k, |u(k)| ¯ u ∀k, |y(k)| H ¯ u Worst-Case Peak Gain:

H

= h1 = |d| +

∞

k=0
cAkb
Signal processing / computer arithmetic – Metalibm workshop
T. Hilaire (thibault.hilaire@lip6.fr)

7/23

SLIDE 14

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Plan

1 Context and problematics 2 Determining the Most Significant Bit 3 Compute the quantized transfer function 4 WCPG

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

8/23

SLIDE 15

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Reliable numerical algorithms

Simulink to SIF Simulink graph SIF H(z)

H(z)

to SIF SIF

formats

SIF Fixed-Point algorithm Code generation C

Software

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

9/23

SLIDE 16

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Reliable numerical algorithms

Simulink to SIF Simulink graph SIF H(z)

H(z)

to SIF SIF

formats

SIF Fixed-Point algorithm Code generation C

Software

Most Significant Bit Least Significant Bit Fixed-Point format

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

9/23

SLIDE 17

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Problem of the format choice

Constraints: wordlength wy is fixed for the output variable y
Goal: no overflows, rigorous error bounds
Bonus: minimize the errors

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

10/23

SLIDE 18

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Problem of the format choice

Constraints: wordlength wy is fixed for the output variable y
Goal: no overflows, rigorous error bounds
Bonus: minimize the errors

H x(k + 1) = Ax(k) + bu(k) y(k) = cx(k) + du(k) Problem: find the smallest MSB position my such that for all k |y(k)|

2my

1 − 2−wy+1

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

10/23

SLIDE 19

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Problem of the format choice

Constraints: wordlength wy is fixed for the output variable y
Goal: no overflows, rigorous error bounds
Bonus: minimize the errors

H x(k + 1) = Ax(k) + bu(k) y(k) = cx(k) + du(k) Problem: find the smallest MSB position my such that for all k

H

¯ u = |y(k)|

2my

1 − 2−wy+1 Mathematical solution: my =

log2 (

H ¯ u) − log2

1 − 21−wy

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

10/23

SLIDE 20

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Problem of the format choice

Constraints: wordlength wy is fixed for the output variable y
Goal: no overflows, rigorous error bounds
Bonus: minimize the errors

H x(k + 1) = Ax(k) + bu(k) y(k) = cx(k) + du(k) Problem: find the smallest MSB position my such that for all k

H

¯ u = |y(k)|

2my

1 − 2−wy+1 Mathematical solution: my =

log2 (

H ¯ u) − log2

1 − 21−wy

Practical solution: control the accuracy of the WCPG such that my − my 1

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

10/23

SLIDE 21

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Propagation of the rounding errors

Exact filter H is: H♦ x♦(k + 1) = Ax♦(k) + bu(k) + εx(k) y (k) = cx (k) + du(k)

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

11/23

SLIDE 22

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Propagation of the rounding errors

Implemented filter H♦ is: H♦ x♦(k + 1) = ♦mx(Ax♦(k) + bu(k)) + εx(k) y♦(k) = ♦my (cx♦(k) + du(k)) where ♦m is an operator guaranteeing the faithful rounding: |♦m(x) − x| 2m−w+1

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

11/23

SLIDE 23

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Propagation of the rounding errors

Implemented filter H♦ is: H♦ x♦(k + 1) = Ax♦(k) + bu(k) + εx(k) y♦(k) = cx♦(k) + du(k) + εy(k) with |εx(k)| 2mx−wx+1 and |εy(k)| 2my−wy+1

H∆

y♦(k) u(k) y(k) ∆(k)

H

µ εx(k) εy(k) ∂ rounding errors

H♦

|∆(k)| H∆

2mx−wx+1

2my−wy+1

Signal processing / computer arithmetic – Metalibm workshop
T. Hilaire (thibault.hilaire@lip6.fr)

11/23

SLIDE 24

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Reliable algorithm

H∆

y♦(k) u(k) y(k) ∆(k)

H

µ εx(k) εy(k) ∂ rounding errors

H♦

1 Initial estimation of the MSB my for the exact filter H 2 Taking into account the errors induced by the initial formats

with the help of the filter H∆, then computation of the new MSB m♦

y

3 If m♦

y = my, return m♦ y

else my ← my + 1 and return to the step 2.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

12/23

SLIDE 25

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Plan

1 Context and problematics 2 Determining the Most Significant Bit 3 Compute the quantized transfer function 4 WCPG

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

13/23

SLIDE 26

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Computing the transfer function

Conversion Implemented filter Transfer function !

✑

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

14/23

SLIDE 27

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Computing the transfer function

Conversion SIF State-space Implemented filter Transfer function Conversion Conversion !

Transfer function of a state-space: H(z) = c(zI − A)−1b + d ✑

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

14/23

SLIDE 28

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Computing the transfer function

Conversion SIF State-space Implemented filter Transfer function Conversion Conversion !

Transfer function of a state-space: H(z) = c(zI − XEX−1)−1b + d = cX   

1 z−λ1

...

1 z−λn

   X−1b + d We compute an approximation H(z) of H(z):

H(z) =
i

biz−i

i

aiz−i ✑ The error

H −

H

(ejω)
may be arbitrarily small.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

14/23

SLIDE 29

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Computing the transfer function

Conversion SIF State-space Implemented filter Transfer function Conversion Conversion !

Transfer function of a state-space: H(z) = c(zI − XEX−1)−1b + d = cX   

1 z−λ1

...

1 z−λn

   X−1b + d We compute an approximation H(z) of H(z):

H(z) =
i

biz−i

i

aiz−i ✑ The error

H −

H

(ejω)
may be arbitrarily small.

Problem:

Signal processing / computer arithmetic – Metalibm workshop
T. Hilaire (thibault.hilaire@lip6.fr)

14/23

SLIDE 30

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H S

?

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 31

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H b H S

?

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 32

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H b H S

?

b S

Transformation from H to S is exact:

A =

      − a1 1 . . . ... . . . 1 − an . . .      

b =

     

b1 −

a1 b0 . . . . . .

bn −

an b0      

c =
1

· · ·

d =

b0

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 33

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H b H S

?

b S

− = ∆S

The difference of two filters is defined as:

S ∆S b S

y(k) b y(k) ∆y(k) u(k)

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 34

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H b H S

?

b S

− = ∆S − = ∆H ?

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 35

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

How to bound the error

Computing the transfer function H(z) of the state-space system S:

H b H S

?

b S

− = ∆S − = ∆H

Relationship between ∆S and ∆H:

H −

H

(ejω)
∆S

, ∀ω ∈ [0, 2π] where ∆S is again the WCPG of the system ∆S.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

15/23

SLIDE 36

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Plan

1 Context and problematics 2 Determining the Most Significant Bit 3 Compute the quantized transfer function 4 WCPG

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

16/23

SLIDE 37

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Reliable evaluation of the WCPG

Reliable WCPG is required for:

determination of the range of variables
analysis of the errors induced by finite-precision computations
bounding the error of the approximation of the transfer

function ✑

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

17/23

SLIDE 38

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Reliable evaluation of the WCPG

Reliable WCPG is required for:

determination of the range of variables
analysis of the errors induced by finite-precision computations
bounding the error of the approximation of the transfer

function

Problem in the case of MIMO filters

Compute S, approximation of the matrix

H

= |D| +

∞

k=0
CAkB
,

such that for an a priori given ε | H − S| < ε ✑ We use the dynamic Multiple-Precision arithmetic.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

17/23

SLIDE 39

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Truncation bound

We determine N such that the truncation error verifies

∞
k=0
CAkB
−

N

k=0
CAkB
ε1

Mathematical solution: N

log2

ε1 Mmin

log2 ρ(A)

where

M = n

l=1 |Rl| 1−|λl| |λl| ρ(A)

(Rl)i,j = Ci,lBl,j Practical solution:

Eigenvalues computed with LAPACK
Interval Arithmetic
Theory of Verified Inclusions (Rump)

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

18/23

SLIDE 40

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Evaluation: powering

N

k=0
CAkB
−

→ We propose: A = XEX−1 thus, Ak = XEkX−1 where X are the eigenvectors E are the eigenvalues

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

19/23

SLIDE 41

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Evaluation: powering

N
k=0
CAkB
−

N

k=0
CVTkV−1B
ε2

We propose: A = XEX−1 thus, Ak = XEkX−1 where X are the eigenvectors E are the eigenvalues In practice: V ≈ X et T = V−1AV + ∆

Multiple precision operations with a

priori given error bound

∆ is controlled to satisfy the bound

ε2 on the propagated error

T2 1 rigorously verified with

the Gershgorin theorem

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

19/23

SLIDE 42

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Evaluation: summation

We proceed by analogy

N

k=0

CVTkV−1B
− N

k=0

C′TkB′
ε3
N

k=0

C′TkB′

− N

k=0 |C′PkB′|

ε4
N

k=0 |C′PkB′| − N k=0 |Lk|

ε5
N

k=0 |Lk| −

SN

ε6

✑ Three basic bricks: XY + Z, X + |Y|, X−1

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

20/23

SLIDE 43

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Behind the scenes

To transform a filter into fixed-point code, we use...

Interval arithmetic
Verified Inclusions for the eigenvalues
Truncation bound for N
Dynamic Multiple-Precision arithmetic
Computation of Gershgorin circles
Evaluation du WCPG
Floating-Point IEEE 754 arithmetic
Computation of eigenvalues with LAPACK

. . . to verify an implementation in Fixed-Point arithmetic!

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

21/23

SLIDE 44

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Conclusion

H(z) H(z) to SIF formats SIF Simulink to SIF Fixed-Point algorithm Code generation Verification against specifications C VHDL

Software ASIC FPGA

Simulink graph SIF SIF

We made filter implementation reliable and automatic
Rigorous use of various arithmetics
Fixed-Point implementations that verify a priori error bounds

in both frequency and time domains

Several reliable bricks use together

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

22/23

SLIDE 45

Context and problematics Determining the Most Significant Bit Compute the quantized transfer function WCPG

Tool development

H(z) H(z) to SIF formats SIF Simulink to SIF Fixed-Point algorithm Code generation Verification against specifications C VHDL

Software ASIC FPGA

Simulink graph SIF SIF

Open-source tool
Implementation in C/C++, Sollya,

Python and Matlab

≈ 15000 lines of code
almost plugged with FloPoCo for FPGA

implementation

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

23/23

SLIDE 46

Frequency specifications

Section 4

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

24/23

SLIDE 47

Frequency specifications

Plan

5 Frequency specifications

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

1/5

SLIDE 48

Frequency specifications

Frequency response: H

ejω

=

H
ejω
magnitude

e

phase

∡H
ejω

Specifications:

βp2 βp1 ωs1 ωp1 ω

H(eiω)
ωp2 ωs2

βs π

stop-band pass-band stop-band

β

H(ejω)
β,

∀ω ∈ [ω1, ω2] ⊆ [0, π]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

2/5

SLIDE 49

Frequency specifications

Verification of an implementation

Rigorous verification against the specifications Implemented filter Boolean Frequency specifications

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

3/5

SLIDE 50

Frequency specifications

Verification of an implementation

Rigorous verification against the specifications Implemented filter Boolean Frequency specifications

Existing approaches:

by simulations
approximation of the

frequency response Our reliable approach:

no simulations, only proofs
rational and interval arithmetic

Goal:

Guarantee upon an implementation
Fast computation of this guarantee

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

3/5

SLIDE 51

Frequency specifications

Verification of a transfer function

We need to show that ∀z = ejω, ω ∈ Ω ⊆ [0, π] β |H(z)| β

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

4/5

SLIDE 52

Frequency specifications

Verification of a transfer function

We need to show that ∀z = ejω, ω ∈ Ω ⊆ [0, π] β2 |H(z)|2 β

2

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

4/5

SLIDE 53

Frequency specifications

Verification of a transfer function

We need to show that ∀z = ejω, ω ∈ Ω ⊆ [0, π] β2 |H(z)|2 β

2

where |H(z)|2 = |b(z)|2 |a(z)|2 = b(z)b(z) a(z)a(z) = b(z)b( 1

z )

a(z)a( 1

z ) =: v(z)

w(z), v(z) and w(z) are polynomials with real coefficients.

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

4/5

SLIDE 54

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

z = ejω ∀ω ∈ Ω ⊆ [0, π]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 55

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

No need to deal with complex numbers Change of variable: t = tan ω

2

z = ejω= 1 − t2 1 + t2 + j 2t 1 + t2 z = ejω ∀ω ∈ Ω ⊆ [0, π]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 56

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

β2 v( 1−t2

1+t2 + j 2t 1+t2 )

w( 1−t2

1+t2 + j 2t 1+t2 )

β

2

z = ejω ∀ω ∈ Ω ⊆ [0, π] ↓ t = tan ω

2

∀ω ∈ Ω ⊆ [0, π] t ∈ [−∞, ∞]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 57

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

β2 r(t) + ja(t) s(t) + jb(t)

∈R since |H(z)|2

β

2

Polynomials r, s, a, b ∈ R[t] z = ejω ∀ω ∈ Ω ⊆ [0, π] ↓ t = tan ω

2

∀ω ∈ Ω ⊆ [0, π] t ∈ [−∞, ∞]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 58

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

β2 r(t) s(t) β

2

Now we work only with reals t = tan ω

2 maps ω on the whole R

Change of variable: ξ = t+2−

√ t2+4 2t

z = ejω ∀ω ∈ Ω ⊆ [0, π] ↓ t = tan ω

2

∀ω ∈ Ω ⊆ [0, π] t ∈ [−∞, ∞]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 59

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

β2 r(t) s(t) β

2

β2 r( 1−2ξ

ξ(1−ξ))

s( 1−2ξ

ξ(1−ξ))

β

2

z = ejω ∀ω ∈ Ω ⊆ [0, π] ↓ t = tan ω

2

∀ω ∈ Ω ⊆ [0, π] t ∈ [−∞, ∞] ↓ ξ = t+2−

√ t2+4 2t

ξ ∈ Ξ ⊆ [0, 1]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5

SLIDE 60

Frequency specifications

Simplifying the problem

β2 v(z) w(z) β

2

β2 r(t) s(t) β

2

β2 p(ξ) q(ξ) β

2

We compute the PGCD(p, q) with a rigorous heuristics by Char et al. z = ejω ∀ω ∈ Ω ⊆ [0, π] ↓ t = tan ω

2

∀ω ∈ Ω ⊆ [0, π] t ∈ [−∞, ∞] ↓ ξ = t+2−

√ t2+4 2t

ξ ∈ Ξ ⊆ [0, 1]

Signal processing / computer arithmetic – Metalibm workshop

T. Hilaire (thibault.hilaire@lip6.fr)

5/5