SICE Annual Conference 2007
Sept. 17-20, 2007, Kagawa University, Japan

Optimal Wavelet Expansion via Sampled-Data H ∞ Control Theory
M. Ogura†1 , M. Nagahara†2 , and Y. Yamamoto†3†

Graduate School of Informatics, Kyoto University, Kyoto, Japan1
ogura@acs.i.kyoto-u.ac.jp
2nagahara@ieee.org
3
yy@i.kyoto-u.ac.jp

Abstract: Wavelet expansion of an L2 signal requires the L2 inner product of the original signal and a scaling function.
In digital signal processing, it is common to use sampled data of continuous-time signals instead of the inner product. This
however causes a large reconstruction error, called “wavelet crime.” We therefore design a causal system which produces
an approximation of the inner product via sampled-data H ∞ control theory. We then make extensions to a multi-rate anda multi-wavelet case. By numerical examples, we show the effectiveness of the proposed method.
Keywords: wavelet, sampled-data control, prefiltering, signal reconstruction.

1. INTRODUCTION
Wavelet [3], [9] is an effective tool for analyzing
continuous-time signals. In wavelet theory, a continuoustime signal is expressed as a linear combination of waves(scaling and wavelet functions) which are localized intime and frequency. This expression, called wavelet
expansion, gives us time-frequency analysis of signals,
which is efficient in analyzing in particular non-stationary
signals. For example, the image format JPEG2000 isbased on wavelet expansion and has higher compression
rate and better quality than the conventional JPEG basedon Fourier expansion [10].
Wavelet expansion of a continuous-time signal begins
with calculating coefficients of each wave contained in
the signal. When the signal is in L2 (the Lebesgue spaceof square integrable functions on R), each coefficient,called a wavelet coefficient, is equal to the L2 inner product of the signal and a scaling function. The computationof the inner product requires all values of the signal on R.
In real situation, however, it is usual that only the sampled
data of the original signal are available. In this case, it is
impossible to obtain the exact values of the coefficients.
Therefore, approximation is essential for computing thecoefficients in real signal processing. The roughest approximation is to use the sampled data themselves as thecoefficients. This is often used in real, by which arises alarge reconstruction error, called “wavelet crime [9].”To prevent this crime, prefiltering methods have been
proposed [4], [5] in the situation where only the sampled data are available. They are based on the sampling theorem by Shannon, which assumes that the original continuous-time signals are fully band-limited up to
the Nyquist frequency. However, this assumption fails
to consider practical signals such as rectangle or triangle
waves.
In recent papers [8], [6], new approaches to this problem has been developed in the framework of sampleddata control theory. By using this theory, we can obtain
the optimal system taking account of the frequency characteristic of the original analog signals.
Motivated by these studies, we consider in the article

the approximation problem via sampled-data control theory. In [6], only the analog characteristic of the input signals is considered. On the other hand, we also consider
the characteristic of the reconstructed analog signals. The
design problem is formulated as a sampled-data H ∞ optimization with a generalized hold. We also make extensions to a multi-rate and a multi-wavelet case. By numerical examples, we show the effectiveness of the proposed
method.

2. WAVELET EXPANSIONIn this section, we discuss wavelet expansion. First,
we introduce multi-resolution analysis [7] which plays anessential role in wavelet theory. Multi-resolution analysis
is defined as an increasing sequence {Vj }j∈Z of closed
subspaces in L2 , and each subspace Vj has a Riesz basis{φj,k }k∈Z where
φj,k := 2j/2 φ(2j · −k).
In this definition, the function φ ∈ L2 is called as a scaling function. Next, consider the orthogonal complement
Wj := Vj+1  Vj . It can be shown that there exists afunction ψ ∈ L2 such that each subspace Wj has a Riesz
basis {ψj,k }k∈Z , whereψj,k := 2j/2 ψ(2j · −k).
Wavelet expansion of a continuous-time signal begins
with specifying a resolution J. Let f be in L2 and PJ
be the orthogonal projection in L2 onto the closed subspace VJ ⊂ L2 . By using the basis {φJ,k }k∈Z in VJ , theprojection of f onto VJ can be expressed as
PJ f =cJ (k)φJ,k .

(1)

k∈Z

On the other hand, corresponding to the direct sum
representation- 1422 -VJ = Vj0 ⊕ Wj0 ⊕ · · · ⊕ WJ−1 ,PR0001/07/0000-1422 ¥400 © 2007 SICE

the projection can be also decomposed as
PJ f =

cj0 (k)φj0 ,k +J−1

dj (k)ψj,k .

- e−ls(2)w-

j=j0 k∈Z

k∈Z

