Efficient Implementations Of 2 D Noncausal Iir Filters Circuits ...

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
549
Efﬁcient Implementations of
2-D Noncausal IIR Filters
Michael M. Daniel, Student Member, IEEE, and Alan S. Willsky, Fellow, IEEE
Abstract— In this paper, we propose a framework for the
in terms of a difference equation. In either case, the difference
efﬁcient implementation of two-dimensional (2-D) noncausal inﬁ-
equation by itself isn’t sufﬁcient to completely specify the
nite impulse response (IIR) ﬁlters, i.e., ﬁlter systems described
ﬁltering algorithm, as one must also specify a set of auxiliary
implicitly by difference equations and boundary conditions. A
number of common 2-D LSI ﬁlter operations, (including low-
conditions. In 1-D, for the most part these are speciﬁed
pass, high-pass, and zero-phase ﬁlters), are efﬁciently realized
as a set of initial conditions, leading to causally-recursive
and implemented in this paper as noncausal IIR ﬁlters. The basic
ﬁltering algorithms with computational load per sample point
concepts involved in our approach include the adaptation of so-
proportional to the order of the difference equation. Moreover,
called direct methods for solving partial differential equations
even for noncausal 1-D ﬁlters like zero-phase IIR ﬁlters,
(PDE’s), and the introduction of an approximation methodology
that is particularly well suited to signal processing applications
implementation results in a per-sample computational load
and leads to very efﬁcient implementations. In particular, for an
proportional to ﬁlter order or, equivalently, to the total number
input and output with N 2 N samples, the algorithm requires
of auxiliary (initial and ﬁnal) conditions, (assuming that the
only O(N 2) storage and computations (yielding a per pixel
1-D noncausal IIR ﬁlters are implemented through the com-
computational load that is independent of image size), and has
a parallel implementation (yielding a per pixel computational
bination of a causal recursion requiring initial conditions and
load that decreases with increasing image size). Also, because
an anticausal recursion requiring “ﬁnal” conditions).
our approach allows for the implementation of ﬁlters with space-
In contrast, the dimension of the required auxiliary condi-
varying coefﬁcients on irregularly shaped domains, it should have
tions in 2-D depends not only on the order of the difference
applications in related areas like linear estimation, geophysical
equation but also on the size of the boundary. Since the size
signal processing, or any ﬁeld requiring approximate solutions to
elliptic PDE’s.
of the boundary is proportional to the dimensions of the 2-
D domain of interest, an apparently signiﬁcant increase in
computational complexity results. In addition, since in most
I. INTRODUCTION
2-D applications there is no natural ordering of the sample
FORtwo-dimensional(2-D)signalprocessingapplications, pointsandnonaturaldirectionforrecursion,thereisnoreason
ﬁnite impulse response (FIR) ﬁlters have been over-
to expect that the auxiliary conditions would separate into
whelmingly preferred to inﬁnite impulse response (IIR) ﬁlters
anything that might resemble “initial” or “ﬁnal” conditions, but
[3], [10], [12]. Among the reasons for this preference are:
rather would more naturally be distributed around the entire
a) FIR ﬁlters can be efﬁciently implemented in both one-
boundary of the 2-D domain, leading to 2-D noncausal IIR
dimensional (1-D) and 2-D through the use of the FFT; and
(2-DNC-IIR) ﬁlters that are not recursively computable.
b) FIR ﬁlters are always stable and do not require any notion
On the other hand, if effective methods of implementation
of recursion or ordering of the sample points (in 1-D or in
for 2-DNC-IIR ﬁlters were available, there would be numerous
2-D) in order to be implemented. In contrast, 1- and 2-D IIR
possibilities for their application. For example, one potential
ﬁlters appear to be dramatically different. First, there is usually
advantage retained in 2-D for IIR ﬁlters is that a given set
no natural ordering of the sample points in 2-D, and 2-D IIR
of frequency response characteristics typically may be met by
ﬁlters are difﬁcult to test for stability. More importantly, for
an IIR ﬁlter of considerably lower order than a corresponding
2-D noncausal IIR ﬁlters the lack of efﬁcient implementations
FIR design. Moreover, 2-DNC-IIR ﬁlters arise naturally in
has both limited the investigation of these ﬁlters and led many
applications such as the modeling of random ﬁelds for image
to argue that they cannot be implemented in practice [3], [10],
processing [1], [2], [11] and computer vision [9]. In this
[12].
paper we present an efﬁcient implementation of 2-DNC-IIR
To understand these issues, as well as our approach to
ﬁlters that overcomes the difﬁculties we have described, thus
dealing with them, consider an IIR ﬁlter, in 1- or 2-D, speciﬁed
offering the possibility of recapturing in 2-D the computational
advantages and ﬂexibility that IIR ﬁlters have in 1-D.
Manuscript received March 24, 1995; revised October 24, 1996. This work
was supported by the Ofﬁce of Naval Research under Grant N00014-91-J-
The key to our approach is the recognition of both the
1004, by the Air Force Ofﬁce of Scientiﬁc Research under Grant F49620-95-
similarities and differences between the implementation of
1-0083, and by the Army Research Ofﬁce under Grant DAAL03-92-G-0115
2-DNC-IIR ﬁlters and the solution of ﬁnite-difference and
(Center for Intelligent Control Systems). This paper was recommended by
Associate Editor W.-S. Lu.
ﬁnite-element approximations of partial differential equations
The authors are with the Laboratory for Information and Decision Systems,
(PDE’s). In particular, the equations resulting from such
Massachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail:
willsky@mit.edu).
methods are of the same form as those for 2-DNC-IIR ﬁlters,
Publisher Item Identiﬁer S 1057-7130(97)03652-5.
and thus the many methods that have been developed for
1057–7130/97$10.00 © 1997 IEEE

