CHAPTER 2

Analysis of Two-Dimensional Signals

and Systems

M a n y physical phenomena are found experimentally to share the basic property that
their response to several stimuli acting simultaneously is identically equal to the sum of
the responses that each component stimulus would produce individually. Such phenom-
ena are called lineal; and the property they share is called linearity. Electrical networks
composed of resistors, capacitors, and inductors are usually linear over a wide range of
inputs. In addition, as we shall soon see, the wave equation describing the propagation
of light through most media leads us naturally to regard optical imaging operations as
linear mappings of "object" light distributions into "image" light distributions.
The single property of linearity leads to a vast simplification in the mathematical
description of such phenomena and represents the foundation of a mathematical struc-
ture which we shall refer to here as linear systems theory. The great advantage afforded
by linearity is the ability to express the response (be it voltage, current, light amplitude,
or light intensity) to a complicated stimulus in terms of the responses to certain "elemen-
tary" stimuli. Thus if a stimulus is decomposed into a linear combination of elementary
stimuli, each of which produces a known response of convenient form, then by virtue
of linearity, the total response can be found as a corresponding linear combination of
the responses to the elementary stimuli.
In this chapter we review some of the mathematical tools that are useful in describ-
ing linear phenomena, and discuss some of the mathematical decompositions that are
often employed in their analysis. Throughout the later chapters we shall be concerned
with stimuli (system inputs) and responses (system outputs) that may be either of two
different physical quantities. If the illumination used in an optical system exhibits a
property called spatial coherence, then we shall find that it is appropriate to describe
the light as a spatial distribution of complex-valued field amplitude. When the illumi-
nation is totally lacking in spatial coherence, it is appropriate to describe the light as a
spatial distribution of real-valued intensity. Attention will be focused here on the anal-
ysis of linear systems with complex-valued inputs; the results for real-valued inputs are
thus included as special cases of the theory.
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 5

2.1
FOURIER ANALYSIS IN TWO DIMENSIONS

A mathematical tool of great utility in the analysis of both linear and nonlinear phenom-
ena is Fourier analysis. This tool is widely used in the study of electrical networks and
communication systems; it is assumed that the reader has encountered Fourier theory
previously, and therefore that he or she is familiar with the analysis of functions of one
independent variable (e.g. time). For a review of the fundamental mathematical con-
cepts, see the books by Papoulis [226], Bracewell [32], and Gray and Goodman [131].
A particularly relevant treatment is by Bracewell [33]. Our purpose here is limited to
extending the reader's familiarity to the analysis of functions of two independent vari-
ables. No attempt at great mathematical rigor will be made, but rather, an operational
approach, characteristic of most engineering treatments of the subject, will be adopted.

2.1.1 Definition and Existence Conditions

The Fourier transform (alternatively the Fourier spectrum or frequency spectrum) of

a (in general, complex-valued) function g of two independent variables x and y will be
represented here by F{g} and is defined b y1

The transform so defined is itself a complex-valued function of two independent vari-

ables fx and fr, which we generally refer to as frequencies. Similarly, the inverse
Fourier transform of a function G(fx, fy) will be represented by F 1 { G } and is de-
fined as

Note that as mathematical operations the transform and inverse transform are very sim-
ilar, differing only in the sign of the exponent appearing in the integrand. The inverse
Fourier transform is sometimes referred to as the Fourier integral representation of a
function g(x, y).
Before discussing the properties of the Fourier transform and its inverse, we must
first decide when (2-1) and (2-2) are in fact meaningful. For certain functions, these
integrals may not exist in the usual mathematical sense, and therefore this discussion
would be incomplete without at least a brief mention of "existence conditions". While
a variety of sets of suficient conditions for the existence of (2-1) are possible, perhaps
the most common set is the following:

'When a single limit of integration appears above or below a double integral, then that limit applies to both
integrations.
6 Introduction to Fourier Optics

1. g must be absolutely integrable over the infinite ( x , y )plane.

2. g must have only a finite number of discontinuities and a finite number of maxima
and minima in any finite rectangle.
3. g must have no infinite discontinuities.
In general, any one of these conditions can be weakened at the price of strengthen-
ing one or both of the companion conditions, but such considerations lead us rather far
afield from our purposes here.
As Bracewell [32] has pointed out, "physical possibility is a valid sufficient condi-
tion for the existence of a transform." However, it is often convenient in the analysis of
systems to represent true physical waveforms by idealized mathematical functions, and
for such functions one or more of the above existence conditions may be violated. For
example, it is common to represent a strong, narrow time pulse by the so-called Dirac
delta function2 often represented by
6(t) = lim Nexp(-N 2 ~t 2), (2-3)
N+m

