CDC01-INV1701

Optimal Wavelet Expansion
via Sampled-Data Control Theory
Kenji Kashima1

Yutaka Yamamoto 2

Abstract
Wavelet theory provides a new type of function expansion and and has found many applications in signal
processing. Discrete wavelet transform of a signal x(t)in L2 (R) is usually computed by the so-called pyramid algorithm. It however requires a proper initialization, i.e., expansion coeﬃcients with respect to thebasis of one of the desirable approximation subspaces.
An interesting question is how we can obtain such coefﬁcients when only sampled values of x(t) are available.
This paper provides a design method for a digital ﬁlter
that optimally gives such coeﬃcients assuming certaina priori knowledge on the frequency characteristic of
the target functions. We then extend the result to thecase of non-orthogonal wavelets. Examples show theeﬀectiveness of the proposed method.

Masaaki Nagahara 3

subsequent approximation accuracy in the pyramid algorithm depends on these coeﬃcients.
In digital signal processing, however, it is more common
practice that only sampled values of signals are available. Then it is not possible to accurately compute
such expansion coeﬃcients in the ﬁrst place. One maythen directly employ the sampled values as such coeﬃcients, with the idea that the intersampling changes are
not large. This is convenient, but obviously not satisfactory from thea theoretical point of view, and hencereferred to as “wavelet crime” by Strang and Nguyen[8].
Pu and Francis [7] gave an optimal ﬁlter to initialize
the expansion from sampled values, and also gave anerror analysis. They have shown that the wavelet crime
can lead to a very large error. On the other hand, theirﬁlter is not causal and inﬁnite-dimensional, and hence
not easily realizable.

1 Introduction
Wavelet theory provides a new type of signal representation: it expands functions in terms of small waves,localized both in time and frequency, in contrast toFourier analysis where basis functions are sinusoids. It
hence applies to many cases where local information
is important, for example, image compression. In particular, discrete wavelet transform ﬁts naturally withthe digital computer with its basis functions deﬁned
by summations but not integrals or derivatives.
Wavelet analysis of a signal begins by approximation by
projecting it onto one of the approximation subspaces,
which is spanned by shifted scaling function, constituting multiresolution analysis. This determines the ﬁnest
resolution, and one is led to ﬁnd expansion coeﬃcientsof coaser orders. The pyramid algorithm [8] can thenbe invoked to obtain such lower scale approximationcoeﬃcients. Furthermore, this can be implemented bya multirate ﬁlter bank [6, 8, 9]. In this procedure, one
has to ﬁrst ﬁnd expansion coeﬃcients with respect tothe basis formed by the shifted scaling functions. The
1 kashima@acs.i.kyoto-u.ac.jp, Department of Applied
Analysis and Complex Dynamical Systems, Graduate School ofInformatics, Kyoto University, Kyoto 606-8501, JAPAN
2 yy@i.kyoto-u.ac.jp, the same address as above; author towhom all correspondence should be addressed.
3 nagahara@acs.i.kyoto-u.ac.jp, the same address as
above.

We instead propose to use a ﬁnite-dimensional andcausal (with ﬁnite-time delay) digital ﬁlter; for its design, we assume a certain frequency domain characteristic, and employ the sampled-data control theory as
per [5, 11, 12, 13] etc. While the performance may
become less optimal compared to that of [7], it is easily realizable and, moreover, as we see later, it can beextended to non-orthogonal wavelets.
The paper is organized as follows: We ﬁrst introduce
discrete wavelet transform and formulate our problem.
We reduce this to an equivalent discrete-time H∞ problem. Then these results are extended to the case ofoversampled signals and non-orthogonal scaling functions. The optimal ﬁlter is compared with the conventional wavelet crime to exhibit the diﬀerence.
2 Problem Formulation
2.1 Discrete Wavelet Transform
We begin with the basic setup in wavelet theory [6, 8].
For the signal space L2 (R), a multiresolution analysis(MRA) is a series {Vj }j∈Z of closed subspaces having
following properties:
1. Vj ⊂ Vj+1 for all j ∈ Z

p. 1

2.

lim Vj = {0} and lim Vj = L2 (R)j→−∞cJ (k)

j→∞

r

✲ HPF

✲ ↓2

dJ−1 (k)
✲

✲ LPF

✲ ↓2

cJ−1 (k)✲3. f (t) ∈ Vj ⇔ f (2t) ∈ Vj+1
4. f (t) ∈ V0 ⇔ f (t − kh) ∈ V0
5. V0 is spanned by {φ(· − kh)}k∈Z

