Applying the Z-transform for the stati analysis of
oating-point numeri al lters

arXiv:0706.0252v1 [cs.PL] 2 Jun 2007

David Monniaux
April 3, 2018

1 Introdu tion
The stati

analysis of

ontrol/ ommand programs, with a view to proving the

absen e of runtime errors, has re ently pi ked up steam, with the in eption
of analyzers

it is nowadays possible to build
to realisti
(sour e

sound

and

industrial situations. A stati

ode or obje t

stati

analyzers s aling up

analyzer takes as input a program

ode) and outputs a series of fa ts, warnings or other

indi ations obtained by automati
A stati

pre ise

apable of s aling up to real industrial programs. In parti ular,

analysis of that program.

analyzer is said to be sound if all the fa ts that it derives from

a program (say, variable

x

is always positive) are always true, regardless

semanti s

of how and on whi h inputs the program is run. Sound stati
based on a
exe utions.

analyzers are

, that is, a mathemati al denition of possible program

It is well-known that any method for program veri ation

annot be

at the same time sound (all results produ ed are truthful), automati
human intervention),

omplete (true results

minating (always produ es a result)

(no

an always be proved) and ter-

1 unless one supposes that program

false

memory is nite and thus that the system is available to model- he king

alarms

te hniques. As a result, sound stati
ality.

analyzers are bound to produ e

pre ise

sometimes; that is, warnings about issues that

One thus wants analyzers that are

annot happen in re-

, that is, model reality so

losely that they seldom produ e false alarms  but also, one wants analyzers that are e ient, taking only reasonable amounts of time and memory
to perform an analysis.
One

ru ial

lass of errors for

overows  say, when
using oating-point

ontrol/ ommand systems is arithmeti

onverting some value to an integer  in programs

omputations.

Su h errors have already proved to be

extremely dangerous, having for instan e
1

The formal version of this result is a

aused the explosion of the Ariane 5

lassi

theorem.

1

of re ursion theory, known as Ri e's

on its maiden ight [16℄. In order to prove the absen e of su h errors, stati
analyzers su h as Astrée

2 [1, 2℄ have to bound all oating-point variables in

the program. It is impossible to do so using simple interval arithmeti ; in
order to bound the output of a numeri al lter, one has to make the analyzer
understand the stability

onditions of the numeri al pro essing implemented

in the appli ation to be analyzed.
In

urrent

ontrol/ ommand designs, it is

utable is obtained by
by automati

ompiling C

ommonpla e that the exe-

ode, or assembly

ode, itself obtained

translation from a high-level spe i ation. This high-level spe -

i ation is typi ally given in a high-level language su h as Simulink

3 Lustre

TM ,4 These languages, in their simplest form, onsider programs

[4℄ or S ade

to be the software

ounterpart of a network of ele troni

integrators, rate limiters...)

onne ted by wires; this is a tually how several

of these languages represents programs graphi ally. Several
grouped into a

ir uits (lters,
ir uits

an be

ompound lter.

One advantage of these high-level languages is that their semanti s is
onsiderably
and
of

leaner than those of low-level languages su h as C. The lter

ompound lter

onstru tions provide natural boundaries for blo ks

omputations that belong together and probably have some interesting

ompositional

modular

and identiable properties. It is thus interesting to be able to analyze these
languages in a

and

fashion; that is, the analysis of

some blo k ( ompound lter) is done independently of that of the rest of the
ode, and the result of that analysis may be plugged in when analyzing
larger programs.

ear

This paper deals with the
lters.

ompositional and modular analysis of

By this, we mean lters that would be linear had they been

implemented over the real eld. Of

ourse, in reality, these lters are imple-

mented over oating-point numbers and none of the
relationships hold.
point

lin-

lassi al mathemati al

We nevertheless provide sound semanti s for oating-

omputations and sound analysis for su h lters.

1.1 Digital ltering
Control/ ommand programs in embedded appli ations often make use of
linear lters (for instan e, low-pass, high-pass, et .). The design prin iples
of these lters over the real numbers are well known; standard basi

designs

(Butterworth, Chebyshev, et .) and standard assembly te hniques (parallel,
serial) are taught in the
2
3

urri ulum of signal pro essing engineers.

http://www.astree.ens.fr
Simulink

TM

is a tool for modelling dynami

e.g. networks of numeri

lters. The

systems and

ontrol appli ations, using

ontrol part may then be

ompiled to hardware or

software.

http://www.mathworks. om/produ ts/simulink/
4

Ample

S ade is a

ommer ial produ t based on LUSTRE.

http://www.esterel-te hnologies. om/produ ts/s ade-suite/
2

literature has been devoted to the design of digital lters implementing some
desirable response, for implementation in sili on or in software, in xed-point
and in oating-point.[12℄
However, dis rete-time lters are often dis ussed assuming
on real numbers. There is still some

omputations

onsiderable literature on the impli a-

tions of xed-point or oating-point numbers, but the vast majority of the
work has fo used on usual

ideal

ase or average

that worst- ase bounds on

ase bounds  it is even argued

lters on real numbers are too pessimisti

and not relevant for lter design [12, 11.3℄. The study of the quantization
and roundo noise generated by xed-point or oating-point implementations has mostly been done from a sto hasti
average

point of view, in order to prove

ase properties.

For our analysis purposes, we need sound worst- ase bounds, and pra tial means for obtaining them with reasonable

omputational resour es. For

these reasons, the point of view of the designers of stati
from that of the lter designers.
A favorite tool of lter designers is the
whi h the overall ideal (i.e.

analyzers is dierent

Z-transform

[12,

hapter 3℄, with

implemented over the real numbers) transfer

fun tion of a lter is summarized in a rational fun tion with real
whose poles and zeroes determine the frequen y response.
we shall show how we

oe ients,

In this paper,

an use this transform to automati ally summarize

networks of linear lters; how this transform allows us to

ompute pre ise

bounds on the out ome of the lter, and to stati ally summarize

omplex

lters; and how to deal with roundo errors arising from oating-point

om-

putations.

1.2 Contributions of the arti le
This arti le gives a sound

abstra t semanti s

for linear numeri al lters im-

plemented in oating-point or xed-point arithmeti s, given as the sum of a
linear part (using the Z-transform) and a nonlinear part (given using ane
bounds); this latter part

omes from the roundo noise (and, possibly, some

optional losses of linear pre ision done for the sake of the speed of the analysis). (Se t. 4 for the ideal, linear part, 7 for the nonlinear part).
In many o

real

(Se t. 2.3) of
the

asions, the

omputed bounds are obtained from the norms

ertain power series. In Se t. 5, we give ee tive methods on

numbers for bounding su h norms.

In Se t.8 we explain how to

implement some of these methods e iently and soundly using integer and
oating-point arithmeti s. In Se t. 9 we study a few

ases.

As with other numeri al domains su h as those developed for Astrée, we
pro eed as follows:

the exa t oating-point

on rete semanti s is overap-

proximated by a mathemati ally simple semanti s on real numbers, whi h
is itself overapproximated by proved bounds, whi h are themselves further
overapproximated by an exe utable semanti s (implemented partly in exa t

3

×α0

+

+
normal
input

delay

×α1

delay

+

×β1

+

delay

feedback input
delay

×α2

×β2

feedback sub−filter

compositional blocks

Figure 1: De omposition of the TF2 lter

+ β1 Sn−1 + β2 Sn−2
hained by serial
tary gates are

Sn = α0 En + α1 En−1 + α2 En−2

into elementary blo ks.

omposition. Inside ea h

The

ompositional blo ks are

ompositional on the left, elemen-

omposed in parallel. On the right hand side, a feedba k loop

is used.

arithmeti s, partly using some variant of interval oating-point
tions). This ensures the soundness of the ee tive

omputa-

omputations.

This paper is an extended version of [18℄.

1.3 Introdu tion to linear lters and Z-transforms
Let us

onsider the following pie e of C

ode, whi h we will use as a running

example ( alled TF2):

Y = A0*I + A1*Ibuf[1℄ + A2*Ibuf[2℄;
O = Y + B1*Obuf[1℄ + B2*Obuf[2℄;
Ibuf[2℄=Ibuf[1℄; Ibuf[1℄=I;
Obuf[2℄=Obuf[1℄; Obuf[1℄=O;
All variables are assumed to be real numbers (we shall explain in later
se tions how to deal with xed- and oating-point values with full generality
and soundness). The program takes

time

onstant
the

oe ients. This pie e of

I as an input and outputs O; A0 et

. are

ode is wrapped inside a (rea tive) loop;

is the number of iterations of that loop. Equivalently, this lter

an

be represented by the blo k diagram in Fig. 1.
onstants and in (resp. yn , on ) the
n. Then, assuming ok = 0 for k < 0, we an
develop the re urren e: on = yn + b1 .on−1 + b2 .on−2 = yn + b1 .(yn−1 + b1 .on−2
+b2 .on−3 )+b2 .(yn−2 +b1 .on−3 +b2 .on−4 ) = yn +b1 .yn−1 +(b2 +b21 b0 ).yn−2 +. . .
where . . . depends solely on yk with k < n − 2. More generally: there
P
exist oe ients c0 , c1 . . . su h that for all n, on =
k=0 ck yn−k . These
Let us note

value of

I

a0 et .
Y, O)

(resp.

the values of the

at time

4

oe ients solely depend on the
for

bk ;

we shall see later some general formulas

omputing them.
But, itself,

oe ients

c′n

yn = a0 .in + a1 .in−1 + a2 .in−2 . It follows that
P there exist
ak and the bk ) su h that on = k=0 c′k in−k .

(depending on the

onvolution kernel

We again nd a similar shape of formula, known as a
The

c′k

sequen e is

alled a

, mapping

MI

Let us now suppose that we know a bound

|in | ≤ MI ;

inequality,

we wish to derive a bound

P

′
k=0 |ck |.MI .
′
onvolution kernel c .

|On | ≤

norm of the

MO

onvolution produ t
i

to

on the input: for all

ompositionally

omputes

n,

on the output. By the triangle

The quantity

P

′
k=0 |ck | is

alled the

What our method does is as follows: from the des ription of a
linear lter, it

.

o.

l1-

omplex

ompa t, nite representations of

onvolution kernels mapping the inputs to the outputs of the sub-blo ks of
the lter, and a
a

urately

omputes the norms of these kernels (or rather,

lose upper bound thereof ).

As a result, one

an obtain bounds on any

variable in the system from a bound on the input.

2 Linear lters
In this se tion, we give a rough outline of what we designate by linear lters
and how their basi

properties allow them to be analyzed.

2.1 Notion of lters

ausality

We deal with numeri al lters that take as inputs and output some (unbounded) dis rete streams of oating-point numbers, with
the output of the lter at time
(times 0 to

t),

t

but not on the future inputs.

mented with a state, and the output at time
time

t

; that is,

depends on the past and present inputs

5 In pra ti e, they are imple-

t

is a fun tion of the input at

and the internal state, whi h is updated.

Su h lters are typi ally

implemented as one pie e of a syn hronous rea tive loop [2, 4℄:

while(true) {
...
(state, output) = filter(state, input);
}
5

There exist non- ausal numeri al ltering te hniques One striking example is Mat-

lab's

filtfilt

fun tion, whi h runs the same

ausal lter in one dire tion, then in the

reverse-time dire tion over the same signal; the overall lter has zero phase shift at all
frequen ies, a very desirable

hara teristi

in some appli ations. Unfortunately, as seen

on this example, non- ausal lters require buering the signal and thus are not usable for
real-time appli ations. They are outside the s ope of this paper.

5

2.2 Linear lters
We are parti ular interested in lters of the following form (or
thereof ): if

(sk )

and

(ek )

the lter, there exist real
for all time

t, s t

ompounds

are respe tively the input and output streams of
oe ients

(the output at time

st =

n
X

α0 , α1 , . . . αn and β1 ,
t) is dened as:

αk et−k +

k=0

m
X

. . . βm su h that

βk st−k

(1)

k=1

or, to make apparent the state variables,




0
1
st−m+1
 ..
..
 ..   .
.
 . =
 0 ···
st
βm · · ·
If the

β

..

.

0
β2









0

···

 st−m
 .
  ..   ..
. .  + 
 0 ···
1
st−1
β1
αn · · ·

0





 et−n
  .. 
. . 
0
en
α0
.
.
.

(2)

nite impulsional response
innite impulsional response

are all null, the lter has ne essarily

(FIR) while in the opposite

ase, it may have

(IIR). The reason for this terminology is the study of the rea tion of the
system to a unit impulse (e0
FIR lter,

n + 1 time

= 1

and

∀k > 0 ek = 0).

In the

ase of a

units after the end of the impulse, the output be omes

permanently null. In the

ase of an IIR lter, the output (when

stable

omputed

ideally in the real numbers) never be omes permanently null, but rather
follows some exponential de ay if the lter is
lter may be unstable.

. A badly designed IIR

Furthermore, it is possible to design lters that

should be stable, assuming the use of real numbers in

omputation, but that

exhibit gross numeri al distortions due to the use of oating-point numbers
in the implementation.
Linear lters are generally noted using their

Z-transform

α0 + α1 z + · · · + αn z n
1 − β1 z − · · · − βm z m
The reasons for this notation will be made
all the ideal

6

(3)

lear in Se t. 4.5. In parti ular,

ompound linear lters expressible with elementary elements

su h as produ ts by

onstants, delays, et .

an be summarized by their Z-

transform (Se t. 4); that is, they are equivalent to a lter whose output is
a linear

ombination of the last

will also be

n

inputs and

m

outputs. The Z-transform

entral in the semanti s of oating-point and xed-point lters

(Se t. 7).
6

An alternate notation [12℄ repla es all o

urren es of

z

by

z −1 .

In su h a formalism,

onditions su h as the poles must have a module greater than 1 are repla ed by the
equivalent for the inverse, e.g. the poles must have a module stri tly less than 1. We
hose polynomials in

z

be ause they allow using normal power series instead of Laurent

series.

6

To summarize some salient points of the following se tions, FIR lters

α's

given by

and de ay

are very easy to deal with for our purposes, while the stability

onditions of IIR lters are determined by the study of the above

zero

rational fun tion and espe ially the module of the zeroes of the

1 − β1 z − · · · −

βm z m polynomial (z0 is a

of

Q

if

Q(z) =

Q(z0 ) = 0).

Those

roots are the inverses of the eigenvalues of the transition matrix. Spe i ally,

the lter is stable if all the zeroes have module greater than 1.

2.3 Bounding the response of the lter
The output streams of a linear lter, as an element of

RN , are linear fun

tions

of the inputs and the initial values of the state variables (internal state
variables).

onvolu-

More pre isely, we shall see later that, negle ting the oating-point errors

tion produ t Q ⋆ E

of the input

(qn )n∈N

a sequen e

E

by some

of reals su h that for any

lter is FIR if this

S

onvolution kernel
P Q

and assuming zero in the initial state variables, the output

is the

: there exists

n
k=0 qk en−k . The

n, sn =

onvolution kernel is null ex ept for the rst few values,

and IIR otherwise. If the initial state values r1 , . . . , rn are nonzero, then
S = Q0 ⋆ E + r1 Q1 + rn Qn where the Qi are onvolution kernels.
Let E : (ek )n∈N be a sequen e of real or omplex numbers. We all ∞
of E , if nite, and note kEk∞ the quantity supk∈N |ek |. Be ause of the

L -

norm

isomorphism between sequen es and formal power series, we shall likewise

note

k

P

k

ak z k k∞ = supk |ak |.

We are interested in bounding the response

of the lter with respe t to the innite norm: i.e. we want to

f

fun tion
the
has

su h that

kSk∞ ≤ f (kEk∞ ).

Said otherwise, if for all the past of

omputation sin e the last reset of the lter,

|s|

has been always less than

f (|M |)

If we do not have initialization

onstru t a

|e|

was less than

M,

then

sin e the last reset.

onditions nor oating-point errors,

f

will

be linear, otherwise it will be ane. Let us pla e ourselves for now in the
former

ase: we are trying to nd a number

innite norm

For any linear fun tion
of

f,

f

g

su h that

kSk∞ ≤ g.kEk∞.

mapping sequen es to sequen es, we

noted,

kf k∞

all

subordinate

supkxk∞ =1 kf (x)k∞ , assuming
g = kE 7→ Q ⋆ Ek∞ . If this quantity is

the quantity

is is nite. We are thus interested in

nite, the lter is stable; if it is not, it is unstable: it is possible to feed an
input sequen e to the lter, nitely bounded, whi h we result in arbitrarily
high outputs at some point in time.

L1-norm

For a sequen e (or formal series)

A, we note kAk1 =

, if nite. Then we have the following

[12, 11.3℄:
Lemma 1.

Proof.

kE 7→ Q ⋆ Ek∞ = kQk1

We shall rst prove that

sequen es

Q

and

E , kQ ⋆ Ek∞

P∞

k=0 |ak |,

alled its

ru ial and well-known result

.

kE 7→ Q ⋆ Ek∞ ≤ kQk1 ; that is, for any
≤ kQk1 .kEk∞ . Let us note C = Q ⋆ E .
7

P
cn = nk=0 qk en−k ,
kek∞ .kQk1 .

therefore

Pn

|cn | ≤

k=0 |qk ||en−k |

≤ kek∞ .

Pn

k=0 |qk |

≤

< kOk1 . Re all that kQk∞ =
P∞We shall then show equality. Let M P
N
|q
|
.
Then
there
exists
N
su
h
that
k=0 |qk | ≥ M . Choose ek = 1
k=0 k
if k ≤ N and qn−k ≥ 0, ek = −1 otherwise. Clearly, kEk∞ = 1, and
Pn
Pn
cn =
k=0 ek qn−k =
k=0 |qn−k | ≥ M , therefore kQ ⋆ Ek∞ ≥ M and
kE 7→ Q ⋆ Ek∞ ≥ M . Sin e this is valid for any M < kQk1 , then the
kE 7→ Q ⋆ Ek∞ = kQk1 equality holds.
Note that most of the dis ussion on numeri al lters found in the signal
pro essing literated is based on the L2-norm
is adapted to energy

kxk2 =


2 1/2 (whi h
k=0 |xk |

P∞

onsiderations)  for instan e, for estimating the fre-

quen y spe trum of the rounding noise. We shall never use this norm in this
arti le.

3 Convolution kernels as formal power series
In the pre eding se tion, we said that the output of the ideal lter is just the

formal power series
rational fun tions

onvolution of the input with some (possibly innite) kernel. In this se tion,
we show how

onvolution, and basi
given as

are a good framework for des ribing this

fa ts about the kernels of the lters we are interested,
.

3.1 Formal power series
We shall rst re all a few denitions and fa ts about formal power series.

K[[X]] over a eld K = R or C is the ve tor
N
spa e of ountably innite sequen es K where the produ t of two sequen es
A : (aP
k )k∈N and B : (bk )k∈N is dened as A.B : (ck )k∈N by, for all n ∈ N,
cn = nk=0 ak bn−k ( onvolution). Remark that for any algebra operation
(addition, subtra tion, multipli ation) and any N , we obtain the same results for the oe ients cn for n ≤ N as if A and B were the oe ients of
polynomials and we were omputing the oe ient cn , the n-th degree oef7
 ient of the polynomial A.B . For this reason, we shall from now on note
P∞
k
A(z) =
k=0 ak z by analogy with the polynomials. Note that for most
The algebra formal power series

of this arti le, we are interested in

formal

possible interpretation as holomorphi
all if the

onvergen e radius of the

the rare o
the needed

asions when we need

7

One

with the

P∞

A

k=0 ak z

ak

are null ex ept for a nite number,

is a polynomial.

K[[X]] as the proje tive limit of the K[X]/X n quotient
K[X]/X n+1 → K[X]/X n morphisms in the ategory of rings.

an therefore see
anoni al

k series is null); we shall note

onvergen e properties (and we shall prove

onvergen es). If all the

the formal series

power series and not with their

fun tions (i.e. it is not a problem at

8

rings

Wherever we have a
lently

onvolution

onsider a produ t

A.B

(ak ) ⋆ (bk )

inverse

of sequen es, we

an equiva-

of formal series.

We shall often wish to take the

of a power series, and the quotient

P

k
k bk b su h that b0 is
(n) as follows: A(0) = A, A(n+1) =
not null. We dene a sequen e of series A
(n)
(n)
(n)
k
qn = an /b0 . Note that for all
n ∈ N, k < n Ak = 0
A − qn ∗ z B where
P
P
n
n
(n+1) +
k
k
n
and A = A
k=0 qk z B ; thus for all n, A ≡
k=0 qk z B (mod X ),
n
whi h may equivalently written as A ≡ Q.B (mod X ). Therefore, A =

A/B

of two series.

Q.B ,

whi h explains why

This is possible for any series

A very important
Another important

Q

an be

S

and

E

quotient

ase for the rest of the paper is

A by B . P
k
1/(1 − z) = ∞
k=0 z .
of

onstatation is that this quotient formula applied to

S = E.
where

alled the

α0 + α1 z + · · · + αn z n
1 − β1 z − · · · − βm z m

(4)

are expressed as formal power series is equivalent to running

the IIR lter dened by the above rational fun tion with

S

E

the inputs and

the output.

3.2 Stability ondition
We manipulate

onvolution kernels expressed as rational fun tions where the

oe ient of degree

0

of the denominator is

1.

We shall identify a rational

fun tion with the asso iated formal power series. Using

omplex analysis, we

shall now prove the following lemma, giving the stability
to lter designers:

kQk1 < ∞

unit dis .

Lemma 2.

|z| ≤ 1

ondition familiar

if and only if all the poles of Q are outside of the

That is: a lter is stable in ideal real arithmeti s if and only if all its
poles have module greater than 1.

Proof.

Consider the poles of the rational fun tion

|z| ≤ 1

unit dis , then the radius of

the meromorphi

fun tion

Q

around 0 has a radius of

greater than 1. This implies that the series
and thus that

kQk1

onverges absolutely within the

|z| ≤ 1

If none are in the
onvergen e stri tly

z=1
kQk1 < ∞ then the series

onverges absolutely for

is nite. On the other hand, if

that dis .

Q.

onvergen e of the power series of

unit dis

and no pole

an be within

4 Compositional semanti s: real eld

exa t

Now, we have a se ond look at the basi
to give a pre ise and
the real numbers.

ompositional

semanti s of linear lters, in order
semanti s of

ompound lters on

We show that any linear lter with one input and one

output is equivalent (on the real numbers) to a lter as dened in 2.2.

9

4.1 Denition
A lter or lter element has

• ni

inputs

I1 ,

• nr

reset state values

...,

In i

( olle tively, ve tor

I ),

ea h of whi h is a

of real numbers;

r1 ,

...,

rnr

R),

( olle tively, ve tor

stream

whi h are the

initial values of the state of the internal state variables of the lter
(inside delay operators) at the beginning of the

• no
If

M

be the

output streams

O1 ,

...,

Ono

omputation;

( olle tively, ve tor

O).

is a matrix (resp. ve tor) of rational fun tions, or series, let

oordinate-wise appli ation of the norm

k·kx

Nx (M )

to ea h rational fun tion,

or series, thereby providing a ve tor (resp. matrix) of nonnegative reals. We
note

mi,j

the element in

R(z)

We note by

M

at line

i

and

olumn

j.

R and by Q[z](z) the
P (z)/(1 − zQ(z)) where P and Q are

the eld of rational fun tions over

ring of rational fun tions of the form

polynomials (that is, the ring of rational fun tions su h that the
term of the denominator is not null).

8 When

L1-norm of the asso iated power series.
When

•
•

omputed upon the real eld, a lter

a matrix
response

F

is

we note

onstant

kF k1

h that

ing) linear response of output stream

F (I, R)

i

di,j

linear

hara terizes the (de ay-

with respe t to reset value

the ve tor of output streams of lter

F

j.

over the reals, on

R.

Then we have

∀I ∈ (RN )ni ∀R ∈ Rnr F (I, R) = T F .I + D F .R

(5)

the ve tor of input streams

I

the

hara terized by:

T F ∈ Mno ,ni (Q[z](z) ) su h that ti,j hara terizes the
of output stream i with respe t to input stream j ;
D F ∈ Mno ,nr (Q[z](z) ) su

a matrix

We note

F ∈ Q[z](z) ,

and the ve tor of reset values

When the number of inputs and outputs is one, and initial values are
assumed to be zero, the

hara terization of the lter is mu h simpler 

all matri es and ve tors are s alars (reals, formal power series or rational
fun tions), and

TD

framework on this

is null. We re ommend that the reader instantiates our
ase for better initial understanding.

4.2 Basi arithmeti blo ks
Plus node implemented in oating point type

8

(z)



T = 1 1 , D = [];

This last ring is the lo alization of the ring

generated by

z,

thus the notation.

10

R[z]

f : ni = no = 1,

of real polynomials at the prime ideal

m
n

F

n

unit delay

Figure 2: A feedba k lter

S ale by

k

 
f : T = k , D = [];
 n
ks): T = z , D = 0;

node implemented in oating point type

n lo k ti
 
 
T = z ,D= 1 ;

Delay without initializer (delay for
Unit delay with initializer :

4.3 Composition
Parallel

omposition



T
0
T = 1
0 T2

omposition through lter


T = T2 .T1 , D = T2 .D1 D2 .

Serial



,



D1 0
D=
0 D2



;

1, then 2:

4.4 Feedba k loops
Let us

onsider a lter

onsisting of a lter

and feedba k loops running the
delays.

(Fig. 2) We split

TO ∈ Mn,n (Q[z](z) )

TF

n

F

m+n inputs and n outputs
n inputs through unit
es TI ∈ Mn,m (Q[z](z) ) and

with

outputs to the last

into sub-matri

representing respe tively the responses to the global

inputs and to the feedba k loop. The system then veries the linear equation
over the ve tors of formal power series:

(Idn −

zTOF )O

D F .R.

O = TIF .I + zTOF .P + D.R, and

thus

+
9
F
By Cor. 4, Idn − zTO is invertible in Mn,n (Q[z](z) ), thus T = (Idn −
F
−1
F
F
−1
F
zTO ) .TI and D = (Idn −zTO ) .D . Se tion 8.2 explains how to perform

su h

=

TIF .I

omputations in pra ti e.

4.5 Examples
A se ond order IIR linear lter is expressed by

α2 .delay2 (E) + β1 .delay1 (S) + β2 .delay2 (S).
(α0 + α1 z + α2 z 2 )E + (β1 z + β2 z 2 )S . This
S = (α0 + α1 z + α2 z 2 )(1 − β1 z − β2 z 2 )−1 .E .
9

S = α0 .E + α1 .delay2 (E) +
S =

This yields an equation

equation is easily solved into

This result is not surprising, be ause the system, by

solutions.

11

onstru tion, must admit

ausal

+

E

TF2

TF2

N1

×−k

Figure 3:

A

S

N2

delay (2)

ompound lter

onsisting of two se ond order lters and a

feedba k loop. Ea h TF2 node is a se ond-order lter whose transfer fun tion
is of the form

(α0 + α1 z + α2 z 2 )(1 − β1 z − β2 z 2 )−1 .

In Fig. 3, we rst analyze the two internal se ond order IIR lters separately and obtain

α0 + α1 z + α2 z 2
1 − β1 z − β2 z 2
a0 + a1 z + a 2 z 2
Q2 =
1 − b1 z − b 2 z 2

Q1 =

(6)

(7)
(8)

The we analyze the feedba k loop and obtain for the whole lter a rational
fun tion with a 6th degree dominator:

S=

Q1 .Q2
.E
1 + kz 2 .Q1 .Q2

where Q1 and Q2 are the transfer fun
α1 z + α2 z 2 )(1 − β1 z − β2 z 2 )−1 ), whi h

(9)

tion of the TF2 lters (form
we

(α0 +

omputed earlier.

5 Bounding the 1-norm of series expansions of rational fun tions
5.1 Inverses of produ ts of ane forms
Let

ξi be

omplex numbers of module stri tly greater than 1. Let

formal power series
The

q (n) ,

n-th

the

degree

Qn

i=1 Qi

where the

oe ient of

oe ient of

zn

qi

is

Qi (z)

(n+1)

−ξi

Q(z) be the
(z − ξi )−1 .

are the power series

, by the easy expansion:

1
−1/ξ
=
z−ξ
1 − z/ξ
in the

Q

(10)

power series, is obtained by su

essive

onvolution produ ts; it is

q (n) =

X

P
∀i,ki ∈N∧ i ki =n
12

Y

(ki )

qi

(11)

We

an therefore bound its module:

q (n) ≤

X

P
∀i,ki ∈N∧ i ki =n

Y

(ki )

qi

(12)

The right hand side of the pre eding inequality is just the

Qn

(n)
(n)
= qi
i=1 Q̇i where q̇i
1
oe ient of the
|ξi |−z series. Sin e |ξ1 | >

q̇ (n) of the series

|ξi |(n+1) is the

=
1,

the

at

z = 1,

that is,

1
|ξi |−1 . We

n-th

order

onvergen e radius of

this last series is stri tly greater than 1; furthermore, all its
nonnegative; therefore, the sum of its

oe ient

oe ients are

oe ients is the value of the fun tion

an therefore give an upper bound:

1
(z − ξ1 ) · · · (z − ξn )

1

≤

1
(|ξ1 | − 1) · · · (|ξn | − 1)

(13)

5.2 Rough and less rough approximation in the general ase
P (z)/Q(z) be a rational fun tion, with P (z) a polynomial of degree
m Q(z) a moni polynomial of degree n. Let zeroes(Q) be the multiset of
zeroes of Q (multiple zeroes are ounted with their multipli ity). P (z) =
P
P
k
k |pk |.kQk1 . Therefore
k pk z Q(z), thus kP k1 ≤
Let

P
Q

This is, however, a very

1

≤Q

kP k1
ξ∈zeroes(Q) (|ξ| − 1)

oarse approximation. Intuitively, the mass of the

onvolution kernel expressed by the

P/Q

series lies in its initial terms. Still,

with the above formula, we totally negle t the
the

(14)

an ellations that happen in

omputation of this initial part of the kernel; i.e. instead of

|a−b|, we bound it by |a|+|b|. The solution is to split kP/Qk1
≥N
and kP/Qk1 . We shall elaborate on this in Se t. 5.5.

into

onsidering

kP/Qk<N
1

5.3 Se ond degree denominators with omplex poles
A

ommon

ase for ltering appli ations is when the denominator is a se ond

degree polynomial
are two
follows:

where

onjugate

Q

of negative dis riminant. In this

omplex numbers

ξ

and

ξ̄

λ = P (ξ)/(ξ − ξ̄).

We shall for now leave
oe ients

ak = −



λ
ξ k+1
13

ak
+

Q

and the de omposition is as

λ
λ̄
P (z)
= P0 (z) +
+
Q(z)
z − ξ z − ξ̄

We are interested in the

ase, the roots of

P0

(15)

out.

of this series:

λ̄
ξ̄ k+1



(16)

Let us write

ak = −

λ = |λ|eiα

and

ξ = |ξ|eiβ ;

then


|λ|  iα −i(k+1)β
−iα i(k+1)β
e
.e
+
e
.e
|ξ|k+1
= −2

|λ|
cos (α − (k + 1)β)
|ξ|k+1

(17)

1/|ξ| mod|λ|/|ξ|. ThereIf Q is moni

To summarize, the sequen e is a de reasing exponential of rate
ulated by a sine wave and multiplied by a

onstant fa tor

omputing |λ| and |ξ| will be of prime importan e.
Q(z) = z 2 + z1 x + z0 , then |ξ|2 = ξ ξ¯ = c0 . In the ase of a rational

fore,

of the form

then
with

α0 + α1 z + α2 z 2
P (z)
=
Q(z)
1 − β1 z − β2 z 2

|ξ| = |β2 |−1/2

and

omplex numbers,

λ = P (ξ)/(ξ − ξ̄).

|λ|2 = λλ̄ =
The numerator is a symmetri
fore

are polynomials in the

P,

and

ξ̄ ,

(19)

Q,
Q ; its

roots of

oe ients of

and thereoe ients

therefore the whole polynomial

oe ients of

P

and

Q.

The denominator



α2 2 + β2 −α1 2 − α0 α1 β1 + α0 2 β2 + α2 α1 β1 + α0 β1 2 + 2 β2
−(β1 2 + 4 β2 )

we

an

Q.

We are now interested in bounding

≤ 1,

ξ

polynomial in

be expressed as a polynomial in the

|λ| =

ompute

P (ξ)(ξ̄ − ξ) + P (ξ̄)(ξ − ξ̄)
(ξ − ξ̄)2

oe ients of

is just the dis riminant of

(18)

Should we prefer not to

an be expressed as a polynomial in the

2

fun tion

|ak |.

If we just use

(20)

| cos (α − (k + 1)β) |

ome ba k to the earlier bounds obtained by totally separating the

series arising from the two poles.
We shall now obtain a better bound using the following
any real

| cos θ| =
using the

√

p
(1 + cos(2θ)) /2 ≤ 2−1/2 (1 + cos(2θ)/2)
√
inequality
1 + x ≤ 1 + x/2. Therefore

cos2 θ =

on avity

|ak | ≤

√

Now, we are interested

r<1

onstatation: for

θ,

∞
X
k=0

|λ|
(1 + cos(2(α − (k + 1)β))/2)
|ξ|k+1
P∞
in bounding
k=N [ak |. For any a and b,

(21)

2

cos(a + kb)r k =

cos a − r cos(a − b)
1 − 2r cos b + r 2

14

(22)

and

0≤
(23)

Let us now see the quality of su h bounds

S1 =

∞
X
k=0

S2 =

S1 ≤ S2 , S1 ≤ S3 :

| cos(a + kb)|r k
∞
X

(24)

1
1−r

rk =

k=0

∞

(25)

1 X
(1 + 2 cos(2(a + kb)))r k
S3 = √
2 k=0


1
1 cos(2a) − r cos(2(a − b))
1
+ .
S3 = √
1 − 2r cos b + r 2
2 1−r 2

(26)

(27)

Note that S3 is not ne essarily better than S2 (for a = 0 and b = 0, S3 /S2 =
√
3/(2 2) ≃ 1.06). However, some moderate gains may be obtained ; for
instan e, for r = 0.7, a = 0 and b = 0.3, S1 ≃ 2.60, S2 ≃ 3.33 and S3 ≃ 2.80.
For pra ti al purposes, the bound obtained using S2 is very su ient and
easy to

ompute. We thus opt for this one.

5.4 Finer bounds using partial fra tion de omposition
Q is their
P prime
with Q the fra tion P/Q admits
P
a
as P/Q = P0 +
i Pi /Qi , where P0 is the Eu lidean
quotient of P by Q and the degree of Pi is stri tly less than that of Qi .
Using the fundamental theorem of algebra, it follows that if the ξi are
the distin t roots of Q and mi their multipli ity, then there exist λi,j ∈ C

It is well known that if

Qi

are pairwise prime polynomials, and

partial de omposition

produ t, then for any polynomial

su h that

P/Q = P0 +

mi
XX
i

Sin e

Q

j=1

λi,j
(z − ξi )j

(28)

ξi and
j , λi′ ,j = λ̄i′ , j .

is a real polynomial, its roots are either real, either pairs of

onjugate

ξi′ = ξ̄ , with

the same multipli ity, and also for all

However, while theoreti ally sound, this result is numeri ally deli ate
when there are multiple roots, or dierent roots very
1.3℄ For instan e, let us
degree polynomial

Q,

lose to ea h other.[13,

onsider a rst-degree polynomial

and a se ond-

then

λ1
λ2
P (z)
=
+
Q(z)
z − ξ1 z − ξ2
and we obtain

P

λ1 = P (ξ1 )/(ξ2 − ξ1 )

(and

(29)

λ2 = P (ξ2 )/(ξ1 − ξ2 )). Both
ξ1 − ξ2 . While it is

numbers will get very large, in inverse proportion of

quite improbable that we should analyze lters where two separate poles
have been intentionally be pla ed very

lose together, it is possible that we

analyze lters with multiple poles (for instan e, the

15

omposition of a lter

with itself ), and, with numeri al
lose poles and thus a dramati
We still

an pro eed with a

instead of fa toring

Q

omputations, we would have two extremely

radius r de omposition
numeri al instability.

into a produ t of

z − ξi

of

P/Q

[13, Def 1.3℄:

fa tors, we fa tor it into

a produ t of Qi su h that for any i, and any roots ξ1 and ξ2 of Qi ,
|ξ1 − ξ2 | < 2r . The same referen e des ribes algorithms for performing

then
su h

de ompositions. We obtain a de omposition of the form

X
P
= P0 +
Pi /Qi
Q

(30)

i

where the roots of ea h
the degree of

Pi .

Qi

P
Q
whi h we
We

are

lose together, the degree of

Pi

is less than

From this we obtain the bound

1

≤ kP0 k1 +

X
i

kPi k.

1
Qi

(31)

1

an bound using the inequalities given in the pre eding subse tions.

an, as before, improve on this bound by splitting the series between an

initial sequen e and a tail.

5.5 Development of rational fun tions and normed bounds
P (z)/Q(z) ∈ Q[z](z) be a rational fun tion representing a power series
by its development (un )n∈N around 0. We wish to bound kuk1 , whi h we
shall note kP/Qk1 . As we said before, most of the mass of the development
of P/Q lies in its initial terms, whereas the tail of the series is negligible
(but must be a ounted for for reasons of soundness). We thus split P/Q
into an initial development of N terms and a tail, and use
Let

+ kP/Qk≥N
kP/Qk1 = kP/Qk<N
1
1
kP/Qk1

is

omputed by

velopment of

P/Q.

performing su h a

omputing expli itly the

N

(32)

rst terms of the de-

We shall see in Se t. 8.3 the di ulties involved in
omputation soundly using interval arithmeti s.

Let dQ be the degree of Q. The development D of P/Q yields an equation
P (z) = D(z).Q(z) + R(z).z N . We have P (z)/Q(z) = D(z) + R(z)/Q(z).z N
and thus

kP/Qk≥N
= kR/Qk1 ≤ kRk∞ .k1/Qk1
1

(33)

The pre eding sub-se tions give a variety of methods for bounding
using the zeroes of

Q(z);

bounds on the absolute values of the zeroes of
expli it
ompute

omputation of

D

k1/Qk1

Se tion 5.2 gives a rough method based on lower

R

Q(z). kRk∞

is bounded by

using interval arithmeti s; as we shall see, we

until the sign of the terms is unknown  that is, when the norm

of the developed signal is on the same order of magnitude as the numeri al

16

error on it, whi h happens, experimentally, when the terms are small in
absolute values.

Therefore,

approximation used

k1/Qk1

kRk∞

is small, and thus the roughness of the

does not matter mu h in pra ti e.

6 Pre ision properties of xed- or oating-point operations
In this se tion, we shall re all a few fa ts on the errors introdu ed by xedand oating-point arithmeti s. They will be su ient for all our reasonings,
without need for further knowledge about numeri al arithmeti s.
Most types of numeri al arithmeti s, in luding the widely used IEEE-754
oating-point arithmeti , implemented in hardware in all

urrent mi ro om-

puters, dene the result of elementary operations as follows: if

f

is the ideal

operation (addition, subtra tion, multipli ation, division et .) over the real
numbers and
where
that

r

is a

f˜ is

roundo
the

fun tion. The roundo fun tion

and is very

•
•
•

hooses a

an be exa tly represented in the used xed- or oating-point data type,
lose to

x; spe

i ally, most systems, in luding all IEEE-754 sys-

tems, provide the following roundo fun tions:

•

f˜ = r ◦ f
value r(x)

orresponding oating-point operation, then

round to
of

x

in the dire tion

+∞: r(x) is the representable real nearest to x in the dire

tion

−∞: r(x) is the representable real nearest to x in the dire

tion

+∞;

round to
of

is the representable real nearest to

0;

round to
of

0: r(x)

10

−∞;

round to nearest (generally, the default mode):
sentable real nearest to

r(x)

is the repre-

x.

not-a-number

In this des ription, we leave out the possible generation of spe ial values
su h as innities (+∞ and

−∞)

ating undened results su h as
10

and

0/0.

(NaN), the latter indi-

We assume as a pre ondition to the

On Intel x86 systems, the des ription of the exa t properties of the oating-point

arithmeti s is

ompli ated by the fa t that, by default, with most operating systems and

languages, the 80287- ompatible oating-point unit performs
ing 80-bit long double pre ision numbers, even when the

omputations internally us-

ompiled program suggests the

use of standard 64-bit double pre ision IEEE numbers. Note that su h usage of supplemental pre ision for intermediate
The nal result of the

omputations is allowed by the C standard, for example.

omputation may therefore depend on the register s heduling and

optimizations performed by the

ompiler. Sin e we reason by

are always sound (albeit pessimisti ) in the fa e of su h
piler and the system do.

17

maximal errors, our bounds

ompli ations, whatever the

om-

numeri al lters that we analyze that they are not fed innities or NaNs  indeed, in some DSP (digital signal pro essor) implementations, the hardware
is in apable of generating or using su h values, and in many other implementations the system is

ongured so that the generation or usage of innities

any

issues an ex eption resulting in bringing the system into a failure mode. Our
framework provides
quantity

x inside

onstru tive methods for bounding

the lters as

kxk∞ ≤ c0 +

Pn

k=1 ck .kek k∞

the input streams of the system; it is quite easy to
does not overow (kxk
su ient

< M );

one

ek

are

he k that the system

an even easily provide some very wide

onditions on the input (kek k∞

not in lude su h

oating-point

where the

P
≤ (M − c0 )/( nk=1 ck )).

We will

onditions in our des ription, for the sake of simpli ity.

x
and the oating-point result x̃ is bounded, in absolute value, by max(εrel |x|, εabs )
11
(the least positive oating-point number)
where εabs is the
−1074
−324
≃ 4.94 · 10
and
and εrel is the
in urred (εabs = 2
−53
−16
εrel = 2
≃ 1.11 · 10
for IEEE double pre ision operations, for the
For any arithmeti

operation, the dis repan y between the ideal result

absolute error
relative error

worst

ase with respe t to rounding modes). We a tually take the

oarser

inequality

|x − x̃| ≤ εrel |x| + εabs

(34)

See [9℄ for more details on oating-point numbers and [17℄ for more about
the ane bound on the error.
In the

ase of xed-point arithmeti s, we have

εrel = 0

and

εabs = δ (δ

is

the smallest positive xed-point number) if the rounding mode is unknown
(round to

+∞, −∞

et .) and

δ/2

is it is the rounding mode is known to be

round-to-nearest.

7 Compositional semanti s: xed- and oating-point
In this se tion, we give and a

ompositional

the oating-point numbers.

abstra t

semanti s of lters on

7.1 Constraint on the errors
Our abstra t semanti s

•
•
11

hara terizes a xed- or oating-point lter

the exa t semanti s of the asso iated lter
an abstra tion of the dis repan y
ideal and oating-point lters.

The absolute error results from the

F

ondition:

an ex eption), underow is generally a benign

ondition.

18

between the

a number

rounded to 0. Contrary to overow (whi h generates innities, or is
relying on relative error bounds if one wants to be sound.

by:

over the real numbers

∆(I) = F̃ (I) − F (I)

underow

F̃

lose to 0 is

ongured to issue

However, it pre ludes merely

We transform
a

F̃ (I)

into the sum of a term that we

urately using algebra and

∆(I)

that we

an bound very

omplex analysis, and a nondeterministi

annot analyze a

urately and soundly without

di ulties, but for whi h bounds are available:

ane almost linear

onsiderable

assuming for the sake of

simpli ity a single input and a single output and no initialization
we obtain an
to

,

onstraint on the

F .
kIk∞ : k∆(I)k∞ ≤ εFrel kIk∞ + εabs

input

k∆(I)k∞

onditions,

with respe t

In short: sin e the lter is linear, the

magnitude of the error is (almost) linear.
We generalize this idea to the
abstra t semanti s
and

ase of multiple inputs and outputs. The

F
∆ is given by matri es εrel
,T ∈ Mno ,ni (R+ )
no
F
ve tor εabs ∈ R+ su h that

hara terizing

F
εrel
,D ∈ Mno ,nr (R+ )

and a

F
kF (I, R) − F̃ (I, R)k∞ ≤ εFrel,T .N∞ (I) + εrel
,D .N∞ (R) + εabs .
where

F̃ (I, R)

is the output on the stream

numbers on input streams

I

omputed upon the

and initial values

I.

(35)

oating-point

7.2 Basi arithmeti blo ks
Plus node implemented in oating point type



f : ni = no = 1, T = 1 1 ,

i
h
f
f
D = 0, εrel,T = εfrel εrel
, εrel,D = 0, εabs = εabs ;
k node implemented in oating point
f
f
, εrel,D = 0; εabs = εabs ;
εrel,T = |k|.εrel

S ale by

Delay without initializer (delay for

n

 
f : T = k , D = 0,

type

lo k ti ks):

εrel,T = 0, εrel,D = 0, εabs = 0
Unit delay with initializer :

εabs = 0

 
T = z n , D = 0,

 
 
T = z , D = 1 , εrel,T = 0, εrel,D = 0,

omposition blo k matri es and ve tors:

Parallel

εrel,T



ε1
= rel,T
0

0
2
εrel
,T



,

εrel,D



ε1
0
= rel,D 2
0
εrel,D



,

εabs



ε1
= 2abs
εabs



.

7.3 Serial omposition
The serial

omposition of two lters is more involved.

the ideal linear transfer fun tions of both lters, and
fun tions implemented over oating-point numbers.
We have
for

G).

We

F̃

F and G be
G̃ the transfer

Let
and

mutatis mutandis

∀I N∞ (F (I) − F̃ (I)) ≤ εFrel .N∞ (I) + εFabs (
are interested in ε = N∞ (F (I) − F̃ (I)): that

is, a ve tor of

positive numbers indexed by the outputs of the system su h that on every
oordinate

k, the dieren

e

δ

between output

19

k

omputed over the reals and

the oating-point numbers over the same input
extend

≤

to real ve tors

I

veries

oordinate-wise.

kδk∞ ≤ εk .

We

The following is easier to understand when ea h lter has a single input
and a single output; then, all ve tors and matri es are s alars (either in
or

R

Nx (v) is simply kvkx .
R of (re)initialization values is split between RF (those onF ) and RG (those on erning G). We split the overall output error

Q[z](z) ,

and

The ve tor
erning

of the system between the part that was introdu ed by the rst lter (and
then amplied or attenuated by the se ond lter) and the part that was
introdu ed by the se ond lter, and use the triangle inequality:

N∞ ((G ◦ F )(I, R) − (G̃ ◦ F̃ )(I, R))
≤ N∞ (G ◦ F (I) − G ◦ F̃ (I)) + N∞ (G ◦ F̃ (I) − G̃ ◦ F̃ (I))
G
G
G
G
≤ N1 (G).(F (I) − F̃ (I)) + εrel
,T .N∞ (F̃ (I)) + εrel,D .N∞ (R ) + εabs
G
G
G
≤ N1 (G).(F (I) − F̃ (I)) + εG
rel,T .(N∞ (F̃ (I)) + N∞ (F̃ (I) − F (I))) + εrel,D .N∞ (R ) + εabs
G
G
G
G
≤ (N1 (G) + εG
rel,T ).N∞ (F̃ (I) − F (I)) + εrel,D .N∞ (R ) + εrel .N∞ (F (I)) + εabs
F
F
F
F
≤ (N1 (G) + εG
rel,T ).(εrel,T .N∞ (I) + εrel,D .N∞ (R ) + εabs )
G
G
G
G
+ εrel
,T .N1 (F ).N∞ (I) + εrel,D .N∞ (R ) + εabs
i
h
F
G
G
≤ (N1 (G) + εrel
,T ).εrel,T + εrel,T .N1 (F ) .N∞ (I)
h
i
h
i
F
F
G
G
+ (N1 (G) + εG
rel ).εrel,D .N∞ (R ) + εrel,D .N∞ (R )
i
h
F
G
+ (N1 (G) + εG
rel ).εabs + εabs

Thus

G◦F
εrel
,D

(36)

G◦F = (N (G) + εG ).εF + εG .N (F ),
εrel
1
1
rel
rel
rel

 ,T
G
F
G◦F = (N (G)+ εG ).εF + εG .
ε
).ε
,
and
εabs
= (N1 (G) + εG
1
rel
rel,D
rel,D
rel
abs
abs

7.4 Feedba k loops
Let us

all

o(n)

the ve tor of outputs of the lter at step

linear fun tion of the

n.

It is, ideally, a

urrent input, the pre eding inputs, and the pre eding

On = L(I≤n , O<N ). Let us all L̃ the asso iated oating-point
tion and Õ the oating-point output of the lter. Let us all ∆ = Õ − O .

outputs.
fun

∆n = L̃(I≤n , Õ<N ) − L(I≤n , O<N )

= L̃(I≤n , Õ<N ) − L(I≤n , Õ<N )) + L(I≤n , Õ<N ) − L(I≤n , O<N )


= L̃(I≤n , Õ<N ) − L(I≤n , Õ<N )) + L(0, ∆<N )

(37)

Cn = L̃(I≤n , Õ<N ) − L(I≤n , Õ<N )) be the sequen e of ve tors of error
F
reations at ea h iteration. Then ∆ veries the equation ∆ = C + zTO .∆.
F
−1
As before, this means ∆ = (Idn − zTO )
.C and thus that N∞ (∆≤n ) ≤

F
−1
.N∞ (C≤n ).
N1 (Idn − zTO )
F
F
F
Let us split εrel,T ∈ Mn,n+m (R+ ) into εrel,I ∈ Mn,m (R+ ) and εrel,O ∈
Let

20

Mn,n (R+ ).

Then

F
F
F
F
N∞ (C≤n ) ≤ εrel
,I .N∞ (I≤N ) + εrel,O .N∞ (Õ<N ) + εrel,D .N∞ (R) + εabs
F
F
≤ εFrel,I .N∞ (I≤N ) + εrel
−O )
,O .N∞ (O<N ) + εrel,O .N∞ (Õ
| N {z <N}
∆<N

F

F
+ εrel,D .N∞ (R) + εabs

But then, noting

A = N1 (Idn − zTOF )−1



(38)

,

F
N∞ (∆≤n ) ≤ A.(εFrel,I .N∞ (I≤N ) + εrel
,O .N1 (T ).N∞ (I≤N )

F
F
F
+ εrel
,O .N∞ (∆<N ) + εrel,D .N∞ (R) + εabs )

Let

F
K1 = A.εrel
,O ∈ Mn,n (R+ )

(39)

and

F
F
K2 (ι, ρ) = A. εFrel,I + εFrel,O .N1 (T )).ι + εrel
,D .ρ + εabs

(40)

N∞ (∆≤n ) ≤ K1 .N∞ (∆<n )+K2 (N∞ (I), N∞ (R)). This means that the
sequen e un = N∞ (∆<n ) veries u0 = 0 and un+1 ≤ K1 .un +K2 (N∞ (I), N∞ (R)).
This implies that for all n, un is less than the least xed point L of v 7→
K1 .v + K2 (N∞ (I), N∞ (R)).
Re all that the spe tral radius of a matrix M of real numbers is the
greatest absolute values of its eigenvalues. If K1 is ontra ting (spe tral
radius less than 1), then v 7→ K1 .v + K2 (N∞ (I), N∞ (R)) has a unique
xed point, by Bana h's xed point theorem; and 1 − K1 is invertible.
−1 K (N (I), N (R)). Let ε
This xed point is v = (1 − K1 )
rel,T = (1 −
2
∞
∞


Then

F
F
K1 )−1 .A. εrel
,I + εrel,O .N1 (T )

F .
K1 )−1 .A.εabs

,

εrel,D = (1 − K1 )−1 .εFrel,D ,

εabs = (1 −

and

N∞ (∆) ≤ εrel,T .N∞ (I) + εrel,D .N∞ (R) + εabs .
K1 = A.εFrel,O ∈ Mn,n (R+ ) where A is the matrix of norms

N1 (Idn − zTOF )−1 ; K1 bounds the amount of oating-point impre ision
that feeds ba k into the system. A is the ampli ation bounding matrix
Then

Re all that

of the lter

onsisting merely of the feedba k loop of the original lter; if

the original lter is stable and well-designed, the

F
moderate. εrel,O measures the

oe ients of

A

reation of impre ision in one iteration of the

internal lter; if the lter is numeri ally well-designed, then its
are very small. On real-world examples,
of

should be

10−15 .

K1

oe ients

was on the order of magnitude

This suggests an ee tive method for bounding from above the various
quantities of the form
(if

y

(1−K1 )−1 .y

that we listed, where

is a matrix, then split it into its

y

is a

olumn ve tor

olumn ve tors).

d∞ = (1 − K1 )−1 .y =
21

∞
X
k=0

K1k .y

(41)

is the unique xpoint of

φ = x 7→ K1 .x + y ,

whi h is monotoni

tra ting. Consider the matrix norm subordinate to

kK1 k = sup
i

This gives a rough bound on

kd∞ k∞ ≤

and

on-

on ve tors:

k1i,j

(42)

j

d∞ :

∞
X
k=0

kK1 kk .kyk∞ =
Pn

dn = (x 7→ K1 .x + y)n (y) =

Let

X

k · k∞

n
k=0 K1 .y .

kd∞ − dn k∞ ≤

kyk∞
.
1 − kK1 k

(43)

d∞ − dn = K1n+1 .d∞ ,

thus

K1n+1
.kyk∞ .
1 − kK1 k

(44)

d∞ :

.kyk∞ .V1

Therefore, the following is an upper bound on



B = dn + −
where

V1

K1n+1
kK1 k − 1

is a ve tor of ones of the same dimension as

(45)

y.

This

may be ee tively performed in oating-point arithmeti

omputation

in order to yield

a sound upper bound by

omputing Eqn. 42 and 45 in round-to-+∞ mode

(x

Remark that we

of

7→ −1/x is monotoni ).
the resulting B̃ by he

king that

an dire tly prove the soundness

K1 .B̃ + y

is less than

B̃

(this

he k-

ing phase, though unne essary assuming a sound implementation, may be
heaply performed for the sake of se urity; while it is possible that the result
should be
and

orre t and the

he k fails, this seems very unlikely in pra ti e,

an be worked around by

hoosing a slightly larger

B̃ ).

7.5 Trading some a ura y for omputation speed; nonlinear
elements
We have split the behavior of the lter into the sum of the

onvolution of

the input signal by the power development of a rational fun tion, representing the exa t behavior, and some error term.
fun tions exa tly over
pensively large.

Q[z](z) ,

then the rational

If we

ompute the rational

oe ients might grow ex-

It seems silly to use high pre ision for the

oe ients of

a system parameterized by oating-point numbers and implemented with
oating-point errors. Indeed, we may redu e the pre ision of the

oe ients

of the rational fun tion at the expense of adding to the margin of error.
An ideal lter of Z-transform the rational fun tion

P (z) =

Pdp

k
k=0 pk z and

Q(z) =

Pdq

is equivalent to a lter with ideal input Z-transform
Z-transform

Q

(Fig 4).

P (z)/(1−Q(z)) where

k
k=1 qk z with non initialization

P

ondition

and ideal feedba k

Su h a lter may be soundly approximated by a

22

P

+

delay

Q

Figure 4: An ideal lter equivalent to a lter of Z-transform

♯

P (z)/(1−Q(z)).

♯

F ♯ with TIF = P ♯ , TOF = Q♯ , εrel,I = kP ♯ − P k1 ,
♯
εrel,I = kQ − Qk1 , εabs = 0, whi h we know how to solve from Se t. 7.4.
♯
More generally: a lter F may be approximated by a lter F with
♯
♯
♯
F
G
F
G
G
F
F
F
transfer fun tion T
= T , εrel,T = εrel,T + εrel,T , εrel,D = εrel,D + εrel
,D ,

non-ideal feedba k lter

♯

F = εF where G is the feedba k lter with internal lter H
εabs
abs
H
H
H
♯
♯
P ♯ , TOH = Q♯ , εrel
,I = kP − P k1 , εrel,I = kQ − Qk1 , εabs = 0.
Note that this gives a generi

ements o

given

TIH =

method for approximating non-linear el-

uring in lters, provided that it is possible to split them into a

linear part and a nonlinear part, the output of whi h

an be bounded by an

ane fun tion of bounds on the absolute value of the inputs.

8 Numeri al onsiderations
We have so far given many mathemati al formulas that are exa t in the

real

eld. In this se tion, we explain how to obtain sound abstra tions for these
formulas using oating-point arithmeti s.

8.1 Interval arithmeti s

upward rounded rounded-

IEEE oating-point arithmeti s [9℄ and good extended pre ision libraries

to-+∞

downward rounded rounded-to-˜−∞

su h as MPFR [7℄ provide fun tions

omputing

(or

f (x1 , . . . , xn )
f˜+ are the asso iated
oating-point downward and upward operations, then f (x1 , . . . , xn ) is guar−
+
anteed to be in the interval [f˜ (x1 , . . . , xn ), f˜ (x1 , . . . , xn )], whi h will guar) and

(or

is the exa t operation on real numbers and

antee the

soundness˜

f˜− (x1 , . . . , xn )

better bounds
will guarantee

of our approa h.

f−

) results: that is, if

and

Furthermore, for many operations,

+
and f (x1 , . . . , xn ) are guaranteed to be optimal; that is, no

optimality

an be provided within the desired oating-point format; this
of

ertain of our elementary operations.

23

8.2 Approximate algebrai

omputations

In many o

P asions, we ideally would like to ompute on real polynomiPP
= nk=1 pk z k but instead we ompute on oating-point polynomials
n
k
P̃ =
k=1 [lk , hk ]z abstra ting the set γ(P̃ ) of polynomials P su h that
∀k pk ∈ [lk , hk ]. In pra ti e, it will often be ne essary that 0 ∈
/ [lk , hk ]
als

in order to avoid un ertainties on the degree of the polynomial.

All the

usual polynomial operations (addition, multipli ation by a s alar, subtra tion, multipli ation) may be abstra ted using interval arithmeti s. We also
in lude a test

ontains0 (P̃ ) whether the null polynomial is in

this stru ture an

abstra t ring

Given a abstra t ring

R,

.

we

onstru t the

abstra t eld

γ(P̃ ).

We

all

of fra tions over

p1 /q1 +p2 /q2 = (p1 q2 +p2 q1 )/(q1 q2 ),
k.(p/q) = (k.p)/q , (p1 /q1 ).(p2 /q2 ) = (p1 .p2 )/(q1 .q2 ), (p1 /q1 )/(p2 /q2 ) = (p1 .q2 )
/(q1 .p2 ), ontains0 (p/q) = ontains0 (p). We an make a simple attempt at
that ring using the following operations:

redu ing the fra tions by

he king that there are no trivial

an ellations

abstra t ring of matri es

between the numerator and denominator in produ ts and quotients.

K,

Given an abstra t ring

we

onstru t the

over that ring with the usual operations: if
if

M = A.B , mi,j =

P

k

ai,k .bk,j .

If

K

is an abstra t eld, we

implement Gaussian elimination in order to
matrix

A

and a matrix

su h that

B.

M = A + B , mi,j = ai,j + bi,j ;
ompute

A−1 .B

an also

given a square

When we look for a pivot, we sele t elements

e

ontains0 (e) is false.

Unfortunately,

omputations on su h approximate stru tures may yield

unfavorable results. In parti ular, the absen e of simpli ation between the
numerator and denominator may yield fra tions
have some

P̃ (z)/Q̃(z)

where

P̃

and

Q̃

ommon zeroes. The spurious poles that are introdu ed not be

that mu h of a problem if we use partial fra tion de omposition (Se t. 5.4),
for they will yield very small
will make the

oe ients in the de omposition; however, they

omputations more

of Se t. 5.1, the results may be
A solution is to perform all
over

Q[z](z) .

Then,

omplex. If using the simple tail bounds

onsiderably worse.
omputations on rational fun tions exa tly

an ellation between a numerator and a denominator

be performed exa tly by division by their greatest

an

ommon divisor, whi h is

obtained from Eu lide's algorithm over the Eu lidean division of polynomials. No spurious poles may be introdu ed. However, on large lter networks,
exa t

omputations may produ e ex eedingly large integer numerators and

denominators.

It is then possible to apply the approximation s heme of

Se t. 7.5 in order to trade speed for potential pre ision. This is the solution
that we implemented in our system: exa t

omputations on rational num-

bers and safe approximations to limit the length of the numbers involved in
the

omputations.

24

8.3 Computation of developments
kP/Qk1 of a series quotient of two polynomials,
N initial terms of development, whi h we ompute

When bounding the norm
we split the series into its

expli itly, and a tail whose norm we bound. The rst idea is to

N

ompute the

rst terms of the series by quotienting the series, as explained in Se t. 3.1

or, equivalently, by running the lter for

1, 0, 0, . . . .

N

iterations on the Dira

input

In order to provide a sound result, one would work using interval

arithmeti s over oating-point numbers. However, as already noted by Feret,
after some number of iterations the sign of the terms be omes unknown and
then the magnitude of the terms in rease fast; it is therefore indi ated to
ompute the development until the rst term of unknown sign is rea hed,

N a ordingly (one may still also enfor e a maximal number of
Nmax ). In order to be able to develop the quotient further with

and assign
iterations

good pre ision, one

an use a library of extended-pre ision oating-point

omputations with sele table rounding dire tion, su h as the MPFR library
now part of GNU MP [7℄.

8.4 Bounding the roots
In order to bound
interval

kP/Qk1 ,

where

P

Q

and

may possibly be given using

oe ients, we have to bound the roots of

Q.

More formally, we

P̃ (z) =
k su h that 0 ∈
/
[l
,
h
]
,
nd
a
family
(
ξ̃
,
ρ
)
(
[l
;
h
]z
n n
k k 1≤k≤n ξk ∈ C
k=1 k k
with ℜξk and ℑξk oating-point numbers, ρk ∈ R+ a oating-point number)
Pn
su h that for any polynomial P =
k=1 pk su h that ∀k pk ∈ [lk , hk ], then,
up to a permutation, the n roots (ξk )1≤k≤n of P are su h that ξk ∈ D(ξ̃k , ρk )
where D(z, r) is the losed dis of enter z and radius r .
Often, what we need is a tually bounds on the |ξk |; this an easily be

have to solve the following problem: given an interval polynomial

Pn

obtained from the pre eding bounds using interval arithmeti

on plus, minus,

multiply and square root.
Our

oe ients are intervals

errors of oating-point

hk − lk

[lk , hk ]

in order to a

omputations. As a

ommodate possible

onsequen e, it is expe ted that

are small. This suggests to us a two-step method for obtaining the

desired bounds:

approximations x

1. Use an e ient and, in pra ti e, very a

j to the roots of

nomial).

Pn

urate algorithm to obtain

lk +hk k
k=1
2 z (the midpoint poly-

2. From those approximations, obtain bounds on the radius of the error
ommitted.
There exist a variety of methods and implementations to perform the rst
point. We used

gsl_poly_ omplex_solve of the GNU S

whi h is based on an eigenvalue de omposition of the

25

ienti

Library [8℄,

ompanion matrix.

For the se ond step, Rump des ribes a variety of bounding methods [22℄
whi h take a polynomial and approximate roots as an input and output
error radii; these methods may be performed using interval arithmeti s. We

ξj

implemented the simplest and roughest one:

xj − p j

and radius

|pj |

is in a

losed dis

of

enter

where

pj =

nP (xj )
Q
,
pn k6=j xj − xk

(46)

whi h is easily implemented using interval arithmeti s (P be omes

P̃

et .).

9 Implementation and ase studies
We implemented the algorithms des ribed above in a simple Obje tive Caml
[15℄ program: lters are represented by a re ord of all their

hara teristi s

(transfer matri es, bounding matri es); fun tions (in the OCaml) sense
stru t lter re ords, or perform
The formal

omposition operations.

omputations on fra tions are performed over

using GNU MP's

mpq

on-

type [7℄. We initially

Q, implemented

onsidered using MPFR [10℄, an

extended pre ision library with sound rounding modes, for interval

omputa-

tions; instead, we simply use the IEEE-754 rounding modes of the hardware
oating-point unit, whi h is mu h faster.

9.1 Composition of TF2 lters
Let us re all the example of Se t. 4.5. It is a

omposition of two TF2 lters

with a feedba k loop around it. The serial

omposition of the lter in Fig. 3

and another TF2 lter, all with realisti

oe ients, is analyzed in about

0.10 s on a re ent PC; the analyzer nds that

εrel ≃

10−12 and

εabs ≃

10−305 .

kSk ≤ gkEk

with

g ≃ 2,

with

The power series developments of rational fun tions (Se t. 8.3) are done

up to around the 27th order.

9.2 Complex nonlinear iterated lter
We now

onsider a nonlinear, iterated lter due to Roozbehani et al. [21℄[5℄.

We rst analyze separately

filter1() (2nd-order linear lter) and filter2()

(2nd-order ane lter). So as to simplify matters, we do not give the transfer fun tions using matri es, matri es inverses et .

but as the solution of

z.

We obtain that system

a system of linear equations over polynomials in

very simply from the program: whenever we see an assignment
turn it into an equation

x=e

x := e,

we

(we assume without loss of generalities that

variables are only assigned on e in a single iteration step), where
original expression where a variable

v

e

is the

that has not yet been assigned in the

26

urrent iteration is repla ed by
value of

iv + z.v , iv

standing for the initialization

v.

void filter1 () {
stati float E[2℄, S[2℄;
if (INIT1) {
S[0℄ = X; P = X;
E[0℄ = X; E[1℄=0; S[1℄=0;
} else {
P =0.5*X-0.7*E[0℄ +0.4*E[1℄ p = 0.5e − 0.7(ie0 + z.e0 )
+1.5*S[0℄-S[1℄*0.7;
+0.4(ie1 + z.e1 )

}

}

E[1℄ = E[0℄;
E[0℄ = X;
S[1℄ = S[0℄;
S[0℄ = P;
X=P/6+S[1℄/5;

+1.5(is0 + z.s0 ) − 0.7(is1 + z.s1 )
e1 = ie0 + z.e0
e0 = e
s1 = is1 + z.e1
s0 = p
x = p/6 + s1 /5

We all e the input value for X. We solve the system and obtain x =
Q.e + Qie0 .ie0 + Qie1 .ie1 + Qis0 .is0 + Qis1 .is1 . The ommon denominator of
2
the Q fra tions is 10 − 15z + 7z , whi h has omplex onjugate roots z su h
that |z| ≃ 1.2. ie1 = is1 = 0 and ie0 = is0 = ι (the last value for input
e su h that INIT1 is true), thus kxk∞ ≤ kQk1 .kek∞ + kQie0 + Qis0 k∞ .kιk.
With a pre ondition kek∞ ≤ 400, this yields kxk∞ < 339. If we take the
oarser inequality kxk∞ ≤ kQk1 .kek∞ + (kQie k∞ + kQis k∞ ).kιk we get
0
0
kxk∞ < 528. Roozbehani et al. nd a bound ≃ 531.

void filter2 () {
stati float E2[2℄, S2[2℄;
if (INIT2) {
S2[0℄ =0.5*X; P = X;
E2[0℄ = 0.8*X; E2[1℄=0; S2[1℄=0;
} else {
P =0.3*X-E2[0℄*0.2+E2[1℄*1.4 p = 0.3e − 0.2(ie0 + z.e0 )
+S2[0℄*0.5-S2[1℄*1.7;
+1.4(ie1 + z.e1 )

}

}

E2[1℄ = 0.5*E2[0℄;
E2[0℄ = 2*X;
S2[1℄ = S2[0℄+10;
S2[0℄ = P/2+S2[1℄/3;
X=P/8+S2[1℄/10;

+0.5(is0 + z.s0 ) + 1.7(is1 + z.s1 )
e1 = 0.5(ie0 + z.e0 )
e0 = 2e
s1 = is0 + z.s0 + τ
s0 = p/2 + s1 /3
x = p/8 + s1 /10

τ = 10/(1 − z) and
x = Q.e + Qie0 .ie0 + Qie1 .ie1 + Qis0 .is0 + Qis1 .is1 + Qc . The ommon

We pro eed similarily (with the introdu tion of
obtain

27

60 + 35z + 51z 2 , with omplex onjugate roots z
su h that |z| ≃ 1.08. Then kxk∞ ≤ kQk1 .kek∞ + k0.8Qie + 0.5Qis k∞ .kιk +
0
0
kQc k∞ . This yields kxk∞ ≤ 1105.
denominator of the

Q

is

The two linear lters are

ombined into an iterated nonlinear lter.

filter1() (resp. filter2()) is run with
(resp. [−800, 800]). We repla e the all
X ∈ [−339, 339] (resp. X ∈ [−1105, 1105]).

a pre- ondition of

X ∈ [−400, 400]

to the lter by its post ondition

void main () {
X = 0;
INIT1 = TRUE; INIT2=TRUE;
while (TRUE) {
X = 0.98 * X + 85;
if (abs(X)<= 400) {
filter1 ();
X=X+100;
INIT1=FALSE;
} else
if (abs(X)<=800) {
filter2();
X=X-50;
INIT2=FALSE;
}
}}
The program then

an be abstra ted into:

while (TRUE) {
X = 0.98 * X + 85;
maybe hoose X in [−1155, 1055];
}
We obtain X ∈ [−1155, 4250.02] by

running Astrée with a large number

of narrowing iterations, whereas Astrée
pre isely and

(x

annot bound

= 0.98x+ 85

X.

In this

has for unique solution

algebrai ally, but in more

annot analyze the original program
ase, the exa t solution

x = 4250)

ase.

4560.

Note that the non-abstra ted program
for

omputed

omplex lters this would not have been the

Roozbehani et al. have a bound of

X ∈ [0, 209].

[−1155, 4250]

ould have been

onverges to a value

≃ 205,

with

However, this very simple program illustrates our methodology

ompositional analysis: nding the optimal solution is possible here be-

ause the program is simple, but would not be possible in pra ti e if we had
added more nonlinear behavior and nondeterministi
rea tive

inputs, as in real-life

ode; whereas by analyzing pre isely ea h linear lter and plugging

the results ba k into a generi

analyzer, we get reasonable results.

28

10 Related works
In the eld of digital signal pro essing, some sizable literature has been devoted to the study of the ee ts of xed-point and oating-point errors on
numeri al lters. In the area of xed-point

omputation, bounds on the sizes

of the various operands are of paramount importan e: operands that leave
the pres ribed range will undergo saturation and the output signal will be
distorted. For these reasons, operands are s aled so as not to produ e digital
saturation; yet, the s ale fa tor should be made large enough that rounding
errors are very small

ompared to the typi al magnitude of the signal. While

the fa t that the l1-norm of the

onvolution kernel is what matters for judg-

ing overow, it is argued that this norm is overly pessimisti  [12, 11.3℄
[11, eq 13℄, not to mention the di ulties in estimating it. In pra ti e, lter
designers have preferred

riteria that indi ate no saturation for most  om-

monpla e inputs, ex luding pathologi al inputs. Our vision is dierent: our
results must be sound in all

ir umstan es, even pathologi al inputs.

The impa t of xed- and oating-point errors in digital lters was

las-

si ally studied from by modeling the errors as random sour es of known
distribution, independent of ea h other and with no temporal
(i.e.

orrelations between su

in reality, false: the
random; the

orrelation

essive values) [3, 20℄. These assumptions are,

omputational pro ess is fully deterministi , and not

omputations are generally interdependent (all

omputations

inside a lter depend on the past of the input variables); and there are temporal

orrelations. However,

ir uit designers are

on erned with the spe tral

distribution of output noise [11℄, and optimization of hardware or software
implementations with respe t to this noise, and these tools are adequate
for this. On the other hand, we merely aim at providing sound bounds for
the outputs of the system, but the bounds that we provide must be sound
without any extra and unfounded suppositions.
J. Feret has proposed an abstra t domain for analyzing programs
prising digital linear lters [6℄.
se ond degree lters. In
works, in a

om-

He provides ee tive bounds for rst and

omparison, we

onsider more

omplex lter net-

ompositional fashion; but we analyze spe i ations, and not C

ode (whi h is usually
loss of stru ture).

ompiled from those spe i ations, with

onsiderable

Another dieren e is that we do not perform abstra t

iterations. Feret's method

urrently

onsiders only se ond-order lters (i.e.

TF2), though it may be possible to adapt it to higher-order lters.
se ond-order lters, the bounds
in this paper are very

On

omputed by Feret's method and the method

lose (sin e both are based on a development of the

onvolution kernel, though they use dierent methods of tail estimation).
Lamb et al. [14℄ have proposed ee tive methods, based on linear algebra,
for

omputing equivalent lters for DSP optimization. They do not

ompute

bounds, nor do they study oating-point errors.
Roozbehani et al. [21℄ nd program invariants by Lagrangian relaxation

29

and semidenite programming, with quadrati

invariants. In order to make

problems tra table, they too apply a blo kwise abstra tion. The

lass of pro-

grams that they may analyze dire tly is potentially larger, but the results
are less pre ise than our method on some linear lters. They do not handle oating-point impre isions (though this

an perhaps be added to their

framework).
One possible appli ation of our method would be to integrate it as a
pre-analysis pass of a tool su h as Astrée [5℄. Astrée

omputes bounds on all

oating-point variables inside the analyzed program, in order to prove the
absen e of errors su h as overow.
reasonably a
by-wire

In order to do so, it needs to

ompute

urate bounds on the behavior of linear lters. A typi al y-

ontroller

ontains dozens of TF2 lters, some of whi h may be

integrated into more

omplex feedba k loops; in some

of the lters may yield too

ases, separate analysis

oarse bounds.

11 Con lusions and future works
We have proposed ee tive methods for providing sound bounds on the outome of

omplex linear lters from their ow-diagram spe i ations, as found

in many appli ations. Computation times are modest; furthermore, the nature of the results of the analysis may be used for modular analyses  the
analysis results of a sub-lter
the sub-lter

an be stored and never be re omputed until

hanges.

The usefulness of these methods is twofold. First, they
implemented in the graphi al user interfa e for designing
then be able to
a

ould be dire tly

ir uits. Users may

ompute gains or to

he k the stability of lters, taking into

ount oating-point errors (whi h

onventional Z-transform te hniques do

not

onsider). Se ond, they

stati

an be used as a way to automati ally obtain

analysis transformers or transfer fun tions: a stati

analysis tool

su h as Astrée may dete t that some program sequen e implements su h or
su h

omplex linear lter, and apply some invariant relation

omputed using

the te hniques in that paper.
In future works, we will examine the
positional, modular analysis.
non-linear lters

•

ase of non-linear lters and

The analysis of a

an be done in two ways or a

ombination of linear and

ombination thereof:

the overall behavior of a nonlinear lter may be
input-output relationship su h as

om-

onstrained by some

kOk∞ ≤ (1 + ǫ)kIk∞

rate limiter), and this input-output relationship

(example of a

an be integrated into

the abstra t semanti s as in Part 7;

•

the overall behavior of a linear lter
this bound information

an be pre isely bounded, and

an be fed into an analysis of a larger nonlinear

30

lter, su h as one based on stati ally

omputed relationships between

intervals [19℄

Referen es
[1℄ B. Blan het, P. Cousot, R. Cousot, J. Feret, L. Mauborgne, A. Miné,
D. Monniaux, and X. Rival. Design and implementation of a spe ial-

The Essen e of Computation: Complexity, Analysis, Trans-

purpose stati

formation

software. In

program analyzer for safety- riti al real-time embedded

, number 2566 in Le ture Notes in Computer S ien e, pages

85108. Springer Verlag, 2002.
[2℄ B. Blan het, P. Cousot, R. Cousot, J. Feret, L. Mauborgne, A. Miné,

PLDI

D. Monniaux, and X. Rival. A stati
software. In

analyzer for large safety- riti al

, pages 196207. ACM, 2003.

[3℄ Bru e W. Bomar, L. Montgomery Smith, and Roy D. Joseph. Round-

IEEE Trans. on Cir uits and Systems

o noise analysis of state-spa e digital lters implemented on oating-

II

point digital signal pro essors.
, 44(11):952955, 1997.

POPL '87:
Pro eedings of the 14th ACM SIGACT-SIGPLAN symposium on Priniples of programming languages

[4℄ Paul Caspi, Daniel Pilaud, Ni olas Halbwa hs, and John A. Plai e. LUSTRE: a de larative language for real-time programming. In

, pages 178188. ACM Press, 1987.

[5℄ Patri k Cousot, Radhia Cousot, Jérme Feret, Laurent Mauborgne, An-

ESOP

toine Miné, David Monniaux, and Xavier Rival. The ASTRÉE analyzer.
In

, number 3444 in Le ture Notes in Computer S ien e, pages

2130, 2005.
[6℄ Jérme Feret. Stati

analysis of digital lters. In

ESOP '04

, number

2986 in Le ture Notes in Computer S ien e. Springer-Verlag, 2004.

library

[7℄ Free Software Foundation.
, 2004.

GMP  GNU multiple pre ision arithmeti

GSL  GNU s ienti library
Standard for Binary Floating-Point Arithmeti
The MPFR Library

[8℄ Free Software Foundation.
[9℄ IEEE.

[10℄ INRIA et al.

, 2004.

. standard 754.

.

The Bell System Te hni al J.

[11℄ Leland B. Ja kson. On the intera tion of roundo noise and dynami
range in digital lters.
February 1970.

, 49(2):159184,

31

[12℄ Leland B. Ja kson.

Digital Filters and Signal Pro essing

. Kluwer, 1989.

Newton iteration for fa torization in
1998.

C[z].

C(z) and simultaneous

J. Complexity

[13℄ Peter Kirrinnis. Partial fra tion de ompostion in

, 14(3):378444,

PLDI '03

[14℄ Andrew A. Lamb, William Thies, and Saman Amarasinghe.
analysis and optimization of stream programs. In
25. ACM, 2003.

guide

[15℄ Xavier Leroy.
. INRIA.

Linear

, pages 12

The Obje tive Caml system, do umentation and user's

[16℄ Ja ques-Louis Lions et al.

Ariane 501: Flight 501 failure.

Te hni al

report, ESA / CNES, 1996. Available on WWW.
[17℄ A. Miné.

ESOP'04

LNCS

Relational abstra t domains for the dete tion of oating-

point run-time errors. In
Springer, 2004.
[18℄ David Monniaux.
meri al lters. In

, volume 2986 of

, pages 317.

Computer-aided veri ation: CAV '05

Compositional analysis of oating-point linear nu, number 3576

in Le ture Notes in Computer S ien e, pages 199212. Springer Verlag,
2005.

Stati analysis symposium (SAS)

[19℄ David Monniaux.

Optimal abstra tion on real-valued programs.
, 2007. To appear.

Trans. on Signal Pro essing

[20℄ Bhaskar D. Rao.

Floating point arithmeti

and digital lters.

, 40(1):8595, January 1992.

[21℄ M. Roozbehani, E. Feron, and A. Megretski.
and

omputation for software veri ation.

In

LNCS, page 606. Springer, 2005.

HSCC

In

IEEE

Modeling, optimization
, number 3414 in

J. of Computational and Applied Math.

[22℄ Siegfried M. Rump. Ten methods to bound multiple roots of polynomials.

, 156(2):403432, 2003.

M , let us note minori,j (M ) the determinant of the matrix
i and olumn j from M . We re all that for any
dimension n
n
X
(−1)n−1 mi,j .minor1,j (M )
det(M ) =
(47)

For any matrix

obtained by removing line
matrix

M

of

j=1

and that the determinant is

n-linear.

Re all that for any matrix

M

of in-

vertible determinant,

M −1 = det(M )−1 .

32



minori,j (M )

t

(48)

Id
Proof.
det(

If

, then there exists B ∈ Q[z] su h that

A ∈ Mn,n (Q[z](z) )
−
zA)
=
1 − zB
n

Lemma 3.

.

Proof by indu tion on

onsider

n.

(z)

The

ase

n = 1

is trivial.

Now let us

n > 1.

det(Idn − zA)
= (1 − za1,1 )minor1,1 (Idn − zA) +
= minor1,1 (Idn − zA) + z

n
X

j=2
n
X
j=1

(−1)n za1,j minor1,j (Idn − zA)

(−1)n za1,j minor1,j (Idn − zA)

(49)

The result follows by the appli ation of the indu tion hypothesis, and the
fa t that

BQ[z](z)

is a ring and thus the determinant of any matrix over that

ring is itself in the ring.

.

Corollary 4.

Mn,n (Q[z](z) )

Proof.

If A ∈ M

, then Id

n,n (Q[z](z) )

n

− zA)

has an inverse in

det(Idn −zA) is of the form 1−zP (z)/Q(z),
where P and Q are polynomials su h that the onstant oe ient of Q is
1, therefore (det(Idn − zA))−1 = Q(z)/(Q(z) − zP (z) is in Q[z](z) . All the
minori,j (Idn − zA) are elements of Q[z](z) , the result follows by applying
By the pre eding lemma,

Equ. 11.

33