550
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
the efﬁcient solution of PDE’s can be used to implement 2-
The order of this difference equation is deﬁned to be
.
DNC-IIR ﬁlters. These methods by themselves, while offering
However, (1) provides only a partial speciﬁcation of a system,
considerable savings in computational complexity and storage,
as it must be accompanied by a set of auxiliary conditions. If
may not reduce these loads enough to make 2-DNC-IIR ﬁlters
the ﬁlter whose input and output satisﬁes (1) is stable, then
attractive. However, by taking advantage of a fundamental
the auxiliary conditions cannot generally be organized as a
difference in objective between solving PDE’s and performing
simple set of “initial” or “ﬁnal” conditions, as considered in
2-D ﬁltering, we can reduce the computational complexity
[3], [12], but must be speciﬁed around the entire boundary of
even further, resulting in implementations with complexity
the 2-D ﬁltering domain; these auxiliary conditions are referred
per 2-D data point independent of domain size—the same
to as boundary conditions (BC’s). If BC’s are speciﬁed, no
attractive feature as in 1-D. While the focus in PDE’s is
simple recursive solution is possible, and all of the output
typically on obtaining numerically very accurate solutions to
values
must in principle be computed simultaneously.
speciﬁc 2-D difference equations, and hence accurate solutions
Algorithms for computing the outputs of such 2-DNC-IIR
to the corresponding physical models, in 2-D ﬁltering the
ﬁlters are discussed in Sections III and IV. The remainder of
difference equation is not the fundamental object. Instead,
this section addresses the some of the issues raised by the
the initial criterion is a set of ﬁlter or frequency response
imposition of boundary conditions on 2-D IIR ﬁlters.
speciﬁcations. A 2-D difference equation is then chosen to
A fundamental issue is the effect of the boundary conditions
meet these speciﬁcations within some tolerance. Consequently,
upon the response of 2-DNC-IIR ﬁlters. The effect of boundary
approximations to the solution of the difference equation are
conditions is an important issue for any system deﬁned on
acceptable as long as they lead to ﬁlters that also meet the
a ﬁnite domain, even in 1-D and for FIR ﬁlters, although
desired tolerances.
this issue is rarely addressed [3], [12]. For IIR ﬁlters, both
In the next section, we introduce the class of 2-DNC-IIR
in 1- and 2-D, the method for limiting the effect of the
ﬁlters and discuss the role of boundary conditions in these
BC’s upon the ﬁlter output is to require that the system be
systems. In Section III, we make the connection between
stable. Not only does stability guarantee that the effect of
implementations of 2-DNC-IIR ﬁlters and methods for solv-
the BC’s will be limited to the boundary regions and decay
ing sparse linear systems of equations, such as those which
with distance from the boundary, stability implies a particular
arise when solving PDE’s. One of these methods involves
choice of BC’s. To illustrate this subtle relationship between
organizing the 2-D data points into 1-D columns, whose
the choice of BC’s and stability, ﬁrst consider 1-D IIR ﬁlters,
dimensions are equal to the linear dimension of the ﬁltering
for which analysis is simpler. For the 1-D IIR ﬁlter whose input
domain. The PDE or ﬁlter solution is then given by processing
and output satisfy
, the imposition
these columns sequentially. While this algorithm is generally
of stability implies that the boundary condition is either an
inefﬁcient, if we view each of these sequential processing
initial rest or a ﬁnal rest condition. A value of
steps as being itself a 1-D processing procedure along the
implies an initial rest condition and
implies a ﬁnal
1-D data set, we are led to the idea of approximating this
rest condition. For higher-order ﬁlters, stability leads to three
step using low-order IIR ﬁltering methods. This idea, which is
possible boundary conditions, where again the exact choice of
developed in Section IV, results in very efﬁcient 2-DNC-IIR
boundary conditions is determined by the ﬁlter coefﬁcients:
ﬁltering procedures applicable to a large class of noncausal,
a) if all the poles are inside the unit circle, then the BC’s
nonseparable ﬁltering applications. In Section V, the efﬁcient
are initial rest conditions; b) if all the poles are outside the
implementation is applied to several 2-DNC-IIR ﬁlters, some
unit circle, then the BC’s are ﬁnal rest conditions; and c) if
of which are zero-phase. Zero-phase ﬁlters are of considerable
the poles are both inside and outside the unit circle, then the
interest in practice, and the apparent difﬁculty in implementing
system is noncausal and the BC’s correspond to both initial
2-D IIR ﬁlters with zero phase has often been cited as one of
and ﬁnal rest conditions. The BC’s for the third possibility
the reasons that FIR ﬁlters are commonly used [10]. We now
can be seen by splitting the IIR ﬁlter into a parallel realization
can implement zero-phase IIR ﬁlters efﬁciently, removing a
of a causal ﬁlter (with poles inside the unit circle) and an
major obstacle to their use in practice.
anticausal ﬁlter (with poles outside the unit circle), where the
causal ﬁlter satisﬁes initial rest and the anticausal ﬁlter satisﬁes
ﬁnal rest. This parallel realization is given by a partial fraction
II. TWO-DIMENSIONAL IIR FILTERS
expansion of the system’s frequency response. The conclusion
AND BOUNDARY CONDITIONS
to be drawn from this discussion is that, for 1-D IIR ﬁlters,
For a rich class of 2-D IIR ﬁlters, the inputs
and
stability both determines the form of the boundary conditions
outputs
satisfy linear constant-coefﬁcient difference
and limits the effect of the boundary condition upon the system
equations (LCCDE’s) of the form
response, i.e., the transients, to the region in which the BC’s
are applied. However, the rate of decay of the transients is a
function of the ﬁlter difference equation, e.g., the magnitude
of
in the ﬁrst-order example.
For 2-D IIR ﬁlters, the link between stability and boundary
conditions is much the same. First of all, given a stable ﬁlter
(1)
satisfying (1), the frequency response follows immediately
from the difference equation. Stability in this case corresponds

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
551
interval, say
. For instance, how would one even deﬁne
the shift of a signal
deﬁned only on
. One option
is to deﬁne
to be equal to zero outside
, and
then to compute
for all
. To guarantee
shift-invariance, the IIR ﬁlter is implemented as a parallel
realization of a causal IIR recursion initialized by initial rest
conditions and an anticausal IIR recursion initialized by ﬁnal
rest conditions. However, while the resulting ﬁlter is shift-
invariant, such parallel realizations do not generally exist for
2-DNC-IIR ﬁlters, unless the ﬁlter is separable. Furthermore,
there are no analogous notions of initial and ﬁnal rest for 2-D
ﬁlters. A notion of shift-invariance which does extend to 2-
Fig. 1.
The elements of N (denoted by ) and @ (1; 1)
N
(denoted by ),
DNC-IIR ﬁlters is given by comparing the following two 1-D
drawn for N = 5.
signals: 1) the response
deﬁned on
to an input
to the usual notion of bounded-input bounded-output stability,
deﬁned on
, and 2) the response
deﬁned
as well as to the concept that the effect due to the BC’s on
on
to the input
which is also
the ﬁlter response near the center of the 2-D domain decays
deﬁned on
. If the system is shift-invariant,
to zero as the boundaries recede to inﬁnity. Again, the form
then
for
. When
of the boundary conditions is completely determined by the
extended to 2-D systems, this notion of shift-invariance applies
imposition of stability, and the width of the annular region near
to 2-DNC-IIR ﬁlters, i.e., ﬁlters satisfying (1) and constrained
the boundaries where the BC’s signiﬁcantly effect the ﬁlter
by boundary conditions.
response is determined by the coefﬁcients of the 2-D difference
III. NONCAUSAL IIR FILTERS AS
equation. However, a major difference between IIR ﬁlters in 1-
LINEAR SYSTEMS OF EQUATIONS
and 2-D is the ability to determine BC’s which lead to stable
systems. As noted earlier, determining such BC’s for 1-D IIR
A. Direct versus Iterative Methods
ﬁlters is straightforward, and follows from a partial fraction
In this section we make precise the connection between the
expansion of the frequency response. For 2-D IIR ﬁlters, there
problem of implementing 2-DNC-IIR ﬁlters and the general
is no general method for determining the BC’s which lead to a
problem of solving large, sparse, sets of linear equations, in
stable system. The lack of such method is due to the inability
particular those arising in the solution of linear PDE’s. The
to factor a 2-D system function and to the complexity of the
methods that result from this connection are quite broadly
boundaries in 2-D, which can cover large regions and have
applicable. For example, our methodology can be used for
complicated geometries. However, as shown in the Section V,
linear difference equations which are not constant-coefﬁcient,
we can often ﬁnd boundary conditions which lead to stable
for regions of support
which are nonsquare and irregu-
ﬁlters. As the algorithms proposed in Sections III and IV
larly sampled, and for various types of boundary conditions.
are motivated by numerical solutions to PDE’s, the boundary
However, for notational simplicity in this and the following
conditions chosen in Section V are discrete equivalents of
sections, we assume that the difference equation is LCCDE,
boundary conditions which lead to stable solutions of PDE’s,
that the ﬁlter domain is square, and that the boundary condi-
where stability again refers to limiting the region in which
tions are of the Dirichlet type. One square domain of
the BC’s signiﬁcantly effect the solution of the PDE. Two
samples is
. For an order
such conditions are Dirichlet and Neumann conditions [15].
ﬁlter, the corresponding Dirichlet conditions are to
For a 2-DNC-IIR ﬁlter of order (1, 1), Dirichlet conditions
set
to some known function
on the annular
correspond to specifying the value of
on the boundary
ring
of the ﬁltering domain. For a ﬁlter of higher order, Dirichlet
. Note that the width of
conditions correspond to setting the value of
on an
this annular ring is determined by the order of the difference
annular ring around the boundary of the ﬁltering domain. As
equation. Both
and
are illustrated in Fig. 1 for a
will be described more clearly in the next section, the width
and a ﬁlter of order (1, 1).
of this annular ring is a function of the ﬁlter order.
For clarity of exposition, we also assume that
Another issue raised by the imposition of boundary con-
in (1), where
is the Kronecker delta function.
ditions is that traditional notions of shift-invariance do not
Since implementing the right-hand side of (1) is equivalent
easily extend to 2-DNC-IIR ﬁlters. To appreciate this subtle
to implementing an FIR ﬁlter, a more complicated right-hand
issue, ﬁrst consider 1-D IIR ﬁlters. Shift-invariance for a 1-
side adds only notational but not conceptual complexity. These
D difference equation requires that a) the locations of the
assumptions lead to 2-DNC-IIR ﬁlters whose inputs and output
boundary conditions are not ﬁxed, but instead adjust to the
satisfy the difference equation
location of the nonzero values of the input, e.g., initial rest
conditions, and b) both the input and output are deﬁned for all
time, i.e.,
. As a consequence, the commonly
deﬁned property of shift-invariance cannot be applied to
(2)
systems with inputs and outputs deﬁned only over a ﬁnite