Figure 1: First stage of hte pyramid algorithm.

where h > 0 is ﬁxed. The function φ in condition
5 is called a scaling function and Vj is spanned by{φj,k }k∈Z where
φj,k (t) := 2j/2 φ(2j t − kh)for all j ∈ Z according to condition 3. Thus for x(t) ∈L2 (R), its orthogonal projection vj (t) onto Vj is viewed
as the approximation of x(t) at scale j. The scale of
the resolution increases as we move up in the nest.
The wavelet subspace W0 is deﬁned to be the orthogonal complement of V0 in V1 .
V1 = V0 ⊕ W0
which extends to
L2 (R) = Vj0 ⊕ Wj0 ⊕ Wj0 +1 ⊕ · · · .

(1)

Since these wavelets reside in the space spanned by thesubsequent space spanned by the shifted and dilated
scaling functions Wj ⊂ Vj+1 follows; i.e.,ψ(t) =
h(k)φ(2t − k), k ∈ Z(2)

LPF denote, respectively, the ideal highpass and lowpass ﬁlters with cutoﬀ frequency π/h[rad/sec], and ↓ 2
is a downsampler deﬁned by↓ 2 : 2 (Z) −→ 2 (Z) : {x(k)}k∈Z → {x(2k)}k∈Z .
In practice, computation of the DWT of x(t) ∈ L2 (R)goes as follows:
1. Select some level J where vJ (t) ∈ VJ represents
x(t) to a desired degree of resolution.
2. If the approximation coeﬃcients at this level{cJ (k)} k∈Z are available, we can obtain{dJ−1 (k)}k∈Z and {cJ−1 (k)}k∈Z by ﬁlter bankabove.
3. Since we can similarly obtain {dj−1 (k)}k∈Z and{cj−1 (k)}k∈Z along with {cj (k)} k∈Z , repeat anumber of times down to the coarsest desired
scale j0 .
This algorithm is called the pyramid algorithm.

kholds. For an appropriate series {h(k)}, the function
generated by (2) gives the mother wavelet ψ(t) suchthat Wj is spanned by {ψj,k }k∈Z where

As shown in step 2, the pyramid algorithm should beinitialized by {cJ (k)} k∈Z . In what follows, we ﬁrst
consider the case of orthogonal wavelets.

ψj,k (t) := 2j/2 ψ(2j t − kh).

Assumption 1 Scaling function φ is orthonormal,i.e.,φ(t), φ(t − kh)L2 (R) = δk,0(3)

We have now constructed a set of functions that could
span all of L2 (R). According to (1), any function
x(t) ∈ L2 (R) could be written as a series expansion
in terms of the scaling function and wavelets:
x(t) =∞
k=−∞cj0 (k)φj0 ,k (t) +

∞
∞



dj (k)ψj,k (t).

k=−∞ j=j0(4)This is called the wavelet expansion and Discrete
Wavelet Transform (DWT) is the map from x(t) ∈L2 (R) to {(cj0 (k), dj0 (k), dj0 +1 (k), . . .)}k∈Z ∈ l2 (Z).
2.2 Computation of DWTFirstly it is known that we can apply the ﬁlter bankshown in Figure 1 to compute {dj−1 (k)}k∈Z and{cj−1 (k)}k∈Z along with {cj (k)}k∈Z [8]. Here HPF andfor all k ∈ Z.
where δij denotes the Kronecker delta.
This assumption makes {φJ,k }k∈Z an orthonormal system of VJ and the projection coeﬃcients are given by
the inner product x, φJ,k L2 (R) . However even in thiscase it will be diﬃcult to compute by means of its definition itself. Since often only sampled data of x(t)are available, common practice is to use these sampled values directly in place of the correct coeﬃcients{cJ (k)} k∈Z . Because this is of course for convenienceand causes errors, it is called the “wavelet crime” in [8].
Our problem here is to compute {cJ (k)} k∈Z directlyfrom these samples via digital ﬁlters. This problem isp. 2

called the optimal initialization problem of the DWT
and is formulated and solved by giving an optimal solution in [7].
2.3 Problem Formulationw(t)✲x(t) ✲
r
W
e−mhs
✲ Sh

✲ S̃
h