This is the wavelet expansion of f at the resolution J, andthe coefficients {cj0 , dj0 , . . . , dJ−1 } are called wavelet
coefficients. Since equations (1) and (2) are different expressions of the same signal PJ f , the wavelet coefficientscan be calculated efficiently from only {cJ (k)}k∈Z by afilter bank. From the above discussion, it can be seen thatthe accuracy of the wavelet expansion at the resolution
J depends on that of the projection onto VJ . In otherwords, it is necessary to accurately compute the coefficients {cJ (k)}k∈Z .
In particular, if {φJ,k }k∈Z is an orthonormal basis of
VJ , the coefficients {cJ (k)}k∈Z are given as follows:cJ (k) = (f, φJ,k )L2
 ∞
J/2
=2
f (t)φ(2J t − k)dt.
−∞

However, as explained in Section 1, this inner productcannot be computed exactly in practice. In order to overcome this difficulty, we propose in the next section a design method of systems, which produce approximated coefficients using only sampled data.

3. COMPUTATION OF COEFFICIENTSBY H ∞ OPTIMIZATIONIn this section, we propose a design method of systems which produce an approximation of wavelet coefficients of a continuous-time signal from only its sampled
data. The problem is formulated as a sampled-data H ∞optimization one. In this section, it is assumed that theresolution J is equal to 0 and that the signal we aim to
expand is in L2 [0, ∞). Moreover, we make an important
assumption about the scaling function.
Assumption 1: There exists m > 0 such that the support of the scaling function φ is in [0, m].

First, we characterize the wavelet coefficient c0 of thesignal f ∈ L2 [0, ∞). From the projection theorem, thecoefficient sequence c0 ∈ 2 is characterized as the solution c ∈ 2 of the following optimization problem:f − Hφ cL2 [0,∞) → min,(3)

where Hφ is defined as
Hφ : 2 → L2 [0, ∞) : c →c(k)φ0,k .

(4)

k≥0

Next, let us consider the error system in Fig. 1. On
the lower path, the exogenous signal w ∈ L2 [0, ∞) isband limited by a transfer function F which has the frequency characteristic of the continuous-time signals we

F

f s S

- K

c̃ - H

φv+ c?e−Fig. 1 Error system
aim to expand, and the signal f is the model signal to
be expanded. This signal f is then sampled by the ideal
sampler
S : F L2 [0, ∞) → 2 ,(Sf )(k) := f (k), k = 0, 1, 2, . . . .
where
F L2 [0, ∞) := {F w : w ∈ L2 [0, ∞)}.
The sampled data are processed by the digital filter K,
and then transformed into a continuous-time signal v =Hφ c̃ ∈ V0 .
Since the wavelet coefficient c0 of the signal f is characterized by (3), the signal c̃ in Fig. 1 approximates c0 if
the error e = f − Hφ c̃ is small. Hence, in order to obtaina digital filter which produces the approximation of thecoefficient c0 , it should be designed to make the error e
to be as small as possible.
Note that the signal f is delayed by e−ls : f → f (·−l)
on the upper path. Without this delay, the reconstruction
process becomes anti-causal one, and hence it is desirable
to take l to be no less than m.
Now, we formulate the design problem as follows.
Problem 1: Given a stable and strictly proper transfer function F , a delay l > 0, a scaling function φ, anda prespecified bound γ > 0, find a digital filter K satisfying Tew ∞ < γ, where Tew ∞ is the L2 [0, ∞)induced norm of the operator Tew from w to e in Fig. 1.
Invoking fast-sample/fast-hold approximation [11], we
can reduce this sampled-data H ∞ optimization problem
to a discrete-time one. This discrete-time problem is a
finite dimensional one by the assumption 1, and can beeffectively solved by MATLAB routines. Hence, by bisection search, we can easily obtain a sub-optimal digital
filter minimizing Tew ∞ .

4. EXTENSIONS
4.1 Higher resolution
In the above section, we have shown an approximation
method to compute wavelet coefficients at the resolution
0. Then, what should we do in order to compute higher
resolution coefficients?
Let us consider the error system in Fig. 2. Comparedwith Fig. 1, the difference is that the sampled data of theoriginal signal f are first upsampled by ↑ 2J . The upsampler ↑ 2J is defined as follows:- 1423 -↑ 2J : {x(k)}∞k=0 → {x(0), 0, 0, . . . , 0, x(1), 0, . . .}.
  
2J −1

- e−ls- e−lsw-

F

f sv + c?e
S - ↑ 2J - K - HφJ −

- K1
w-

Fig. 2 Error system: higher resolution case

F

f s-

S

s-

k≥0