552
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
conditioned conjugate gradient or multigrid, might be just as
or more effective for some applications, (especially for 3-D
problems), but we show here that direct methods allow for
very efﬁcient implementations of classes of 2-D ﬁlters.
The LU factorization
yields a unit lower-triangular
matrix
and an upper-triangular matrix
. Given
and
,
the solution to
can be found very efﬁciently by sequen-
tially solving the following two triangular systems:
and
. Solving for
is called forward-substitution, while
Fig. 2.
The output mask for the 9-point NNM difference equation.
solving for
from
is called back substitution. Assuming
has dimension
, solving for
by explicitly computing
for all
and satisfy the Dirichlet conditions
and then computing
requires
computa-
for all
.
tions (measured in terms of ﬂoating point adds and multiplies).
Equation (2) can be cast in matrix form as
If
is dense, the LU factorization approach yields at most
minimal computational savings, as the LU factorization alone
requires
computations, while the substitutions require
(3)
only
additional computations. However, if
is sparse, as
where the nonzero elements of
are the ﬁlter coefﬁcients
for 2-DNC-IIR ﬁlter systems, the savings in both storage and
. Vectors
and
contain the ﬁlter input
and output
computations can be tremendous, especially if proper orderings
, respectively, in
, and
contains the contribution of
are used to minimize the amount of ﬁll-in (loss of sparsity)
the Dirichlet conditions entering through the ﬁlter difference
that occurs during the factorization (see [4]). Amortizing the
equation. The order in which the variables
appear
costs of the factorization over a large number of ﬁlter inputs
in
is the ordering of
, or the ordering of
. For
further decreases the effective computation requirements for
direct methods (see below), this ordering can drastically alter
the LU approach. In our application, however,
is equal to
apparent complexity.
the number of 2-D data points,
, and thus computations of
Note that
has dimension
. A nice property
order greater than linear in
can still make this approach
of IIR ﬁlters is that they generally require a small number
prohibitive. Fortunately, as we will see, in the context of 2-D
of coefﬁcients, so that
and
. In other
ﬁltering there are natural and very accurate approximations to
words,
will be very sparse. This obviously suggests the
the LU factorization approach that do result in total complexity
use of numerical methods developed for solving large sparse
that is linear in
.
systems which take advantage of this sparsity to minimize
B. Columnwise Orderings
computational and storage requirements. In particular, there
are two distinct classes of methods for calculating the output
To make the following discussion explicit we focus here on
in (3)—iterative and direct methods. Iterative methods begin
a common 2-D LCCDE of order (1, 1), the 9-point nearest
with an estimate
of , and produce at each step an estimate
neighbor model (NNM). This difference equation also arises
which theoretically converges as
; however,
quite frequently in engineering applications [6], [9], [11],
in practice the series must converge within a tolerable error in
[15], most notably as the ﬁrst-order and second-order ﬁnite-
a ﬁnite number of steps. Direct methods consist of variants of
difference and ﬁnite-element approximations to elliptic PDE’s.
the LU factorization [6], [7], and produce the exact solution
The constant-coefﬁcient form of the 9-point NNM is given by
(disregarding numerical errors) in a ﬁnite number of steps.
For signal processing applications, the same ﬁlter is typi-
cally applied to a large number of inputs, and thus
must
be factored only once for a direct method. This factorization
can be done off-line, i.e., the factorization costs can either
(4)
be considered part of the ﬁlter design process or amortized
over the large number of inputs. This property of direct
The output mask of this difference equation is illustrated in
methods motivates us to focus here on direct implementations
Fig. 2. Note that the LSI system characterized by the frequency
of 2-DNC-IIR ﬁlter systems. Iterative methods, such as pre-
response from difference equation (4) is zero-phase if
,
. .
.
.
.
.
.
.
. .
. .
.
.
.
.
.
.
(5)
. .
.
.
. .

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
553
,
, and
. Consider a 2-DNC-IIR
Although
,
, and
are very sparse, the matrices
and
ﬁlter which satisﬁes (4) on
and Dirichlet conditions on
are generally full, with the exception of
. Storing
. If
and
are ordered columnwise into
each of these matrices requires
storage elements and
dimensional vectors
computing each
requires
computations, leading to
and
, respectively, (3)
a total of
storage elements and
computations
becomes1 (5), as shown at the bottom of the previous page.
for the entire factorization.
The structure of
in (5) is block tridiagonal, and the
The lack of sparsity in the blocks of (8) also implies a large
dimensional blocks
,
, and
are tridiagonal. Note that
computational burden for the on-line solution. First note that,
(5) allows for a space-varying NNM difference equation, but
since the factorization
is needed to compute
for the constant-coefﬁcient difference equation the subscripts
at each step of the recursion, these factors can be stored
on the blocks of
can be dropped. In this case, the nonzero
in place of
in (8). The solution to (5) is then given by
elements of
,
, and
are given by
forward-substitution, (initialized by
),
(9)
followed by back-substitution, (initialized by
)
(10)
Since
and
are generally full, the solution of (9) and
(10) requires
computations per step and hence
total computations. Thus, not only is the off-line computational
load
, but the on-line storage and computations are both
, signiﬁcantly greater than the
goal.
For ﬁlters of order
, a columnwise ordering of
However, if an approximate solution can be tolerated, the
leads to a matrix
which has block bandwidth
, while each
block LU factorization based on the columnwise ordering
of the blocks has bandwidth
. (A matrix
has bandwidth
leads to an efﬁcient approximation strategy which achieves
if element
is nonzero only for
.)
the goal of
storage elements and
computations
A simple recursive algorithm can be invoked to com-
for both the off-line factorization and the on-line solution.
pute the factorization of
in (5). This algorithm factors
In particular, note that (9) requires ﬁrst solving the lower
“block-by-block,” and is thus referred to as the block LU
triangular equations
factorization [4], [7]. This algorithm recursively computes the
(11)
matrices
and
using the following
equations:2
followed by solving the upper triangular system
“compute
”
(6)
(12)
“compute
”
(7)
Recall that we have organized our variables into 1-D columns,
The recursion is initialized with
. Solving (6) at each
and thus the solution of the lower triangular system (11) can be
step is performed by an LU factorization on
, and, as we
thought as a causal 1-D recursion, beginning at the bottom of
discuss next, this factorization
is needed on-
the column (
) and proceeding recursively to the top of the
line. For these recursions to be well-posed, the matrices
column (
). The upper triangular system (12) corresponds
must be invertible for all
. Conditions which
to an anticausal recursion proceeding from top to bottom. The
guarantee this are discussed in [7]. For the examples presented
back-substitution ﬁltering, (10), requires implementing an FIR
in Section V, the matrices
are invertible.
ﬁlter, in the form of a matrix multiplication, along a single
The block LU factorization resulting from the procedure
column of
.
just described yields
Thus we can view (11) and (12) as 1-D recursive ﬁltering
operations, albeit shift-varying recursions, since in general
.
and
will not be Toeplitz. If
and
are full,
.
.
.
. .
.
then the order of these recursive ﬁlters equals the length
.
.
(8)
.
. .
of the column, and it is the need to determine (off-
line) and then implement (on-line) these high-order recursions
While the recursion given by (6) and (7) is conceptually
that leads to the severe computational burden. However, if
straightforward, its computational load can be overwhelming.
these recursive 1-D ﬁlters can be approximated by lower-
This is a direct result of the columnwise ordering, which leads
order recursions—e.g., if
and hence
and
can
to a destruction of the sparsity of
during the factorization.
be approximated by banded matrices, then both the storage
and computational requirements for the forward-substitution
1 For any matrix in this paper, such as A in (5), block entries not indicated
phase of the on-line solution can be reduced to
. For
are zero.
2 The validity of the recursions (6) and (7) can be veriﬁed directly by
the back-substitution, the computational burden is governed
equating A in (5) with the expression in (8).
at each step by multiplication of the matrix
with
.