c0 (k)✲ Kd+✲ ❞❄e(k)✲−Here kˆd (ejω ), ŵ(jω) and φ̂(jω) denote discrete-timeor continuous-time Fourier transform of Kd , W and φrespectively.
Unfortunately, this ﬁlter in Theorem 1 is neither causal
nor ﬁnite-dimensional, hence diﬃcult to implement.
We thus reformulate the problem as follows: First notethat
(S˜h x)(k)We formulate the design problem of a digital ﬁlter for
optimal initialization. Let J = 0 without loss of generality. Consider the block diagram Figure 2. We ﬁrst
take m = 0. It is appropriate to consider a class offunctions x(t) for such expansion, and as in H∞ control
theory, we attempt to minimize the L2 -induced norm
from w to the error e. W is an LTI weighting low-passﬁlter that makes this problem well-posed. This governs
the frequency domain characteristic of the target functions. For a given scaling function φ, S̃h is the operator
deﬁned byS̃h : L2 (R) −→ 2 (Z) : x(t) → {c0 (k)}k∈Z .
Sh denotes the sampler with sampling period h. Kd
is the digital ﬁlter to be designed. The upper pathgives the idealized wavelet expansion coeﬃcients for x
whereas the lower path gives an approximant resulting
from a sampler and a digital ﬁlter. The objective here
is to minimize the L2 -induced norm of the operator
2

Tew : L (R) −→  (Z) : x(t) → {e(k)}k∈Z .

=Sh (Fx)(k).

=Figure 2: Error system2=

x, φ0,k L2 (R) ∞x(t)φ(t − kh)dt−∞
 ∞x(t)φ̌(kh − t)dt=−∞

Here F denotes a continuous-time linear time-invariant
system with impulse response φ̌(t) := φ(−t). Thismeans that Problem 1 is an optimal discretization
problem of non-causal continuous-time system F. In
view of this, we make yet another assumption.
Assumption 2 There exists an integer l > 0 such thatthe scaling function φ(t) has compact support in [0, lh],i.e.,φ(t) = 0, t ∈/ [0, lh].
(6)
Many well-known wavelets satisfy this assumption, forexample, Haar scaling function φH and 2nd order Bspline (triangle) spline φT deﬁned byφH (t) =(5)

In other words, the worst case norm of this block diagram against signals in L2 (R).
Problem 1 For m = 0, given stable W (s) and orthonormal scaling function φ, find a digital filter Kd (z)which minimizes Tew  in (5).
Pu and Francis [7] have shown the following theorem
for this problem.
Theorem 1 Problem 1 has the (not necessarily
unique) solution|ŵ(jωh + j2πk)|2 φ̂(−jωh − j2πk)kˆd (ejωh ) = k2k |ŵ(jωh + j2πk)|1,0,t ∈ [0, h]otherwise,




and
φT (t) =t,t ∈ [0, h]
2h − t, t ∈ [h, 2h] .
0,otehrwise,(7)

With this assumption, for m ≥ l allowing an m step delay, the error system showed in Figure 2 become causal.
Our design problem then becomes as follows:
Problem 2 Given an integer m > 0, stable W (s)and scaling function φ satisfying Assumptions 1 and
2, find a digital filter Kd (z) that minimizes Tew  in(5), where (A, B, C) denotes a minimal realization ofW.
3 Reduction to finite dimensional problem

Suppose l = 1 for simplicity; the same method applies
and the corresponding worst case error norm isto the case of l > 1.
1/2

2

 

First
we apply the lifting technique [1, 10] to.
sup kˆd (ejωh ) − φ̂(−jωh − j2πk) ŵ(jωh + j2πk)
ω
continuous-time
signals in Figure 2, and x̃[k](θ) denoteskp. 3

lifted signal, i.e., x̃[k](θ) := x(kh + θ). By lifting x andw, we can obtain a lifted representation of system W
X(k + 1)x̃[k](θ)

where A ∈ Rn×n and X(k) ∈ Rn . We then have, from
(6), that (k+1)h(S̃h x)(k)

=

=

Adt
Cdt



h
BeA(h−τ ) w̃[k](τ )dτ0θ
A(θ−τ )
B w̃[k](τ )dτ.
0 Ce

eAh X(k) +CeAθ X(k) +

=
=


Gdt

:=

φ(t − kh)x(t)dt

Bdt
Ddt

Ad
 Cd1

 0

 ..
 .

 0

 0
C



··· ······ ···
0 ···
.. ..
.
.
··· 0
1
0 ··· 00 ··· ···
0
0
1
..
.

0
0
0
..
.

Bd1
Dd11
0
..
.

0
1
0

0
0
0

kh
 h

=

φ(τ)x(kh + τ )dτ

✛

0

 h
=



0
0
0
..
.








0 

−1 
0
(8)

✛
Gdt

φ(τ)x̃[k](τ )dτ.

✛0

The error system in Figure 2 for m = l − 1 can be
rewritten as that in Figure 3 where✲ Kd