where the limit operation provides a convenient mental construct but is not meant to be
taken literally. See Appendix A for more details. Similarly, an idealized point source of
light is often represented by the two-dimensional equivalent,
6 ( x , y) = ~lim
- + N~
m ~ X ~ [ - N ~ T
y2)].
( ~+
Such "functions", being infinite at the origin and zero elsewhere, have an infinite dis-
continuity and therefore fail to satisfy existence condition 3. Qther important examples
are readily found; for example, the functions

both fail to satisfy existence condition 1.

If the majority of functions of interest are to be included within the framework of
Fourier analysis, some generalization of the definition (2-1) is required. Fortunately, it
is often possible to find a meaningful transform of functions that do not strictly satisfy
the existence conditions, provided those functions can be defined as the limit of a se-
quence of functions that are transformable. By transforming each member function of
the defining sequence, a corresponding sequence of transforms is generated, and we
call the limit of this new sequence the generalized Fourier transform of the original
function. Generalized transforms can be manipulated in the same manner as conven-
tional transforms, and the distinction between the two cases can generally be ignored,
it being understood that when a function fails to satisfy the existence conditions and
yet is said to have a transform, then the generalized transform is actually meant. For a
more detailed discussion of this generalization of Fourier analysis the reader is referred
to the book by Lighthill [194].
To illustrate the calculation of a generalized transform, consider the Dirac delta
function, which has been seen to violate existence condition 3. Note that each member
function of the defining sequence (2-4) does satisfy the existence requirements and that
each, in fact, has a Fourier transform given by (see Table 2.1)

2Fora more detailed discussion of the delta function, including definitions, see Appendix A.
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 7

F{N2 exp[-N2.sr(x2 [
+ y2)]) = exp - -'f$ f3].

Accordingly the generalized transform of S(x, y) is found to be

Note that the spectrum of a delta function extends uniformly over the entire frequency
domain.
For other examples of generalized transforms, see Table 2.1.

2.1.2 The Fourier Transform as a Decomposition

As mentioned previously, when dealing with linear systems it is often useful to decom-
pose a complicated input into a number of more simple inputs, to calculate the response
of the system to each of these "elementary" functions, and to superimpose the individ-
ual responses to find the total response. Fourier analysis provides the basic means of
performing such a decomposition. Consider the familiar inverse transform relationship

expressing the time function g in terms of its frequency spectrum. We may regard this
expression as a decomposition of the function g(t) into a linear combination (in this
case an integral) of elementary functions, each with a specific form exp(j2.srf t). From
this it is clear that the complex number G(f ) is simply a weighting factor that must
be applied to the elementary function of frequency f in order to synthesize the desired
g(t).
In a similar fashion, we may regard the two-dimensional Fourier transform as a de-
composition of a function g(x, y) into a linear combination of elementary functions of
+
the form e x p [ j 2 ~fxx
( fry)]. Such functions have a number of interesting properties.
Note that for any particular frequency pair ( fx, fy) the corresponding elementary func-
tion has a phase that is zero or an integer multiple of 2.sr radians along lines described
by the equation

where n is an integer. Thus, as indicated in Fig. 2.1, this elementary function may be
regarded as being "directed" in the (x,y) plane at an angle 8 (with respect to the x axis)
given by

In addition, the spatial period (i.e. the distance between zero-phase lines) is given by
8 Introduction to Fourier Optics

\ \ FIGURE
- - - - 2.1
-.-

Lines of zero phase for the fur

exp[j2.rr(fxn + fry)].

In conclusion, then, we may again regard the inverse Fourier transform as providing a
means for decomposing mathematical functions. The Fourier spectrum G of a function g
is simply a description of the weighting factors that must be applied to each elementary
function in order to synthesize the desired g. The real advantage obtained from using
this decomposition will not be fully evident until our later discussion of invariant linear
systems.

2.1.3 Fourier Transform Theorems

The basic definition (2-1) of the Fourier transform leads to a rich mathematical
structure associated with the transform operation. We now consider a few of the
basic mathematical properties of the transform, properties that will find wide use in
later material. These properties are presented as mathematical theorems, followed
by brief statements of their physical significance. Since these theorems are direct
extensions of the analogous one-dimensional statements, the proofs are deferred to
Appendix A.
1. Linearity theorem. F { u g + ph} = uF{g) + PF{h}; that is, the transform of a
weighted sum of two (or more) functions is simply the identically weighted sum of
their individual transforms.
2. Similarity theorem. If F{g(x, y)) = G(fx, fv), then

F{g(ax, by)) =

that is, a "stretch" of the coordinates in the space domain (x, y) results in a contrac-
tion of the coordinates in the frequency domain ( fx, fy), plus a change in the overall
amplitude of the spectrum.
3. Shift theorem. If F{g(x, y)} = G(fx, fy), then

that is, translation in the space domain introduces a linear phase shift in the fre-
quency domain.
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 9

4. Rayleigh's theorem (Parseval's theorem). If F{g(x, y)} = G(fx, fy), then

The integral on the left-hand side of this theorem can be interpreted as the energy
contained in the waveform g(x, y). This in turn leads us to the idea that the quantity
IG(fx, fy)I2 can be interpreted as an energy density in the frequency domain.
5. Convolution theorem. If 3{g(x, y)) = G( fx, fy) and F{h(x, y)} = H ( fx, fy), then

The convolution of two functions in the space domain (an operation that will be
found to arise frequently in the theory of linear systems) is entirely equivalent to
the more simple operation of multiplying their individual transforms and inverse
transforming.
6. Autocorrelation theorem. If 3{g(x, y)} = G(fx, fy), then

Similarly,

This theorem may be regarded as a special case of the convolution theorem in which
we convolve g(x, y) with g*(-x, -y).
7. Fourier integral theorem. At each point of continuity of g,

At each point of discontinuity of g, the two successive transforms yield the angular
average of the values of g in a small neighborhood of that point. That is, the suc-
cessive transformation and inverse transformation of a function yields that function
again, except at points of discontinuity.

The above transform theorems are of far more than just theoretical interest. They
will be used frequently, since they provide the basic tools for the manipulation of Fourier
transforms and can save enormous amounts of work in the solution of Fourier analysis
problems.
10 Introduction to Fourier Optics

2.1.4 Separable Functions

A function of two independent variables is called separable with respect to a specific

coordinate system if it can be written as a product of two functions, each of which
depends on only one of the independent variables. Thus the function g is separable in
rectangular coordinates (x, y) if

while it is separable in polar coordinates (r, 8) if

Separable functions are often more convenient to deal with than more general
functions, for separability often allows complicated two-dimensional manipulations to
be reduced to more simple one-dimensional manipulations. For example, a function
separable in rectangular coordinates has the particularly simple property that its two-
dimensional Fourier transform can be found as a product of one-dimensional Fourier
transforms, as evidenced by the following relation:

Thus the transform of g is itself separable into a product of two factors, one a function
of fx only and the second a function of f y only, and the process of two-dimensional
transformation simplifies to a succession of more familiar one-dimensional manipula-
tions.
Functions separable in polar coordinates are not so easily handled as those sep-
arable in rectangular coordinates, but it is still generally possible to demonstrate that
two-dimensional manipulations can be performed by a series of one-dimensional ma-
nipulations. For example, the reader is asked to verify in the problems that the Fourier
transform of a general function separable in polar coordinates can be expressed as an
infinite sum of weighted Hankel transforms

where

and X k { }is the Hankel transform operator of order k, defined by

 CHAFTER 2 Analysis of Two-Dimensional Signals and Systems 11

Rk{g~(r))= 2~
I,"r g ~ ( rJk(2m-p)
) dr.