554
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
Since
is generally full and not Toeplitz, this operation
This leads to the following approximation to (6) and (7).
will require
computations per step. However, if we
Suppose that
is an approximation to
which is -banded
similarly approximate
with a lower-order FIR ﬁlter, i.e.,
(i.e., banded with bandwidth ). Note that
is exactly
by approximating
with a banded matrix, the total com-
banded, so we set
. We then compute a
-banded
putational and storage requirements for the on-line solution
approximation of
:
reduce to the
goal.
Note, however, to reduce the off-line computational load
(13)
to
, it is not sufﬁcient that
and
are well
otherwise
approximated by matrices with narrow bandwidth; a method
where, from the corollary,
requires
must exist for determining these approximate matrices in
computations. Assuming that
is a good approx-
computations per stage. In the next section, we describe
imation to
and that
is a good approximation
such an approximation procedure in detail.
to
, then
. Since
is
IV. EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
tridiagonal, the product
requires
computations. Note, however, that this product has bandwidth
A. Development of the Approximate Block LU Algorithm
, which will in turn require a growing bandwidth at each
The approximate implementation of 2-DNC-IIR ﬁlters de-
step. However, under the key assumption that the matrices of
scribed in this section is motivated by the fact that, for
interest are approximately
-banded, we can neglect elements
many ﬁlters, a small number of elements in the blocks
,
outside the
-bandwidth. Thus, deﬁne the truncation operator
, and
dominate the rest of the elements. An efﬁcient
as
approximation to the on-line solutions follows by setting to
zero the insigniﬁcant elements of
,
, and
. Recursions
(14)
(9) and (10) then can be implemented very efﬁciently if one
otherwise
takes care to avoid operating on the zero elements.
which yields the following approximation to (6):
However, the approach of simply discarding the insignif-
(15)
icant elements of (8) is not enough. First, the off-line fac-
torization will still require
computations. Secondly,
Similarly, substituting
into (7) in lieu of
yields
searching for the signiﬁcant elements of each block in (8) and
, requiring
calculations. Once
storing them in data structures for efﬁcient implementation
again, this yields a matrix that is (
)-banded, and applying
can be a costly procedure. Many ﬁlters, however, have a
our assumption of -bandedness, we obtain our approximation
property which allows us to overcome these difﬁculties. For
to (7):
these ﬁlters, all of the matrices of interest—
,
, and
—can be well approximated by banded matrices of some
(16)
bandwidth
. Thus, we know a priori what elements of
Thus, (15) and (16) can be repeated iteratively, each stage
these matrices must be stored. Since
is not a function of
,
requiring
computations. For
(15)
each of these matrices has
nonzero elements.3 Since
and (16) form an approximation to (8) requiring a total of
there are
such matrices, the total required storage is
computations and
storage elements. This
, as desired. Furthermore, as we now describe, we can
procedure results in the following approximation of (8):
compute these approximations with an overall computational
load of
.
.
.
The key assumption required for these approximations to
. .
. .
. .
. .
(17)
yield good results is that, for
, the blocks
.
.
are approximately banded, i.e., well approximated by setting
to zero all the elements which do not fall within a small
with
-banded
and
stored in place of
, where
bandwidth of the main diagonal. If
is approximately
is the LU factorization of
. The approximate on-
banded, then from the recursions (6) and (7), it is apparent that
line solution is given by substituting
and
into (9) for
the blocks
and
will generally be approximately
and
, and
into (10) for
.
banded as well. Furthermore, as the following corollary states,
Note that this approximation method extends equally well
if we have a banded approximation to
, we can efﬁciently
to higher-order ﬁlters, with the elements in
are again
compute a banded approximation of its inverse.
ordered columnwise. Each of the blocks in the block LU
Corollary of [5] (See the Appendix for Proof): If
is an
factorization is analogously approximated by a banded matrix.
matrix with bandwidth
, the elements of
which
This extension of the approximation algorithm is demonstrated
lie within the bandwidth
can be computed (exactly) in
for an order (2, 2) ﬁlter in Section V.
operations.
What follows is an example intending both to suggest the
class of ﬁlters for which this approximate implementation will
offer signiﬁcant savings and to justify why in such cases we
3 By O(Np) computations we mean that the number of computations
divided by N p tends to a constant as N ! 1; writing O(qNp) denotes
expect to be able to choose a small value of
independent of
that this constant is a polynomial in of order q.
the size
of the image domain.

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
555
B. Analysis of the Block LU Approximation
the tighter the clustering of signiﬁcant values about the main
Consider a 2-DNC-IIR ﬁlter satisfying (5). To make sub-
diagonal of
, and thus the smaller the bandwidth
needed
sequent analysis simpler, we make a slight deviation from
to approximate
at a desired level of accuracy.
the constant-coefﬁcient model of (4). Namely, assume that
While this example is particularly simple, it does serve
and
for all
and
to demonstrate the intuition and plausibility of our approxi-
. For
, we assume
and
mation. Moreover, a general guideline veriﬁed by extensive
. For the point we wish to make, the values of
numerical simulations, some of which are given in Section V,
the other NNM coefﬁcients are of no consequence. The ﬁrst
is that the approximate factorization applies to order (1,1)
block row of (5) follows as:
ﬁlters for which
and becomes more accurate
(for ﬁxed
) as the ratio
decreases. These
observations are consistent with the preceding example, since
as
. Also note that
.
is a measure of the “degree” of diagonal dominance of the
.
.
.
.
.
.
.
. .
. .
. .
. .
.
.
elements in the blocks
. In fact, for 2-DNC-IIR ﬁlters of
any order, the diagonal dominance of the blocks
seems
to be a useful guideline for determining which ﬁlters can be
implemented by the algorithm in Section IV-A. Determining
(18)
the exact class of ﬁlters which can be approximated by this
algorithm is beyond the scope of this paper.
The factorization of
in (5) begins by factoring
:4
C. A Parallel Approximate Implementation
.
In this section, we brieﬂy discuss a straightforward par-
.
.
.
.
.
.
allelization of the algorithm discussed in Section IV-A. In
.
.
.
.
.
.
. .
particular, a variant of the serial block LU factorization is
.
.
.
.
.
.
.
.
.
.
. .
. .
. .
.
.
.
.
. .
. .
cyclic block reduction [4], which is easily implemented in
parallel. The block LU factorization proceeds by eliminating
columns
sequentially from
to
. However,
it is possible to eliminate columns in the interior of
inde-
(19)
pendently. For example, consider again the block tridiagonal
The next step of the factorization is to compute
. Note
matrix
given in (5). Assume, for simplicity, that
is odd.
that, even though
and
are banded, this operation
If we order the even columns last and the odd columns ﬁrst,
will require
computations. Also note that
(5) takes the form
, where
(21)
. .
.
.
.
.
Because of the coupling implied by the NNM difference
. ..
(20)
equation, the elimination of the odd columns of
for the
.
.
.
.
.
block factorization of (21) can be performed in parallel.
.
.
.
. .
. .
Upon completing this step, the second block equation of (21)
becomes (22), shown at the bottom of the page, where the
The simple form of (20) leads to the expression
superscript
is used to relabel the variables after the -stage
. From this expression, we see that
of the cyclic block reduction. The blocks of
are given by
for
the values of
decay geometrically away
from the main diagonal. Thus, the smaller the value of
,
(equivalently, the greater the diagonal dominance of
),
4 The Toeplitz structure of L11 and U11 is a result of the particular choice
of difference equation, but is not a general property even of ﬁlters described
(23)
by constant-coefﬁcient difference equations.
. .
.
.
.
.
.
. .
. .
.
.
.
.
(22)