Ad := eAh ∈ Rn×n

Figure 4: Equivalent discrete time system

2

B1 : L [0, h] −→ R : x̃(τ ) →n

 h
Cd1 :=D110 heA(h−τ ) B x̃(τ)dτ

0

CeAθ φ(θ)dθ ∈ R1×n

: L2 [0, h] −→ Rh θ
: w̃(τ ) → 0 0 CeA(θ−τ ) Bφ(θ)w̃(τ )dτdθ.



e(k)
✛

Ad
 Cd1
C

B1
D11
0



✛
0
−1 
✛
0

w̃[k]

w(t)
✲ W

Figure 3: Error system

This can be converted ([2, 3]) to an equivalent discretetime system by deﬁning Bd1 , Dd11 as


B1
D11





B1∗

∗
D11




,

Bd1
Dd11

We now give some generalizations, including the caseof non-orthogonal wavelets.
4.1 Projection onto higher resolution subspace
We have given a procedure to optimally initialize thewavelet expansion, i.e., optimally projecting onto V0
which has the same scale of resolution as the samplingrate. It is however possible to project onto a subspacewith a higher scale resolution than V0 , by making useof a given frequency distribution (i.e., weighting W ) of
the original signals. This leads us to the design problem
depicted in Figure 5.

✲ Kd

LLT :=For l > 1, we can apply the same method by lifting the scaling function. This problem is solvable(sub)optimally and the dimension of the obtained ﬁlter
does not exceed n + m.
4 Extensions


:= L.

Therefore we can obtain the following theorem.
Theorem 2 Problem 2 is equivalent to an 2 induced
norm minimization problem of discrete-time system
in Figure 4. Here Gdt is given by (8) and Adt ∈R(n+m)×(n+m) .

x(t)
q
✲ Sh

✲S̃h/2JcJ (k)✲ ↑ 2J✲ Kd+ e(k)✲❄❞ ✲−Figure 5: Error system
The diﬀerences with Problem 2 are S̃h/2J and ↑ 2J .
First deﬁne S̃h/2J :S̃h/2J : L2 (R) −→ 2 (Z) : x(t) → {cJ (k)}k∈Z ,and the output of upper path is cJ (k). In accordancewith this, the sampled values are also upsampled byp. 4

↑ 2J to match with the upper path in the sampling
rate:
↑ 2J

: 2 (Z) −→ 2 (Z): v(k) → u(k) :=v(k)0, k = 2J l, l = 0, 1, . . .
,otherwiseand processed by Kd (z). This gives an approximation
problem in the space VJ instead of V0 .
A diﬃculty here is that the overall error system becomes multirate due to the upsampler, but one can reduce this problem to a single-rate problem via discretetime lifting [4, 3, 11, 12].
4.2 Non-orthogonal scaling function
When Assumption 1 is not satisﬁed, taking the innerproduct with the scaling function does not directly give
projection coeﬃcients. In such a case we can employ
the dual basis {ξk (t)}k∈Z such that
φ0,p , ξq L2 (R) = δp,q .
According to MRA, we can takeg(k)φ(t − k), k ∈ Z.
ξk (t) =

3

2

1

0

−1

−2

−3

−4
30

32

34

36

38

40

42

44

46

48

50

time[sec](10)

Figure 6: Time response: source signal (dotted), wavelet
crime(dash) and proposed (solid)

for the dual basis. Here desirable coeﬃcients given byx, ξk L2 (R) are represented asc0 (k) =

4

(9)

k



We ﬁrst compare the present method with that by the“wavelet crime.” Note that the Haar wavelet satisﬁes
Assumption 1. Figures 6, 7 show the time responses of
the two methods against the source signal x(t) = 4 sin t.
Figure 6 shows the projection onto V0 whereas Figure 7
shows the projection onto V1 , per the method given in
Section 4.2. The latter exhibits more ﬁdelity. Figure 8
exhibits the corresponding squared error. The present
ﬁlter shows much advantage over the wavelet crime.

g(l)c0 (k + l), k ∈ Z

3

2

l

where c0 (k) = x, φJ,q L2 (R) . Truncating the sum
above suitably to approximate the expansion, we canobtain c0 (k) along with c0 (k) via FIR digital ﬁlter
KL (z) =

4

lf


1

0

−1

−2

g(l)z

l−lf

.

−3

l=ls
−4
30

32

34

For c0 (k) we can apply the method in the previous
section.
5 A Design Example

36

38

40

42

44

46

48

50