which acts with it’s sampling period 1/2J .
The optimal digital filter K is obtained by minimizing
L2 [0, ∞) induced norm of the operator from w to e inFig. 2. This minimization is also a multi-rate sampleddata H ∞ optimization. The multi-rate system is equivalently transformed into a single-rate one by discrete-time
lifting [1] or polyphase representation [9]. Then, by fastsample/fast-hold approximation, the problem is reduced
to a discrete-time one.
4.2 Multi-wavelet
Multi-wavelet is an extension of the conventional
wavelet and uses a wavelet basis generated by several
scaling functions called multi-scaling functions. By thisextension, there arises great flexibility in the design of awavelet, and desirable properties such as symmetry or orthonormality of scaling functions can be realized simultaneously [2].
In multi-wavelet, the multi-resolution analysis is constructed by multi-scaling functions φ(1) , . . . , φ(r) ∈ L2 ,(1)(r)
and the subspace VJ has a Riesz basis {φJ,k1 , . . . , φJ,kr :k1 , . . . , kr ∈ Z}, where(l)

φJ,k := 2J/2 φ(l) (2J · −k).
Then, the orthogonal projection of a signal f ∈ L2 [0, ∞)onto the closed subspace VJ ⊂ L2 has the form
PJ f =

r 

(l)
(l)
cJ (k)φJ,k (k).
l=1 k∈Z(1)(r)

As in the scalar case, these coefficients cJ , . . . , cJ
are called as wavelet coefficients of the signal f at theresolution J, and can be characterized as the solutionc(1) , . . . , c(r) ∈ 2 of the following optimization problem:


r




(l) 
Hφ(l) c 
→ min,
(6)
f −


l=1

L2 [0,∞)

φ

..
.

- Kr
The upsampling operation is implemented by inserting
2J − 1 equidistant zero-valued samples between two consecutive samples of x(k) before the sampling rate is multiplied by the factor 2J . Then the signal is processed by
the digital filter K, which estimates the data between theoriginal samples. Then the estimated (virtually) higher
resolution signal is transformed to a continuous signal by
the following generalized hold:c(k)φJ,k ,(5)
HφJ : 2 → L2 [0, ∞) : c →

c̃(1)- H (1).. - c?
v+ c?e.
−6c̃(r)- H (r)φFig. 3 Error system: multi-wavelet casewhere Hφ(l) is defined in the same way as (4) and (5).
Now, let us consider the error system in Fig. 3. As
in Section 3, the resolution J is assumed to be equal to
0. Since the wavelet coefficients are characterized as thesolution of the optimization problem (6) similar to thescalar case one (3), we can formulate the design problem
in the same way as problem 1.
Since H ∞ optimization is executed in the state spacein which MIMO (multi-input multi-output) systems canbe treated in the same way as SISO (single-input singleoutput) systems, the formulated problem can be treated asa sampled-data H ∞ optimization one (in this case, the filter to be designed is a single-input and r-output system),and so can be reduced to a discrete-time finite dimensional one. And hence sub-optimal filters K1 , . . . , Kr
which compute the approximation of the wavelet coefficients c(1) , . . . , c(r) respectively are obtained.

5. NUMERICAL EXAMPLESIn this section, we present numerical examples.
Throughout this section, the frequency characteristic F
is taken as follows:
F (s) =

1.
(7s + 1)(0.7s + 1)

5.1 Scalar wavelet
As an example of a scalar case, we consider the orthogonal Daubechies-2 wavelet [3]. We make comparisons between the wavelet crime, the proposed methodwith J = 0, and one with J = 2. For the proposed
method, the optimal filters are designed with l = 5 inboth cases.
First, we make a comparison by observing the frequency responses of the error systems shown in Fig. 1(J = 0) or Fig. 2 (J = 2). Fig. 4 shows the frequency responses. In this figure, the error frequency response of the
wavelet crime reconstruction (i.e., K = 1) is optimized
by varying the delay l so that the H ∞ norm is minimized.
It can be seen that the errors by the proposed method areboth smaller than the wavelet crime in a wide frequency
range. In low frequency range, the wavelet crime reconstruction shows a better response. It follows that the error
in wavelet crime appear particularly in a middle and highfrequency range.
We can also observe that, by taking the resolution J
higher than 0, we can obtain more accurate reconstruction. This shows effectiveness of the use of the interpo-- 1424 -Fig.
4 Frequency responses of the error systems:
wavelet crime (dotted), proposed with J = 0(dashed), and J = 2 (solid)

Fig. 5 Time responses of expansions: Source signal(dash-dotted), sampled data (circle), wavelet crime(dotted), proposed with J = 0 (dashed), and J = 2(solid)lator. Note that the sampling period in the sampler is notchanged in the case of J = 0 and J = 2.
Next, we execute the wavelet expansions of a chirp
signal,
3
t
+ .
(7)
9 10
This signal is a typical non-stationary sinusoid signal
whose frequency increases linearly with respect to timet. The reconstructed signals are shown in Fig. 5. The
wavelet crime reconstruction shows a peaked response.
In Fig. 6, we illustrate the squared reconstruction errors.
In this figure, our reconstruction shows smaller error than
the wavelet crime. It should be noticed that the error by
the wavelet crime reconstruction increases as time passes.
This justifies the error frequency response Fig. 4.
We can conclude from these simulations that• wavelet crime shows in particular middle- and highfrequency errors,∞• sampled-data H
optimization is effective insampled-data wavelet expansion,• upsampling is also effective.
f (t) = sin(ω(t) · t),ω(t) =Fig.
6 Time responses of squared expansion errors: wavelet crime (dotted), proposed with J = 0(dashed), and J = 2 (solid)