556
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
Since
is again block tridiagonal, the odd-even reordering
V. EXAMPLES
can continue recursively, where roughly one half of the
In this section, 2-DNC-IIR ﬁlters are implemented5 using
remaining columns are eliminated in parallel at each stage
the approximation algorithm of Section IV-A. The 2-DNC-
of the algorithm. Rather than the
stages required for the
IIR ﬁlters are assumed to satisfy (2) on
and homogeneous
block LU algorithm (and corresponding on-line solution),
Dirichlet conditions on
. The system functions of
approximately
stages are required for block cyclic
these ﬁlters have the form
reduction, and each of the columns in each stage can be
operated on in parallel. If more coarse-grained parallelism
is required,
can be partitioned into
regions, where
the columns in each region can be eliminated independent of
the other regions [6]. (In the case of cyclic block reduction,
.
.)
.
.
An approximate block cyclic reduction algorithm follows
for any of these parallel structures by noting the strong
similarities between implementing (23) and (6), (7). Namely,
if a processor is allocated for each of the odd columns of
,
.
.
.
.
.
.
.
.
.
the ﬁrst stage of the (parallelized) approximate cyclic block
reduction is
(a)
factor
and compute
.
(24)
on processors
.
.
(b)
compute
on processors
where
is the frequency in the -direction and
is the
compute
frequency in the
-direction. Adding a polynomial
to the numerator of (24) would obviously allow one to
on processors
sharpen the ﬁlter frequency responses, such as by narrowing
(c)
compute
the transition regions or placing zeros in the stopbands [8],
on processors
[16]; however, the purpose of the following examples is
only to demonstrate the utility of the approximation given in
compute
Section IV-A and the viability of 2-DNC-IIR ﬁlters. Thus it
on processors
makes sense to implement only the recursive portion of the
(d)
compute
difference equation.
Some example ﬁlters are given by the following coefﬁcient
on processors
matrices:
compute
on processors
(e)
compute
on processors
Note that Step (e) requires communication between the
processors, since
and
are computed on different,
but “neighboring,” processors. For the ﬁrst stage, the compu-
(25)
tational load for each processor is
, where the asymp-
totic complexity is determined primarily by Step (a). The
computational load for each processor will remain constant
for subsequent stages of the algorithm. Ignoring interprocessor
communication costs after each stage of the recursion, the
total factorization will require
computations.
Since there are
pixels, the per pixel computation time for
the fully parallel implementation is
, which
(26)
decreases with increasing image size. This algorithm extends
to higher-order ﬁlters just as cyclic reduction can be extended
5 All of the Matlab ﬁles used to generate the examples in the section are
to matrices
with larger block bandwidth.
available by anonymous ftp at lids.mit.edu in the directory pub/ssg/outgoing.

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
557
(a)
(b)
(c)
(d)
Fig. 3.
64 2 64 point sampling of Hi(ej! ; ej! ) for four ﬁlters. (a) Low-pass ﬁlter H1 given by J1. (b) High-pass ﬁlter H2 given by J2. (c) Fan ﬁlter
H3 given by J3. (d) Low-pass ﬁlter H4 given by J4. Note that the only frequency with nonzero phase is H3, in which case the magnitude of H3 is plotted.
where
. The coefﬁcients of each
In the following examples, the four 2-DNC-IIR ﬁlters are
ﬁlter are normalized such that
. The frequency
implemented exactly (to within round-off error) using the
responses
of all four difference equations are
nested dissection algorithm [6] and approximately using the
illustrated in Fig. 3. Filters
and
are low-pass ﬁlters,
algorithm of Section IV-A. For analyzing the errors introduced
is the edge-enhancing ﬁlter given in [13], and
corresponds
by the approximate implementation, deﬁne the error signal
to a primitive (low-order) fan ﬁlter [8]. Note that ﬁlters
,
to be
, where
is the exact
, and
have zero phase.
ﬁlter output and
is the output obtained with an
Note that, while ﬁlters
,
, and
have relatively
approximation bandwidth of
. Two of the error measures
nonsharp transition bands, they have many practical applica-
used are
tions. Edge enhancement is one application requiring ﬁlters
with nonsharp transition bands. In fact, all of the edge-
enhancing ﬁlters given in [12], [13], including system
,
and
have frequency responses with nonsharp transition bands.
When enhancing edges, the width of the transition band is
(27)
determined by the model of the smooth edges which are to be
enhanced. This model usually implies a frequency response
where
is a vector containing
for all
,
with a smooth transition band. The smooth transition band
and
is the standard
norm.
can also be used to avoid “over-enhancing” the image, or to
Example 1. The Response of 2-DNC-IIR Filters to Sinusoidal
make the enhanced image less sensitive to noisy data. Another
Inputs: In this example, we consider the response of the
application requiring ﬁlters with smooth frequency responses
2-DNC-IIR example ﬁlters to sinusoidal inputs of the form
is optimal linear estimation. The low-pass ﬁlters given by
and
are of exactly the same order and structure as those
commonly arising in the estimation of Markov random ﬁelds
from noisy measurements [1], [14].
(28)