Here the function Jk is the kth-order Bessel function of the first kind.

2.1.5 Functions with Circular Symmetry: Fourier-Bessel lkansforms

Perhaps the simplest class of functions separable in polar coordinates is composed of

those possessing circular symmetry. The function g is said to be circularly symmetric
if it can be written as a function of r alone, that is,

Such functions play an important role in the problems of interest here, since most optical
systems have precisely this type of symmetry. We accordingly devote special attention
to the problem of Fourier transforming a circularly symmetric function.
The Fourier transform of g in a system of rectangular coordinates is, of course,
given by
m

-w

To fully exploit the circular symmetry of g, we make a transformation to polar coordi-

nates in both the (x, y) and the ( fx, fy) planes as follows:

r= Jm x= rcos~

0 =
(3
arctan -

For the present we write the transform as a function of both radius and angle,3

Applying the coordinate transformations (2-26) to Eq. (2-25), the Fourier transform
of g can be written

G0(p 4 ) = l ' 1;
o
2
d0 dr rga(r) e x p l j2nrp(cos 0 cos 4 + sin 0 sin +)I (2-28)

or equivalently,

3Note the subscript in Go is added simply because the functional form of the expression for the transform
in polar coordinates is in general different than the functional form for the same transform in rectangular
coordinates.
12 Introduction to Fourier Optics

Finally, we use the Bessel function identity

where Jo is a Bessel function of the first kind, zero order, to simplify the expression for
the transform. Substituting (2-30) in (2-29), the dependence of the transform on angle
4 is seen to disappear, leaving Go as the following function of radius p,

Thus the Fourier transform of a circularly symmetric function is itself circularly

symmetric and can be found by performing the one-dimensional manipulation of (2-
3 1). This particular form of the Fourier transform occurs frequently enough to warrant
a special designation; it is accordingly referred to as the Fourier-Bessel transfomz, or
alternatively as the Hankel transform of zero order (cf. Eq. (2-23)). For brevity, we
adopt the former terminology.
By means of arguments identical with those used above, the inverse Fourier trans-
form of a circularly symmetric spectrum Go@) can be expressed as

Thus for circularly symmetric functions there is no difference between the transform
and the inverse-transform operations.
Using the notation a { ) to represent the Fourier-Bessel transform operation, it fol-
lows directly from the Fourier integral theorem that

at each value of r where gR(r) is continuous. In addition, the similarity theorem can be
straightforwardly applied (see Prob. 2-6c) to show that

