Engineering Perspective, Part Discrete Wavelet Techniques': Tutorial Wavelets From An Electrical 1

A Tutorial on Wavelets from an Electrical
Engineering Perspective, Part 1:

Discrete Wavelet Techniques’
T.K. Sarkar’, C. Su’, R. Adve’, M. Salazar-Palma’, L.
Garcia-Castillo’, Rafael R. B o g
‘Department of Electrical and Computer Engineering

Syracuse University
121 Link Hall
Syracuse, New York 13244-1240 USA
Tel: +1 (315) 443-3775
Fax; +1 (315) 443-4441
E-mail: tl<[email protected]
Web: http://web.syr.edu/-tksarkar
2Dept. Senales and Sistemas

Polytechnique University of Madrid
28040 Madrid
Spain
E-mail: [email protected]
3Dept. Electronica y Electromagnetism0

Universidad de Sevilla
41012 - Sevilla
Spain
Keywords: Wavelets, wavelet transforms, filters, matrix equations, 2. Introduction

operators
Many books and numerous papers have been published
describing wavelets. It is not possible for us to include all the ref-
1. Abstract erences. Selected references [ 1-31 have been chosen to illustrate
where additional materials are available. No attempt has been made
The objective of this paper is to present the subject of wave- to provide the earliest reference material.
lets from a filter-theory perspective, which is quite familiar to
electrical engineers. Such a presentation provides both physical Wavelets are a set of functions that can be used effectively in
and mathematical insights into the problem. It is shown that taking a number of situations, to represent natural, highly transient phe-
the discrete wavelet transform of a function is equivalent to filter- nomena that result from a dilation and shift of the original wave-
ing it by a bank of constant-Q filters, the non-overlapping band- form. For example, when a pulse propagates through a layered
widths of which differ by an octave. The discrete wavelets are pre- medium, due to dispersion and for different electrical properties of
sented, and a recipe is provided for generating such entities. One of the layers, the pulse is dilated and delayed, due to the finite veloc-
the goals of this tutorial is to illustrate how the wavelet decompo- ity of propagation. The application of wavelets (which literally
sition is carried out, starting from the fundamentals, and how the translates from ondellets in French into English as small waves)
scaling functions and wavelets are generated from the filter-theory was first made in the area of geophysics [4], in 1980, by the French
perspective. Examples are presented to illustrate the class of prob- geophysicist J. Morlet, of Elf-Aquitane. A good history fi-om the
lems for which the discrete wavelet techniques are ideally suited. It mathematical perspective is available in the special issue of the
is interesting to note that it is not necessary to generate the wave- IEEE Proceedings [5].
lets or the scaling functions in order to implement the discrete
wavelet transform. Finally, it is shown how wavelet techniques can In electrical engineering [6-81, however, wavelets have been
be used to solve operatodmatrix equations. It is shown that the popular for some time, under the various names of multirate sam-
“orthogonal-transform property” of the discrete wavelet techniques pling, quadrature-mirror filters, and so on. Since the majority of the
does not hold in numerical computations. readers of this article are assumed to have an electrical-engineering
background, it will perhaps be useful to describe the methodology
in terms of filter theory. One of the objectives is to illustrate that
performing the discrete wavelet transform is equivalent to filtering
a signal by a band of constant-Q filters, the non-overlapping band-
widths of which differ by an octave. It is hoped that this mode of
‘This is part 1 of a two-part article. Part 2, which treats the presentation will make wavelets easier to visualize, conceptualize,
continuous case, will appear in the December issue. and apply to the problem at hand-if the wavelet theory is relevant!
IEEE Antennas and Propagation Magazine, Vol. 40, No. 5, October 1998 IEEE
1045-9243/98/$10.0001998 49
One of the goals of this paper is to illustrate how one can generate
the scaling functions and the wavelets, specially tailored to one's
needs.
Consider a signal, x ( t ) , the Fourier transform of which is

X ( w ) . In this paper, we deal only with discrete wavelet tech-
niques. Discrete wavelet techniques are quite suitable for discrete
signal processing, for example, in speech and image processing. In Transmitter Receiver
particular, their applications are very desirable in data compres-
sion. Since a complex matrix is a two-dimensional system, the Figure 1. The principles ok sub-band filtering.
solution of a matrix equation may be posed as an image-analysis
problem. As we shall observe, discrete wavelet techniques may be of two, these signals can be decimated by a factor of two without
suitable for the solution of large complex matrix equations. any aliasing. This is equivalent to reducing the sampling rate.
Decimation or down-sampling by a factor of two implies that
Continuous techniques, on the other hand, may be suitable altemate samples are dropped, and the data are compressed, as
for time-domain processing, where the wavelet transform can be shown in Appendix 1 (Section 8). The purpose of decimation is to
interpreted as a windowed Fourier transform. This we shall deal reduce the sampling rate and, thereby, the bandwidth of the signal.
with in the second part of the paper. This down-sampling i s possible because and ~ ' ( nhave
.I(.) ) an
effective bandwidth of f = 1I 4 or w = z1 2 , because they have
been filtered.
3. Development of the discrete wavelet methodology from
filter-theory concepts
Next, both ~ ' ( n and
) ~ ' ( n are
) down-sampled by a factor of
3.1 Preliminaries two, resulting in u(n) and v(n). This sub-sampling can continue
further, as we shall see later on. Here, we will restrict ourselves to
Consider the signal ~ ( tthat
) is discrete, so that it is repre- the two stages of filtering of x ( n ) by h ( n ) and g ( n ) , for illustra-
sented by the sequence tion purposes. Since both u(n) and v(n) have a smaller band-
width, they can be easily sampled, quantized, coded, and transmit-
x( .): Iz = 0,1,2,... (1) ted. In the receiver, the quantized, sampled, and coded )(. and
v(n) are received. The down-sampled versions, .(a) and ~ ( n ) ,
Then, its Fourier transform is best handled by the z transform (the can be transmitted at a much lower bit rate than the original signal,
lowercase letters are for functions in the original domain, and the without any loss of information, as they have a smaller bandwidth
uppercase letters are used for the z transform): than the original signal. The problem is how to reconstruct the sig-
nal x(n) back again, given the sub-sampled versions u(n) and
v(n) of x ( n ) ,without any aliasing or distortion.
Using this type of sub-band splitting has many advantages:

where w = 2 $ is the angular frequency, and z = e J o = e J 2 ~
Consider the sampled signal x(n) to be bandlimited, and assume it 1. This methodology results in a "maximally decimated" sig-
is sampled (say) at f = 1Hz. Thus, the sampling interval is nal (Le., some of the sample values of the signal have been deleted
At = 1 sec. From the Nyquist sampling criterion, it is then neces- or set to zero), namely, the sampling rate can be reduced, without
sary for the signal x ( n ) to be bandlimited to 112 Hz, so that it can any loss of information. [Note that v(n) of x(n) have a lower
be sampled at two times its bandwidth at the Nyquist frequency bandwidth than the original signal x ( n ) and, hence, the sample
without aliasing. The bandwidth of the signal (in angular fie- rate has been reduced by a factor of two]. This is equivalent to
qUency), mband > is saying that u(n) and v(n) have been decimated by a factor of two.
2. Even if .(E) and v(n) are quantized in a very rough fash-

ion (say, quantized into two bits, rather than the conventional eight
and the sampling angular frequency, w,a,,p, is bits or 16 bits), then the reconstruction is really remarkable, even
though such large quantization errors are introduced in the coding
=2z. (4) of u(n) and v(n) [3-81. This has been shown, at least in the area of
image processing. We will also demonstrate that this happens in
Now, let us filter the signal .(a) by a low-pass filter, H ( z ) , of the compression of matrices arising in electromagnetic-field prob-
bandwidth 7112 [Le., 0 5 w 5 z12 1, and by a high-pass filter, G(z) , lems. Our interest in this is due to the fact that a matrix is essen-
tially a discretized version of an image. So, can we apply this
of bandwidth w = n 12 to z.Let u'(n) be the low-pass-filtered methodology to the efficient solution of operator equations? We
signal [Le., ~ ' ( n has
) been obtained by passing x(n) through the will address this issue later on.
low-pass filter h ( n )1, and let ~ ' ( nbe
) the high-pass-filtered signal
3. Mother Nature, in many cases, performs processing in such
[i.e., ~ ' ( n has
) been obtained by passing x ( n ) through the high-
a way. For example, human ears and eyes, at least in the first stage
pass filter g ( n ) ] .This is shown in Figure 1. Now, since the band- of decoding sound and sight, perform processing by constant-Q
width of the signals ~ ' ( nand
) ~ ' ( nhas
) been reduced by a factor filters that are very closely related to wavelets.
50 /€€E Antennas and Propagation Magazine, Vol. 40, No. 5 , October 1998
Next, the signal is reconstructed from the decimated trans- H ( z ) = h(O)+ h(l)z-'+ ...+~ ( N ) z - ~ . (14)
mitted signals (.) and .(.). The signals are up-sampled by a
factor of two to produce ~ " ( n and) ~ " ( n )The
. principle of up- Without loss of generality and for convenience we choose
sampling is shown in Appendix 2 (Section 9). Then, they are
filtered by two receiving filters, g'(n) and h'(n) . The outputs are El'(.)= z-"H(z-') so that h ' ( n ) = h ( N - n ) . (15)
combined to form ."(n). Now, let us see how x(n) is related to its
estimate X(n) , and then the methodology to extract x(n) will be The factor Z-" is used to guarantee causality of the filters [Le.,
obvious. A'(.) = 0 for n < 01. The high-pass filter g'(n) is chosen in such a
way that
Please observe that
G ' ( z )= z - ~ G ( z - ' )with g ' ( n ) = g ( N - a ) . (16)

U ' ( 2 )= H ( z ) X ( z ) , (5)
V ' ( z )= G ( z ) X ( z ) . (6)
In addition, we define the high-pass filter g ' ( n ) in terms of the
From Appendix 1 (since ~ ( n and
) v(n) have been decimated by a low-pass filter coefficients h'(n) by choosing
factor of two)
g'(n)=-(-l)"h'(N -.)=-(-l)"h(n),
1
U ( Z ) = -[ut(&)
2 + V'(-&)], (7)
g(.) = (-l)"h(N- n) = (-l)%'(n),
V ( z ) = '2[ V ~ ( & ) + V ~ ( - & ) ] . so that
By using the results of Appendix 2 (where ~ " ( n and

) ~ " ( nhave
)
been up-sampled by a factor of two),
U"(z)=u ( z 2 ) = 21 [ u ~ ( z ) + z i ' ( - z ) ] ,
(9)
Substitution of Equations (18) into Equation (13) demonstrates that
all four equations are consistent, and the aliased component due to
V " ( z )= V ( z 2 )= ,[V'(z)+
1 V+)] (10) X(-z) is zero.
Furthermore, by utilizing Equations (15) and (16), Equation