time[sec]
Figure 7: Time response: source signal (dotted) and pro-We present a design example for the case of the Haarand 2nd order B-spline wavelets. Put the sampling
period h = 1, delay step m = 2 and weighting function
W (s) be1
W (s) =(T s + 1)(10T s + 1)
where T := 22.05/π. Note here that the sampling rate
for audio CD (and MD) is 44.1kHz, and the above
time constants correspond to the frequecies 1kHz and
10kHz when normalized to h = 1. The weighting W
above simulates a fairly wide-range orchestral musicfrequency energy distribution observed in many commercial CD’s.

posed (projection onto V1 ) (solid)

Now consider the 2nd order B-splines. The scaling
function is φT given by (7). Note however that thisdoes not satisfy Assumption 1 and its corresponding
subspaces consist of piecewise linear functions. For thisscaling function the dual basis is given by (10) where√ √g(l) = 3( 3 − 2)|l| .
Here we take
KL (z) =3
g(k)z l−3

l=−3

p. 5

15

0.08

0.07

0.06
10
0.05

0.04

0.03
5
0.02

0.01

0
30

32

34

36

38

40

42

44

46

48

50

0
30

31

32

33

34

35

36

37

38

39

40

time[sec]time[sec]

Figure 8: Time responses of squared errors: wavelet crime

Figure 10: Time responses of squared errors: wavelet(dash), projection onto V0 (dotted) and projection onto V1 (solid)

and the total delay is taken to be 5 steps. The resultsfor x(t) = sin 0.3t + 0.5 sin t are shown in Figure 9, 10.
1.5

1

0.5

0

−0.5

−1

−1.5
30

31

32

33

34

35

36

37

38

39

40

time[sec]
Figure 9: Time response: source signal (solid), wavelet
crime(dash) and projection onto V0 (dotted)

Both examples exhibits quite an admirable performance of the proposed method.

6 Concluding Remarks
We have presented a new method, based on sampleddata control theory, for designing a digital ﬁlter in the
optimal initialization problem of the DWT. The obtained ﬁlters are ﬁnite-dimensional and show advantages to a naive representation of using sampled values directly. We have also extended the results tothe higher order expansion and to the non-orthogonal
wavelets.

References[1]
B. Bamieh, J. B. Pearson, B. A. Francis and A. Tannenbaum, A lifting technique for linear periodic systems

crime (dash) and projection onto V0 (solid)

with applications to sampled-data control systems, Syst.
Control Lett., 17: 79–88, 1991.
[2]
T. Chen and B. A. Francis, “On the L2 -induced normof a sampled system,” Systems and Control Letters, vol. 15,
pp. 211–219, 1990.
[3]
T. Chen and B. A. Francis, Optimal Sampled-data
Control Systems, Springer, 1995.
[4]
T. Chen and B. A. Francis, “Design of multirate ﬁlter
banks by H∞ optimization,” IEEE Trans. Signal Processing, SP-43, pp. 2822–2830, 1995.
[5]
P. P. Khargonekar and Y. Yamamoto, “Delayed signal reconstruction using sampled-data control,” Proc. of
35th Conf. on Decision and Control, pp. 1259–1263, 1996.
[6]
S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, 1998.
[7]
PU Qian and B. A. Francis, “Optimal Initialization of
the Discrete Wavelet Transform,” submitted for publication
May 5, 1999.
[8]
G. Strang and T. Nguyen, Wavelet and Filter Banks,
Wellesley-Cambridge, 1997.
[9]
P. P. Vidyanathan, Multirate Systems and Filter
Banks, Prentice Hall, Englewood Cliﬀs, 1993.
[10] Y. Yamamoto, “New approach to sampled-data systems: a function space method,” Proc. 29th Conf. on Decision and Control, pp. 1881–1887, 1990.
[11] Y. Yamamoto, H. Fujioka and P. P. Khargonekar,“Signal reconstruction via sampled-data control with multirate ﬁlter banks,” Proc. 36th Conf. on Decision and Control, pp. 3395–3400, 1997.
[12] Y. Yamamoto and P. P. Khargonekar, “Fromsampled-data control to signal processing,” in Learning,
Control and Hybrid Systems, Springer Lecture Notes inControl and Information Sciences, vol. 241, pp. 108–126,
1998.
[13] Y. Yamamoto, M. Nagahara and H. Fujioka, “Multirate signal reconstruction and ﬁlter design via sampled-data
H ∞ control,” MTNS2000, 2000.
[14] G. Zelniker and F. J. Taylor, Advanced Digital Signal
Processing, Dekker, 1994.

p. 6