B{g~(ar))= a+o
(s)
When using the expression (2-31) for the Fourier-Bessel transform, the reader should re-
member that it is no more than a special case of the two-dimensional Fourier transform,
and therefore any familiar property of the Fourier transform has an entirely equivalent
counterpart in the terminology of Fourier-Bessel transforms.

2.1.6 Some Frequently Used Functions and Some Useful Fourier Transform
Pairs

A number of mathematical functions will find such extensive use in later material that
considerable time and effort can be saved by assigning them special notations of their
own. Accordingly, we adopt the following definitions of some frequently used func-
tions:
 CHAPTER 2. Analysis of Two-Dimensional Signals and Systems 13

FIGURE 2.2
Special functions.

Rectangle function
0 otherwise
sin(~x)
Sinc function sinc(x) =
TX

1 x > o
Signum function
- 1 x < o

Triangle function A(x) =

otherwise

Comb function cornb(x) = 6(x - n)

14 Introduction to Fourier Optics

T A B L E 2.1
Transform pairs for some functions separable in
rectangular coordinates.
Function Transform

e x - ( a +2
2
y 2 ?- exp
lab1
[ (2+ 2 )]
-n - -

I
6(ax,by) -
lab1
exp[j.rr(ax+ by)] 6(fx - al2, f~ - bl2)
ab
---1 1
sgn(ax) sgn(by)
lab1 jnfx jnfy
1
comb(ax) comb(by) -comb(fxla) comb(fylb)
lab1

exp[jn(a 2 x 2+ b2y2)] -

1 2 2
exp[-(alxl + blyl)l lab1 1 + ( 2 nf ~ l a )1 ~+ (27rfylb)2

Circle function circ( -)/, = -/, = 1

( 0 otherwise.
The first five of these functions, depicted in Fig. 2.2, are all functions of only one in-
dependent variable; however, a variety of separable functions can be formed in two
dimensions by means of products of these functions. The circle function is, of course,
unique to the case of two-dimensional variables; see Fig. 2.3 for an illustration of its
structure.
We conclude our discussion of Fourier analysis by presenting some specific two-
dimensional transform pairs. Table 2.1 lists a number of transforms of functions sep-
arable in rectangular coordinates. For the convenience of the reader, the functions are
presented with arbitrary scaling constants. Since the transforms of such functions can
be found directly from products of familiar one-dimensional transforms, the proofs of
these relations are left to the reader (cf. Prob. 2-2).
On the other hand, with a few exceptions (e.g. exp[-v(x2 + y2)], which is both
separable in rectangular coordinates and circularly symmetric), transforms of most
circularly symmetric functions cannot be found simply from a knowledge of one-
dimensional transforms. The most frequently encountered function with circular sym-
metry is:

circ(r) = 3
"'
r= 1 .
( 0 otherwise
 (a)
FIGURE 2.3
(a) The circle function and (b) its transform.
16 Introduction to Fourier Optics

Accordingly, some effort is now devoted to finding the transform of this function. Using
the Fourier-Bessel transform expression (2-3 I), the transform of the circle function can
be written
1
B{circ(r)) = 2 7 ~ r J o ( 2 ~ r pdr.
)
I0
Using a change of variables r' = 2 ~ r and
p the identity

we rewrite the transform as

where J 1 is a Bessel function of the first kind, order 1 . Figure 2.3 illustrates the circle
function and its transform. Note that the transform is circularly symmetric, as expected,
and consists of a central lobe and a series of concentric rings of diminishing amplitude.
Its value at the origin is T . AS a matter of curiosity we note that the zeros of this trans-
form are not equally spaced in radius. A convenient normalized version of this function,
with value unity at the origin, is 2 F . This particular function is called the "besinc"
function, or the "jinc" function.
For a number of additional Fourier-Bessel transform pairs, the reader is referred to
the problems (see Prob. 2-6).

2.2
LOCAL SPATIAL FREQUENCY
AND SPACE-FREQUENCY LOCALIZATION

Each Fourier component of a function is a complex exponential of a unique spatial fre-

quency. As such, every frequency component extends over the entire (x, y) domain.
Therefore it is not possible to associate a spatial location with a particular spatial fre-
quency. Nonetheless, we know that in practice certain portions of an image could con-
tain parallel grid lines at a certain fixed spacing, and we are tempted to say that the
particular frequency or frequencies represented by these grid lines are localized to cer-
tain spatial regions of the image. In this section we introduce the idea of local spatial
frequencies and their relation to Fourier components.
For the purpose of this discussion, we consider the general case of complex-valued
functions, which we will later see represent the amplitude and phase distributions of
monochromatic optical waves. For now, they are just complex functions. Any such func-
tion can be represented in the form

where a(x, y) is a real and nonnegative amplitude distribution, while 4(x, y) is a real
phase distribution. For this discussion we assume that the amplitude distribution a(x, y)
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 17

is a slowly varying function of (x, y), so that we can concentrate on the behavior of the
phase function +(x, y).
We define the local spatialfrequency of the function g as a frequency pair (fix, fiy)
given by
1 d 1 d
fix = Gz4(x' Y) fiy = T;; &4(* Y).
In addition, both fix and fiYare defined to be zero in regions where the function g(x, y)
vanishes.
Consider the result of applying these definitions to the particular complex function

representing a simple linear-phase exponential of frequencies (fx, fy). We obtain

Thus we see that for the case of a single Fourier component, the local frequencies do
indeed reduce to the frequencies of that component, and those frequencies are constant
over the entire (x, y) plane.
Next consider a space-limited version of a quadratic-phase exponential function:
which we call a "finite chirp" f ~ n c t i o n , ~

g(x, y) = (x2 + y2)l rect

exp[j7.r~ (&) rect(&).
Performing the differentiations called for by the definitions of local frequencies, we find
that they can be expressed as

We see that in this case the local spatial frequencies do depend on location in the
(x, y) plane; within a rectangle of dimensions 2Lx X 2Ly, f i x varies linearly with the
x-coordinate while fiYvaries linearly with the y-coordinate. Thus for this function (and
for most others) there is a dependence of local spatial frequency on position in the (x, y)
plane.6
Since the local spatial frequencies are bounded to covering a rectangle of dimen-
sions 2Lx X 2Ly, it would be tempting to conclude that the Fourier spectrum of g(x, y)
is also limited to the same rectangular region. In fact this is approximately true, but not
exactly so. The Fourier transform of this function is given by the expression

4Fora tutorial discussion of the importanceof quadratic-phase functions in various fields of optics, see [229].
5Thename "chirp function", without the finite length qualifier, will be used for the infinite-length quadratic
phase exponential, exp[j.rrp (xZ+ y 2 ) ] .
6From the definition (2-37) the dimensions of fix and fiY are both cycles per meter; in spite of what might
appear to be a contrary implication of Eq. (2-39). The dimensions of P are meters-2.
18 Introduction to Fourier Optics

This expression is separable in rectangular coordinates, so it suffices to find the one-

dimensional spectrum
Lx
Gx(fx)= ,jrPxZej2rfxx dX.

Completing the square in the exponent and making a change of variables of integration
from x to t = ,/@ (x - $)yields

This integral can be expressed in terms of tabulated functions, the Fresnel integrals,
which are defined by

The spectrum Gxcan then be expressed as

The expression for G y is of course identical, except the Y subscript replaces the X
subscript. Figure 2.4 shows a plot of IGx(fx)l vs. fx for the particular case of Lx = 10
and p = 1. As can be seen, the spectrum is almost flat over the region (-Lx,Lx)and

0.8.

0.6.

0.4 -

0.2.
FIGURE 2.4
-15 -10 -5 5 10 15 The spectrum of the finite chirp function,
fx Lx=1O,/3=1.
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 19

almost zero outside that region. We conclude that local spatial frequency has provided
a good (but not exact) indication of where the significant values of the Fourier spectrum
will occur. However, local spatial frequencies are not the same entity as the frequency
components of the Fourier spectrum. Examples can be found for which the local spatial
frequency distribution and the Fourier spectrum are not in as good agreement as found in
the above example. Good agreement can be expected only when the variations of +(x, y)
are sufficiently "slow" in the (x, y) plane to allow +(x, y) to be well approximated by
only three terms of its Taylor series expansion about any point (x, y), i.e. a constant
term and two first-partial-derivative terms.
Local spatial frequencies are of special physical significance in optics. When the
local spatial frequencies of the complex amplitude of a coherent optical wavefront are
found, they correspond to the ray directions of the geometrical optics description of that
wavefront. However, we are getting ahead of ourselves; we will return to this idea in
later chapters and particularly in Appendix B.

2.3
LINEAR SYSTEMS

For the purposes of discussion here, we seek to define the word system in a way suf-
ficiently general to include both the familiar case of electrical networks and the less-
familiar case of optical imaging systems. Accordingly, a system is defined to be a map-
ping of a set of input functions into a set of output functions. For the case of electrical
networks, the inputs and outputs are real-valued functions (voltages or currents) of a
one-dimensional independent variable (time); for the case of imaging systems, the in-
puts and outputs can be real-valued functions (intensity) or complex-valued functions
(field amplitude) of a two-dimensional independent variable (space). As mentioned pre-
viously, the question of whether intensity or field amplitude should be considered the
relevant quantity will be treated at a later time.
If attention is restricted to deterministic (nonrandom) systems, then a specified in-
put must map to a unique output. It is not necessary, however, that each output corre-
spond to a unique input, for as we shall see, a variety of input functions can produce no
output. Thus we restrict attention at the outset to systems characterized by many-to-one
mappings.
A convenient representation of a system is a mathematical operator, S O , which we
imagine to operate on input functions to produce output functions. Thus if the function
gl (xl, yl) represents the input to a system, and g2(x2,y2) represents the corresponding
output, then by the definition of S {I, the two functions are related through
g 2 ( ~ 2~, 2 =) S{gl (XI,YI )I. (2-41)
Without specifying more detailed properties of the operator SO, it is difficult to state
more specific properties of the general system than those expressed by Eq. (2-41). In
the material that follows, we shall be concerned primarily, though not exclusively, with
a restricted class of systems that are said to be lineal: The assumption of linearity will
be found to yield simple and physically meaningful representations of such systems; it
will also allow useful relations between inputs and outputs to be developed.
20 Introduction to Fourier Optics

2.3.1 Linearity and the Superposition Integral

A system is said to be linear if the following superposition property is obeyed for all
input functions p and q and all complex constants a and b:

As mentioned previously, the great advantage afforded by linearity is the ability

to express the response of a system to an arbitrary input in terms of the responses to
certain "elementary" functions into which the input has been decomposed. It is most
important, then, to find a simple and convenient means of decomposing the input. Such
a decomposition is offered by the so-called sifting property of the 6 function (cf. Section
1 of Appendix A), which states that

This equation may be regarded as expressing gl as a linear combination of weighted

and displaced 6 functions; the elementary functions of the decomposition are, of course,
just these 6 functions.
To find the response of the system to the input gl, substitute (2-43) in (2-41):

Now, regarding the number gl(5, q ) as simply a weighting factor applied to the ele-
mentary function 6(x1 - 5, yl - q), the linearity property (2-42) is invoked to allow
S { )to operate on the individual elementary functions; thus the operator S{}is brought
within the integral, yielding

As a final step we let the symbol h(x2, y2; 6, q ) denote the response of the system at
point (x2, y2) of the output space to a 6 function input at coordinates (6, q ) of the input
space; that is,

The function h is called the impulse response (or in optics, the point-spreadfunction)
of the system. The system input and output can now be related by the simple equation

This fundamental expression, known as the superposition integral, demonstrates

the very important fact that a linear system is completely characterized by its responses
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 21

to unit impulses. To completely specify the output, the responses must in general be
known for impulses located at all possible points in the input plane. For the case of
a linear imaging system, this result has the interesting physical interpretation that the
effects of imaging elements (lenses, stops, etc.) can be fully described by specifying
the (possibly complex-valued) images of point sources located throughout the object
field.

2.3.2 Invariant Linear Systems: Transfer Functions

Having examined the input-output relations for a general linear system, we turn now to
an important subclass of linear systems, namely invariant linear systems. An electrical
network is said to be time-invariant if its impulse response h(t; 7)(that is, its response at
time t to a unit impulse excitation applied at time 7)depends only on the time difference
(t - 7 ) . Electrical networks composed of fixed resistors, capacitors, and inductors are
time-invariant since their characteristics do not change with time.
In a similar fashion, a linear imaging system is space-invariant (or equivalently,
isoplanatic) if its impulse response h(x2, y2; 5 , q ) depends only on the distances (x2 -6)
and (y2 - q ) (i.e. the x and y distances between the excitation point and the response
point). For such a system we can, of course, write
h(x2, ~ 25,; q ) = h(x2 - 6, y2 - q).
Thus an imaging system is space-invariant if the image of a point source object changes
only in location, not in functional form, as the point source explores the object field. In
practice, imaging systems are seldom isoplanatic over their entire object field, but it
is usually possible to divide that field into small regions (isoplanatic patches), within
which the system is approximately invariant. To completely describe the imaging sys-
tem, the impulse response appropriate for each isoplanatic patch should be specified;
but if the particular portion of the object field of interest is sufficiently small, it often
suffices to consider only the isoplanatic patch on the optical axis of the system. Note
that for an invariant system the superposition integral (2-47) takes on the particularly
simple form
m

(2-49)

which we recognize as a two-dimensional convolution of the object function with the

impulse response of the system. In the future it will be convenient to have a short-
hand notation for a convolution relation such as (2-49), and accordingly this equation
is written symbolically as

where a @ symbol between any two functions indicates that those functions are to be
convolved.
The class of invariant linear systems has associated with it a far more detailed math-
ematical structure than the more general class of all linear systems, and it is precisely
22 Introduction to Fourier Optics

because of this structure that invariant systems are so easily dealt with. The simplicity of
invariant systems begins to be evident when we note that the convolution relation (2-49)
takes a particularly simple form after Fourier transformation. Specifically, transforming
both sides of (2-49) and invoking the convolution theorem, the spectra G2(fx, fy) and
Gl(fxJ fy) of the system output and input are seen to be related by the simple equation

where H is the Fourier transform of the impulse response

The function H, called the transfer function of the system, indicates the effects of the
system in the "frequency domain". Note that the relatively tedious convolution opera-
tion of (2-49) required to find the system output is replaced in (2-50) by the often more
simple sequence of Fourier transformation, multiplication of transforms, and inverse
Fourier transformation.
From another point of view, we may regard the relations (2-50) and (2-51) as indi-
cating that, for a linear invariant system, the input can be decomposed into elementary
functions that are more convenient than the 6 functions of Eq. (2-43). These alternative
elementary functions are, of course, the complex-exponential functions of the Fourier
integral representation. By transforming gl we are simply decomposing the input into
complex-exponential functions of various spatial frequencies ( fxJ fy). Multiplication
of the input spectrum G1 by the transfer function H then takes into account the effects
of the system on each elementary function. Note that these effects are limited to an
amplitude change and a phase shift, as evidenced by the fact that we simply multiply
the input spectrum by a complex number H ( fx, fy) at each (fx, fy). Inverse transfor-
mation of the output spectrum G2 synthesizes the output g2 by adding up the modified
elementary functions.
The mathematical term eigenfunction is used for a function that retains its original
form (up to a multiplicative complex constant) after passage through a system. Thus
we see that the complex-exponential functions are the eigenfunctions of lineal; invari-
ant systems. The weighting applied by the system to an eigenfunction input is called
the eigenvalue corresponding to that input. Hence the transfer function describes the
continuum of eigenvalues of the system.
Finally, it should be strongly emphasized that the simplifications afforded by
transfer-function theory are only applicable for invariant linear systems. For applica-
tions of Fourier theory in the analysis of time-varying electrical networks, the reader
may consult Ref. [158]; applications of Fourier analysis to space-variant imaging
systems can be found in Ref. [199].

2.4
TWO-DIMENSIONAL SAMPLING THEORY

It is often convenient, both for data processing and for mathematical analysis pur-
poses, to represent a function g(x, y) by an array of its sampled values taken on a
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 23

discrete set of points in the (x, y) plane. Intuitively, it is clear that if these samples
are taken sufficiently close to each other, the sampled data are an accurate representa-
tion of the original function, in the sense that g can be reconstructed with considerable
accuracy by simple interpolation. It is a less obvious fact that for a particular class of
functions (known as bandlimited functions) the reconstruction can be accomplished ex-
actly, provided only that the interval between samples is not greater than a certain limit.
This result was originally pointed out by Whittaker [298] and was later popularized by
Shannon [259] in his studies of information theory.
The sampling theorem applies to the class of bandlimited functions, by which we
mean functions with Fourier transforms that are nonzero over only a finite region R
of the frequency space. We consider first a form of this theorem that is directly analo-
gous to the one-dimensional theorem used by Shannon. Later we very briefly indicate
improvements of the theorem that can be made in some two-dimensional cases.

2.4.1 The Whittaker-Shannon Sampling Theorem

To derive what is perhaps the simplest version of the sampling theorem, we consider a
rectangular lattice of samples of the function g, as defined by

gdx, y) =
(4 (3
comb - comb - g(x, y).

The sampled function g, thus consists of an array of S functions, spaced at intervals of

width X in the x direction and width Y in the y direction, as illustrated in Fig. 2.5. The
area under each 6 function is proportional to the value of the function g at that particular
point in the rectangular sampling lattice. As implied by the convolution theorem, the
spectrum Gsof g, can be found by convolving the transform of comb(x1X) comb(y1Y)
with the transform of g, or

G,y(
fx, f ~ =
)
{ (3 I);(
F comb - comb - 8 G(fx, f ~ )

FIGURE 2.5
The sampled function.
24 Introduction to Fourier Optics

FIGURE 2.6a
Spectrum of the original function.

where the @ again indicates that a two-dimensional convolution is to be performed.

Now using Table 2.1 we have

1 (4I):(
F comb - comb - = XY comb(Xfx)comb(Yfy)

while from the results of Prob. 2-1 b,

FIGURE 2.6b
Spectrum of the sampled data (only three periods
are shown in each direction for this infinitely
periodic function).
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 25

It follows that

Evidently the spectrum of g, can be found simply by erecting the spectrum of g about
each point (nlX, mlY) in the (fx, fY) plane as shown in Fig. 2.6b.
Since the function g is assumed to be bandlimited, its spectrum G is nonzero over
only a finite region R of the frequency space. As implied by Eq. (2-53), the region over
which the spectrum of the sampled function is nonzero can be found by constructing
the region R about each point (nlX, mlY) in the frequency plane. Now it becomes clear
that if X and Y are sufficiently small (i.e. the samples are sufficiently close together),
then the separations 11X and 1/Y of the various spectral islands will be great enough
to assure that the adjacent regions do not overlap (see Fig. 2.6b). Thus the recovery of
the original spectrum G from Gs can be accomplished exactly by passing the sampled
function g, through a linear invariant filter that transmits the term (n = 0, m = 0) of
Eq. (2-53) without distortion, while perfectly excluding all other terms. Thus, at the
output of this filter we find an exact replica of the original data g(x, y).
As stated in the above discussion, to successfully recover the original data it is
necessary to take samples close enough together to enable separation of the various
spectral regions of Gs. To determine the maximum allowable separation between sam-
ples, let 2Bx and 2By represent the widths in the fx and fy directions, respectively, of
the smallest rectangle7 that completely encloses the region R . Since the various terms
in the spectrum (2-53) of the sampled data are separated by distances 11X and IIY in
the fx and fy directions, respectively, separation of the spectral regions is assured if
1 1
X I - and YI -
2Bx 2By'
The maximum spacings of the sampling lattice for exact recovery of the original func-
tion are thus (2Bx)-' and BY)-'.
Having determined the maximum allowable distances between samples, it remains
to specify the exact transfer function of the filter through which the data should be
passed. In many cases there is considerable latitude of choice here, since for many
possible shapes of the region R there are a multitude of transfer functions that will pass
the (n = 0, m = 0) term of Gs and exclude all other terms. For our purposes, however,
it suffices to note that if the relations (2-54) are satisfied, there is one transfer function
that will always yield the desired result regardless of the shape of R, namely

The exact recovery of G from Gs is seen by noting that the spectrum of the output of
such a filter is

'For simplicity we assume that this rectangle is centered on the origin. If this is not the case, the arguments
can be modified in a straightforward manner to yield a somewhat more efficient sampling theorem.
26 Introduction to Fourier Optics

The equivalent identity in the space domain is

where h is the impulse response of the filter,

h(x,y) = { (kX)
3-' rect - rect(&-)} = 4BxBy sinc(2Bxx) sinc(2Byy).

Noting that

Eq. (2-56) becomes

Finally, when the sampling intervals X and Y are taken to have their maximum allow-
able values, the identity becomes

Equation (2-57) represents a fundamental result which we shall refer to as the

Whittaker-Shannon sampling theorem, It implies that exact recovery of a bandlimited
function can be achieved from an appropriately spaced rectangular array of its sampled
values; the recovery is accomplished by injecting, at each sampling point, an interpola-
tion function consisting of a product of sinc functions, where each interpolation function
is weighted according to the sampled value of g at the corresponding point.
The above result is by no means the only possible sampling theorem. Two rather
arbitrary choices were made in the analysis, and alternative choices at these two points
will yield alternative sampling theorems. The first arbitrary choice, appearing early
in the analysis, was the use of a rectangular sampling lattice. The second, somewhat
later in the analysis, was the choice of the particular filter transfer function (2-55).
Alternative theorems derived by making different choices at these two points are no less
valid than Eq. (2-57); in fact, in some cases alternative theorems are more "efficient" in
the sense that fewer samples per unit area are required to assure complete recovery. The
reader interested in pursuing this extra richness of multidimensional sampling theory
is referred to the works of Bracewell [31] and of Peterson and Middleton [230]. A
more modern treatment of multidimensional sampling theory is found in Dudgeon and
Mersereau [85].
 CHAPTER 2 Analysis of Two-Dimensional Signals and Systems 27

2.4.2 Space-Bandwidth Product

It is possible to show that no function that is bandlimited can be perfectly space-limited

as well. That is, if the spectrum G of a function g is nonzero over only a limited re-
gion R in the (fx, fy) plane, then it is not possible for g to be nonzero over only a
finite region in the (x, y) plane simultaneously. Nonetheless, in practice most functions
do eventually fall to very small values, and therefore from a practical point-of-view
it is usually possible to say that g has signijicant values only in some finite region.
Exceptions are functions that do not have Fourier transforms in the usual sense, and
have to be dealt with in terms of generalized Fourier transforms (e.g. g(x, y) = 1,
g(x, Y) = cosP.rr(fxx + fry)], etc.).
If g(x, y) is bandlimited and indeed has significant value over only a finite region
of the (x, y) plane, then it is possible to represent g with good accuracy by aJinite number
of samples. If g is of significant value only in the region - Lx 5 x < Lx, -Ly Iy <
Ly, and if g is sampled, in accord with the sampling theorem, on a rectangular lattice
with spacings (~Bx)-', BY)-' in the x and y directions, respectively, then the total
number of significant samples required to represent g(x, y) is seen to be

which we call the space-bandwidth product of the function g. The space-bandwidth

product can be regarded as the number of degrees of freedom of the given function.
The concept of space-bandwidth product is also useful for many functions that
are not strictly bandlimited. If the function is approximately space-limited and ap-
proximately bandlimited, then a rectangle (size 2Bx X 2By) within which most of the
spectrum is contained can be defined in the frequency domain, and a rectangle (size
2Lx x 2Ly) within which most of the function is contained can be defined in the space
domain. The space-bandwidth product of the function is then approximately given by
Eq. (2-58).
The space-bandwidth product of a function is a measure of its complexity. The
ability of an optical system to accurately handle inputs and outputs having large space-
bandwidth products is a measure of performance, and is directly related to the quality
of the system.

PROBLEMS-CHAPTER 2

2-1. Prove the following properties of 6 functions:

(a) @x, by) = &S(x, y).

m m

(b) comb(ax) comb(by) = ,& 1 1 8 (x z, y

n=-mm=-m
- -
m
%).

2-2. Prove the following Fourier transform relations:

(a) F{rect(x) rect(y)) = sinc(fx) sinc(fy).