558
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 4.
The DFT magnitude of the exact responses and the approximation errors for a passband input (left column) and a stopband input (right column) to
the 2-DNC-IIR low-pass ﬁlter given by H1(ej! ; ej! ). The inputs are given by (28) for (k1; k2) = (3, 2) and (k1; k2) = (25, 20). Note the scalings
of the vertical axes. (a) Response to pass-band signal. (b) Response to stop-band signal. (c) Error in pass-band response for = 2. (d) Error in stop-band
response for = 2. (e) Error in pass-band response for = 4. (f) Error in stop-band response for = 4.
The
term serves only to normalize the discrete Fourier
Fig. 4(a) and (b). The Fourier transform
is
transform (DFT) of
. Assume for this example that
shown in lieu of
in order to demonstrate the 9-to-1
64. Consider ﬁrst the low-pass ﬁlter
. For
(3,
(passband to stopband) selection ratio of the low-pass ﬁlter.
2), the sinusoidal input lies in the ﬁlter’s passband, while
The small ridges in the DFT’s of the ﬁlter responses are
for
(25, 20) the input lies in the stopband.
due to the transients introduced by the Dirichlet boundary
The exact responses to these two inputs are illustrated by
conditions. These ridges are barely noticeable for the response

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
559
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 5.
(a) The impulse response of the fan ﬁlter H3. (b) The impulse response of the low-pass ﬁlter H4. (c) The error in the fan ﬁlter impulse response
for = 2. (d) The error in the low-pass ﬁlter impulse response for = 2. (e) The error in the fan ﬁlter impulse response for = 4. (f) The error
in the low-pass ﬁlter impulse response for = 4.
to the stopband input. The effect of these transient signals is
for the response to the passband input.) A similar geometric
negligible near the center of the ﬁlter domain
, and in fact
decrease in the error signal is obtained for larger values of
.
becomes zero as the boundaries recede to inﬁnity (
).
As veriﬁed by extensive numerical simulation, these results
The accuracy of the algorithm of Section IV-A is also
generalize to every sinusoidal signal input to each of the
demonstrated in Fig. 4 for approximation bandwidths
four example ﬁlters. Namely, for every sinusoidal input to
2 and 4. The DFT of the approximation error
is
the 2-DNC-IIR example ﬁlters, the ﬁlter response is given
illustrated for both inputs and both values of
. Note that
by a weighting of the input by the frequency response, plus
the approximation errors are small even when
2, and they
some transients whose effect is limited to the boundaries of
decrease by an order of magnitude when
increases from 2 to
the ﬁltering domain. The effect of these transients on the ﬁlter
4. (In terms of the error metrics,
16 and
22
output near the center of
decreases to zero as the boundary

560
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
Fig. 6.
Approximation errors " versus , for N = 64 and a unit impulse
Fig. 7.
Approximation errors " versus N . The approximation bandwidth is
input applied at (i; j) = (32, 32).
ﬁxed at = 4 and the input is a unit impulse applied at (i; j) = (N=2; N=2).
recedes to inﬁnity (
). Also, for each sinusoidal input,
approximation errors to decrease geometrically with increasing
the approximate response is accurate even for small values
approximation bandwidth.]
of
. More importantly, the approximation errors decrease
Example 3. The Independence of the Approximation Accu-
geometrically in magnitude with
. This geometric decrease
racy upon
: In Section IV-A, the computational and storage
is further demonstrated in the following examples.
loads of the approximation algorithm were shown to be
Example 2. The Impulse Response of 2-DNC-IIR Filters:
and
, respectively, meaning that the per
Because LSI ﬁlters are completely characterized by their
pixel computational and storage loads converge to polynomials
impulse responses, we now consider the impulse responses
in
of order two and one, respectively. However, if the
of both the exact and the approximate implementations of the
per pixel computational and storage loads are to be truly
example 2-DNC-IIR ﬁlters. The impulse response
is the
constant, the approximation bandwidth needed for a desired
ﬁlter output in response to a unit impulse applied at the center
approximation accuracy must not be an increasing function
of
. For
64, this input is
. The
of
. For a unit impulse applied at the center of
, Fig. 7
impulse response of both the fan ﬁlter
and the order (2,2)
shows that the approximation error metric
remains constant
low-pass ﬁlter
are plotted at the top of Fig. 5. Because
over a wide range of
for a ﬁxed value of
4. Identical
the impulse responses are essentially zero outside the region
results are obtained for other values of
. The independence
{
48}, the responses are not plotted outside
of
upon
for a desired solution accuracy is consistent
this region. The error signals plotted in Fig. 5 correspond to
with the analysis of Section IV-B, where the rate of decrease
the difference between the true impulse responses and those
with distance from the diagonal of the elements in the blocks
of the approximate implementations for
2 and
4.
and
was shown to be independent of
. Thus, the
Like the approximate responses to the sinusoidal inputs given
approximation algorithm has constant per pixel computational
in Example 1, the approximate impulse responses are accurate
and storage loads.
even for
2, and the errors the decrease geometrically
Example 4. Edge Enhancement of a Square Pulse: For this
with increasing
. Note especially the accuracy of the
-
example, we examine the response of the edge-enhancing ﬁlter
banded approximation for the order (2,2) low-pass ﬁlter. While
to a square pulse. The square pulse input and the response
the frequency response of this system is very sensitive to
of the edge-enhancing ﬁlter are illustrated at the top of Fig. 8.
the values of the coefﬁcients stored in
, the approximate
Note that the response to the square input is analogous to
implementation is quite accurate even for small
. For
4,
the step response of the ﬁlter. The approximate responses for
the maximum value of
is on the order of
.
2 and 4, and the corresponding approximation errors,
To illustrate that the errors decrease geometrically as
are also illustrated in Fig. 8. Again, the approximation errors
increases beyond four,
is plotted versus
in Fig. 6 for
are small and visually unrecognizable even for these small
ranging from two to ten. Fig. 6 also demonstrates that the
approximation bandwidths, and the errors decrease by an order
geometric decrease in the error with increasing
applies to
of magnitude as
increases from 2 to 4.
the other three example ﬁlters. [Recall that in Section IV-
B we argued that the approximation errors will be small
VI. CONCLUSION
when
. For these order (1,1) ﬁlters, the
In this paper we describe an approach to the efﬁcient imple-
elements in the blocks
and
decreased geometrically
mentations of 2-DNC-IIR ﬁlters. In addition to efﬁciency, we
in magnitude with distance from the diagonal. Thus, for
-
were also motivated to consider ﬁlters speciﬁed by boundary
banded approximations of these blocks, one would expect the
rather than initial conditions, as the former are frequently the

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
561
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 8.
(a) Square pulse input signal. (b) Exact response of edge-enhancing ﬁlter. (c)–(d) Filter output and error signals for = 2. (e)–(f) Filter output
and error signal for = 4. Note the scalings of the vertical axes.
natural choice and are required, for example, if zero-phase
application of concepts from the direct solution of PDE’s to the
ﬁltering is desired. Indeed, some methodologies now exist for
calculation of the solution of a 2-D difference equation; and
designing 2-D difference equations to meet desired frequency-
b) the development of new approximations, motivated by and
selective speciﬁcations, and we demonstrated that imposing
appropriate for ﬁltering applications, that reduce complexity
boundary conditions upon these difference equations can lead
to desired levels. In particular, the algorithms resulting from
to the desired frequency selectivity.
our procedure have constant computational complexity per
The approach we developed for efﬁciently implementing 2-
pixel and, if implemented in maximally parallel form, have
DNC-IIR ﬁlters involves a combination of two things: a) the
total computation time per pixel that decreases as image size