Therefore,
(12) simplifies to
'[
k ( z ) = - z- N G ( ~ ) G ( ~ - z' )- "+H ( z ) H ( z - ' ) ] X ( z )
2
1
= -2- { [ G ( z ) X ( z +
) G(-z)X(-z)]G'(z) (19)
= -[H(-z-')H(-z)
1 + H(z)H(z-')]z-"X(z).
+[ H ( z ) x (2) + N(- z ) X ( -z)]H'( z ) } (1 1 ) 2
1 If the filter H ( z ) is chosen in such a way that

= -{[ ~+(H ( ~ ) H ~ ( ~.
G ( ~ ) G.) ))
]x(
2
+ [ G ( - z ) G ~ ( z+)H ( - ~ ) H ~ ( ~ ) ] x ( - ~ ) } , (12) H(-z)H(-z-') +H(z)H(z-') = 2,
where k(z), G ' ( z ) , H ' ( z ) , G ( z ) , H ( z ) , and X ( z ) are the z or
transforms of Z ( n ) , g ' ( n ) , h ' ( n ) , g(n), h ( n ) , and A(.). The z
transform has been defined by Equation (2).
The estimated signal X ( z ) contains the original signal

(which is given by the first term) and an aliased part (which is then one would have a perfect reconstruction property:
given by the second term of Equation (12)). Now, to remove the
aliasing effect, the second term must be zero, Le., qz) = z-NX(z); X(.) = x(. -N ); (22)
and the estimated signal is delayed by N samples.
Let H ( z ) be a FIR (finite-impulse-response) filter of order N + 1. The filters H ( z ) and G ( z ) are called quadrature mirror fil-
Then, h(n) will have N +1 terms. We consider N to always be ters (QMF), because they have symmetry around the point 7d2, as
odd. Then, shown in Figure 2. Ideally, we would prefer H ( z ) and G ( z ) to be
IEEE Antennas and Propagation Magazine, Vol. 40, No. 5,October 1998 51
So, if we solve only for H ( z ) , the other three filters of Figure 1
will then be given by Equations (24),(25), and (26).
Now, we show how to solve for H ( z ) . In this example, we

have N = 3 , and the order of the filter, L, is given by
L =order of filter= N +1= 4 . (27)
Utilizing Equation (20), we have
[h(O)- h(1)z-l + h(2)z-Z - h(3)z"I

Figure 2. The quadrature mirror filters. x[h(O)- h(1)z + h(2)z2 - h(3)z3] (28)
+[h(O)+ h(1)z-I + h(2)z-Z + h(3)Z-q
x [ ~ ( o )+ ~ ( I )+zh(2)z2+ h(3)z3]= 2 .
Equating all the coefficients related to the individual powers of z in

Equation (28) leads to
d n 2 7t '
0 h(O)h(2)+ h(l)h(3)= 0 (for z 2 and z - ~ ) . (30)
Figure 3. Non-overlapping ideal filters.

We need two more equations to solve for H ( z ) . Since G ( z ) is a
high-pass filter, then at w = 0,
non-overlapping, as shown in Figure 3. However, due to
realizability conditions, they are like those presented in Figure 2.
Moreover, we would like to have N ( z ) and G ( z ) be FIR filters
(finite impulse response), as opposed to IIR (infinite impulse
response), for ease of numerical computation. To illustrate the From Equation (23),this leads to
nature of the various FIR filters given by the various equations,
consider as an example N = 3 [9]. Then, the various filters H ( z ) , -h(3)+h(2)-h(l)+h(O)= 0
G ( z ), N ' ( z ) , and G ' ( z ) will have four nonzero coefficients in
their expansion. Or, equivalently, h(n), g ( n ), h'(n) , and g'(n) In addition, from Equations (21) and (3 1) at w = 0 we have
will have four entries. They will have the following form. If we
choose
N
H ( z )= z h ( n ) z - " + h(1)z-' + h ( 2 ) Y 2+ h ( 3 ) ~ ,- ~(23)
= h(0)
From Equations (23) and (33) we get
n=O
then from Equation (18), H(e") = H(1) = h(0)+ h(1) + h(2)+ h(3)= A . (34)
G ( ~=)- z - 3 ~ ( - z - l ) From Equations (32) and (34) we see that

= -z-3[q~)
- h(qz+I + h(2)z+2- h(3)z+3] 1
h(0)+ h(2)= -,
z/z (35)
= +q3) -~ ( Z ) ~ - I +h ( ~ )-~h -( ~~ ) ~ - ~ ,
and from Equation (15), 1

h(1)+ h(3) = -
a' (36)
H ' ( z ) = z-"H(z-')
In addition one has Equations (29) and (30). We need one more
= h(3) + h(2)z-1 + h ( ~ ) +~h -( ~~ ) ~ - ~ , equation, as these four equations [Equations (29), (30), (35), and
(36)] are linearly dependent, since Equations (35) and (36), in
and finally from Equation ( 1 8), conjunction with Equation (30), lead to Equation (29). The
question is how to find the fourth equation. Without the fourth
G ' ( z )= z-NG(z-*)= -H(-z) equation it is not possible to get the complete solution. Here, the
various methodologies differ, and different researchers have come
= -h(O) + h(1)z-I - h(2)z-Z + h(3)z-3. up with different procedures.
52 IEEE Antennas and Propagation Magazine, Vol. 40,No. 5, October 1998

For example, Daubechies [11 originally developed this meth-
odology by constraining the filter N ( z ) to be smooth, by enforcing
all derivatives to be zero at w = 0 , up to order p , or, equivalently,
G‘P(l) = 0 , where the superscript p denotes the pth derivative of
G’ . The actual value for p is determined from the number of equa-
tions needed to solve for the values of h(n), n = OJ,. ..,A’. This
leads to taking the various moments of g’(n) and setting them I
equal to zero. Daubechies chose this procedure because enforcing I
the above conditions guarantees smooth wavelets, which we will

define later. However, for the present case, where the number of
discrete wavelet coefficients is finite, the smoothness of the wave-
lets is a moot point. This is true since for the discrete case, the
wavelet transform can be implemented (for the examples that we
-150 -
are interested in) without explicitly specifying the wavelets. The
smoothness of the wavelets at this point, then, is of no concern to
-200 I
us! Another approach, for example, is given in [lo, 111. 0 0.5 1 1.5 2 2.5 3 3.5
-w
In mathematical terms, setting the derivative of G’(1) to zero Figure 4b. The phase responses of the fourth-order filters.
leads to
-Oh(O) - Ih(1) + 2h(2) - 3h(3) = 0 . (37)

center is proportional to f i . The transition from jH(z)\= 0.98&
Solution of Equations (30), (35), (36), and (37) leads to to jH(z)l= 0.02& is over an interval of length 4/&.
I+& Instead of having a two-stage decomposition of the signal

h‘(0)= h(3) = __
4Jz ’ x ( n ) into ~ ’ ( nand
) ~ ’ ( n, one
) can perform a multistage decom-
position by applying the filters successively to each stage of the
3+& decomposition, as shown in Figure 5. The electrical-filter-theory
h‘(l) = h(2)= _ _
4Jz ’ equivalent is shown in Figure 6. This is the decomposition part (or
the part labeled “transmitter” in Figure I). The coefficients at the
3-A output of Figure 5 are then thresholded. This is equivalent to
h’(2)= h ( l )= _ _ keeping only those coefficients that are bigger in magnitude than
4Jz ’ some constant E . The values smaller than E are set equal to zero,
resulting in the approximate filter output. Now, the interesting part
1-43 is that these “filtered” thresholded coefficients can now be used to
h‘(3)= h(0) = __
44‘5 recover the original signal with an accuracy better than E (!). (We
will see this feature later on). The reconstruction algorithm is
The magnitude and the phase responses of H ( z ) , H‘(z) , G(z) , depicted in Figure 7 , where the approximated coefficients are up-
G‘(z) are shown in Figures 4a and 4b, respectively. Note that sampled, and then filtered the way that is depicted in the receiver
of Figure 1.
these filters have no ripples, with many zeros at n.The slope at the
The type of decomposition outlined in Figure 5 is identical to
a wavelet decomposition, as the next section will illustrate. For
example, filtering a function by h(n) is equivalent to fitting a
scaling function at a certain scale, and filtering by g ( n ) is equiva-
lent to curve fitting x ( n ) by wavelets at the same scale as the
scaling function. The mathematical connection is now established
between wavelet theory and filter theory.
3.2 Connection between the filter theory and the mathematical

theory of wavelets
To establish a connection between the filter theory and the

mathematical theory of wavelets, we will deal with functions of a
real continuous variable. Discrete samples then will be viewed as a
“sample and hold” of the these functions of a continuous variable.
0 0.5 1 1.5 2 2.5

.- 3 3.5
Such an approach is absolutely necessary, because the wavelets
cannot be expressed in terms of discrete samples. Moreover, the
dilation equation, which is at the heart of wavelet theory, has no
-w
solutions for the discrete samples. So, the discrete wavelet trans-
Figure 4a. The magnitude responses of the fourth-order filters. form implies that the functions we are dealing with are functions of
/€€€Antennas and Propagation Magazine, Vol. 40, No. 5, October 1998 53

Figure 5. A multistage decomposition of the signal .(a). Figure 7. The reconstruction of the original signal, x ( n ) ,from
the approximate coefficients of the discrete wavelet transform.
a continuous variable. They can have integer shifts. In addition,

they can be scaled up and down by integer multiples.
We consider the original signal, x ( t ) , with bandwidth from 0

to n.Also, consider some functions #(t). We assume that the inte-
ger shifts m of &), namely #(t - m ) , are a Riesz basis for the
original space [38]. Let us denote
I 7-1
Hence, one can represent a dilated version of &)-namely,
&2)-by a combination of the functions 4(t - m ) [or dm ( t ) ]with
some coefficients h'(nz), resulting in the dilation equation
4 4
or, equivalently, [with h'(n) = h(N - m) , as defined by Equa-

tion (15 ) ]
N N
4(t)=~ C h f ( m ) 4 ( 2 1 - m ) = 1 / 2 C h ( N - m ) ~ ( 2 t - m )(41)
.
m=O nz=o
It is interesting to note that the coefficients h'(m) tum out to be the

same coefficients that we have described earlier in the receiver part
of Figure 1 and Figure 7.
As an example, consider 4(t) to be a pulse function of mag-

nitude 1, located between 0 and 1, as shown in Figure 8. Then,
4(t/2) has support from 0 to 2, or is a dilated version of $ ( t ) .
Please note that 4(t) can be represented by 4(2t) and 4(2t - I),
x(n) where each of these functions is defined between 0 to 112 and 112
to 1, respectively. In this case, h'(0)= h'(1) = I/&. Equation (41)
essentially reflects the basis of the dilation equation: Le., the func-
tion can be approximated by a weighted sum [through h'(m)] of
the shifted and dilated versions of the same function.
The fimction d ( t ) , which solves the dilation equation for a

particular h ' ( m ) , i s called the scaling function, and is also called
the father of wavelets. Also, note that the functions #(t - m) are
assumed to be normalized, Le.,
Figure 6. A filter-theory representation of the discrete wavelet I 4(t - m)dt = 1

transform. I
54 IEEE Antennas and Propagation Magazine, Vol. 40, No. 5, October 1998
Here, " d t "can be represented as the incremental interval length.
Therefore, from Equation (41), utilizing Equation (42), we get 4(t)= f i [ * 2 - ' / 2 h ' [ 2 k t ) ]where
, * denotes a convolution. Note
k=l
that the sequence of convolutions is carried out by various com-
j C ( t - k)dt = 1 =
I
Jic m=O
hQ)j
1
4(2t-2k-m)
2
JiN
d(2t) = - h'(m)
2 m=O
pressed versions of the same signal. When the sequence of convo-
lutions converges, it yields the function # ( t ) .
or, equivalently, When N = 0 , the dilation equation becomes, for hb = &,
nz=O and the scaling function becomes a delta function. If we choose

L = N + 1 = 2 and
Hence, the coefficients h'(m), satisfying the dilation equation,
must satisfy Equation (44). This is the same equation as Equa-
tion (34), restricted to N = 3 terms. The corresponding wavelets
are given by (compare with Equation (17))
Hence, in this case, the dilation equation is &) = 4(2t) + 4(2t - 1).
One possible solution of this dilation equation is the pulse func-
m1=0 tion:
(45)
N
= Jz C-(-l)'%(N - m)4(2t - m) , 1ifO<t<l
nz=O 0 otherwise
where N + 1 is the order of the filter (which is always even). Here, The wavelet is generated from Equation (45), and is given by
the filter coefficients h'(m) and g'("z) have been assumed to be
, a given g ' ( m ), is called the mother of
real. The function ~ ( t )for -1 f o r O < t < l / 2
wavelets. (This might be viewed as the traditional Judeo Christian
concept of mother: where the mother is generated from the father!)
0 otherwise
So, from an electrical-engineering perspective, if we have the
filter coefficients h ' ( m ) , then the scaling function 4(t) can be since y/(t) = -4(2t)+ 4(2t - 1). This is shown in Figure 8. This is
obtained by iteratively solving the dilation Equation (41). This is the pulse doublet, and it is related to the Haar wavelet. Please note
carried out by starting with letting #(t) be a pulse function, and that the scaling function @(t)is in this case orthogonal with respect
then filtering by h ' ( m ) , and continuing until convergence is to its own translates, i.e.,
reached. Once the scaling function is known, the wavelets are a)
given by Equation (45).
p(t)+(t - m)dt = 0. (52)
-m
The scaling function can also be solved for in the transform
domain. By taking the Fourier transform of Equation (41), one can
This is also true for the wavelets:
write the dilation equation in the w domain as
m
A
= --N
$(a>
I
) $(w/2),
-JW/2
(e (46) j ~ ( t ) l y (-
t m)dt = 0 (53)
2 -m
where i ( w ) is the Fourier transform of &), and H ' ( z ) is the z For N = 3 , we obtain the case presented earlier. Note that for this
transform of h'(n) (see Equation (15)). One can repeatedly apply case, the filter coefficients hi are given by Equation (38). Hence,
Equation (46) to obtain
Owing to Equation (42), $(O) = 1. Therefore, in the limit n + ca,

Equation (47) becomes
The solution for $(a)is point-wise convergent, provided the infi-

nite series converges. In the original domain, this is equivalent to Figure 8. An example of the dilation equation for A4 = 2 .
the scaling function can be obtained from the solution of the The confusing part here is what to use as the starting point for
dilation equation, a scientific endeavor to carry out a numerical analysis utilizing
wavelets? There are two choices:
&)= $[(1+.11)~(21)+(3+.11)~(2t-l)
1. Do we start with + ( t ) ,construct h ' ( n ) ,and then generate y ( t )
+ (3 - &)&t - 2)+ (1 -.11)4(2t - 3)]. (54a) [from Equation (45)]?; or
The wavelets are generated in an analogous fashion by using 2 . Do we start with h ' ( n ) , and then create $(t) and ~ ( t? )
Equations (45) and (38), resulting in
For the discrete case that we arc dealing with, the answer is
v(t)=$[(1-&)4(2t)-(3-&j$(2t-l) straightforward: that is, we design h'(n), and then obtain # ( t ) and
~ ( t. This
) is also much simpler in practice. However, for the dis-
+ (3 + &)4(2t - 2)- (1 + &)4(2t - 3)]. (54b) crete wavelet transform, as we shall see, 4(t) and ~ ( t are
) not at
all required in the numerical computation!
We are not going to delve further into the solution of the dilation
equation, since for the discrete wavelet transform and for the
In summary, the mathematical basis of the wavelets has been
problems that we are interested in electromagnetics-namely, solu-
presented from a filter-theory perspective. How to construct scal-
tions of large matrix equations-the scaling functions and the
ing €unctions 4 and wavelets y has been shown, starting from the
wavelets are really not necessary. This is because the discrete
wavelet representation can be carried out from the knowledge of filters h'(m) , and utilizing the perfect-reconstruction argument
only h(m)! However, we present some other wavelets for illustra- presented with subband-filtering techniques. Once h'(m) i s known,
tion purposes. 4 can be generated from Equation (41), and v , from Equa-
tion (45).
For example, consider the Shannon wavelet, which is the
dual of the Haar wavelet. The scaling function is given by
4. Approximation of a function by wavelets
Consider a function ~ ( t. The

) objective is to approximate it
by the wavelets ~ ~ , ~Le.,
( t ) ,
and its transform is given by
i(4= { 1 for o</w(<z

0 otherwise
(55)
where we define the wavelets by

The wavelet can be generated from Equation (45), and is given by
yn&) = 2-"'2y(2-"t - k ) . (56)
y ( t ) = s- ipn-a, 3nt
2
Note that v[2-"t) represents a dilated version of ~ ( t )For
.
and its transform is n = -1, we say the scale is the finest. This is because for n > -1,
the function gets dilated and becomes wider. ~ ( t - k represents
)
the shifts. Therefore, we approximate the function ~ ( tby
) a dilated
and shifted version of the same function ( ~ ( .t We
) assume that we
are dealing with orthogonal wavelets, hence
Lagrange half-band filters are used when one needs a non-negative
frequency response of the filters h ( n ) . A closed-form expression
was given by Ansari [12], starting from
H ' ( z )= - c
1 + ' h;(2n - 1)[z-**+' + 2 - 1 1 , The scale factor 2?12, appearing in Equation (56), is there so as to
make the functions ~ ~ , ~orthonormal,
( t ) Le.,
2 n=l
with 4i - 1 coefficients. The coefficients hi are determined using

the Lagrange interpolation formula:
It is interesting to note that ~ ~ , ~has ) DC value
( tzero
21
(-l)n+J-'n(i- k + 1/ 2 ) ( = 0) = 0 ,whereas x ( t ) may not! Therefore, if we have to
h;(2n - 1) = k=l talk about convergence of the series in Equation (55), we can only
(i - n)!(i - 1+ n)!(2n- 1) '
talk about mean-square convergence. This dichotomy, however,
does not arise when the sum in Equation (55) is finite. Utilizing
These filters are very regular, as they have a 2i-fold zero at z = -1 . Equations (58) and (55), we observe that
and these are the discrete wavelet coefficients of the function ~ ( t. ) and from Equation (45),
These are the same values shown in Figure 5. Our objective in this
section is to establish that isomorphism.
Now, if we have to carry out the inner products in Equa-

tion (59), it will be extremely time consuming, as the inner prod- It is interesting to note that a byproduct of Equations (66) and (67)
ucts have to be carried out for all values of n and k. Here, the is that
strength of the wavelet techniques comes in, as they provide a fast
and accurate way to recursively evaluate the inner products. This is
accomplished through the introduction of the scaling functions
) 2-"24(2-'t - j ) .
@ l , j ( t= (60)
We further assumeiutilize the orthogonality relationships between which can be generalized to

the scaling functions and the wavelets, and between the scaling
functions themselves, i.e.,
and Next, observe that
m m
2 ( n )= j p ( t ) @ (-
t n)dt = j p ( t ) ChF(t)&@(2t - 2n - k)dt
-m -m k
We now define the coefficients q nthrough
(68)
= Ch'(k)a(-')(2n+k)=Co(-l)(k)h'(k-2n),
k k
and
It is now shown how the dk,n are evaluated recursively through
the C k , n .
bo(.)= Ca(-l)(k)g'(k-2?2).
CI
We have, from Equations (58), (61), and (62), that the fol-
lowing orthonormal set Given the existence of relationships like Equations (68) and (69),
and drawing the isomorphism between
represents a basis, as the functions involved in the set are orthogo-

nal. From the dilation equation, Equation (41), and from Equation
(45), we note that
we can generalize the expressions of Equations (68) and (69) to
(&5(2t - k)Im
k=-m
~ E+~ ~, , ~ h ' ( n - - 2 k )C
c ~ ,= = c n , , h ( N + 2 k - n ) for j 2 - 1 , (70a)
n n
also form an orthonormal set for the same space. This is for scale
n = -1, as defined by Equation (56). Therefore, we can expand any
function p ( t ) by
These results have been derived utilizing Equations (1 5) and (16).

The above recursive relations show that we need to compute the
inner product of Equation (63) at the highest scale, n = -1, only
once (instead of using Equation (59)), and then the wavelet coeffi-
cients dk,n of the XDwr(n,k ) ( n = 0,1,2,. ..) are computed recur-
sively from Equations (70a) and (70b). From a filter-theory point
Here, u(-')(k) are the coefficients used to represent the function
of view, Equations (70a) and (70b) show that the c , , ~need to be
p ( t ) by @(2t- k ) . The superscriptsj on the coefficients d ( k ) and convolved with h(n) and g(n) , then down-sampled by a factor of
b j ( k ) represent the value of the scale j at which they are two. Through Figure 5 and from the above development, it is clear
represented. We have from Equation (41), that for the computation of the discrete wavelet transform (DWT),
/€E€ Antennas and Propagation Magazine, Vol. 40, No. 5,October 1998 57
it is not necessary to even know what the scaling functions and
wavelets are, as one can directly use Equations (70a) and (70b)
without going through the mathematical derivations, as Figures 5
and 6 illustrate. One starts with c ~ , - which
~ , is the coefficient gen-
erated by correlating 4(2t) with the function ~ ( ,t and
) then recur-
sively computes the discrete wavelet transform mathematically
through Equations (70a) and (70b), and graphically using Figure 5 ,
which is easier to visualize from a filter-theory perspective. The
methodology is the same. The process described so far is similar to
the transmitter part as labeled in Figure 1. If the 4(2t) are the
impulse scaling functions, then q - 1 will be equivalent to the
sampled version of ~ ( t namely
), x(n).
There is another subtle point that one should introduce now!

So far, in the approximation in Equation (55), the limits are infin-
ity. This is good from a mathematical perspective. However, from
a practical reality, the limits have to be finite. Hence, in addition to
dk,n, we also have c ~ , ,The
~ . approximation of ~ ( tfrom) a practi-
cal standpoint is done by
Figure 9b. A compressed version of the picture of “Lena,”

utilizing QMF compression (40:l compression).
Hence, we have both wavelets and scaling functions in the
approximation. This very important feature is generally not clearly
delineated in many presentations. Now, observe that Equation (71)
is the representation of ~ ( t. )
In summary, the evaluation of the discrete wavelet coeffi-

cients in Equation ( 5 5 ) is equivalent to filtering the coefficient
c ~ , -obtained
~ from ~ ( t(or,
) equivalently, the sampled values of
~ ( t for
) a certain class of scaling functions) by a cascade of mutu-
ally orthogonal filters, as shown in Figure 5. The bandwidth of the
filters is reduced by a factor of two as one goes towards the DC
value. For the infinite sum in Equation ( 5 9 , the scaling function
does not enter into the final sum, except in the intermediate com-
Figure 9c. A compressed version of the picture of “Lena,” util-

izing JPEG compression (40:l compression).
putations. However, if the sum is finite in Equation ( 5 5 ) , then one

obtains Equation (71), and the scaling functions are needed in the
summation.
In practice, once the coefficients d,c,n and q m have been

obtained, those coefficients with magnitudes below a certain
threshold value of E ( e g , E = are set equal to zero. At this
point, the transmitter then sends out the approximate coefficients
-
dk,n and c k , M , which are the values obtained after thresholding
dk,n and q mby E . Now the problem is, how does one recover
Figure 9a. The original picture of “Lena.” ~ ( tfrom
) these approximate coefficients in a fast efficient way?
By utilizing Equations (66) and (67) in Equation (65b), we bits are utilized to code the compressed image, the reconstruction
is still better than the conventional JPEG algorithm [21].
As an example, consider the original image of “Lena,” shown

in Figure 9a. The images in Figures 9b and 9c show the result of a
40: 1 compression of “Lena,” using wavelet techniques utilizing
both signal-dependent QMF filters and JPEG compression. (JPEG
is the current standard for image compression, and the coefficients
of the filters h(m) are determined from the signal). Here,
(72)
compression refers to the total number of bits required to store the
image, as opposed to the total number of bits of the original image.
Even though there is noticeable degradation in both of the images,
the two methods perform reasonably in the reconstruction of the
image. This is due to the fact that the image is compressed by
blocks. The image compressed utilizing the signal-dependent QMF
decomposition looses detail in the local (i.e., high-frequency) edge
information. In particular, notice that the sharpness in the eyes and
By equating Equation (72) to Equation (65a), we find
the detail in the feather cap are blurred. Unlike the JPEG-com-
pressed image, however, there are no objectionable artifacts, such
as the “blockiness” mentioned earlier.
(73) As another example, consider the application to halftone

images. Unlike continuous-tone images, such as photographs, a
= C[h ( N + 2n -k)#)(n) + g ( N + 2n - k)h(o)(n)] halftone image is digital in nature: all pixels are either on or off.
These digital pixels are densely placed on the paper (as many as
/I
2400 per inch, to emulate high-quality magazine halftones). Half-
tone images are produced by photographing a continuous-tone
Hence, we start with the approximate coefficients Ck,M and d”k,,, ,
image through a very fine screen, typically having 100-133 lines
and we then recursively generate, through the use of a generaliza- per inch. They are also often simulated-with lower resolution and
tion of Equation (73), the coefficients quality-on common digital printers, including laser, ink-jet, and
dot-matrix printers. It has become common for such images to be
scanned into a computer using an 8-30 bitsipixel optical-input
scanner, with the resulting image stored in the computer. Such a
digital simulation of a typical halftone image is shown in Fig-
Note that the are the estimates of the discrete values of x ( n ) ure 10a (with a magnification of three for demonstration purposes).
for the impulse scaling functions. Also, observe that Equation (74) Notice the dot pattern of the background, which is not present in
is equivalent to Figure 7, from a filter-theory perspective, and is most continuous-tone images. This image was printed by a 400
identical to the receiver point in Figure 1. The most remarkable dots-per-inch printer. Figures 10b and 10c show the result of a 30:l
point here is that even though the original coefficients have been compression of Figure 1Oa, using both the signal-dependent QMF
thresholded by c to produce the approximate coefficients, filters and the JPEG compression. Notice that the QMF compres-
sion attenuates the dot pattern associated with the original. This is
and CT,‘,,,
, the original sampled function x(n) can be reconstructed because the high detail content of the dot pattern is lost in the QMF
through Y ( n ) with an accuracy better than E . This we will illus- compression. In this instance, attenuation has become an asset,
trate through numerical examples in the next section. actually improving the appearance of the image. The JPEG-com-
pressed image seems to produce a noise pattern in the output. This
I n summary, in order to implement aiid carry out the discrete is indicative of the problems faced when using JPEG compression
wavelet transform, it is not even necessary to introduce the concept on halftone images.
of scaling functions and wavelets. The filter-theory approach
essentially provides the same methodology, in a simpler, practical By utilizing a wavelet technique, it is possible to quantize
fashion. images with bit rates as low as 0.4 bits per pixel, while maintaining
a sufficiently high quality of reconstruction [9].
As an illustration of how to utilize the discrete wavelet trans-
form, we consider the compression of an image. For the case of
images, one is dealing with the two-dimensional discrete wavelet 5. Relevance of discrete-wavelet techniques in computational
transform, which is a generalization of the one-dimensional case. electromagnetics
In the following examples, we consider the image to be a There are essentially two different ways to solve operator
two-dimensional array. We take the discrete wavelet transform in equations utilizing wavelet techniques. The first possibility is to
two dimensions by utilizing a recursive relationship similar to use the wavelets as basis functions, in a conventional Method-of-
Equations (70a) and (70b). We utilize an eighth-order filter that has Moments formulation. The second approach is to use traditional
been designed to match the signal [21]. Once the discrete wavelet sub-sectional basis functions, and obtain a dense complex matrix.
coefficients are obtained, they are thresholded, and then the origi- The wavelet techniques are then used to compress the elements of
nal image is reconstructed utilizing the recursive relation of Equa- the matrix, and either a direct sparse solver or an iterative method
tion (74). The objective is to illustrate that even though only a few is used to obtain the solution of the sparse-matrix system. The
IEEE Antennas and Propagation Magazine, Vol. 40,No. 5,October 1998 59

basic philosophy of the two methodologies is the same: namely, to
obtain a sparse complex matrix instead of a full one, as results
utilizing MOM and sub-sectional basis functions. However, the
advantage of forming a sparse matrix is offset by the problem that
the condition number of the transformed matrix may be worse.
This is particularly significant for the second method, where the
transformation utilizing a wavelet-like basis, from a complex-full
matrix to a sparse-complex matrix, is orthogonal, from a strictly
mathematical point of view. However, from a purely numerical
perspective, as will be shown later, the results do not confirm this
orthogonal transformation, as the condition number of the system
changes and, in some cases, the change is by an order of magni-
tude.
5. 1 Application of wavelet basis for the solution of operator

equations
In [22], a Galerkin method, utilizing a wavelet basis, has

been used for the analysis of radiation for TM-scattering problems.
There, a method is also proposed to compress the impedance
matrix, utilizing wavelet techniques. This is accomplished through
a compression process in which only the significant terms in the
expansion of the (yet unknown) current are retained and subse- Figure lob. A compressed version of the simulated halftone
quently derived. picture of Figure loa, utilizing QMF 30:l compression.
Wavelet bases have also been used to solve the Fredholm

integral equation of the first kind, which leads to the TM-scattering
from conducting cylinders [23]. Here, the compactly supported
semi-orthogonal bases have been used. The wavelets have been
specially constructed for the bounded interval for solving first-kind
integral equations. It has been observed that the use of cubic-spline
wavelets almost diagonalizes the matrix. Explicit closed-form
polynomial representations for the scaling functions and wavelets
are given.
In [24], the wavelet-expansion method, in combination with

the boundary-element method (BEM), has been used to solve the
integral equation for the surface currents. The unknown current is
Figure 10c. A compressed version of the simulated halftone

picture of Figure loa, utilizing JPEG 30:l compression.
expanded in terms of a basis derived from a periodic, orthogonal

wavelet in a finite interval. Because the geometrical representation
of the BEM is employed to establish the mapping between the
curved computational domain and the interval [0,1], it would be
very interesting to find out how this can be extended to three-
dimensional problems, where an arbitrary surface needs to be dis-
cretized.
An adaptive multi-scale Moment Method is presented in [25], for

the solution of the Fredholm integral equation of the first kind. An
Figure loa. The original, simulated halftone picture printed on adaptive procedure is outlined, which refines the unknown solution
a 400 dpi digital printer. in regions utilizing multi-scale wavelet-like functions. Even though
68 IEEEAntennas and Propagation Magazine, Vol. 40,No. 5,October 1998

the matrix is highly sparse, the deterioration in the condition num- ent, or even orthogonal! This point was illustrated by [13-171,
ber of the matrix-as opposed to the condition number of the origi- through the solution of the following differential equation:
nal matrix utilizing a conventional sub-sectional basis-is clearly
visible. One disturbing fact is that the results are not consistent
d2Y
-=sinx+2 for 0 1 x < 2 z ,
with the increase of the order of the filter, nor with the level of (78)
dx2
thresholding used. The methodology can also be used for solution
of the differential form of Maxwell’s equations, utilizing finite
elements [26].Similar wavelet-like functions can also be used to which has a solution y = -sinx + x(x - 2 z ) when the boundary
compress the elements of the impedance matrix, utilizing a thresh- conditions are y ( x = 0) = 0 and y(x = 2 n ) = 0 .
olding operation.
A particular choice of the basis functions e, for the unknown
In [27], the problem of electromagnetic scattering from per- y could be the Fourier series, i.e.,
fectly conducting strips, coated with thin dielectric material, is
analyzed, utilizing an adaptive multi-scale Moment Method. The
method of a non-uniform grid and the multi-scale technique, which y = a0 + cm
n=l
(a, cos nx + b, sin nx) , (79)
generates a locally finer grid, are usually used when the solutions
of the integral equations or the differential equations are known to
vary widely in different domains. By non-uniform gridding, one in which the Fourier bases are linearly independent and orthogonal.
can reduce the size of the problem and improve the accuracy. The This yields a solution y = -sinx . Hence, it does not provide the
multi-level or the multi-grid technique has been widely used in d2Y formed
solving differential equations and integral equations [28-321. correction solution. This is because the series for -,
dx
Kalbasi and Demarest [33, 361 applied multilevel concepts to solve after the double differentiation of Equation (79), is not complete,
the integral equation by the Moment Method on different levels, as the a. term is missing. So, Ae, does not form a complete set
which has been called the Multilevel Moment Method. No matter [13-171. Hence, Equation (79) does not and cannot provide the cor-
what the multi-grid technique is, the basis functions for an rect solution. To rectify the error, one should choose the basis
improved approximation have to be reconstructed again. By using functions as [ 15,161
the multi-scale technique in one dimension, the basis functions for
the new scale have to be reconstructed. The new approximation
grids formed by the multi-scale technique are the same as those for
the multilevel technique; however, the constructions for the func-
y = co + c p + c2x2 + c
m
U=1
(a, cos nx + b,, sin nx) . (80)
tions are different.
These basis functions are not even linearly independent in the
In addition to these numerical concems, there are some philo- interval 0 to 2z,as both x and x2 can be approximated by the
sophical concems that need to be addressed in the selection of an remainder of the h c t i o n s representing the Fourier series. How-
appropriate set of basis functions. Consider the solution of an
operator equation ever, if we form __,
d2Y we obtain
dx2
AX=Y, (75) m
do + (d, cos nx + e, sin nx) ,

where, in general, A is a known integro-differential operator, and n=l
X is the unknown to be solved for a given excitation, Y. For the
solution of boundary-value problems, the most widely used tech- (with do = 2c2, d, = -n2a,, e, = n2b,,), which, indeed, form a
niques, like MOM, FEM, and BEM, convert the functional equa-
complete set. Utilizing Equation (76), it is then possible to get the
tion to a matrix equation. The matrix equation is then solved for
exact solution. This example clearly illustrates that the choice of
some unknown constants, instead of obtaining the solution in terms
the bases is important. Moreover, it is not necessary for the basis
of unknown functions. The solution procedure starts by expanding
functions e, to be orthogonal or even complete. What is necessary
the unknown, X, in terms of known basis functions e , ( t ) , with
is that Ae, must form an orthogonal set [13-171. Hence, it is nec-
some unknown constant multipliers a, in front, i.e.,
essary to address the question, how does the set Ae, form a com-
plete set when one chooses a wavelet-like basis?!
N
~ ( tz )C a l e i ( t ) .
,=I For the wavelet methodology to be successful, the wavelet
techniques have to ensure that the wavelet decomposition is carried
We then substitute Equation (76) into Equation (75), and form the out for Ae, , and not for e , . This is perhaps extremely difficult to
error or the residual (R), defined by do. Also, it needs to be rigorously shown that all the nice proper-
ties of perfect reconstruction, using the wavelets described earlier,
remains valid after a couple of differentiations of the series, pre-
sumably carried out term by term. Another hurdle is how to extend
the one-dimensional technique to two-dimensional and three-
dimensional problems. These are open problems.
Then, the objective is to make the residuals zero with respect to
some weighting functions, W j . From Equation (77), it is quite
In this first Section 5.1, we described how to use the wavelets
clear that Aej must approximate Y in some sense. Hence, it is as basis functions in the conventional MOM problems leading to
required that Aei must he linearly independent, and must form a sparse matrices. Also, in [26],how to use this approach for imple-
complete set. It is immaterial whether the ei are linearly independ- mentation of the differential form of Maxwell’s equations utilizing
IEEE Antennas and Propagation Magazine, Vol. 40, No. 5, October 1998 61
finite elements has been illustrated. Initial results illustrate how to g‘(n) (see Equation (17) for the relationship between h‘(n) and
use the wavelet techniques to generate a sparse matrix, utilizing a g‘(n)). Therefore, multiplying a vector (say [Y]) by [PI is
wavelet-like basis. A similar methodology can also be used in equivalent to filtering (in the transmitted portion) in Figure 1, fol-
making sparse matrices-arising in the implementation of finite lowed by a sub-sampling of two. The matrix-vector product yields
elements in the solution of differential forms of Maxwell’s equa- u(n) and v(n) . The first four elements are equivalent to u(n) and
tions-still sparser.
the last four are v(n). The matnx-vector product has already
In the next Section 5.2, how to use the discrete wavelet tech- incorporated the sub-sampling by a factor of two. The sub-sam-
niques to transform a complex dense matrix, arising in conven- pling by a factor of two is accomplished by the shift between the
tional MOM, into a sparse matrix, utilizing a set of orthogonal elements of each row of the matrix [PI. Now, for [PI to be an
transformations, is illustrated. Some work [35-381 has already been orthogonal matrix, it is necessary that the following three equa-
done to address this topic. This is described next. tions, similar to Equation (29) [the normalization constant is set to
unity] and Equation (30) [the filter is orthogonal to its two-shifted
version] hold:
5.2 Solution of large matrix equations by the discrete wavelet
transform C h 2 ( i )= 1,
1
Consider the solution of a matrix equation [A][X]=[Y],

where [A] is a known Q x Q matrix, [Y] is a known Q x 1 vector, h(O)h(2)+ h(l)h(3)+ h(2)h(4)+ h(3)h(5)= 0 , (83)
and [XI is the unknown to be solved for. Note that Q has to be an
integer power of two for the wavelet techniques to be applicable. If h(O)h(4)+ h(l)h(5)= 0 . (84)
the original matrix [A] is not of size 2m, then the matrix [A] can
Finally, from the boundary conditions for the filter,
be augmented by a diagonal identity matrix to make it 2m.First, H ( z = I ) = & and G(z =1)= 0, we have
we illustrate how to obtain the wavelet transform of a vector [Y],
and then we will illustrate how to take the wavelet transforms of a 1
h(0)+ h(2)+ h(4)= -= h(1)+ h(3)+ h(5) . (85)
matrix [A]. We first select h‘(n).We then use h’(n) and g‘(n) to Jz
find the discrete wavelet coefficients, utilizing Equations (70a) and
(70b). What is going to be different is that we are going to express Equations (82)-(85) provide four independent equations, and one
Equations (70a) and (70b) as circular correlations with respect to needs two more to uniquely solve for the h ( i ) s . If one follows
h’(n) and gf(n), or as circular convolutions with respect to h(n) Daubechies’ procedure for making the wavelets smooth, one would
and g ( n ) , respectively. This is a computationally efficient way to need the derivatives of G’P(z= 1) = 0 for p = 1 and p = 2 , lead-
carry out convolution by utilizing the FFT. As an example, we first ing to (from Equation (18))
create the following 8 x 8 orthogonal matrix [PI as
-0-h(O)+lh(l)+2h(2)+3h(3)-4h(4)+5h(5)=0 , (86)
-O*h(O)+lh(l)-4h(2)+9h(3)-16h(4)+25h(5)= 0 . (87)
The setting of the higher-order moments of g’(n) to zero guaran-

tees the smoothness of the wavelets. This, in tun, provides a recipe
so that the discrete-wavelet coefficients of the transform drop off
rapidly, as one goes to the dilated scales from a fine scale. The
number of zeros of G’P(z) at z = 1 tells us how many basis func-
tions are needed in Equation (55) for approximating ~ ( t )The .
smoother the function and the higher the order of zeros, the faster
the expansion coefficients go to zero, and the fewer coefficients we
need to keep. For piecewise functions, a wavelet basis is better.
These piecewise functions may have jumps. They may be smooth
and then suddenly rough. We keep more Coefficients in the rough
neighborhoods by going to a smaller scale T J .The mesh adapts
to ~ ( tin) a way that Fourier methodology finds difficult.
I J
0 0 -h, +h, -h2 +h3 4 4 +h5
-h4 +A5 0 0 -h, +h, -h, +h3 If x ( t ) has p derivatives, its wavelet coefficients decay like
-h2 +h3 -h4 +hs 0 0 -h, +hi 2-np [38]:
[In Equation (81), we have used hk to represent h(k) , to conserve

space]. This equation is the matrix form of the discrete wavelet
transform of Equations (70a) and (70b). Inside [PI we have the
filter coefficients, six in number. Please note that the first four where J is a constant, and &’)(t) represents the pth derivative of
rows are due to the filters h’(n), and the last four rows are due to 4) ’
However, since we are dealing with a finite number of terms,
the drop-off rate of the wavelet coefficients is of little significance.
The solution of the above equations can be obtained analytically,
and has been given by Daubechies [ 11 as
l + f i + J X
h(5) = h’(0) = -
16&
[p2116x16[y~~16x1 [y~116x1
wavelet coefficients
5 +f i +3 & Z T
h(4) = h’(1)= -
16&
[[y,’l1+ scaling function contribution
10 - 2 4 1 3 + 2 4 5 + 2&0
h(3) = h’(2) = -
16&
Figure 11. The principles of the wavelet transform, applied to a

10 - 2 f i - 2 J G G matrix.
h(2) = h’(3) = - (884
16&
where
Ck,2 for k = 1,...,4

1+&&6Z7iT dk.4,2 fork = 5 ,...,8
h(0)= h’(5)= - T ( k )= (89)
16& dk-8,1 f o r k = 9,..., 16 ‘
dk-lG,O fork = 17,. ..,32

Once the filters are available, one can form the orthogonal matrix
[PI by substituting in the above values of h. Let us say we have a The coefficients c k , 2 are due to the scaling functions & k , and the
vector [Yl], of length Q = 25 = 32. We now know how to extend remainder of the coefficients, dk,, ( i = 0,1,2) are due to wavelets
the wavelet transform to [ Y ] ] We
. create a [PI] matrix, which is V I & Note
. that in this case, even though we have carried out a dis-
32 x 3 2 , and only six elements of any of its rows are populated by crete wavelet transform, we do not need to h o w anything about
the elements of Equation (88); and a matrix of size 32 x 32, similar the scaling and the wavelet functions 42,kand vl,k.The discrete
to Equation (81), is also formed. 26 elements per row of the matrix wavelet transform of [Y,] is depicted in Figure 5 . The choice of
are zero. Note that the first 16 rows of [PI] are formed by h’(m), the filter h’(m) completely defines the entire procedure. Also, note
and the last sixteen rows, by g’(m) . We pre-multiply [Yl] by [I?,], that for an infinite sum in Equation ( 5 5 ) , we never talk about the
and obtain a vector [Y2]. The last 16 elements of [Y2] are dk,o, scaling function in the summation. However, for the finite-sum
case, the wavelets by themselves are not complete. One needs the
and they are fixed. This is the result of filtering [Y,] by g(n) , in
contribution of the scaling functions. Unfortunately, this very
Figure 5 . We pre-multiply the first 16 elements of [Y2],Le., [Y;] important feature is not spelled out explicitly in the literature.
as shown in Figure 11, by an orthogonal matrix [P2], which is
To compute the inverse of the discrete wavelet transform of
16 x 16. Again, only six of the elements of any row of matrix [P2]
the resultant vector, one simply reverses the procedure. In this
are nonzero. The result will be a vector [Y3]of 16 elements. The case, we start with the smallest level of the hierarchy, and work our
last eight elements are fixed, as they are dk,, . The first eight ele- way through. This is expressed by Equation (74). The inverse
ments of [Y3], namely CY;], are again pre-multiplied by the matrix of all the [PIS in this case is simply their transpose, as [PIS
orthogonal matrix [PI of Equation (81), as shown in Figure 11, are real orthogonal matrices. To compute the two-dimensional
wavelet transform, one follows the rules of the FFT. One deals
producing c k , 2 and d k , 2 . The final result is then the wavelet
with all the rows and then all the columns.
decomposition of the vector Yl . The resultant composite vector of
32 elements includes the wavelet coefficients resulting from the Next, we consider the solution of the matrix equation
discrete wavelet transform. It is interesting to note that when deal- [A][ X] = [Y] , utilizing the discrete wavelet transform. The discrete
ing with a finite-length vector, the widely presented formula of wavelet transform actually does not solve any matrix equations.
Equation ( 5 5 ) is no longer applicable. Here we use Equation (71).
What the wavelet transform does is to preprocess the matrix [A]
Instead, we have a series containing both the scaling functions and
wavelets. Specifically, for the case we have just considered, we and make it sparse, if some thresholding is applied [31].
obtain a transformed vector of 32 elements, [TI, which is given by
The basic principle, as outlined by Beylkin, Rokhlin, and
Coifman [ 181, based on [191, is as follows. Let the matrix elements
be generated from a kemel, such that the magnitude of the ele-
ments of the matrix [A] decay from the diagonal as -,where
li- jl“
i a n d j may he the row or the column number. Then, if the two-
dimensional discrete wavelet transform is applied to the system
matrix [A], the resulting system matrix will be sparse if all the
elements below a threshold are set to zero. Typically, one would
have only 10Qloglo

I:( elements in the sparse system, where Q is
the size of the matrix and E is the truncation level, i.e. elements of
the resultant matrix whose absolute value is less than E will be Order YOof Error in Cond[B,]
discarded. So, for a 2048x 2048 matrix, the resultant system of Filter Nonzero Reconstruction
matrix would be sparse by a factor of about 30 [20]. ~ + 1Elements 6
4 7.58% ,58 10-4 4.66 x 106
Consider a real matrix [A], of dimension Q x Q , where Q is
8 6.28% ,54 10-4 5.56 x lo6
large. Let the elements of [A] be defined by
16 6.41% ,52x10-4 5.74 107
32 6.54% ,5iX10-4 1.57 107
( 1 if i = j
We now apply a wavelet transform to the matrix [A]. This is

equivalent to pre- and post-multiplying [A] by a number of
orthogonal matrices. Let [SI be the product of the orthogonal
matrices [Pi], as explained earlier for the one-dimensional discrete
wavelet transform explained by Figure 11: Order YOof Error in Cond[B,]
of Filter Nonzero Reconstruction
Even though the sizes of the various [Pl]s are not the same, we
make them the same by supplementing, say, [P2] by a diagonal
unity matrix, [I], to make it the same size as [P,]. Thus,
rows and columns, which is similar to carrying out a two-dimen-

(92)
sional FFT.
Now, we consider [A] to be of the form in Equation (go),

Since the product of all orthogonal matrices is an orthogonal
and choose various order filters (i.e., M = N + 1) for h(n), with 4,
matrix, [SI is an orthogonal matrix. When the wavelet transform is
applied to the system of equations [A][X] = [ Y ] one
, obtains 8, 16, or 32 terms. Next, we compute [B] = [S][A][SIT ([SI being
given by Equation (91)), and then apply a threshold to the elements
of [B] to obtain matrix [B,]. The matrix [B,] is an extremely
(93)
sparse matrix. For example, if the threshold is set at and the
where T denotes the transpose of a matrix. Since [SI is an size of [A] is Q = 5 12, and if we then apply a fourth-order filter
orthogonal matrix, we have (Le., M = N+1= 4 ) h(n), we find that only 7.58% of the ele-
ments of [Ba] are nonzero (see Table 1).
Next, take that sparse matrix [Ba], and try to reconstruct [A]
where [I] is the identity matrix. So, we take the wavelet transform by computing
of [XI and [Y]to form [X'] and [Y'], and we take the wavelet
transform of [A] to form [B] . Hence, Equation (93) reduces to
[B][X'] = [Y'] with [B] = [S][A][S]T. (95) Define an average value, 6, of the error between the elements of
[ A a l and [AI by
The unknown [XI is solved for fi-om this by
[XI = [sIT[xr]. (96)
[B] is the two-dimensional wavelet transform of [A], and has From Tables 1-6, it is seen that the matrix [A] can be recovered
been computed by a series of one-dimensional transforms to its from [B,] to provide [A,] with an average error that is lower than
64 IEEEAntennas and Propagation Magazine, Vol. 40, No. 5, October 1998

Table 3. The sparseness, error in reconstruction, and condition the value of the threshold used to eliminate the wavelet coefficients
number of the two-dimensional wavelet transform of a matrix of [A], (or, equivalently, the elements of matrix [B] ). The number
[A] of the form in Equation (go), as a function of the order of of [P,] matrices (see Equation (92)) used for obtaining the results
the filter, with a threshold of 0.001, Q = 2048, and the condi-
tion number of [A] given by Cond[A]= 6.80 x 10’. [ [21)]
appearing in Tables 1-6 is INT log2 __ (LNT stands for “the
integer part of’). Please note that the wavelet transform of a real
o:-[ 1 Error in
Reconstruction 1 Cond[B,] function is always real.
The other interesting property to note is that utilizing a

higher-order filter does not necessarily produce larger sparsity, as
shown in Tables 2,4, and 5.
The computational time in the computation of the wavelet

transforms in the compression of a matrix is now investigated. We
consider a filter of length L = N + 1 , and the data matrix is of
Table 4. The sparseness, error in reconstruction, and condition length Q. Then the one-dimensional wavelet transform of [Y] may
number of the two-dimensional wavelet transform of a matrix be done (we carry out the initial product at the highest stage of
[A] of the form in Equation (90), as a function of the order of resolution and then down-sample) by the following number of
mathematical operations:
tion number of [A] given by Cond[A] = 4.45 x 10‘.
G:)
QL 1+-+-+ ... =2QL (99)
To carry out the two-dimensional wavelet transform of [A], we

Elements
require (2QL)* operations. Therefore, to produce the sparse sys-
tem of Equation (95), we require 4Q2L2+ 2QL operations, result-
ing in a matrix [B] that contains, at most, of the order of O(Q)
elements. If we now apply the conjugate-gradient method to solve
the sparse system, per iteration we will require O(2Q) multiplica-
tions to carry out two matrix-vector products. Even if the conju-
Table 5. The sparseness, error in reconstruction, and condition gate-gradient method converges in at most Q steps (where Q is the
number of the two-dimensional wavelet transform of a matrix number of unknowns), then we have solved Equation (95) in an
[A] of the form in Equation (90), as a function of the order of operation count of Q 0 ( 2 Q ) , in addition to 4Q2L2+ 2QL opera-
tions. Observe that this O(Q2) is significantly lower than the con-
tion number of [A] given by Cond[A] = 3.02 x lo7.
ventional Q3/3 operations typically required in a solution of a
oi=
: I Error in
Reconstruction 1 Cond[B,] matrix equation of size Q. This i s essentially the contribution of
Beylkin, Coifman, and Rokhlin [18].
However, it is interesting to note that the nature of the varia-

tion ____ may also be the result from a convolution. In that
li - jl“
case, the FFT would be much faster than the discrete wavelet trans-
form, as the FFT essentially diagonalizes the operatorimatrix.
However, for other cases, even when the variation is not due to a
Table 6. The sparseness, error in reconstruction, and condition convolution, the wavelet result still holds.
number of the two-dimensional wavelet transform of a matrix
[A] of the form in Equation (90), as a function of the order of The most disturbing fact about the discrete wavelet transform
the filter, with a threshold of 0.0001, Q = 2048, and the condi- is that the condition number of the matrix changes after the trans-
form. It has been shown earlier that the wavelet transform of the
tion number of [A] given by Cond[A] = 6.80 x lo7.
matrix [A], which is [B] , has been formed through a series of
orthogonal transformations. Therefore, by definition, the condition
Error in Cond[B,] number of the matrix should not change as one goes from [A] to
Elements [B] . However, as is clear from the six tables, there is a change in
the condition number of the transformed matrix when one uses a
different-order filter. Also, the condition number is dependent on
the threshold used to truncate the elements. There appears to be a
lack of systematic change in the results. In this case, the matrix
[A] is real.
As a final example, consider the electromagnetic scattering
from an array of wires, randomly spaced. We considered 56 thin-
wire antennas. Six were of length 2.7A and radius 0.0051. The
remaining 50 were 32 long and of the same radius. The 56 wires
were located inside a parallelepiped of dimensions
271xX5/2x21/2. The usual MOM application led to a
2096 x 2096 matrix. The matrix was compressed, utilizing a filter
h ( n ) of length 16. The compression for the real part of the imped-
ance matrix was 17.8%, Le. 749289 of the elements of the matrix
c x (4 YD(4
Figure 12. The symbol used to represent decimation by a factor
of two.
were above a threshold of For the imaginary part of the

impedance matrix, only 3.28% of the elements were nonzero. This
sparse impedance matrix was then used in a conjugate-gradient
routine, to solve the transformed Equation (93). Convergence in
the residuals of was obtained in 95 iterations. This simple
example demonstrates the potential for the solution of large matrix
equations.
There are a few points that are worth mentioning. First of all,
if the compression was not done on the real and the imaginary
4 3
>
parts of the matrix separately, then the degree of compression was
merely 35%, as opposed to 17.8%. This is significant. Secondly,
the size of the impedance matrix has to be a power of two. Thirdly,
the conjugate-gradient method takes the same number of iterations
(95) to converge to the solution when applied to the original dense
matrix, or to the sparse matrix, as the transformation is presumably
orthogonal, from a strictly theoretical point of view. However, now
as the entire compressed matrix is in memory, the number of page
faults is small, and so the result can be obtained quite efficiently.
Note that one of the disadvantages with this procedure is that

the size of the matrix Q has to be a power of two for efficient 0 1
implementation of the wavelet transform.
6. Conclusion Figure 13. The waveform obtained in the sampled domain

from decimation by a factor of two.
The discrete wavelet transform has been presented from first
principles, utilizing the basic concepts of filter theory. It has been
shown how to construct the filters h(m) that produce the wavelets P"'
and the scaling functions. However, for the discrete case, the intro-
duction of wavelets and scaling functions are not at all necessary.
Finally, it has been shown how to apply this technique to the solu-
tion of large matrix equations. This was accomplished by com-
pressing a large matrix by means of the discrete wavelet transform.
The disturbing point is that the condition number before and after Figure 14. The spectrum of the original signal.
the transform and thresholding is quite different as a function of
the order of the filter.
8. Acknowledgment
The contributions of the Editor-in-Chief in helping to prepare

this paper for publication are gratefully acknowledged.
Figure 15. The spectrum of the down-sampled signal.
9. Appendix 1. The principle of decimation by a factor of two
Pictorially, decimation by a factor of two is represented by

the symbol shown in Figure 12, where the decimated signal,
y D ( n ) ,has been generated from the original signal, x ( n ) . In the or in the z transform domain,
sampled domain, this is equivalent to obtaining the waveform
shown in Figure 13. Note that altemate sample values have been Yo(,)= - y y D ( n ) z - n= CX(M)Z-"2
dropped. Therefore, from Figure 13, in the sampled domain n form even
= ‘[x(S)
2
+ Xj-&)I.
If we observe the spectrum, then one observes that the original sig- Figure 18. The spectrum of the original signal.
nal has the spectrum given in Figure 14. Once the signal is down-
sampled, the spectrum is as given in Figure 15. Hence, the spec-
trum of YD(e’”) is aliased, and it is the sum of
Figure 19. The spectrum of the up-sampled signal.

10. Appendix 2. The principle of expansion by a factor of two
Pictorially, up-sampling can be represented by the symbol or, in the transform domain
shown in Figure 16. In the sampled domain, this is equivalent to
inserting a zero between the sampled signals, as shown in Figure
17. Mathematically, this is equivalent to
If we observe the spectrum of Y L ( z ) ,then we observe that the

original signal in Figure 18 is transformed into the spectrum of Y2 ,
as shown in Figure 19.
11. References
1. Ingrid Daubechies, Ten Lectures on Wavelets, CBMS-NSF

Regional Conference Series in Applied Mathematics, Philadelphia,
SIAM, 1992.
Figure 16. The symbol used to represent up-sampling by a fac-
tor of two. 2 . C. K. Chui, An Introduction to Wavelets, New York, Academic
Press, 1992.
P 3. P. P. Vaidyanathan, Multirute Systems and Filter Banks,

Englewood Cliffs, NJ, Prentice-Hall, 1993.
4. J. Morlet, G. Arens, I. Fourgeau and D. Giard, “Wave Propaga-

tion and Sampling Theory,” Geophysics, 47, 1982, pp. 203-236.
5 . Special issue ofIEEE Proceedings, 84,4, April 1996.
6. D. Esteban and C. Galland, “Application of Quadrature Mirror

Filters to Split-Band Voice Coding Schemes,” Proceedings of the
IEEE International Conference on Acoustics, Speech, and Signal
Processing, Hartford, Connecticut, 1977, pp. 191-195.
7. M. J. T. Smith and T. P. Barnwell, 111, “Exact Reconstruction

Techniques for True Structured Subband Coders,” IEEE Transac-
tions on Acoustics, Speech, and Signal Processing, ASSP-39,
1986, pp. 434-441.
8. M. Vettereli, “Multidimensional Subband Coding: Some Theory

and Algorithms,” Signal Processing, 6 , 1984, pp. 97-1 12.
9. T. H. Koomwinder (ed.), Wavelets: An Elementary Treatment of

Theory and Applications, New Jersey, World Scientific, 1993.
Figure 17. Up-sampling by a factor of two is equivalent to
inserting a zero between the sampled signals in the sampled 10. S . Schweid and T. K. Sarkar, “A Sufficiency Criteria for
domain. Orthogonal QMF Filters to Ensure Smooth Wavelet Decomposi-
/E€€ Antennas and Propagation Magazine, Vol. 40, No. 5,October 1998 67
tion,” Applied and Computational Harmonic Analysis, 2, 1995, pp. 25. C. Su and T. K. Sarkar, “A Multiscale Moment Method for
61-67. Solving Fredholm Integral Equation of the First Kind,” in J. A.
E=mg (ed.), Reports on Progress in Electromagnetic Research,
11. S. Schweid and T. K. Sarkar, “Iterative Calculating and Fac- Volume 17, Cambridge, Massachusetts, EMW Publishing, 1997,
torization of the Autocorrelation Function of Orthogonal Wavelets pp. 237-264.
with Maximal Vanishing Moments,” IEEE Transactions on Cir-
cuits and Systems, CAS-42, 11, November 1995, pp. 694-701.
26. T. K. Sarkar, R. S. Adve, L. Castillo and M. Salazar, “Utiliza-
tion of Wavelet Concepts in Finite Elements for an Efficient Solu-
12. R. Ansari, C. Guillemot and J. F. Kaiser, “Wavelet Construc-
tion of Maxwell’s Equations,” Radio Science, 29, July 1994, pp.
tion Using Lagrange Half-Band Filter,” IEEE Transactions on Cir-
cuits andSystems, CAS-38, 1991, pp. 1116-1118. 965-977.
13. M. A. Krasnoselskii, G. N. Vainikko, P. P. Zabreiko, Ya B. 27. C. Su and T. K. Sarkar, “Electromagnetic Scattering from
Rutitskii and V. Ya Stetsenko, Approximate Solution of Operator Coated Strips Utilizing the Adaptive Multiscale Moment Method,”
Equations, Groningen, The Netherlands, Walters-Noordhoff, 1972. in J. A. King (ed.), Reports on Progress in Electromagnetic
Research, Volume 18, Cambridge, Massachusetts, EMW Publish-
14. S. G. Mikhlin, The Numerical Performance of Variational ing, 1998, pp. 173-208.
Methods, Groningen, The Netherlands; Walters Noordhoff, 1971.
28. A. Brandt, “Multi-level Adaptive Solutions to Boundary Value
15. K. Rektorys, VariationalMethods in Mathematics, Science and Problems,” Mathematics of Computation, 31, 1997, pp. 330-390.
Engineering, Dordrecht, Holland, D. Reidel Publishing Company,
1975. 29. W. Hackbusch, Multigrid Methods and Applications, New
York, Springer-Verlag, 1985.
16. T. K. Sarkar, “A Note on the Choice of Weighting Functional
in the Method of Moments,” IEEE Transactions on Antennas and 30. S. F. McCormick, Multigrid Methods: Theoly, Applications
Propagation, AP-33,4, April 1985, pp. 436-441. and Supercomputing, New York, Marcel Dekker, 1988.
17. T. K. Sarkar, A. R. Djordjevic and E. h a s , “On the Choice of 3 1. J. Mandel, “On Multilevel Iterative Methods for Integral Equa-
Expansion and Weighting Functions in the Numerical Solution of tions of the Second Kind and Related Problems,” Numer. Math.,
Operator Equation,” IEEE Transactions on Antennas and Propa- 46, 1985, pp. 147-157.
gation, AP-33, 9, September 1985, pp. 988-996.
32. P.W. Hemker and H. Schippers, “Multiple Grid Methods for
IS. G. Beylkin, R. Coifman and V. Rokhlin, “Fast Wave Trans- the Solution of Fredholm Integral Equations of the Second Kind,”
forms and Numerical Algorithms,” Communications on Pure and Mathematics of Computation, 36, 153, 1981.
AppliedMathematics, 44, pp. 141-183.
33. K. Kalbasi and K. R. Demarest, “A Multilevel Enhancement of
19. A. P. Calderon, “Intermediate Spaces and Interpolation, the the Method of Moments,” in Seventh Annual Review of Progress in
Complex Method,”Stud. Math., 24, 1964, pp. 113-190. Applied Computational Electromagnetics, March 199 1, Monterey,
California, Naval Postgraduate School, pp. 254-263.
20. W. H. Press, S, Teukolsky, W. T. Vetterling and B. P. Flemery,
Numerical Recipes-The Art of Scientijic Computing, Cambridge, 34. K. Kalbasi and K. R. Demarest, “A Multilevel Formulation of
Cambridge University Press, 1992. the Method of Moments,” IEEE Transactions on Antennas Propa-
gation, AP-41, 5, May 1993, pp. 589-599.
21. S. C. Schweid, “Projection Method Minimization Techniques
for Smooth Multistage QMF Filter Decompositions,” PhD thesis, 35. Z. Baharav, Y . Leviatan, “Impedance Matrix Compressing
Syracuse University, May 1994. Using Adaptively Constructed Basis Functions,” IEEE Transac-
tions on Antennas and Propagation, AP-44, 9, 1996, pp. 1231-
22. B. Z.Steinberg and Y. Leviatan, “On the Use of Wavelet 1238.
Expansions in the Method of Moments,” IEEE Transactions on
Antennas and Propagation, AP-41, 5, 1993, pp. 610-619. 36. R. L. Wagner and W. C. Chew, “A Study of Wavelets for the
Solution of Electromagnetic Integral Equations,” IEEE Transac-
23. J. C. Goswami, A. K. Chan and C. K. Chui, “On Solving First- tions on Antennas and Propagation, AP-43, 8, 1995, pp. 802-810.
Kind Integral Equations Using Wavelets 011 a Bounded Interval,”
IEEE Transactions on Antennas and Propagation, AP-43, 6, June 37. F. X. Canning, J. F. School, “Diagonal Preconditioners for the
1995, pp. 614-622. EFIE Using a Wavelet Basis,” IEEE Transactions on Antennas and
Propagation, AP-44, 9, 1996, pp. 1239-1246.
24. G. Wang, “A Hybrid Wavelet Expansion and Boundary Ele-
ment Analysis of Electromagnetic Scattering from Conducting 38. G. Strang and T. Nguyen, Wavelet and Filter Banks, Cam-
Object,” IEEE Transactions on Antennas and Propagation, AP-43, bridge, Wellesley-Cambridge Press, 1996.
2, February 1995, pp. 170-178.
68 IEEE Antennas and Propagation Magazine, Vol. 40, No. 5,October 1998
Introducing Feature Article Authors PhD degree at Syracuse University. His research interests include
numerical electromagnetics and the use of signal-processing tech-
niques in numerical electromagnetics.
Tapan Kumar Sarkar received the B. Tech. degree from the

Indian Institute of Technology, Kharagpur, India, the MScE degree Magdalena Salazar-Palma was born in Granada, Spain. She
from the University of New Brunswick, Fredericton, Canada, and received the PhD degree in Ingeniero de Telecomunicacion from
the MS and PhD degrees from Syracuse University, Syracuse, New the Universidad Politecnica de Madrid (Spain), where she is a Pro-
York in 1969, 1971, and 1975, respectively. fesor Titular of the Departamento de Seiiales, Sistemas y Radio-
comunicaciones (Signals, Systems and Radiocommunications
From 1975 to 1976, he was with the TACO Division of the Department) at the Escuela Tecnica Superior de Ingenieros de
General Instruments Corporation. From 1976 to 1985, he was with Telecomunicacion of the same university. She has taught courses
the Rochester Institute of Technology, Rochester, NY. From 1977 on electromagnetic field theory, microwave and antenna theory,
to 1978, he was a Research Fellow at the Gordon McKay Labora- circuit networks and filter theory, analog and digital communica-
tory, Harvard University, Cambridge, MA. He is now a Professor tion systems theory, and numerical methods for electromagnetic
in the Department of Electrical and Computer Engineering, Syra- field problems, as well as related laboratories. Her research within
cuse University; Syracuse, NY. He has authored or co-authored the Grupo de Microondas y Radar (Microwave and Radar Group)
more than 170 journal articles, and has written chapters in ten is in the areas of electromagnetic field theory, and computational
books. His current research interests deal with numerical solutions and numerical methods for microwave structures, passive compo-
of operator equations arising in electromagnetics and signal proc- nents, and antenna analysis; design, simulation, optimization,
essing, with application to system design. implementation, and measurements of hybrid and monolithic
microwave integrated circuits; and network and filter theory and
Dr. Sarkar is a registered Professional Engineer in the State design.
of New York. He was an Associate Editor for Feature Articles of
the IEEE Antennas and Propagation Society Newsletter, and he She has authored a total of 10 contributions for chapters and
was the Technical Program Chairman for the 1988 IEEE Antennas articles in books, 15 papers in international journals, and 75 papers
and Propagation Society International Symposium and URSI in international symposiums and workshops, plus a number of
Radio Science Meeting. He was the Chairman of the Intercommis- national publications and reports. She has lectured in several short
sion Working Group of URSI on Time-Domain Metrology. He is a courses, some of them in the framework of European Community
member of Sigma Xi and USNC/URSI Commissions A and B. He Programs. She has participated in 19 projects and contracts,
received one of the “best solution” awards in May, 1977, at the financed by international, European, and national institutions and
Rome Air Development Center (RADC) Spectral Estimation companies. She has been a member of the Technical Program
Workshop. He received the Best Paper Award of the IEEE Trans- Committee of several international symposiums, and has acted as
actions on Electromagnetic Compatibility in 1979, and at the 1997 reviewer for international scientific journals. She has assisted the
National Radar Conference. He is a Fellow of the IEEE. He Comision Interministerial de Ciencia y Tecnologia (National
received the degree of Docteur Honoris Causa from the Universite Board of Research) in the evaluation of projects. She has also
Blaise Pascal, Clermont-Ferrand, France, in 1998. served in several evaluation panels of the Commission of the
European Communities. She has acted in the past and is currently
Chaowei Su was born in Jiangsu, People’s Republic of acting as topical Editor for the disk of references of the Review of
China, on April 2 3 , 1961. He received the BS and MSc degrees Radio Science. She is a member of the editorial board of two sci-
from Northwestern Polytechnical University (NPU), both in the entific journals. She has served as Vice Chair and Chair of the
Department of Applied Mathematics, in 1981 and 1986, respec- IEEE MTTIAP-S Spanish joint Chapter, and is currently serving as
tively. He joined the Department of Applied Mathematics at NPU Chair of the Spain Section of the IEEE. She has received two indi-
in 1982. Mr. Su was appointed as an Associate Professor and a full vidual research awards from national institutions.
Professor at NPU in 1993 and 1996, respectively. He was a Visit-
ing Scholar at SUNY, Stony Brook, New York, USA, between
December, 1990, and June, 1992, supported by a Grumman Fel-
lowship. Since August, 1996, he has been a Visiting Professor at
Syracuse University, IJSA. He is an author of Numerical Methods
in Inverse Problems of Partial Differential Equations and Their
Applications (published by NPU Press). His main research inter-
ests are in the area of numerical methods of electromagnetic scat-
tering and inverse scattering, and inverse problems of partial dif-
ferential equations.
Raviraj S. Adve was born in Bombay, India. He received the

BTech degree in electrical engineering from the Indian Institute of Luis-Emilio Garcia-Castillo was born in 1967 in Madrid,
Technology, Bombay, in 1990. He is currently working toward the Spain. In 1992. he received the degree of Ingeniero de Telecomu-
nicaci6n from the Universidad PolitCcnica de Madrid (Spain). In
1993, he became a Research Assistant in the Departaniento de
Seiiales, Sistemas y Radiocomunicaciones (Signals, Systems and
Radiocommunications Department) in the Escuela TBcnica Supe-
rior de Ingenieros de Telecomunicaci6n of the same university.
Since 1997, he has been an Associate Professor of the Departa-
mento de Ingenieria Audiovisual y de Comunicaciones (Audiovis-
ual and Communications Engineering Department) in the Escuela
Univesitaria de Ingenieria TCcnica de Telecomunicaci6n of the
same university, where he teaches microwave theory and the
related laboratory. His research activities and interests are focused
on the application of numerical methods, mainly finite elements, to
electromagnetic problems, including research on curl-conforming
finite elements, the characterization of multiconductor and Rafael Rodriguez Boix received the BSc, MSc, and PhD
waveguide structures, analysis of scattering and radiation prob- degrees in physics from the University of Seville, Spain, in 1985,
lems, and the use of wavelet theory in computational electromag- 1986, and 1990, respectively. Since 1985, he has been with the
netics. Other research areas of his interest are network theory and Electronics and Electromagnetics Department at the University of
filter design. Seville, where he became an Associate Professor in 1994. Durmg
the summers of 1991 and 1992, he was at the Electrical Engineer-
He has authored four contributions for chapters and articles ing Department of UCLA as a Visiting Researcher Also, d u n g
in books, six papers in international joumals, and 29 papers in the summer of 1996, he was at the Computer and Electrical Engi-
intemational symposiums, plus a number of national publications neering Department of Syracuse University as a Visiting
and reports. He has participated in seven projects and contracts Researcher. His current research activities are focused on the
financed by international, European, and national institutions and numerical electromagnetic analysis of planar microwave circuits
companies. and printed circuit antennas. %’
August 13-21,1999, University of Toronto, Canada
irst Call for Papers Now Available
The first announcement booklet and call for papers for the XXVIth General Assembly of the
Intemational Union of Radio Science is now available. It includes the schedule for the sessions
of the 10 URSI Commissions, as well as the instructions and format for submitting papers. There
is also information on the Young Scientists Program and a Canadian Student Competition.
The information in this first announcement booklet is essential for anyone wishing to submit an
abstract of a paper to be presented at the General Assembly. The booklet can be obtained by
sending a request with full address and contact information to
URSI GA ‘99 Management Office

National Research Council Canada
Montreal Road, Building M-19
Ottawa, Ontario, Canada K1A OR6
Tel: (613) 993-7271; Fax: (613) 993-7250
E-mail: [email protected]
Those interested can also obtain the information, and request being placed on the mailing list, at
The deadline for receipt of abstracts is January 15,1999.

Engineering Perspective, Part Discrete Wavelet Techniques': Tutorial Wavelets From An Electrical 1

Uploaded by

Copyright:

Available Formats

Engineering Perspective, Part Discrete Wavelet Techniques': Tutorial Wavelets From An Electrical 1

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Engineering Perspective, Part Discrete Wavelet Techniques': Tutorial Wavelets From An Electrical 1

Uploaded by

Copyright:

Available Formats

A Tutorial on Wavelets from an Electrical

Engineering Perspective, Part 1:

‘Department of Electrical and Computer Engineering

2Dept. Senales and Sistemas

3Dept. Electronica y Electromagnetism0

Keywords: Wavelets, wavelet transforms, filters, matrix equations, 2. Introduction

Consider a signal, x ( t ) , the Fourier transform of which is

Using this type of sub-band splitting has many advantages:

2. Even if .(E) and v(n) are quantized in a very rough fash-

G ' ( z )= z - ~ G ( z - ' )with g ' ( n ) = g ( N - a ) . (16)

V ( z ) = '2[ V ~ ( & ) + V ~ ( - & ) ] . so that

By using the results of Appendix 2 (where ~ " ( n and

Furthermore, by utilizing Equations (15) and (16), Equation

1 If the filter H ( z ) is chosen in such a way that

The estimated signal X ( z ) contains the original signal

and the estimated signal is delayed by N samples.

Now, we show how to solve for H ( z ) . In this example, we

L =order of filter= N +1= 4 . (27)

Utilizing Equation (20), we have

[h(O)- h(1)z-l + h(2)z-Z - h(3)z"I

Equating all the coefficients related to the individual powers of z in

Figure 3. Non-overlapping ideal filters.

G ( ~=)- z - 3 ~ ( - z - l ) From Equations (32) and (34) we see that

and from Equation (15), 1

52 IEEE Antennas and Propagation Magazine, Vol. 40,No. 5, October 1998

equal to zero. Daubechies chose this procedure because enforcing I

the above conditions guarantees smooth wavelets, which we will

-Oh(O) - Ih(1) + 2h(2) - 3h(3) = 0 . (37)

I+& Instead of having a two-stage decomposition of the signal

3.2 Connection between the filter theory and the mathematical

To establish a connection between the filter theory and the

0 0.5 1 1.5 2 2.5

/€€€Antennas and Propagation Magazine, Vol. 40, No. 5, October 1998 53

a continuous variable. They can have integer shifts. In addition,

We consider the original signal, x ( t ) , with bandwidth from 0

or, equivalently, [with h'(n) = h(N - m) , as defined by Equa-

It is interesting to note that the coefficients h'(m) tum out to be the

As an example, consider 4(t) to be a pulse function of mag-

The fimction d ( t ) , which solves the dilation equation for a

Figure 6. A filter-theory representation of the discrete wavelet I 4(t - m)dt = 1

or, equivalently, When N = 0 , the dilation equation becomes, for hb = &,

nz=O and the scaling function becomes a delta function. If we choose

Owing to Equation (42), $(O) = 1. Therefore, in the limit n + ca,

The solution for $(a)is point-wise convergent, provided the infi-

Consider a function ~ ( t. The

i(4= { 1 for o</w(<z

where we define the wavelets by

with 4i - 1 coefficients. The coefficients hi are determined using

Now, if we have to carry out the inner products in Equa-

We further assumeiutilize the orthogonality relationships between which can be generalized to

and Next, observe that

represents a basis, as the functions involved in the set are orthogo-

These results have been derived utilizing Equations (1 5) and (16).

There is another subtle point that one should introduce now!

Figure 9b. A compressed version of the picture of “Lena,”

In summary, the evaluation of the discrete wavelet coeffi-

Figure 9c. A compressed version of the picture of “Lena,” util-

putations. However, if the sum is finite in Equation ( 5 5 ) , then one