5.2 Multi-wavelet
In this subsection, we show an example of a multiwavelet case using the symmetric and orthogonal multiwavelet with approximation order 2 proposed in [2]. This
multi-wavelet has the multi-scaling functions shown in
Fig. 7. With this multi-wavelet, we compare the proposed
method with one proposed in [5].
The design parameters are as follows. The reconstruction delay l = 4, the number of scaling functions is r = 2,
and the subspace is V0 (i.e., J = 0).
Fig. 8 shows the frequency responses of the error system shown in Fig. 3. Note that in the conventional case,the reconstruction delay l is equal to 3/2. In this figure,the proposed method shows the better error frequency response than that of the conventional method, in particularin a middle- and high-frequency range. The maximum
difference of the response is about 7 [dB] at around 0.5[rad/sec].

- 1425 -(a) φ(1)(b) φ(2)Fig. 7 Multi-scaling functions of Chui-Lian wavelet

Fig. 8 Frequency responses of the error systems: conventional (dashed) and proposed (solid)Fig.

10 Time responses of squared expansion errors:
conventional (dashed) and proposed (solid)

the discrete-time lifting and an MIMO filter design, respectively. Finally, by numerical examples, the effectiveness of the proposed method have been shown.

REFERENCES

Fig. 9 Time responses of expansions: source signal (dotted), sampled data (circle), conventional (dashed),
and proposed (solid)
Then, we simulate the reconstruction with the chirp
signal in (7). Fig. 9 shows the time response against thechirp signal. The conventional reconstruction shows agood response at the sampling instants, while in the intersamples it shows large gaps. On the other hand, the response of the proposed method shows the better responsethan the conventional one. Fig. 10 shows the squared
errors. This figure shows that the error of the proposed
method is smaller than the conventional one. From thesesimulation results, the sampled-data H ∞ optimization iseffective also in the multi-wavelet case.

6. CONCLUSION
In this article, a new design method of digital filters
which compute approximated wavelet coefficients has
been proposed, in the case that only the sampled-dataof the original signal are available. The design problem
is formulated as an sampled-data H ∞ optimization one.
We have shown the problem can be reduced to a discretetime one by the fast-sample/fast-hold method. The obtained discrete-time H ∞ optimization can be effectively
solved by MATLAB routines. This method can be extended to a multi-rate and a multi-wavelet case by using[1] T. Chen and B. Francis, Optimal Sampled-Data
Control Systems, Springer, 1994.
[2] C. K. Chui and J. Lian, A study of orthonormal
multi-wavelets, Appl. Numer. Math., 20, pp. 273–
298, 1996.
[3] I. Daubechies, Ten Lectures on Wavelets, SIAM,
1992.
[4] S. Ericsson and N. Grip, Efficient wavelet prefilterswith optimal time-shifts, IEEE Trans. Signal Processing, 53, No. 7, pp. 2451–2461, 2005.
[5] T.-C. Hsung, M.-C. Sun, D.P.-K. Lun and W.-C.
Siu, Symmetric prefilters for multiwavelets, IEE
Proc.-Vis. Image Signal Processing, 150, No. 1, pp.
59–68, 2003.
[6] K. Kashima, Y. Yamamoto and M. Nagahara, Optimal wavelet expansion via sampled-data control
theory, IEEE Signal Processing Lett., 11, No. 2, pp.
79–82, 2004.
[7] S. G. Mallat, A theory of multiresolution signal
decomposition: the wavelet representation, IEEE
Trans. Patt. Anal. Machine Intell., 11, pp. 674–693,
1989.
[8] P. Qian and B. A. Francis, Optimal initialization of
the discrete wavelet transform, Proc. Workshop Recent Advances in Control, 1998.
[9] G. Strang and T. Nguyen, Wavelet and Filter Banks,
Wellesley-Cambridge Press, 1997.
[10] D. S. Taubman and M. W. Marcellin, JPEG2000:
standard for interactive imaging, Proceedings of theIEEE, 90, No. 8, pp. 1336–1357, 2002.
[11] Y. Yamamoto, A. G. Madievski and B. D. O. Anderson, Approximation of frequency response forsampled-data control systems, Automatica, 35, No.
4, pp. 729–734, 1999.

- 1426 -