562
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—II: ANALOG AND DIGITAL SIGNAL PROCESSING, VOL. 44, NO. 7, JULY 1997
increases. In particular, our approximation is based on the
For certain matrices, such as those discussed in Section IV-A,
columnwise ordering of data points and the block LU fac-
the elements of
which fall within the bandwidth
can be
torization of the linear system that results from this ordering.
seen as a reasonable approximation to
.
While exact factorization is still complex computationally, the
To implement this algorithm, note that
must be
observation that each successive block of computation could
the ﬁrst element computed of
. Also, to compute any
be viewed as a 1-D ﬁltering operation along a column of the
element
, all elements
such that
,
image led to the idea of a reduced-order approximation of
, and
must already have been computed.
each of these 1-D columnwise ﬁlters. In matrix terms, this
With these restrictions placed on the recursion, the number
corresponds to a banded approximation to each of the blocks
of computations to compute
within a bandwidth
is
in the block LU factorization, with bandwidth (and 1-D ﬁlter
bounded above closely by
, where the ﬁrst
order)
. The resulting algorithm was shown both to achieve
term represents the computations for the diagonal elements
the computational levels mentioned previously and to yield
only. The computational complexity of the algorithm is thus
excellent results using small values of
for a number of
.
low-order frequency-selective ﬁlters.
The approach that we have described is, in principle, appli-
REFERENCES
cable to a broad range of ﬁltering problems, e.g., higher-order
[1] R. Chellappa and S. Chatterjee, “Classiﬁcation of textures using gaussian
or nonconstant-coefﬁcient difference equations and irregularly
Markov random ﬁelds,” IEEE Trans. Acoust., Speech, Signal Processing,
shaped regions. Indeed, the success we have demonstrated here
vol. ASSP-33, pp. 959–963, 1985.
together with the guidelines we have described for situations
[2] H. Derin and P. A. Kelly, “Discrete-index Markov-type random pro-
cesses,” Proc. IEEE, vol. 77, Oct. 1989.
in which our approximation should work well provide ample
[3] D. E. Dudgeon and R. M. Mersereau, Multidimensional Digital Signal
motivation for the application of this methodology and for a
Processing.
Englewood Cliffs, NJ: Prentice-Hall, 1984.
[4] I. S. Duff, A. M. Erisman, and J. K. Reid, Direct Methods for Sparse
complete investigation of general conditions on the difference
Matrices.
Oxford, England: Oxford Univ. Press, 1987.
equation coefﬁcients under which our approach is guaranteed
[5] A. M. Erisman and W. F. Tinney, “On computing certain elements of
to provide accurate answers for small values of
. However,
the inverse of a sparse matrix,” Commun. ACM, vol. ACM-18, no. 3,
1975.
the implementations of 2-DNC-IIR ﬁlters are certainly not
[6] A. George and J. W. Liu, Computer Solution of Large and Sparse
limited to direct methods. Many of the efﬁcient iterative
Positive Deﬁnite Systems.
Englewood Cliffs, NJ: Prentice-Hall, 1981.
algorithms developed for solving PDE’s will likely provide
[7] G. H. Golub and C. F. Van Loan, Matrix Computations.
Baltimore,
MD: John Hopkins Univ. Press, 1990.
efﬁcient solutions to 2-DNC-IIR ﬁlters, and perhaps to a
[8] Q. Gu and M. N. S. Swamy, “On the design of a broad class of 2D
different class of ﬁlters than can be implemented by the
recursive digitial ﬁlters with fan, diamond, and elliptically-symmetric
responses,” IEEE Trans. Circuits Syst. II, vol. 41, pp. 603–614, 1994.
algorithm of Section IV-A.
[9] B. K. Horn and B. G. Schunck, “Determining optical ﬂow,” Artif. Intell.,
vol. 17, pp. 185–203, 1981.
APPENDIX
[10] Mathworks Inc., Matlab Image Processing Toolbox, 4.1 ed.
Natick,
MA, June 1993.
COMPUTING CERTAIN ELEMENTS OF
[11] B. C. Levy, M. B. Adams, and A. S. Willsky, “Solution and linear
THE INVERSE OF A BANDED MATRIX
estimation of 2-D nearest-neighbor models,” Proc. IEEE, vol. 78,
1990.
A special application of the results in [5] allows for the
[12] J. S. Lim, Two-Dimensional Signal and Image Processing.
Englewood
efﬁcient computation of certain elements of the inverse of
Cliffs, NJ: Prentice-Hall, 1990.
[13] W. Lu and A. Antoniou, Two-dimensional Digital Fitlers.
New York:
a banded matrix. Namely, if
is an
dimensional
Marcel Dekker, 1992.
matrix with bandwidth
and
, then the elements
[14] M. R. Luettgen, W. C. Karl, A. S. Willsky, and R. R. Tenney,
“Multiscale representations of Markov random ﬁelds,” IEEE Trans.
of
which lie within a bandwidth
of the diagonal can
Signal Processing, vol. 41, p. 3377, Dec. 1993.
be computed in
computations. The results of [5]
[15] A. R. Mitchell and D. F. Grifﬁths, The Finite Difference Method in
are based on the following observation: given
Partial Differential Equations.
New York: Wiley, 1980.
[16] N. A. Pendergrass, S. K. Mitra, and E. I. Jury, “Spectral transformations
(where
and
are unit lower-triangular and upper-triangular,
for two-dimensional digital ﬁlters,” IEEE Trans. Circuits Syst., vol.
respectively), then
CAS-23, pp. 26–35, 1976.
(29)
(30)
From these relations,
for
is given by (31),
Michael M. Daniel (S’93) received the B.S. degree
which does not depend upon computing any elements in
from the University of California, Berkeley, in
1990 and the Ph.D. degree from the Massachusetts
and
.
Institute of Technology (MIT), Cambridge, in
1997.
He has been with Schlumberger–Doll Research,
Ridgeﬁeld, CT, and the Cambex Corporation,
Waltham, MA, and has been a teaching assistant for
the Circuits, Signals, and Systems and the Discrete-
(31)
Time Signal Processing courses at MIT. He recently
coauthored the undergraduate-level book Computer
Explorations for Signals and Systems Using MATLAB. His research interests
include signal and image processing, stochastic processes and uncertainty
analysis, and geophysical inverse problems.
Dr. Daniel is a member of SIAM, Sigma Xi, Tau Beta Pi, and Eta Kappa Nu.

