418 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 52, NO. 2, FEBRUARY 2004
ForWaRD: Fourier-Wavelet Regularized
Deconvolution for Ill-Conditioned Systems
Ramesh Neelamani, Member, IEEE, Hyeokho Choi, Member, IEEE, and Richard Baraniuk, Fellow, IEEE
Abstract—We propose an efficient, hybrid Fourier-wavelet reg-
ularized deconvolution (ForWaRD) algorithm that performs noise
regularization via scalar shrinkage in both the Fourier and wavelet
domains. The Fourier shrinkage exploits the Fourier transform’s
economical representation of the colored noise inherent in de-
convolution, whereas the wavelet shrinkage exploits the wavelet
domain’s economical representation of piecewise smooth signals
and images. We derive the optimal balance between the amount of
Fourier and wavelet regularization by optimizing an approximate
mean-squared error (MSE) metric and find that signals with
more economical wavelet representations require less Fourier
shrinkage. ForWaRD is applicable to all ill-conditioned deconvolu-
tion problems, unlike the purely wavelet-based wavelet-vaguelette
deconvolution (WVD); moreover, its estimate features minimal
ringing, unlike the purely Fourier-based Wiener deconvolution.
Even in problems for which the WVD was designed, we prove
that ForWaRD’s MSE decays with the optimal WVD rate as the
number of samples increases. Further, we demonstrate that over
a wide range of practical sample-lengths, ForWaRD improves on
WVD’s performance.
Index Terms—Deblurring, deconvolution, restoration, wavelet-
vaguelette, wavelets.
I. INTRODUCTION
D
ECONVOLUTION is a recurring theme in a wide variety
of signal and image processing problems. For example,
practical satellite images are often blurred due to limitations
such as aperture effects of the camera, camera motion, or atmo-
spheric turbulence [1]. Deconvolution becomes necessary when
we wish a crisp deblurred image for viewing or further pro-
cessing.
A. Problem Statement
In this paper, we treat the classical discrete-time deconvolu-
tion problem. The problem setup and solutions are described in
one dimension (1-D), but everything extends directly to higher
Manuscript received October 21, 2002; revised March 27, 2003. This work
was supported by the National Science Foundation under Grant CCR-99-73188,
the Air Force Office of Scientific Research under Grant F49620-01-1-0378,
the Office of Naval Research under Grant N00014-02-1-0353, the Defense
Advanced Research Projects Agency under Grant F30602-00-2-0557, and
the Texas Instruments Leadership Universities Program. The associate editor
coordinating the review of this paper and approving it for publication was Dr.
Chong-Yung Chi.
R. Neelamani was with the Department of Electrical and Computer Engi-
neering, Rice University, Houston, TX 77005-1892 USA. He is now with the
Upstream Research Company, ExxonMobil, Houston, TX 77252 USA.
The authors are with the Department of Electrical and Computer Engineering,
Rice University,Houston, TX 77005-1892 USA (e-mail: choi@ece.rice.edu;
richb@ece.rice.edu).
Digital Object Identifier 10.1109/TSP.2003.821103
Fig. 1. Convolution model setup. The observation
y
consists of the desired
signal
x
first degraded by the linear time-invariant (LTI) convolution system
H
and then corrupted by zero-mean additive white Gaussian noise (AWGN)
.
dimensions as well. The observed samples consist of un-
known desired signal samples
first degraded by circular
convolution (denoted by
) with a known impulse response
from a linear time-invariant (LTI) system and then cor-
rupted by zero-mean additive white Gaussian noise (AWGN)
with variance (see Fig. 1)
(1)
Given
and , we seek to estimate .
A naive deconvolution estimate
is obtained using the oper-
ator inverse
as
1
(2)
Unfortunately, the variance of the colored noise
in
is large when is ill conditioned. In such a case, the mean-
squared error (MSE) between
and is large, making an un-
satisfactory deconvolution estimate.
In general, deconvolution algorithms can be interpreted as es-
timating
from the noisy signal in (2). In this paper, we focus
on simple and fast estimation based on scalar shrinkage of indi-
vidual components in a suitable transform domain. Such a focus
is not restrictive because transform-domain scalar shrinkage lies
at the core of many traditional [3], [4] and modern [2], [5] de-
convolution approaches.
B. Transform-Domain Shrinkage
Given an orthonormal basis
for , the naive esti-
mate
from (2) can be expressed as
(3)
1
For noninvertible
H
, we replace
H
by its pseudo-inverse and
x
by its
orthogonal projection onto the range of
H
in (2) [2]. The estimate
x
in (2) con-
tinues to retain all of the information that
y
contains about
x
.
1053-587X/04$20.00 © 2004 IEEE
NEELAMANI et al.: ForWaRD: FOURIER-WAVELET REGULARIZED DECONVOLUTION FOR ILL-CONDITIONED SYSTEMS 419
(a) (b)
Fig. 2. Economy of Fourier versus wavelet representations. (a) Energies in decibels (dB) of the Fourier and wavelet components of the noise
H
colored by
the pseudo-inverse of a 2-D 9
2
9 box-car smoothing operator. The components are sorted in descending order of energy from left to right. The colored noise
energy is concentrated in fewer Fourier components than wavelet components. (b) Energies of the Fourier and wavelet components of the
Cameraman image
x
.
The signal energy is concentrated in fewer wavelet components than Fourier components.
An improved estimate can be obtained by simply shrinking
the
th component in (3) with a scalar , [6]:
(4)
(5)
The
denotes the retained part of
the signal
that the shrinkage preserves from (2), whereas
denotes the leaked part
of the colored noise
that the shrinkage fails to at-
tenuate. Clearly, we should set
if the variance
of the th colored noise component
is small relative to the energy
of the corresponding
signal component and set
otherwise. The shrinkage by
can also be interpreted as a form of regularization for the
deconvolution inverse problem [4].
The tradeoff associated with the choice of
is easily under-
stood: If
, then most of the th colored noise component
leaks into
with the corresponding signal component; the re-
sult is a distortion-free but noisy estimate. In contrast, if
,
then most of the
th signal component is lost with the corre-
sponding colored noise component; the result is a noise-free but
distorted estimate. Since the variance of the leaked noise
in (5) and the energy of the lost signal comprise the MSE
of the shrunk estimate
, judicious choices of the ’s help
lower the estimate’s MSE.
However, an important fact is that for a given transform do-
main, even with the best possible
’s, the estimate ’s MSE
is lower bounded by [5], [7], [8]
(6)
From (6),
has small MSE only when most of the signal en-
ergy (
) and colored noise energy ( )is
captured by just a few transform-domain coefficients—we term
such a representation economical—and when the energy-cap-
turing coefficients for the signal and noise are different. Oth-
erwise, the
is either excessively noisy due to leaked noise
components or distorted due to lost signal components.
Traditionally, the Fourier domain (with sinusoidal
’s) is
used to estimate
from . For example, the LTI Wiener decon-
volution filter corresponds to (4) with each
determined by
the
th component’s signal-to-noise ratio [3], [4]. The strength
of the Fourier basis is that it most economically represents the
colored noise
[see Fig. 2(a) and Section III-B for details].
However, the weakness of the Fourier domain is that it does not
economically represent signals
with singularities such as im-
ages with edges [see Fig. 2(b)]. Consequently, as dictated by the
MSE bound in (6), any estimate obtained via Fourier shrinkage
is unsatisfactory with a large MSE; the estimate is either noisy
or distorted for signals
with singularities [see Fig. 4(c), for
example].
Recently, the wavelet domain (with shifts and dilates of a
mother wavelet function as
’s) has been exploited to estimate
from , for example, Donoho’s wavelet-vaguelette deconvo-
lution (WVD) [5]. The strength of the wavelet domain is that it
economically represents classes of signals containing singular-
ities that satisfy a wide variety of local smoothness constraints,
including piecewise smoothness and Besov space smoothness
[see Fig. 2(b) and Section V-B for details]. However, the weak-
ness of the wavelet domain is that it typically does not econom-
ically represent the colored noise
[see Fig. 2(a)]. Con-
sequently, as dictated by the MSE bound (6), any estimate ob-
tained via wavelet shrinkage is unsatisfactory with a large MSE;
the estimate is either noisy or distorted for many types of
.
Unfortunately, no single transform domain can economically
represent both the noise colored by a general
and sig-
nals from a general smoothness class [5]. Hence, deconvolution
techniques employing shrinkage in a single transform domain
cannot yield adequate estimates in many deconvolution prob-
lems of interest.
C. Fourier-Wavelet Regularized Deconvolution (ForWaRD)
In this paper, we propose a deconvolution scheme that re-
lies on tandem scalar processing in both the Fourier domain,
420 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 52, NO. 2, FEBRUARY 2004
Fig. 3. Fourier-wavelet regularized deconvolution ( ForWaRD ). ForWaRD employs a small amount of Fourier shrinkage (most
1
) to partially attenuate
the noise amplified during operator inversion. Subsequent wavelet shrinkage (determined by
) effectively attenuates the residual noise.
which economically represents the colored noise , and
the wavelet domain, which economically represents signals
from a wide variety of smoothness classes. Our hybrid Fourier-
Wavelet Regularized Deconvolution ( ForWaRD ) technique es-
timates
from by first employing a small amount of scalar
Fourier shrinkage
and then attenuating the leaked noise with
scalar wavelet shrinkage
(see Fig. 3) [9], [10].
Here is how it works: During operator inversion, some Fourier
coefficients of the noise
are significantly amplified; just a
small amount of Fourier shrinkage (most
) is sufficient
to attenuate these amplified Fourier noise coefficients with min-
imal loss of signal components. The leaked noise
that
Fourier shrinkage
fails to attenuate [see (5)] has significantly
reduced energy in all wavelet coefficients, but the signal part
that Fourier shrinkage retains continues to be economically
represented in the wavelet domain. Hence, subsequent wavelet
shrinkage effectively extracts the retained signal
from the
leaked noise
and provides a robust estimate.
For an idealized ForWaRD system, we will derive the optimal
balance between the amount of Fourier shrinkage and wavelet
shrinkage by optimizing over an approximate MSE metric. We
will find that signals with more economical wavelet representa-
tions require less Fourier shrinkage.
Fig. 4 illustrates the superior overall visual quality and lower
MSE of the ForWaRD estimate as compared with the LTI
Wiener filter estimate [3], [4] for the 2-D box-car blur operator,
which models rectangular scanning aperture effects [1], with
impulse response
for and
0 otherwise (see Section VIII for details). For this operator,
the WVD approach returns an esentially zero estimate; scalar
wavelet shrinkage cannot salvage the signal components since
nearly all wavelet coefficients are corrupted with high-variance
noise.
Indeed, even in problems for which the WVD was designed,
we will prove that the ForWaRD MSE also decays with the same
optimal WVD rate as the number of samples increases. Further,
for such problems, we will experimentally demonstrate For-
WaRD ’s superior MSE performance compared with the WVD
over a wide range of practical sample sizes [see Fig. 6(a)].
D. Related Work
Kalifa and Mallat have proposed a mirror-wavelet basis ap-
proach that is similar to the WVD but employs scalar shrinkage
in a mirror-wavelet domain adapted to the colored noise
instead of shrinkage in the conventional wavelet domain [2]. Al-
though the adapted basis improves on the WVD performance
in some “hyperbolic” deconvolution problems, similarly to the
WVD, it provides inadequate estimates for arbitrary convolu-
tion operators. For example, for the ubiquitous box-car blur
,
again, most signal components are lost during scalar shrinkage
due to high-variance noise. Fig. 7(b) illustrates that ForWaRD
is competitive with the mirror-wavelet approach, even for a hy-
perbolic deconvolution problem.
Similar to ForWaRD, Nowak and Thul [11] have first em-
ployed an under-regularized system inverse and subsequently
used wavelet-domain signal estimation. However, they do not
address the issue of optimal regularization or asymptotic per-
formance.
Banham and Katsaggelos have applied a multiscale Kalman
filter to the deconvolution problem [12]. Their approach em-
ploys an under-regularized, constrained-least-squares prefilter
to reduce the support of the state vectors in the wavelet do-
main, thereby improving computational efficiency. The amount
of regularization chosen for each wavelet scale is the lower
bound that allows for reliable edge classification. While similar
in spirit to the multiscale Kalman filter approach, ForWaRD em-
ploys simple Wiener or Tikhonov regularization in the Fourier
domain to optimize the MSE performance. In addition, For-
WaRD employs simple scalar shrinkage on the wavelet coeffi-
cients in contrast to more complicated prediction on edge and
nonedge quad-trees [12]. Consequently, as discussed in Sec-
tion VI-D, ForWaRD demonstrates excellent MSE performance
as the number of samples tends to infinity and is, in fact, asymp-
totically optimal in certain cases. Further, as demonstrated in
Section VIII, ForWaRD yields better estimates than the multi-
scale Kalman filter approach.
There exists a vast literature on iterative deconvolution
techniques; see [4], [13]–[15], and the references therein. In
this paper, we focus exclusively on noniterative techniques for
the sake of implementation speed and simplicity. Nevertheless,
many iterative techniques could exploit the ForWaRD estimate
as a seed to initialize their iterations; for example, see [16].
E. Paper Organization
We begin by providing a more precise definition of the con-
volution setup (1) in Section II. We then discuss techniques that
employ scalar Fourier shrinkage in Section III. After briefly re-
viewing wavelet theory in Section IV, we introduce the WVD
technique in Section V. We present the hybrid ForWaRD scheme
in Section VI and discuss its practical implementation in Sec-
tion VII. Illustrative examples lie in Section VIII. We conclude
and sketch future directions in Section IX. A short WVD review
in Appendix A and technical proofs in Appendices B–D com-
plete the paper.
II. S
AMPLING AND DECONVOLUTION
Most real-life deconvolution problems originate in contin-
uous time and are then sampled. In this section, we sketch the
relationship between such a sampled continuous-time setup and
the setup with discrete-time circular convolution considered in
this paper [see (1)].
NEELAMANI et al.: ForWaRD: FOURIER-WAVELET REGULARIZED DECONVOLUTION FOR ILL-CONDITIONED SYSTEMS 421
(a) (b)
(c) (d)
Fig. 4. (a) Desired Cameraman image
x
(256
2
256 samples). (b) Observed image
y
:
x
smoothed by a two–dimensional (2-D) 9
2
9 box-car blur plus white
Gaussian noise with variance such that the BSNR
=
40 dB. (c) LTI Wiener filter estimate (SNR
=
20.8 dB, ISNR
=
5.6 dB). (d) ForWaRD estimate (SNR
=
22.5
dB, ISNR
=
7.3 dB). See Section VIII for further details.
Consider the following sampled continuous-time deconvo-
lution setup: An unknown finite-energy desired signal
is
blurred by linear convolution (denoted by
) with the known fi-
nite-energy impulse response
of an LTI system and then
corrupted by an additive Gaussian process
to form the ob-
servation
. For finite-support and
, the finite-support can be obtained using circular
convolution with a sufficiently large period. For infinite-support
and , the approximation of using circular
convolution can be made arbitrarily precise by increasing the
period. Hence, we assume that the observation
over a nor-
malized unit interval can be obtained using circular convolu-
tion with a unit period, that is,
with
. Deconvolution aims to estimate from the sam-
ples
of the continuous-time observation . For example,
can be obtained by averaging over uniformly spaced
intervals of length
(7)
Other sampling kernels can also be used in (7);
2
see [17] and
[18] for excellent tutorials on sampling. Such a setup encapsu-
lates many real-life deconvolution problems [1].
The observation samples
from (7) can be closely ap-
proximated by the observation
from setup (1) [1]; that is
(8)
if the continuous-time variables
, , and com-
prising
are judiciously related to the discrete variables
2
For example, impulse sampling samples at uniformly spaced time instants
t
=
n=N
to yield
z
(
n
)=(
t
)
.
422 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 52, NO. 2, FEBRUARY 2004
, , and comprising . We choose to define
. The so defined can be
assumed to be AWGN samples with nonzero variance
;
the large bandwidths of noise processes such as thermal
noise justify the whiteness assumption [1]. Let
and
denote signals obtained by first making and
periodic and then bandlimiting the resulting signals’ Fourier
series to the frequency
Hz (for antialiasing). We define
and define as uniformly
spaced (over
) impulse samples of .
With these definitions, we can easily show that the error
.
For all finite-energy
and , both
and decay to zero with increasing
because they represent the norm of the aliasing components of
and , respectively. Consequently, soon
becomes negligible with respect to the noise variance
and
can be ignored. Hence, solutions to estimate
from the
in (1)—the focus of this paper—can be directly applied
to estimate
from . For a wide range of Besov space
signals, the estimate of
can then be interpolated with
minimal error to yield a continuous-time estimate of
,as
sought in (7) [19], and [20].
3
In Sections V and VI, we will analyze the MSE decay rate
(in terms of
) of the WVD and ForWaRD solutions to the
setup (1) as the number of samples
. At each ,we
assume that the corresponding
and in (1) originate
from an underlying continuous-time
and , as defined
above. Further, we assume that the corrupting
in (1) are
AWGN samples with variance
that is invariant with .
III. F
OURIER-BASED REGULARIZED
DECONVOLUTION (FORD)
A. Framework
The Fourier domain is the traditional choice for deconvolu-
tion [4] because convolution simplifies to scalar Fourier opera-
tions. That is, (1) can be rewritten as
(9)
where
, , , and are the respective length- discrete
Fourier transforms (DFTs) of
, , , and , and ,
(assuming is even) are the nor-
malized DFT frequencies. Rewriting the pseudo-inversion op-
eration [see (2)] in the Fourier domain
if
otherwise
(10)
where
is the DFT of , clearly demonstrates that noise com-
ponents where
are particularly amplified during
operator inversion.
3
The Besov space range is dictated by the smoothness of the sampling kernel.
Let
x
(
t
)
2
Besov space
B
(see Section IV-B for the notation). Then, if the
sampling kernel of (7) is employed, then the interpolation error is negligible
with respect to the estimation error for the range
s>
1
=p
0
1
=
2
; the range
decreases to
s>
1
=p
if impulse sampling is employed [19], [20].
Deconvolution via Fourier shrinkage, which we call Fourier-
based Regularized Deconvolution (FoRD), attenuates the am-
plified noise in
with shrinkage
(11)
The
, commonly referred to as regularization terms
[4], [21], control the amount of shrinkage. The DFT components
of the FoRD estimate
are given by
(12)
The
and comprising denote the respective
DFTs of the retained signal
and leaked noise
components that comprise the FoRD estimate [see (5)].
Typically, the operator inversion in (10) and shrinkage in (12)
are performed simultaneously to avoid numerical instabilities.
Different FoRD techniques, such as LTI Wiener deconvolu-
tion [3], [4] and Tikhonov-regularized deconvolution [21], differ
in their choice of shrinkage
in (12). LTI Wiener deconvolu-
tion sets
(13)
with regularization parameter to shrink more (that is,
) at frequencies where the signal power is small
[3], [4]. Tikhonov-regularized deconvolution, which is similar
to LTI Wiener deconvolution assuming a flat signal spectrum
, sets
(14)
with
[21]. Later, in Section VI, we will put both of these
shrinkage techniques to good use.
B. Strengths of FoRD
The Fourier domain provides the most economical represen-
tation of the colored noise
in (2) because the Fourier
transform acts as the Karhunen-Loeve transform [22] and
decorrelates the noise
. Consequently, among all linear
transformations, the Fourier transform captures the maximum
colored noise energy using a fixed number of coefficients
[23]. This economical noise representation enhances FoRD
performance because the total FoRD MSE is lower bounded by
[5].
4
The best
possible FoRD MSE
is achieved using the LTI Wiener deconvolution shrinkage of
(13) in (12) [7]. When the signal
in (2) also enjoys an econom-
4
The factor
N
arises because
j
X
(
f
)
j
=
N
j
x
(
k
)
j
for any signal
x
.
![Fig. 7. MSE performance of ForWaRD compared to the mirror-wavelet basis approach [2] at different BSNRs.](/figures/fig-7-mse-performance-of-forward-compared-to-the-mirror-351jhst8.webp)
![Fig. 6. MSE performance of ForWaRD compared with the WVD [5] for different N . The MSE incurred by ForWaRD’s wavelet shrinkage step decays much faster than the WVD’s MSE with increasing N , whereas the Fourier distortion error saturates to a constant.](/figures/fig-6-mse-performance-of-forward-compared-with-the-wvd-5-for-8n808m4u.webp)