DANIEL AND WILLSKY: EFFICIENT IMPLEMENTATIONS OF 2-DNC-IIR FILTERS
563
Alan
S.
Willsky
(S’70–M’73–SM’82–F’86)
received the S.B. and Ph.D. degrees from the
Massachusetts
Institute
of
Technology
(MIT),
Cambridge, in 1969 and 1973, respectively.
He joined the MIT faculty in 1973 and is
currently a Professor of electrical engineering.
From 1974 to 1981, he served as Assistant Director
of the Laboratory for Information and Decision
Systems, MIT. He is also a founder and member
of the Board of Directors for Alphatech, Inc.,
Burlington, MA. He is the author of the research
monograph Digital Signal Processing and Control and Estimation Theory
and is coauthor of the undergraduate text Signals and Systems. He has held
visiting positions at Imperial College, London, U.K., l’Universit´e de Paris-Sud,
and INRIA, France. His present research interests are in problems involving
multidimensional and multiresolution signal processing and imaging, discrete-
event systems, and the asymptotic analysis of control and estimation systems.
Dr. Willsky was program chair for the 17th IEEE Conference on Decision
and Control, and has been an Associate Editor for several journals. He has
served as a member of the Board of Governors and was Vice President for
Technical Affairs of the IEEE Control Systems Society. He was program
chairman for the 1981 Bilateral Seminar on Control Systems held in the
People’s Republic of China, and was a Special Guest Editor in 1992 for
the IEEE TRANSACTIONS ON INFORMATION THEORY. In 1988, he was made
a Distinguished Member of the IEEE Control Systems Society. He has
given several plenary lectures at major scientiﬁc meetings, including the
1992 Inaugural Workshop for the National Center for Robust and Adaptive
Systems, Canberra, Australia. In 1975, he received the Donald P. Eckman
Award, in 1979 the Alfred Nobel Prize, and in 1980 the Broder S. Thompson
Memorial Prize recognizing a paper excerpted from his monograph.