GEOMETRIC DESIGN TECHNIQUE FOR
TWO-DIMENSIONAL ANALOG AND DIGITAL FILTERS
Cosk&n Zumrutkaya
Thesis
Z8
NAVAL POSTGRADUATE SCHOOL
Monterey, California
THESIS
TWO-
GEOMETRIC DESIGN TECHNIQUE FOR
DIMENSIONAL ANALOG AND DIGITAL FILTERS
by
Coskun Zumrutkaya
June 19 7 9
Thesis
Advisor: S. R.
Parker
Approved for public release; distribution unlimited
UNCLASSIFIED
SECURITY CLASSIFICATION OF THIS PAGE (Whan Data Entered)
REPORT DOCUMENTATION PAGE
I a.£»OAT SUMII*
READ INSTRUCTIONS
BEFORE COMPLETING FORM
2. GOVT ACCESSION NO.
1. RECIPIENT'S CATALOG NUMBER
4. T|TlE fend Subtttl*)
Geometric Design Technique for Two-
Dimensional Analog and Digital Filters
5- TYPE OF REPORT a PERIOD COVERED
Engineer's Thesis;
June 197 9
« PERFORMING ORG. REPORT NUMBER
7. AUTHORS
Coskun Zumrutkaya
•. CONTRACT OR GRANT NUMBERS.)
>. PERFORMING ORGANIZATION NAME AND ADDRESS
Naval Postgraduate School
Monterey, California 93940
10. PROGRAM ELEMENT. PROJECT, TASK
AREA « WORK UNIT NUMBERS
II. CONTROLLING OFFICE NAME ANO AOORESS
Naval Postgraduate School
Monterey, California
12. REPORT DATE
June 19 7 9
IS. NUMBER OF PAGES
153
U. MONITORING AGENCY NAME ft AOORESSf different from Controlling Ottlee)
IS. SECURITY CLASS, (of thio report)
Unclassified
<Sa. OCCLASSIFI CATION/ DOWNGRADING
SCHEDULE
16. DISTRIBUTION STATEMENT for thit Report)
Approved for public release; distribution unlimited.
17. DISTRIBUTION STATEMENT (o: thm omotrmct entered In Block 30, II different from Report)
IS. SUPPLEMENTARY NOTES
IS. KEY WORDS 'Continue on rereree tide it neceeeerr and Identity or block nimtber)
Analog Filter
Digital Filter
Extended Bode approach
20. ABSTRACT (Continue an referee eide It I
•ary and Identity ey bleek
neer)
A technique for the use of semi-infinite planes to approximate
the log magnitude versus log frequency characteristics of
two-dimensional filters is presented (Extended Bode approach) .
This technique is applied to quarter plane filters and works
well for separable transfer functions. Other non-separable
canonic (basic) transfer functions are also studied in terms
of their planar approximation. The general technique is shown
DO , r ™ n 1473
(Page 1)
EDITION OF 1 MOV «» IS OBSOLETE
S/N 0103-014-6601 |
UNCLASSIFIED
SECURITY CLASSIFICATION OF THIS PAGE fWhan Dote Entered)
UNCLASSIFIED
f*eu*wTv claudication q» Tmt mtQtrw*,^ n»i« i.<, K<
(20. ABSTRACT, Continued)
to be useful for the insight it provides as well as being a
simple approach to design. The double bilinear z-transform
is studied for use with the two-dimensional analog transfer
functions to convert them to the digital domain.
DD Form 1473
„, 1 Jan 73
S/N 0102-014-6601
1 Jan 73 UNCLASSIFIED
N 0102-014-6601 2 iieuAW classification o* ▼»••• natr**™ o«« *«#•*•*)
Approved for public release; distribution unlimited.
Geometric Design Technique for
Two-Dimensional Analog and Digital Filters
by
Coskun Zumrutkaya
Lieutenant /''Turkish Navy
B.S.E.E., Naval Postgraduate School, 1977
M.S.E.E., Naval Postgraduate School, 1979
Submitted in partial fulfillment of the
requirements for the degree of
ELECTRICAL ENGINEER
from the
Naval Postgraduate School
June 1979
ABSTRACT
A technique for the use of semi-infinite planes to
approximate the log magnitude versus log frequency charac-
teristics of two-dimensional filters is presented (Extended
Bode approach) . This technique is applied to quarter plane
filters and works well for separable transfer functions.
Other non-separable canonic (basic) transfer functions are
also studied in terms of their planar approximation. The
general technique is shown to be useful for the insight it
provides as well as being a simple approach to design. The
double bilinear z-transform is studied for use with the
two-dimensional analog transfer functions to convert them
to the digital domain.
TABLE OF CONTENTS
I. INTRODUCTION 8
A. DIGITAL SIGNAL PROCESSING 8
B. AREA OF INVESTIGATION AND OBJECTIVES 18
C. PREVIEW OF RESULTS 20
II. FUNDAMENTALS OF TWO-DIMENSIONAL
RECURSIVE FILTERING 21
A. INTRODUCTION 21
B. DEFINITIONS : 22
1. Two-Dimensional Signals 22
2. Two-Dimensional Systems 28
3. Causality, Separability, Stability 31
4. Two-Dimensional Difference Equations - 36
5. Two-Dimensional z-Tranform 39
6. Two-Dimensional DFT 45
C. STRUCTURES OF TWO-DIMENSIONAL FILTERS 4 7
1. One-Dimensional Digital Filter
Realization 47
2 . Two-Dimensional Digital Filter
Realization 53
D. STABILITY 61
III. MATHEMATICAL BACKGROUND FROM ANALYTIC
PLANAR GEOMETRY 69
IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY
DOMAIN 87
A. SEPARABLE FACTORS 87
B. NON-SEPARABLE FACTORS 10 8
C. CONCLUSIONS 118
V. TWO-DIMENSIONAL DIGITAL FILTER DESIGN
TECHNIQUES 122
A. INTRODUCTION 122
B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN
PROCEDURE BY DOUBLE Z-TRANSFORM 125
1. Double z-Transform 125
2. Design Procedure 130
C. DIGITAL CANONIC FACTORS AND STABILITY
CONSIDERATIONS 132
1. Introduction 132
2. Separable Factors 132
3. Non-Separable Factors 133
4. Conclusions 145
VI. CONCLUSIONS 148
LIST OF REFERENCES 150
INITIAL DISTRIBUTION LIST 153
ACKNOWLEDGMENT
I am fortunate to have this opportunity to thank the
many people who have helped in the preparation of this
thesis. First of all, I am most indebted to my advisor.
Dr. Sidney R. Parker, for his encouragement, patience,
guidance and instruction.
I want to thank the Turkish Navy for the generous
opportunity to pursue this research.
I want to thank Mrs. Rosemary Lande for her excellent
typing.
Last, but not least, I would like to thank my wife,
Sumru and my son Devrim, for their patience and understanding
I. INTRODUCTION
A. DIGITAL SIGNAL PROCESSING
Digital signal processing has its roots in 16th century
mathematics, especially in the fields of astronomy and the
compilation of mathematical tables. Today it has become
a powerful tool in a multitude of diverse fields of science
and technology. The applications of digital signal pro-
cessing varies from low-frequency spectrum seismology
through spectral analysis of speech and sonar into the
video spectrum of radar systems [1] .
In their book, Theory and Applications of Digital
Signal Processing [2] , Rabiner and Gold have combined digi-
tal signal processing theory, with a variety of applications
ranging from sonar, radar, communication, music, seismic
and medical signal processing, with digital component
technology. This technology is the main driving force for
progress in this field as well as the general area of com-
puter design. It was also observed by them, that, although
the formulation of engineering problems is often as vague
as those of the "softer" sciences (such as anthropology,
psychology, etc.), the execution of these problems appears
to depend on greater and greater accuracy and reproducibility
The capability of digital systems to achieve a guaranteed
accuracy and essentially perfect reproducibility is very
appealing to engineers.
A digital filter is defined in [3] to be a computa-
tional process in which a sequence of numbers acting as
input, is transformed into a second sequence of numbers,
or output digital signal; where the term "digital" implies
that both time (the independent variable) and amplitude are
quantized.
The field of one-dimensional digital filtering, which
is outlined in figure 1.1, encompasses recursive, non-
recursive and Fast-Fourier Transform processing. The terms
recursive and nonrecursive, instead of IIR (Infinite Impulse
Response) and FIR (Finite Impulse Response) are preferred
in this thesis. It is shown that recursive processing is
much more efficient than nonrecursive processing. Stockham's
method [4] to perform fast convolution, which later became
known as the FFT method, improved the efficiency of non-
recursive techniques, so that comparisons in one-dimension
are no longer strongly biased toward recursive methods [2] .
One-dimensional recursive digital filter design depends
strongly on the effective and well developed continuous
filter design theory. Moreover, the stability analysis of
higher order one-dimensional recursive realizations can be
solved by investigation of the root distribution of its
factored form representation. It is important to realize
that continuous domain design techniques and factorization
property (the fundamental theorem of algebra) exists only
in one-dimension.
a
CO Eh
>H fa
en Q
i— i
Q
I
8
w
Q
B
g
£
us
I
H
fa fa
o „
z fa
Eh
CO fa
fa
, M
O fa
as
II
fa
CD
9
fa
H
fa
O fa
S3
H Z
en
I
H
fa
388
S l_T 2
O H H Q
u S R ^
las§
O H fa
< i< H O H O M
£ .3 2 H fa U Q
03
2
fa
Eh
02
>H
cn
fa
Eh
H
U
H
Q
fa
Z
O
H
Z
fa
H
a
i
m
z
o
fa
o
fa
Hi
>
fa
fa
>
o
fa
fa
D
10
Areas of application of one-dimensional digital signal
processing are listed in figure 1.1. There are important
military and specifically naval applications in addition
to those listed, such as missile and torpedo guidance,
launcher control system, fire control, combat information
collection, dissemination and display.
Digital filtering in several dimensions, which is over-
viewed in figure 1.2, has gained considerable importance,
especially for the two-dimensional case. Figure 1.3 presents
an example of two-dimensional low-pass and high-pass filtering
Until 196 6, two-dimensional digital filtering was implemented
by nonrecursive or convolutional techniques. In this method
the output is the weighted sum of unit sample responses of
all past input values. The serious disadvantage of the
convolutional method is the requirements of a very large
number of arithmetic operations.
The development of the Fast Fourier Transform (FFT) in
1966, reduced the number of arithmetic operations consider-
ably and is used extensively today. Filtering via FFT is
accomplished by computing the transform of the input func-
tion, multiplying by the frequency response of the filter,
and inverse transforming the result. The recursive
algorithm, which has in general an infinite impulse response
[3] , constitutes another technique for the realization of
two-dimensional digital filters.
Hall [5] compares the amount of computation required for
the three filtering techniques. He shows that, in general,
the FFT and recursive algorithms are preferable to the
11
02
T
*—*
co
o
S
1-1
s
CO
En
3
en
W
>H
2
co
H
Q
J
1
<
Z
Eh
H
s
O
H
M
Q
•r 1
K
fa
O
O
fa
CO
>H
2
05
O
£
fa
Eh
K
1
Eh
CS3
i
1
fa
OS
D
Eh
U
D
«
Eh
CO
CO
H
CO
>•
<
z
<
fa
Eh
<
Eh
CO
05
a
ij
H
& b
§£
H CO
a
z;
o
z
eg
fa
Eh
CO
>H
CO
fa
<
Eh
H
U
fa
<:
z
o
H
CO
fa
S
H
Q
I
Z
fa
o
fa
HI
>
fa
>
o
CM
fa
fa
CJ
H
fa
12
Fig. 1.3. a. An example of two-dimensional filtering
(An original photograph) [2]
13
Fig. 1.3.b.
An example of two-dimensional filtering
(Low pass filtered version of Fig. 1.3a)
[2]
14
.'.-^■■'/r.-a"' ,\'^y' : ^*-"' '■' ■ \*
■ * S '
Fig. 1 . 3c.
An example of two-dimensional filtering
(High pass filtered version of Fig. 1.3a)
[2]
15
nonrecursive one in number of computations and storage
requirements. Hall demonstrates that the recursive filter
algorithm constitutes the best method for large data, i.e.,
it is the fastest and cheapest.
When processing two-dimensional data by recursive filters
a fundamental problem exists due to inherent feedback, namely,
the problem of numerical stability. Since, in general, two
variable polynomials cannot be factored into a product of
first and second order real coefficient polynomials in each
of the variables, it is difficult to solve the stability
problem for two-dimensional recursive filtering. Conse-
quently, the majority of papers published in two-dimensional
digital filtering deal with the design of nonrecursive
filters which are inherently stable [6-10].
There are several papers discussing two-dimensional
recursive digital filters. For example [11] and [12]
formulate two-dimensional recursive filters by the z-
transform and linear difference equations, and although
they investigate problem areas related to stability and
realization, the majority of problems remain to be solved.
The most important applications of two-dimensional
digital filtering of digital data are in picture processing
and geophysical data analysis. Picture processing can be
categorized into:
a) Digital image restoration and enhancement,
b) Computer pictorial pattern recognition.
16
The methods of image restoration and enhancement are
applied to invert degradations, such as aberrations,
atmospheric effect, scanning, motion, and to manipulate
images to improve viewing phenomena experienced by the human
eye. Image restoration and enhancement has been used with
great success in biomedicine, i.e., in extraction of quan-
titative information from x-ray films, chromosome counting,
measuring the extent of arteriosclerosis from arteriograms,
in forensic sciences application, i.e., fingerprint image
enhancement for automatic classification, and in astronomy,
i.e., removal of turbulence from astronomical photography
[13].
Pictoral pattern recognition by digital computer has
gained great importance in satellite surveillance and mili-
tary reconnaissance. The application of two-dimensional
digital filters permits the separation of different horizontal
scans on magnetic and gravity maps in geophysics and can
be used equally well for structural and topographic maps
or for any other type of factor which is available in the
format of a planar grid [14] . Although digital filtering
in several dimensions is gaining importance in medicine,
i.e., heart volume measurements, in fire control problems,
and to analyze complex electronic circuitry, there exists
no general theory concerning structures, analysis and
design.
These considerations can be summarized as follows.
Since one-dimensional recursive digital analysis and design
17
techniques depend strongly on the fundamental theorem of
algebra and on extensive utilization of continuous design
theory, it represents a special, i.e., a non-generalizable
field in the area of N-dimensional recursive digital filtering.
Multi-dimensional digital filtering is well developed in two-
dimensions but is based, due to the absence of recursive theory
analysis and design tools, predominantly on nonrecursive
filtering schemes. Because of their inherent advantages, much
of the current research is directed towards recursive filtering
algorithms and techniques.
B. AREA OF INVESTIGATION AND OBJECTIVES
One area of investigation of this thesis is to see how the
semi-infinite straight line approximation technique of one-
dimensional filters (based upon the log modulus/log frequency
or decibel vs. log- frequency plots) can be extended to two
(or higher) dimensional recursive filters. The advantages of
log-modulus approach is that a transfer function in factored
form becomes a linear sum of terms when the logarithm of its
magnitude is taken, and each of these terms can be approxi-
mated by semi-infinite straight line segments in the one-
dimensional case, and semi-infinite planes in the multi-
dimensional case, when expressed in terms of log-frequency.
The following are specific objectives:
1. To postulate typical (canonic) factors for
two-dimensional filters and to investigate their properties
in the log-modulus/log frequency domain.
18
2. To develop an approximation technique for the
design of two-dimensional recursive filters to meet given
specifications in the frequency domain.
The frequency response of two dimensional filters is
inherently a volume in the space of magnitude, log cj, and
log CO2 • This volume is the given or desired specification
of the filter, which in the separable case can be approxi-
mated by the combination of semi-infinite planes.
In this thesis we investigate how we can approximate
this volume by semi-infinite planes as an extension of
the approximation used in one dimension, i.e., by semi-
infinite lines. After determining the planar equations
which perform the given specification of the filter, we
investigate ways to get a transfer function; first, in the
continuous frequency domain and second for two-dimensional
recursive equations in the discrete variable domain.
The results of this study will classify those charac-
teristics that can be achieved by cascading simple canonic
factors, and indicate limitations that are imposed by
trying to achieve a particular frequency characteristic
(such as the fan filter) in as simple form as possible.
The approach taken here is to investigate the problem
in the continuous domain first, and then study the trans-
formation from continuous frequency (differential equation)
domain to discrete time (sampled data) domain.
19
C. PREVIEW OF RESULTS
We have proposed a geometric method to design two-
dimensional recursive digital filters. This is a new
method which is relatively simple and easy with respect to
other design methods. We have investigated properties of
some typical (canonic) factors and considered their stability
properties. In Chapter II we go through the basic properties
of two-dimensional filtering. In Chapter III we provide
background in terms of planar geometry for the discussion
in subsequent chapters. In Chapter IV we investigate the
properties of typical factors and their stability in continu-
ous variable domain. The new design technique, including
stability considerations, starting with the continuous
frequency domain and extending to two-dimensional recursive
filters in the discrete domain, is presented in Chapter V.
The results are summarized and some comments for further
study are made in Chapter VI. It should be noted that the
results presented apply to quarter plane filters, as defined
in Chapter V.
20
II. FUNDAMENTALS OF TWO-DIMENSIONAL RECURSIVE FILTERING
A. INTRODUCTION
There are many signals that are inherently two-dimensional
for which two-dimensional signal processing techniques are
required. The most common examples are, of course, images
in which the two variables are the spatial coordinates.
Image processing [15], [16] plays an important role in many
areas of scientific and technical research. Some examples
are satellite imagery radiography, radar, and biomedical
images such as acoustical holograms, medical x-rays, and
electron micrographs [2]. Not all two-dimensional data
come from an image. In prospecting for gas and oil, a
common procedure is to monitor the seismic signals produced
by detonating an explosive charge [17] . An array of geo-
phones situated along a line radiating in a vertical or
horizontal direction from the point of explosion is used
to record these signals. In this case, the two variables
are distance (from the point of explosion measured along
the line) and time (the time at which the seismic signal
reaches a given geophone or, equivalently , a given distance
from the explosion) .
Although two-dimensional signals may be processed by
one-dimensional systems, it is generally preferable to con-
sider using two-dimensional systems. Many of the basic
ideas of one-dimensional signal processing may be readily
21
extended to the two-dimensional case. There are some
very important concepts of one-dimensional systems, however,
that are not directly extentable to two dimensions. It
is the goal of this chapter to discuss the basic ideas and
techniques of two-dimensional signal processing and to
illustrate them into the context of two-dimensional filter
design. A more detailed discussion of this theory plus
an overview of recent work in two-dimensional filtering is
contained in the excellent survey paper by Mersereau and
Dudgeon [18] , and H. Chang and J.K. Aggarwal [19] . Other
useful references are the first few chapters of [1] and
chapter 7 of [2] .
B. DEFINITIONS
1. Two-Dimensional Signals
One-dimensional signals are functions of a single
variable and two-dimensional signals are functions of two
integer variables. Consider the two-dimensional sequence
x(m,n) where m and n are integer variables, shown in figure
2.1. As in the one-dimensional case, the notation x(m,n)
is often short hand for a sampled version of a continuous
two-dimensional signal x(s,t): i.e.,
x(m,n) = x(mT,,nT 2 ) = x(s,t)
(2.D
s=mT 1 ,t=nT 2
One interpretation of x is a function which assigns a
(generally complex) number to each integer ordered pair
22
- 7 y- it/a M
*
Fig. 2.1. Pictorial representation for two-
dimensional discrete array
X P £$ flit \-$l> J
23
(m,n) . The function x is not defined when either or both
of its arguments (m,n) is not an integer. Such a signal
will be interchangeably referred to as either a two-dimensional
array, a two-dimensional sequence or simply as a two-dimensional
signal .
Some useful two-dimensional sequences are defined
below and shown in figures 2.2 to 2.4. These include:
1. Digital Impulse or Unit Sample
1 m = n =
u (m,n) = \ (2.2)
elsewhere
2. Digital Step
(
I 1 m,n >_
u ,(m,n) = \ (2.3)
elsewhere
3. Exponential
m n . n
a, a~ m,n
x(m,n) = X (2.4)
otherwise
4. Sinusoid (Complex)
j ((jo^m+oj^n)
x(m,n) = e -°° < m,n < +=° (2.5)
As seen above, the two-dimensional step is related to the
two-dimensional impulse by
24
/
■+
+
+-H-4
-t
■f- f- -h-h
■/ -h-h -h +- /-■/- +-/-++-
-h- -h ■/- f- -H f 4 A -h A 4 +■
f +■ 4 + f 4 4 4 -h + +
4- 4- 4- -hif- 4 4- 4 -h-f- -h
-4 -h -i-4- f- -f- + -h-h-h 4 i-
-h -h 4- -h +- -f-f- 4 1-+- j-4-
+ j-4 4- f -h- 4- ■/- 4 -h -h
i--l-4--/-4- -/-4-4--h+-
-h4—h 4-/- 4- 4 +■ 4 4
I / / > rn.
Fig. 2.2. Digital Impulse
25
- I * * -f 4
f +■ + -t-_zf- X XXX X
-h-h + + f XXXXXX
-h + + 4-4- XX xxxx
-h-h-h h-h.XX XXXX,
-h-t- 4- 4- -4 .XXX XXXX
+.+ + + i .XiXXxxx
/ / / / / / y ) } y } ) > "
4--h4-4-4-4-4-4-i--h4- 4-
-h -h 4- 4- f -4 4 h^4 4
4--A 4- 4 4 4 4 4 + 4--h
-h 4 -h-j-h- 4- 4- -h 4- -h
+ + -h -/ 4 -L- -h 4 4 4
Fig. 2.3. Digital Step
26
-t- -h H- -h -j- -J- -f- ■f- f- -f-
-*- + + ++• f--*- ff f
-/- -f- ■+ 4-+ f 4- 4-f- *■ *■
-f- + + -h -h f 1 * f- *-*
■t- -t + -h HfA- J- j- -f- f- +
m
4 4 -h -h 4 4 4- + -f- -h -h +■
4- + 4 + + 4- 4- + -h +- -h -t-
+-4-4 -h 4 4- 4 4- 4 4 f-
-4 4- 7-/ -h -h A -h +~ 4-
Fig. 2.4. Exponential
27
m n
u_ 1 (m,n) =11 u o (m 1 ,n 1 ) (2.5a)
m, =-<» m, =-oo
2 . Two- Dimensional Systems
A two-dimensional discrete system is defined as
the unique transformation T that maps an input sequence
(x(m,n)} and the initial condition sequence {s(m,n)} into
an output sequence {y(m,n)}
y(m,n) = T[{x(m,n) },{s(m,n) }] (2.6)
In signal processing applications, we assume zero initial
conditions. With this assumption, the above relation may
be shortened to
y(m,n) = T[{x(m,n)}] (2.7)
The system characterized by T is
(a) linear if
T[{ax 1 (m,n) + bx 2 (m,n) }] = aT [ {x 1 (m,n) } ] + bT [ {x 2 (m,n) }]
(2.8)
for arbitrary constants a and b.
(b) shift invariant if
y(m-m 1 ,n-n 1 ) = T [ {x (m-ir^ , n-n 1 ) } ] (2.9)
for all m, and n, .
28
A hardware and software system which is equivalent to the
transformation T of (2.6) is called a digital filter. The
class of complex-valued sequences which are square summable
will play the role of input and output signals of the
digital filter. Linearity and shif t-invariance are indepen-
dent properties of a system in that neither property implies
the other. Thus the system
T[x(m,n)] = h(m,n) x(m,n)
is linear but not shift invariant and the system
2
T[x(m,n)] = x (m,n)
is shift invariant but not linear
For two-dimensional discrete linear systems, the
impulse response is defined as the response of the system
at (m,n) to a two-dimensional digital impulse input at
(m, /n,) with zero initial conditions, i.e.,
h(m,n;m,,n, ) = T [ {u q (m-m 1 ,n-n 1 ) } ] (2.10)
where u (m,n) denotes the digital impulse. In particular,
the impulse response of a two-dimensional discrete LSI
(Linear Shift Invariant) system is
29
h (m,n;m, ,n, ) = h (m-m, ,n-n, )
or
h(m,n) = T[{u o (m,n) }] (2.11)
The response of a two-dimensional LSI* system with zero
initial conditions is completely characterized by its
impulse response. In fact using the identity
x(m,n) = I x(m, ,n , ) u (m-m, ,n-n 1 )
l'"l'"o % V "1'
m 1 ,n 1 =-«
I x(m-m 1 ,n-n 1 )u (m 1 ,n 1 )
m^n^- 00
we get
y(m,n) = T[{ x (m, ,n 1 ) u (m-m, ,n-n, ) } ]
ra r n r
l x(m 1 ,n 1 )T[{u (m-m 1 ,n-n 1 ) }]
m 1 /n 1 =- c
x (m, ,n, ) h (m-m, ,n-n^)
m l' n l = ~°° '
*LSI = Linear Shift Invariance
30
y(m,n) = £ x (m-m^n-r^) h (m, ,n,)
m 1/ n 1 =-<
= x(m,n) * h(m,n) (2.12)
where "*" denotes two-dimensional convolution. In other words,
for LSI systems the basic convolution theorem is valid.
3. Causality, Separability, Stability
In discussing various features of two-dimensional
discrete systems, it is convenient to have spectral nota-
tion [19] to represent certain regions of the spatial domain.
1. Let S(a,3) denote a sector with the angular
interval (a, 3). In polar form S(a,3) may be described by
S(a,3) = ( (r,6) |r > 0,a < 9 < 3) (2.13)
where (a, 3) are measured in radians. In particular, we
give the following notation for the sectors frequently
encountered:
S++ = S[0,tt/2]; S-+ = S[Tr/2,ir]
S— = S[-tt,-tt/2] ; S+- = S[-tt/2,0] (2.14)
S*+ = S[0,tt) ; S*- = S[-tt,0)
where "[" and "]" imply including the boundary and ")" implies
X
not including the boundary.
31
Notice that S++, S-+, S — , and S+- represent first, second,
third, and fourth quadrants, respectively, and S*+ and S*-
represent half-planes which are symmetric to each other with
respect to the origin. (See fig. 2.5.)
2. For a two-dimensional sequence (x(m,n)}, the
set
E v = { (m,n) |x(m,n) ? 0} (2.15)
is called the support of (x(m,n)}.
3. Let {h(m,n)} be the impulse response of a two-
dimensional discrete LSI system, and let E, be the support
of (h(m,n) }.
Then, a two-dimensional LSI system is said to be
causal if E, is the set S++, and semi-causal if E, is the
set S* + . They are shown in Fig. 2.5.
A two-dimensional discrete LSI system is said
to be separable if its impulse response can be factored into
a product of one-dimensional responses (Fig. 2.6); i.e.,
h(m,n) = h x (m) h 2 (n) (2.16)
If the equation (2.16) is not satisfied, the filter is said
to be nonseparable. The advantages of separable filters is
that two dimensional convolution (2.12) may be carried out
as a sequence of one-dimensional convolutions. This can be
seen by rewriting (2.16) as follows.
32
fcm
Fig. 2.5. Support Area of Impulse Response
(a) Causal
(b) Semi-causal S.
(c) Semi-causal S*_
33
X(«r«.)
H,(-0
H 2 («.)
)
1' 2
H(z) = H 1 (z 1 )H 2 (z 2 )
Fig. 2.6. Separable Digital Systems
34
y(m,n) = j» h 1 (m^) h 2 (r^) x (m-m-^n-n-^ (2.17a)
m, ,n =-<»
I h 1 (m 1 ) [ I h 2 (n 1 )x(m-m 1/ n-n 1 ) ] (2.17b)
m =-oo n,=-oo
I h 1 (m 1 ) a(m-m.,n) (2.17c)
m l =_c
where a (m-m, , n) is a sequence of one-dimensional convolutions
(obtained by evaluating the terms inside the bracket of
(2.17b) for each fixed value of m, ) . Equation (2.17c) shows
that y(m,n) may be obtained by a second sequence one-dimensional
convolution.
If both the input sequence x(m,n) and the filter
impulse response h(m,n) are separable, then it is readily
seen that the output sequence y(m,n) is also separable.
In this case we obtain the result [using (2.12) and (2.16)]
00 00
y(m,n) = [ I h ± Cm^b (m-m 1 ) J I I h 2 (n 1 ) c (n-i^) ] (2.18a)
m =-« n =-oo
= a(m)3(n) (2.18b)
where
x(m,n) = b(m) c(n)
35
V
A two-dimensional LSI system is said to be (BIBO)
stable if every bounded input sequence produces a bounded
output sequence. This condition is satisfied if and only
if the two-dimensional impulse response of the system is
absolutely summable, i.e.,
00
I |h(m,n) | < c° (2.19)
m, n=-°°
As in the one-dimensional case, (2.19) can be shown to be a <e-
necessary and sufficient condition for stability [20]. One
problem with (2.19) is that it can be quite difficult to
evaluate for an arbitrary h(m,n).
4 . Two-dimensional Difference Equations
As in the one-dimensional case, two-dimensional LSI
systems can often be described by a two-dimensional linear
constant-coefficient difference equation relating the output
y(m,n) of the filter to its input x(m,n). The most general
form for such a difference equation is shown by
I a(k 1 ,k 2 )x(m-k 1 ,n-k 2 ) - I b ( ^^ y (m-Jl^n-j^) =
(k 1 ,k 2 )eR a (i 1 ,l 2 )eR h
(2.20)
wh
ere R and R are finite sets of spatial grid points,
called the input and output mask, respectively, and {a(k 1 ,k 2 )}
and (bUwJU)} are the set of constant coefficients that
36
characterize the particular filter. Generally, a difference
equation has a family of solutions as is the case with a
differential equation. The difference equation is said to
be recursive if there exists a sequence of computations which
yields the output sequence serially from the input sequence
and that portion of the output sequence which is already
computed/ including the initial condition. An implicit
ordering of grid points in the spatial domain is assumed for
the serial computation of the output sequence.
There are two cases of (2.20):
1. R , K in the S++ (first quadrant), b(0,0) ^
and y(m,n) = for outside of S++ (zero initial conditions),
then (2.20) represents a causal system.
2. If R , R. in S*+, b(0,0) ? and y(m,n) = for
a o
outside of S*+ (zero initial conditions), then (2.20)
represents a semi-causal system.
In both cases we may assume b(0,0) = 1 without loss of
generality and write (2.20) as
y(m,n) = £ a (k-^k^ x (m-k-^n-k^
(k r k 2 )£R a
I b(l 1 ,!L 2 )Y(m-l 1 ,n-l 2 )
(2.21)
(l lf l 2 ) eR b -(0,0)
It is clear that (2.21) is a recursive equation, and the
implicit ordering, for example may be given as
37
{(0,0); (1,0), (0,1), (1,1); (2,0) ,(2,1) , (0,2) ,(1,2); (2,2);...}
for a causal system, and
{(0,0), (0,1) , (0,2) ,...; (-1,1) , (0,1) , (1,1) ;...;...
... (-1,2) , (0,2) , (1,2) ,.. .; . ..}
for a semi-causal system.
If {b(£ 1 ,£ 2 )} = {u Q (£ 1 ,£ 2 )} the equation (2.21)
becomes
y(m,n) = I a(k 1 ,k 2 )x(m-k 1 ,n-k 2 ) (2.22)
(k 1 ,k 2 ) £ R a
so that the output sequence is a weighted moving average of
the input sequence. The impulse response of (2.22) is
unique and is obviously {a(k,,k 2 )}. Such filters are called
finite extent impulse filters (FIR) because the impulse
response has only finitely many nonzero terms. FIR filters
obviously are always BIBO stable because a(k,,k 2 ) is always
absolutely summable. (It has only finitely many nonzero
terms . )
If, on the other hand, the sequence (b(£,,£ 2 )}
has several nonzero terms, then the impulse response will
generally have infinitely many terms. We will say that such
filters are infinite-extent impulse response filters (IIR) .
38
Such filters may be BIBO unstable, and so it is necessary
to devise a stability test for such filters.
FIR filters have the advantages that they are
easier to design, that controlling the phase response is
easier (in particular, it is easy to produce filters which
have a linear phase response) , and that since the input-
output relationship is a finite-extent convolution, Fast
Fourier Transform (FFT) techniques can be used to speed the
computation of the output sequence. On the other hand,
recursive (IIR) filters generally require less high speed
storage and usually require less computation time than FIR
filters. Comparison of FIR and IIR filters is contained
in [5].
5 . Two-dimensional z-tranform
A useful mathematical tool for representation of a
sequence x(m,n) is its two-dimensional z-transform. Given
a sequence x(m,n), the z-transform of this sequence is
defined to be:
<
CO oo
Z[{x(m,n)}] 4 X(z 1 ,z 2 ) = I £ x (m,n) z^ (2.23)
m=-°° n=-°°
where z, and z 2 are independent complex variables. The most
useful property of the z-transform is the fact that the z-
transform of the convolution of two sequences is the product
of their z-transforms , i.e.:
39
Z[{x(m,n) }*{h(m,n) }] = X(z 1 ,z 2 )H(z 1 ,z 2 ) (2.24)
Therefore, the input-output relationship of an LSI filter
may be expressed in the z-transform domain as:
Y(z 1 ,z 2 ) = H(z 1 ,z 2 )X(z 1 ,z 2 ) (2.25)
where H(z-,,z 2 ), the z-transform of the impulse response
is called the transfer function.
A two-dimensional transfer function evaluated on
the unit circle is called the frequency response of the
system, i.e.,
jw, jco 9 ju>, jo> 9 jui, joo 9
Y(e ,e *) = H(e x ,e )X(e x ,e z ) (2.26)
when a two-dimensional LSI system is describable by the
difference equation of (2.20), and we take a z-transform >'
of this, it follows that
-k -k 2
[ I a(k 1 ,k 2 )z 1 z 2 ]X(.z 1 ,z 2 )
(k x ,k 2 ) eR a
- [ I b(£ 1 ,£ 2 )z 1 z 2 Z ]Y(z L ,z 2 ) =
(i lf l 2 ) el^
40
or
H(z 1 ,z 2 ) = Y(z 1 ,z 2 )/X(z 1 ,z 2 )
' k l ~ k 2
I a(k 1 ,k 2 )z 1 z 2
(k l' k 2 )eR a
-A, -l
I b(£ 1 ,£ 2 )z 1 x z 2 z
( *■ -i / & ^ ) s^j-.
= A(z 1 ,z 2 )/B(z 1 ,z 2 ) (2.27)
where A and B are the z-transf orms of (a(k,,k 2 )} and
{b (£■•* & 2 ) } respectively. In particular when (2.27) repre-
sents a causal or semi-causal system, i.e., (2.27) is the
z-transform of (2.21), we can give explicit descriptions for
R and R, in each case. For the causal case we define
a d
(Fig. 2.5a)
R a = * ( k i' k 2> I ° i k l i M l' °i k 2 - M 2 }
(2.28)
Rj^ = ( (l lf i 2 ) | <_l 1 <_N lf 0£^ 2 - N 2^
For the semi-causal case we define (Fig. 2.5b)
R a = { (k i' k 2 } 1° i k l- r V k 2 = ° ; " M ct i k li M 6'° < k 2 - M 2 }
(2.29)
R b = {(£ 1 ,£ 2 )|0<^ ll N 3 ,£ 2 = 0;-N al £ 1 <N a ,Q<£ 2l N 2 }
41
HI
A
+)
•H
+j
3
a
e
o
u
en
■H
a;
o
c
<d
cr
w
4-1
3
a
4->
a M
O 0)
4J
(U H
4-1 <44
5 4J
C
144 13
O <T5
3
s cr
ft
^4 4J
n3 ^4
•H -H
Q 4-J
(T3
CM
en
•H
fa
42
C+
2
—
i
-4
1
jC
■P
•H
T3
Q)
-U
a,
g
o
o
CO
•H
<D
O
c
CD
a 1
CO
-u
o
a
4-> M
o
CD
O 4-1
<-<
w
3
e to
(C U
U I
03 g
•H CD
a co
IT)
CN
CT>
•H
fa
43
and in both cases b(0,0) = 1 is assumed. The inverse
z-transform is defined as
x(m,n) = j j / X(z,,z ) z™ 1 z^ dz,dz„ (2.30)
(2njr n „ 12 12 12
C l c 2
where c^ and c 2 are suitable closed contours in the z, and
z 2 planes.
The two-dimensional Fourier transform of a sequence
is defined as
DO)., 3u ?
X(e ,e *) = F[x(m,n)]
X(z 1 ,z 2 )
2
(z 1 ,z 2 ) eF
- j ( co-j m+cu^n)
£ x(m,n)e x * (2.32)
m r n=-°°
For any given sequence the set of (z,,z 2 ) for which the
z-transform converges is called the region of convergence.
Uniform convergence of the z-transform requires that the
sequence be absolutely summable, i.e.
oo
I |x(m,n) | iz 1 |" m |z 2 |" n < =o (2.32a)
m, n=-°°
44
Similarly, the Fourier transform X(e ,e ) converges
uniformly to a continuous function of co, and w_ if
I |x(m,n) | < oo . (2.32)
m,n=-°°
A power series of the form of (2.23) is known as a two-
dimensional Laurent series in the theory of complex functions
in two variables .
6 . Two-dimensional DFT
Let {x(m,n)} be a finite support sequence of the
size MxN. Then the two-dimensional discrete Fourier
transform (DFT) of (x(m,n)} is defined as
M ~ 1N ~ X -^ ( R k i + £ k 2>
X(k 1# k 2 ) =11 x(m,n)e M L N z (2.33)
m=0 n=0
for <_ k. < M-l and <_ k p < N-l. Similarly, the inverse
discrete Fourier transform (IDFT) is given by
45
X(m,n)
M-l
- I
MN
k,=0
N-l
I
k 2 =o
x(k lf k 2 )e
k, k 2
+:2TT( ir m + ir n)
(2.34)
Note that the DFT corresponds to sampling the Fourier trans-
form of MxN points, i.e.,
X(k x ,k 2 )
= x (e
Ua>i 300.
/e
k l k 2
a) l =27T lT' a) 2 =27T -N-
(2.35)
]W.
30).
More generally, if X(e ,e ) is the Fourier transform
of (x(m,n)} whose size is greater than MxN, then the
IDFT of {X(k 1 ,k )} will result in aliasing. In other words,
the IDFT of (X(k,,k 2 )} will yield one period of (x(m,n)}
where
x(m,n) =
m 1 ,n 1 =-<»
x (m+m,M,n+n,N)
(2.36)
46
C. STRUCTURES OF TWO-DIMENSIONAL FILTERS
1. One-dimensional Digital Filter Realization
One-dimensional digital filters, in general, can
be shown in the z-transform domain as follows:
y
1=0 Y(z i )
H < z l> = ~ = xuJt (2 ' 37)
1 - I a m Z l
m=l
or , in the time domain
M L
y(n) I a y(n-m) + J b x(n-l) (2.38)
m
m=l 1=0
The realization implementing equation (2.38) directly, is
known as the direct form 1 and shown in Figure 2.7.
Equation (2.37) can be rewritten in a slightly
different form by introducing a new variable, W(z 1 ) , such
that
W(z,) Y(z,)
H / z ) = L_ . . x s (2.39)
H(Z 1 J X(z x ) W z 1
The corresponding set of difference equations consists of:
47
I
N
J t— 1
2 t 3
I
H
as
kl
M
r
e
o
t^
+j
o
cu
u
•H
Q
CN
0^
•H
fa
o
5<
48
M
w(n) = x(.n) + I a m w(n-m) (2.40)
m=l
and
L
y(n) I b £ w(n-£) (2.41)
1=0
Equations (2.40) and (2.41) are realized in the
direct form 2 as shown in Figure 2.8 which can be visualized
as a cascade realization of two digital filters, realizing
the denominator and numerator polynomials, respectively.
Figure 2.8 can be redrawn as shown in Figure 2.9. The
resulting realization is called direct form 3 or canonic
form, since it represents a structure with the number of
delays equal to the order.
At this point, an additional realization will be
introduced, the direct form 4, which combines the charac-
teristics of the direct forms 1 and 2 having only two
entries in each summer, which corresponds directly to
hardware implementation and by realizing the numerator and
denominator polynomials separately in cascaded form
(Figure 2. 10) .
There are several other methods to realize a one-
dimensional M order digital filter. For example, the
cascade form of first and second order sections H(z,), as
shown in Figure 2.11, where
49
H(z X )
2 b £ Z l
£=0
-£
M
1 " I a z.
L ml
m=l
-m
X(n)
W(n)
I
«— 2
H
«, r
H
rf 1
^ 3 V
7-1
Fig. 2.8. Direct Form 2
50
H(z 1 )
£=0 * L
M
i - y a z"
L ml
m=0
-m
Y(-)
Fig. 2.9. Direct Form 3
(Canonic)
51
I
X t~o II
J ^
N
WH
l^H-
u-jK
JO
cv^~~
■ ^
a
CjW-
3
— 4
_ —
-±
KT
c
X
£
O
fa
-P
o
CD
S-i
•H
Q
o
H
•H
fa
52
H(z 1 ) = a Q n H i (z 1 ), * = ^ (2.42)
i=l
Each component H^(z^) can be realized in one of the above
outlined forms.
Partial fraction expansion methods applied to
equation (2.37) lead to a parallel form realization of
first and second order sections (Figure 2.12), i.e.,
k
H(z 1 ) - c+ J H i (z 1 ), k= [^+1] (2.43)
i=l
Other forms of realization include hybrid structures,
i.e., parallel-cascade forms, the transpose configuration,
which can be obtained for all previously outlined structures
by reversing the direction of signal flow and by interchanging
all branch and summing modes, wave-digital and continued
fraction expansion filters.
The direct forms discussed above and the series/parallel
arrangement of lower order sections are by far the most
often used ones in computer simulation and in digital
hardware.
2. Two-dimensional Digital Filter Realization
Two-dimensional Quarter plane digital filters, in
general form, can be shown as
53
N
ii
n"
x
s
u
fa
0)
o
en
u
CM
•H
fa
54
HW = C+2H } U,)
Fig. 2.12. Parallel Form
55
£ =0 £ =0 X Z
H(z,,z ) = — - (2.44)
M l M 2
-m.. -m
1 ~ I I a mm z, z
L *< m.iru 1 2
m l =0 m 2 = °
The numerator of the two-dimensional transfer function (2.44)
can be written as a one-variable polynomial in z~, weighted
by coefficients which are functions in z,, i.e.,
L 2
L l
I
( I
«,,=0
J..=0
L 2
' S V
£ 1 _il 2
z 2 (2.45)
" £ 2
:„ (2.46)
The same reasoning applied to the denominator polynomials
leads to
M M
-m, -itu
D ( z r z 2 } = 1 - ^ ( l ^2 Zl )Z2
m,=0 m.=0
(2.47)
M 2
= 1-1 (V. < z l»» Z 2 m (2 - 48)
m 2 =0
56
To obtain the direct form 3 representation, a new function,
W(z 1 z 2 ) is introduced:
W(z,,z 9 ) Y(z,,z 9 )
H(z i z 2> - xi^Tzf) ■ w(^fr (2 - 49)
where each factor can be written using equation (2.46) as
W(z x ,z 2 ) =
*2
£ 2 -0
/ N I 1 £ 2 (Z 1 ))Z 2 *<V*2 J
(2.50)
and equation (2.48) as
M 2
-m
2
Y(z 1 ,z 2 ) = W(z 1 ,z 2 ) + I (D M m (z 1 ))z 2 Z Y(z 1 z 2 )
m 2 =0 ^ 2
(2.51)
The realization of equations (2.50) and (2.51) is shown in
Figure 2.13 using compact notation, and Figure 2.14 where
N L l ( z l^ anc ^ D M m ^ z l^ are ^ m Pl emented f° r each (L-,^ 2 )
and (nunu) by one-dimensional direct form 4 realizations.
The general form of two-dimensional causal digital
filter (2.44) can be rewritten by introducing a new variable,
W(z-,z 2 ) or in the previous function, then,
W ( z , z « ) Y ( z , , z 9 )
H <V Z 2> " X( Zl ,z 2 ) • W(z',z 2 ) (2 " 52 '
57
e
o
•H
-p
+J
2
-u
o
(d
a
e
o
u
a
e
o
+J
u
aj
u
■H
Q
m
H
•H
58
Q
CO
e
S-l
Pn
-P
U
<U
u
•H
Q
■H
fa
59
Each factor can be written as :
L l L 2
-I -I
Y(l l'»2 J = I I W Z l 1 z 2 2 WCz r z 2 )
£ 1 =0 £ 2 =0
(2.53)
and
M M
-m -m«
W(z r z 2 ) =- XC V z 2 ) + I I a z x z 2 W(z r z 2 )
m 1 =0 m 2 =0 - 1 -
(m 1+ m 2 )^0) (2<54)
The equations (2.53) and (2.54) can be written as a differ-
ence equation as follows,
L l L 2
y(m,n) = I I b z z W(m-2. 1 ,n-£ 2 ) (2.55)
l l =0 £ 2 =0 1 2
and
M l M 2
w(m,n) = x(m,n) + I I a m m w ^ m " m 2 ' n ~ m 2 * (2.56)
m 1 =0 m 2 =0
m 1 ,m 2 ^
60
Without loss of generality, L. is chosen equal to M. for
all i. The structure realizing equation (2.55) and (2.56)
is shown in Figure 2.15. It is noted that the number of
delays in Figure 2.15 equals the order of H(z..,z ). The
structure is, therefore, by definition, canonic.
The direct form 4 (Fig. 2.10) realization of H(z.,z 2 )
is shown in Figures 2.16 and 2.17, using compact notation and
one-dimensional direct form 4 (Fig. 2.10) implementation of
each N « and D , respectively.
L 1^2 M l m 2
D. STABILITY
Two major problems in the design of a recursive filter
are approximation and stability. Since the output values
are used through feedback by the recursive algorithm, it is
possible for the output values to become arbitrarily large
independent of the input. A filter of this kind is said to
be unstable, which is an undesirable condition. Thus we
need to know what constraints to put on the recursive filter
coefficients, so 'that the filtering operation will be stable.
This problem is similar to the stability problem for the one-
dimensional case, except that the added dimension inherently
increases the complexity of the analysis. The major diffi-
culty is that the concept of zeros and poles does not hold
in multi-dimensional systems so that simple algebraic systems
extrapolation of one-dimensional results is not obvious.
• Stability Theorems :
Theorem 1 (Shank's theorem): A causal recursive filter
with the z-transform H(ZwZ 2 ) = A(z 1 , z 2 ) /B (z 1 , z 2 ) , where A
61
(Q*— I — hQ
m
(hh— I — ►©
M
Qk
4
©
t-j
TCNJ
c
c
U
Q
CN
CO
e
o
-U
u
CD
U
•H
Q
U")
CN
•H
fa
62
rc\j
Tsl
og
— @* —
M
rsj +
i
x
c
o
■H
+J
id
-u
c
+J
o
td
cu
2
o
Q
u
o
4J
o
0)
S-l
•H
Q
^o
CN
Cn
■H
63
rsj
H
%
J5'
Q
CM
6
S-i
o
fa
u
0)
CM
t7>
■iH
fa
and B are polynomials in z, and z 2 , is stable if and only
if there are no values of z, and z such that B(z ,z„) =
for all z. _> 1 and z« >^ 1. (In general we assume
B(Z 1' Z 2 } = b 00 + ^O 2 !" 1 + b 01 Z 2" 1 + ^l 2 !* 1 ^" 1 ••"
that is an expression in increasing powers of z. and z 9 .)
In other words, the theorem says that if there are any
values (real or complex) of z, and z 2 for which B(z,,z ? )
is zero for z. and z 2 simultaneously greater than or equal
to unity in magnitude, then the filter will be unstable.
It is much more difficult to determine the stability of
two-dimensional filters than of one-dimensional filters.
In the one-dimensional case, it is only necessary to locate
finite set of roots in the z-plane. In order to mechanize
the theorem, we must, in general, find the values of z. and
z_ for which B(z,,z~) = 0. There is no technique for locating
the zeros of a general two-dimensional polynomial in (z,,z_).
Shanks, Tretial and Justice, in their paper [11] proposed a
technique to apply theorem 1 to determine the stability of
a two-dimensional recursive filter, H(z,,z 2 ), by which the
unit disk in the z, plane must be mapped into the z^ plane
by solving the implicit two-variable denominator polynomial
of H(z,,z 2 ). A necessary and sufficient condition for
stability is determined if the map of the unit disk
d, = (z-,z_ _> 1) does not overlap the unit disk
d 2 E (z ? ,z > 1) in the z plane. This method requires an
infinite number of mappings and, therefore, can not be
applied exactly.
65
Huang [20] simplified Shank's stability theorem
considerably by stating the stability theorem as follows
Theorem 2: A causal filter with a z-transform
A(z, ,z.J
H(z w z n ) =
'l'"2' B(z 1 ,z 2 )
where A and B are two-variable polynomials in z,,z , is
stable if and only if:
1) the map of the unit circle 6d = (z , | z, | =1 ) ,
according to B(z ,z 2 ) = 0, lies outside of the unit disk
d_ = (z ,|z |_>1) in the z -plane, and
2) no point in the unit disk d.. = ( z ^ , | z , | _> 1 )
maps into the point z = by the relation B(z..,z 2 ) = 0.
In order to apply this theorem as a stability test, we
have to map the unit disk 5d, = (z , |z 1 |=1) into the z ?
plane according to B(z, ,z 2 ) = or z 2 = f (z,) | i i-, and
to see whether the contour in the z^ plane lies inside the
unit disk d 2 = (z- , | z 2 |<_1) . Also, it is necessary to solve
B(z,,0) = to see whether there are any roots with magnitude
greater than 1.
Although theorem 2 is much simpler than using Shank's
original theorem, the procedure is still infinite in the
sense that, in principle, we have to map the unit circle
5d, into the z ? -plane.
Ansell [21] , whose main contribution is to couple the use
of a Hermite test for checking stability with a series of
66
Sturm tests checking positivity, has reduced theorem 2.
Although Ans ell's results enable us to test stabilty in a
finite number of steps, it still is, unfortunately, very
tedious.
Make the change of variables:
and,
1 - Z 2
1 + z 2
and let
E(p, /P?)
H(ZwZ ) =
l'~2
fcp 1 ,p 2 ;
where E and F are polynomials in p, and p~. We can restate
Theorem 2 as follows.
Theorem 3 (Ansell's theorem): The causal recursive
filter H(z,,z 2 ) is stable if and only if:
1) for all real finite w , the complex polynomial
in p„, F(joj,p 2 ) has no zeros in ReP 2 > 0, and
2) the real polynomial in p-,, F(p,,l) has no
zeros in Rep, 0.
This theorem is essentially the same as Theorem 2, but
an advantage of Theorem 3 is that condition 1) can be
tested using standard techniques of circuit theory.
67
Anderson and Jury [22] modified Huang's method by out-
lining a procedure which replaces the bilinear transforma-
tions by the construction of a Schur-Cohn matrix and checking
for positivity of a set of self-inversive polynomials. It
also replaces the Hermite test by a Schur-Cohn matrix test
and requires a series of Sturm tests or equivalently , a
system of tests establishing the roots distribution of a
polynomial. These modifications represent a substantial
reduction in computations as compared to Huang's method. We
will discuss stability considerations which are related to
our design technique in Chapter IV and Chapter V in detail.
In this chapter, basic properties of two-dimensional
digital filters have been reviewed and summarized including
semi-causality to establish a fundamental theory of two-
dimensional digital filters.
68
III. MATHEMATICAL BACKGROUND FROM ANALYTIC PLANAR GEOMETRY
In this study, we propose using semi-infinite planes to
approximate a given frequency response of a two-dimensional
analog filters by extending one-dimensional semi-infinite
straight line approximation to two-dimensional cases. For
example, as we will see, the frequency response (dB vs. log-
frequency) of the separable canonic factor
H(ja) 1 ,j(u 2 ) = cl+jUiTi ) (l+ja) 2 T 2 )
can be approximated by a collection of planes. This is
discussed in Chapter IV in detail.
To prepare for the discussions in the next chapters,
it is worthwhile to review and develop useful equations
from analytic planar geometry which are going to be used
throughout this study with some examples of how to use them.
The properties of the analytic planar geometry can be
summarized as follows:
1. The general equations for a plane are given by
Ax+By+Cz+D=0 (3.1)
where A, B and C are the direction numbers of the normal
to the plane. The equation of a plane passing through the
point (xwy^z,) and perpendicular to the line which has
direction numbers of A, B and C is given by:
69
A(x-x 1 ) + B(y-y 1 ) + CU-z^ = (3.2)
2. The direction cosines of the normal to the plane
are given by-
cos a = — ^ = T-TTo (3.3a)
(A Z +B Z +C ) 1/Z
cos 6 = — ? ^ — TT7T (3.3b)
(A z +B +C ) ±/Z
cos y =
{A 2 + B 2 + C 2 ) l/2
where a, 8/ and y are the angles that the normal makes with
the x, y, z axes respectively.
The relationship between direction cosines is given by
2 2
Cos a + Cos 3 + Cos y = 1 (3.3d)
3. For a line joining points (x, ,y^/Z-^) and (x 2 ,y2> z 2)
the equation is given by
Cos a _ Cos 3 _ Cos y ,^ ^,
x 2 -x x y 2 -y! z 2 -z 1
4. The equation of a straight line with the direction
numbers a, b, c and passing through the point (x^y^z^ is
given by
x- Xl y-y 1 z- Zl
a ~ b = n — (3.5)
5. The angle 6, between two intersecting lines with
direction numbers a, b, c and a', b f , c' respectively is
given by:
Cos = Cos a Cos a 1 + Cos 6 Cos 3' + Cos y Cos y 1
(3.6a)
or
Cos 9 = aa ' + bb ' + cc ' - (3.6b)
Va 2 +b 2 + c 2 Va ,2 +b ,2 + c'2
6. For parallel lines, the relationship between
direction numbers is given by
a' b^ c* U * /J
7. For perpendicular lines:
aa' + bb' + cc ' = (3.8)
8. The direction numbers of the line of the intersection
between two given intersecting planes, with direction numbers
A, B, C and A' , B', C respectively, is given by
71
a =
b =
B
B* C
A
C A'
(3.9a)
(3.9b)
c =
B
A' B'
(3.9c)
The angle 6 between these two planes is given by
Cos
AA' + BB 1 + CC*
V
A 2 +B'
+ C V A '
(3.10)
2 2 2
9. The equation for a plane in terms of its slope can
be calculated as follows:
As seen from Fig. 3.1, the slope, m, of the normal
line is given by
m = tan
h
g
^hFW
(3.11)
m = slope of the plane
M = tan (90° - 9)
= - 2. = . 1
h m
(3.12)
72
normal to olane
normal to nlane
plane
Fig. 3.1. A General Plane Configuration
73
The equation of the plane is given by one of the
following:
A(x-x 1 ) +By + Cz = (3.13a)
Ax + B(y-y,) + Cz = (3.13b)
and
Ax + By + C(z-z 1 ) = (3.13c)
A*l = By ± = Cz 1 = -D (3.13d)
The intersection of the plane with the z-x plane
is given by
A(x -x 1 ) + Cz = (3.14)
and the slope, nu, of this line is
»i = a§ = -§ (3 - 14a)
Similarly, the intersection of the plane with the
z-y plane is given by
B(y -y, ) + Cz = (3.15)
74
and the slope of this line, nu , is
m 2 ~ " dy~ " " C (3.15a,
Substituting (3.15a) and (3.14a) into (3.13)
yields the following alternate equations for the plane.
-m(x-x-j) - m 2 y + z = (3.16a)
-m,x - m 2 (y-y 2 ) + z = (3.16b)
-itux - m 2 y + Cz - Zj) = (3.16c)
The direction cosines of the normal in terms of
m-, / m 2 are
-m,
Cos a = (3.17a)
Vm^
2 2
-m~
Cos 6 = (3.17b)
v^7
m 2 2 + l
Cos Y =
Vm-L 2 + m 2 2 +1
(3.17c)
Substituting (3.14a) and (3.15a) into (3.11) and
(3.12) yields
m = i- ( 3 - 18a)
m 2
75
and
M = -V m i +m 2 (3.18b)
2
The angle, 8, between two intersecting planes from
(3.10), (3.14a) and (3.15a), is given by
m,m, * + irunu ' + 1 ,- ,„,
Cos 9 = LJ: Li (3.19)
ymj 2 + m 2 2 + 1 Vmp+mp + 1
The direction numbers of the line of intersection,
from (3.9), (3.14a) and (3.15a) are given by
a = (m£ - m 2 ) CC (3.20a)
b = (m£ - m 1 ) CC* (3.20b)
c = (m|m 2 - m^ir^) CC * (3.20c)
The direction cosines of this intersection line
are given by
m 2 ~ m 2 (3.21a)
Cos a =
V2 2
(ml - nu) + (m| - m, ) + (m-jnu - mlm^)
2
m' - m,
Cos 6 = (3.21b)
V2 2
(mi -nu) + ^ m i ~ m i) + (m^m 2 -m^ir^)
76
m l m 2 " m 2 m l ,, ,. ,
Cos y = (3.21c)
V^ 7
2 2 2
m 2 ) + (m* -m,) + (m|m - m'm, )
10. When two planes are added in one dimension to form
a new plane, the following results. Given two planes
-m,x - nuy + z - z, =
-m,'x - mly + z 1 - z' =
Their sum in the z-direction with z" = z* + z is
given by
-(m 1 +m')x - (m 2 +m 2 ')y + z" - ^ + z£) =
The slopes of the intercepts with the zx and zy planes
m 1 " = -(m 1 + m|)
m 2 " = -(m 2 + m 2 ')
The slope of the normal is given by
2 , ,_ , .,2,-1/2
m
" = [ (m, +m[) + (m 2 +m 2 ) ]
(3.22
77
To aid in comprehension, two examples follow:
Example (1) : Given m, (slope in z,x plane), m,
(slope in x,y plane) , and x, (the point at which x-axis goes
through the plane). Find the equation of the plane.
The general plane equation (3.1) is given by
Ax + By + Cz + D =
When z =
B D
x A Y " A
m 3 = If = " I (given)
and
When y =
and
x, = - j (given)
A D
" " C x " C
■l - E - "I (given)
D
Z l C
78
If we substitute these into the general equation, the result is
or
A
Ax - Am-,y - — z - Ax, =
3 J m, 1
x - nuy - — z - x, =
3 m, 1
or
z = nux - nunuy - x i m i (3.23)
As shown here, we can derive the plane equation easily
from the given specification. After finding this equation,
we can find the canonic factor by simply letting
x = log go.
y = log a)-
i i 2
z = 10 log |T|
Alternately, we can go from the canonic factor to equations.
Thus equation (3.23) can be written as
2
10 log |T(joo 1 , joo 2 ) | = m 1 log oj 1 - m^ log oj 2 - 10 log K
(3.24)
where
10 log K = ^i x i
79
This result assumes that
i i2 w l
|T(ju>,,ju> 2 ) | = (3.24a;
Taking the logarithm of (3.24a) yields
2
10 log |T| = 10 p log u^ - 10 q log oj 2 - 10 log K
(3.25)
By equating (3.25) and (3.24), we get
10 p = m 1 (3.26a)
10 q = m 1 m 3 (3.26b]
It follows from (3.26a) and (3.26b) that
m 3 = q/p .
Thus the transfer function of (3.24a) has a logarithmic
characteristic as sketched in Fig. 3.2.
The previous result can be extended to the following
transfer function
(1 + j Wl ) p/2
T(j<D, / ja> 9 ) = —t-t? 7TT? (3.27)
1 Z K ±/Z (l + jw 2 ) q/
By substituting (3.26a) and (3.26c) into (3.27), then
80
Z=io 1^1T|
■s loovJ
*=: looi*/|
Fig. 3.2. The logarithmic characteristic of
P
T(jgo 1 , jo) 2 )
jj-
Kcu.
20 log |t| = m, log co, - m, m log aj_ - 10 log K
81
nu/20
(1 + jco,) L
T( >l' jw 2 ) ■ — Hu/20 (3 * 28)
K 1/Z (l + ja) 2 ) Z
When w, << 1 and oj 2 << 1 , region I
T(jo Jl /jco 2 ) = K~ 1/2
and
20 log |T| = - 10 log K
When co, << 1 and oj- >> 1* region II
T(j Ul ,ja> 2 ) = — ^-720
K 1/2 (ja J2 ) A
and
20 log |T| = - 10 log K - nu log cj.
When oj 2 << 1 and oj, >> 1, region III
111,720
(j«2.)
T(ju) 1 ,jo) 2 ) = -j/2
K
and
20 log |T| = - 10 log K + m, log
0).
82
When u>, >> 1 and oo- >> 1, region IV
nu/20
j OJ-
T(ju), ,jw~) = ,
1 2. . nu/20
K X / 2 H 1 2 '
and
20 log |T| = - 10 logK + m 1 log oj, - m, log to,.
The resulting planar approximation for the logarithmic
transfer function is shown in Figure 3.3
Example 2 ;
In this example, we derive planar equations from the
given transfer function.
Given a separable two-dimensional transfer function:
T(jw, , ju 9 ) = ■= (3.29)
(1 + JU) 1 T 1 ) P (1+ JU) 2 T 2 ) H
The magnitude square of the transfer function (3.29) is
given by
K 2
|T(ja) 1 ,ja) 2 )| = 2 2~p" 2 2~q < 3 - 3 °)
1 A (1 + o Ji Z t 1 Z ) P (1 + ^ 2 T 2 }
Take the logarithm of both sides of (3.30),
10 log |T| 2 = -10 p log(l+co 1 2 x 1 2 ) - 10 q log(l+u> 2 x 2 ) + 20 log K
(3.31)
when w,!, << 1 and co 2 x 2 << 1, the equation (3.31) becomes
83
A fcrtoUjITl
t 7^ta>W
IE
Br
3- l°j^
> TUlaj^l
Fig. 3.3. The logarithmic characteristic of
(1 + jo)) Mi/20
X Z K ±/Z (l + jco 2 ) il 2/ zu
20 log |t I = -10 log K + m, log u>, - m- log go-
84
2
10 log |T| = 20 log K = C
where C is some constant.
When w^t^ >> 1 and oj 2 t 2 >> 1,
or
10 log |T| = -20p log u^ - 20q log a, T
+ 20 log K
i i 2
10 log |T| = -20p log o^ - 20q log a> 2 - 20p log t.
Let
-20q log x 2 + 20 log
K
2
z = 10 log |t|
x = log to.
y = log a).
and we get the planar equation,
Z = -
20px - 20qy - 20p log t, - 20q log t 2 + 20 log K
or
85
p q
T T
(-20p)x + (-20q)y - z = 20 log (-± ? )
K
or
n^x + m 2 y - z - v
where
m^ = slope of the line (x,z) plane
= -20 p db/decade.
m 2 = slope of the line (y,z) plane
= -2 0q db/decade.
In a similar manner the regions (co, t, << 1;
CJ 2 T 2 >> "^ an< ^ ^l 1 ! >> ^' (J °2 T 2 <<c ^ can ^ e i nvest: '- < ? ate( i'
The planar approximation for the logarithmic transfer function
is sketched in Fig. 4.1 where q = p = 1.
In this chapter we have developed some useful planar
geometry formulas and shown how they can be used to determine
a planar approximation for the two-dimensional logarithmic
transfer function for separable cases. Some non-separable
factors are discussed in the next chapter in which experimental
verification of this approximation technique is also presented.
86
IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY DOMAIN
In Chapter III, we discussed some equations of planar geo-
metry and showed how to use them for approximating separable two-
dimensional transfer functions. In order to improve the scope
of the proposed design technique, we now consider the behavior
of other separable and non-separable factors insofar as
their planar geometry approximation is concerned. This will
provide insight into the frequency response characteristics
of two-dimensional systems and indicate what can be expected
by cascading these canonic factors, knowing their individual
behavior. The logarithmic characteristic of cascaded factor
is simply the sum of their individual logarithmic character-
istics. This approach parallels the semi-infinite line
approximation (dB vs. log w ) technique which has been used
so successfully in one-dimensional design.
A. SEPARABLE FACTORS
The following several separable canonic factors are now
considered.
(1) (1 + jaj 1 T 1 ) P (l + jt 2 u3 2 ) q (4. a)
i
(2) (1 + jaJ 1 T 1 ) P (l + ju) 2 T 2 ) q (l + joo 2 T 3 ) q (4.b)
(3) [1 + (ja Jl T 1 ) P ] a • (4.c)
where p, q and a are integers.
87
(1) First, as an example of the technique, consider
the following analog separable transfer function, the first
separable canonic factor with p = q = -1;
H(ju> ,ju> ) = ,-, 7 ■ w-i ■ • r (4.1)
1 2 (1 + joj 1 t 1 ) (1 + joo 2 t 2 )
The logarithm of magnitude-square of this function
is given by
*> *j _
10 log |H(co 1 ,oo 2 ) | = -10 logd+u^ t^ 2 ) - 10 log(l+co 2 t 2 )
(4.2)
Let us define the following for simplicity of notation.
z = 10 log |H(oo 1 ,o) 2 ) I = 20 log |H(o) 1 ,a) 2 ) | (dB)
(4.3a)
x = 20 log oo-l (4.3b)
y = 20 log w 2
(4.3c)
Consider the regions in the (gj, ,u) 2 ) plane (as seen
in Figure 4.2) as follows:
Region I, when oo,x, < 1 and ^ 2 x 2 < 1,
then z = db (4.4a)
Region II, when co^t -^ > 1 and co 2 t 2 < 1
z = -20 log (w 1 t 1 ) db (4.4b)
88
Region III, when ^^l < 1 and W 2 T 2 > 1
z = -20 log(co 2 T 2 ) db (4.4c)
and
Region IV, when ou^ > 1 and u t~ > 1
z =
= -20 log (oj 1 t 1 ) - 20 log (oj 2 t 2 ) (4.4d)
The above equations (4.4a), (4.4b), (4.4c) and (4.4d)
can be expressed as
Region I, z = (4.5a)
Region II, z = -x - 20 log t-, (4.5b)
Region III, z = -y - 20 log t 2 (4.5c)
Region IV, z = -x - y - 20 log t 2 - 20 log x,
(4.5d)
These equations describe four intersecting planes in
x,y,z space as shown in Figure 4.1.
Plane I is the horizontal plane at zero db. Plane II
and Plane III have slopes of -20 db/decade passing through
the corner lines as indicated in Figure 4.1. These lines
are analogous to the break point frequencies of one-dimensional
89
Z = 20 I o g |H( w r w 2 ) I
og w
slope —-20 d B
- 20dB
Fig. 4.1. Planar Approximation for
1
H (co, ,u>~)
1+J0J 1 T 1 ) (1 + J00 2 T 2 )
2 2
20 log I H (oi 1 ,co 2 ) | = -10 log (1+w, T,
-10 log (1+co 2 2 t 2 2 )
90
*2* 1 /^
Y = I O
g w
IT
W
I
X =
°g w.
Fig. 4.2. The regions of
H(j<VJ0) 2 ) " (1 + J03 1 T 1 ) (1 + J0) 2 T 2 )
2 2
20 log |H(o) 1 ,o) 2 ) | = -10 log (1+u^ Tj )
2
2 T 2
2 2
-10 log (1+oj„ t 9 )
in (log u), , log co 9 ) plane
91
log-modulus plots. Several properties of this plane are
summarized. Pertinent equations from analytic geometry are
given in Chapter III.
The actual frequency response of this separable
canonic factor with t, = t, = 1 is shown in Figure 4.3a
(dB vs. log w,, log oj 2 ) which fits the planar approximation ,
In Fig. 4.3b/ the contour plot of the frequency response is
given, which verifies the regions in log uj-, / log u>- plane
shown in Fig. 4.2. The logarithmic plots versus linear
frequency are given in Fig. 4.3c and Fig. 4.3d. It is
apparent that the approximation regions are well defined
when logarithmic frequency is used.
The frequency response of this separable canonic
factor with t, = 1, t- = 2 is given in Figures 4.4. These
results again confirm the logarithmic approximation.
(2) Next consider the canonic form
H(j Ml ,ja> 2 ) = (1 + j^tjj (1 + ju) 2 x 2 ) (1 + ja> 2 T 3 )
where we have two separable factors of oj 2 . The frequency
response of this factor with t, * 1, x 2 = 2 and x^ = 4 is
shown in Figures 4.5. As seen from this figure, we have
-40 dB/decade slope in the direction of u 2 , as expected
from the one-dimensional case. The regions are shown in
the contour plot with dotted line (Fig. 4.5c) .
92
In a similar manner we can analyze these regions as
follows: The logarithm of magnitude square of this function
is given by
10 log |H(u) 1 ,u) 2 ) | 2 = - 10 log(l + oj 1 2 t 1 2 ) - 10 log (1 + w 2 2 x 2 2 )
- 10 log(l + w 2 2 t 2 )
Consider the following regions in (to, , w 2 ) plane (as
seen in Fig. 4.5a) as follows:
Region I, when u-,t-, << 1 and co 2 x 2 < ^ 2 t 3 << 1
then z = dB.
Region II, when u,t, >> 1 and o) 2 t 2 < ^t-d <k 1
z = - 20 log (oj- l t 1 ) dB
Region II, when co-,t-, << 1 and oj 2 t 2 < ^2 T 3 <<; 1
z = -20 log (co 2 x 2 ) - 20 log(oj 2 x 3 ) dB
Region IV, when m,T, >> 1, ^2 T 3 > C0 2 T 2 >:> ^
z = -20 log u,t, - 20 log ^ 2 T 2 " 20 l0g W 2 T 3
93
Region V, when co^ << 1, ^^^ >> 1, ^ T >>o) 2 t 2 >> 1
Z = -20 logU 2 T 3 ) - 20 log(o, 2 T 2 )
Region VI, when ^^l >> ^' W 2 T 3 > T 2 (J °2 >:> 1
Z = - 20 log a) 1 T t - 20 log oo 2 t 2 - 20 log oj 3 t 3
where
Z = 20 log |H(u) 1 ,t 2 ) | .
(3) Consider the separable canonic factor now,
H(a) 1 ,o) 2 ) = [1 + (jw 1 T 1 ) P ] a
We can show that the following regions occur by
using the notation of (4.3a) and (4.3b) for simplicity:
94
Fig. 4.3a. The frequency response of
. 1
H(jo) 1 ,:o) 2 ) - (i +jWi ) (l+jo) 2 )
2
z = 2Q log|H(o) 1/ o) 2 ) | = -10 logd+u^ )
(dB vs. log bu,log o> 2 )
-10 log (l+w 2 )
95
X-SCflLE=l .OOE+00 UNITS INCH
T-SCRLE=I . OOE+00 UNITS INCH
•ig. 4.3b. The contour plot of H^,^) = (1+j ) ( i +j(0 )
z = 20 log |H| = -10 log (l+o^ 2 ) - 10 log (l+io 2 2 )
96
Fig. 4.3c.
The frequency response of
H(joJ 1 ,jco 2 )
(l+jai 1 ) (l+jo) 2 )
z = 20 log|H| = -10 log(l+oj 1 2 ) -10 logd+o^ 2 )
(dB vs. 03, ,0)-) (linear frequency)
97
10.
>",
a-,
3.
I
T
Fig. 4.3d. The contour plot of H(oo..,go_) = . . , . , . = , . r
r 1 2 (1+]U) 1 ) (l+:co 2 )
z = 20 log|H| = -10 logd+u^ 2 ) - 10 log(l+w 2 2 )
(linear frequency)
98
Fig. 4.4a. The frequency response of
H(do) 1 ,do) 2 )
(l+joa 1 ) (l+j2co 2 )
2 2
= 20 log|H| = -10 log(l+oj 1 ) - 10 log(l+4co 2 )
(dB vs. log co, , log co 2 )
99
X-SCRLE-1 .OOE+00 UNITS INCH.
T-SCflLE=l . OOE+00 UNITS INCH.
Fig. 4.4b.
The contour plot of the frequency response of
H(o3 1 w 2 ) =
(l+ja) 1 ) (l+j2u) 2 )
z = 20 log|H| = -10 log(l+co 1 2 ) - 10 log(l+4to 2 2 )
in (log to, , log to 2 ) plane
100
Fig. 4.4c. The frequency response of
1
H ( j co , , j co -, ) = ,-, , • -
J 1 2' (l+jco,) (1 + j2co 2 )
z = 10 log|H| = -10 logd+co, 2 ) - 10 log(l+4co 2 2
(dB vs. co, ,co_) (linear frequency ccale)
101
Fig. 4.4d. The contour plot of the frequency response of
H(jo) 1 ,joj 2 ) =
(1+w, ) (l+4oo 2 )
z = 20 log|H| = -10 log(l+oo 1 2 ) - 10 log(l+4aj 2 )
(dB vs. o) 1 ,o) 2 ) (linear frequency scale)
102
^ ' o g w 2
«H
"2
YE
nr
12
0)i*Tx
i
I
UH-i
-► !og w 1
Fig. 4.5a.
The region of
H(j(D 1 ,ju) 2 ) =
(l+jo) 1 T 1 ) (1+joo t ) (l+jco 2 x 3 )
in (log a), /log co„) plane
103
Fig. 4.5b. The frequency response of
1
H(jo) 1 ,jo) 2 ) - (i+j 0Jl ) (l+j2co 2 ) (l+34w 2 ;
z = 20 log|H| = -10 log(l+03 1 2 ) - 10 log(l+4w 2 )
2
-10 log(l+16u> 2 )
(dB vs. log o) lf log w 2 )
104
X-SCRLE=1 . OOE+00 UNITS INCH.
Y-SCRLE=I . OOE+00 UNITS INCH.
Fig. 4.5c. The contour plot of the frequency response of
H(ju lf ju 2 ) = ( i+j Ul ) (l+j2oo 2 ) (l+34co 2 ;
2 2
z = 20 log|H| = -10 log(l+a) 1 ) - 10 log (1+4^ )
2
-10 log(l+16oj 2 )
(log oo-L vs. log co 2 )
105
Fig. 4.5d.
z =
The frequency response of
H(joJ 1 ,jOi 2 ) = (1+j(Ui) (i +j 2co 2 ) (l+j4oo 2 )
20 log|H| = -10 log(l+ca 1 2 ) - 10 log(l + 2aj 2 2 )
-10 log(l+16oo 2 2 )
(dB vs. (ji-.,uiy) (linear frequency scale)
106
10,
0. -
X-
/«u
Fig. 4.5e. The contour plot of the frequency response of
. 1
H(D(o 1 ,30) 2 J - (1+jc0i ) (l+j2oa 2 ) (l + j4co 2 )
? 2
z = 20 log|H| = -10 log (1+w^) - 10 log(l + 4aj 1 )
-10 log (l+16oo 2 2 )
(linear frequency scale)
107
Region I, when go, > 1/t, ,
z = apx + 20ap log t-, dB
Region II, when to, < 1/t-,
z = dB
These separable canonic factors are identical to
the canonic factors of a one-dimensional transfer function.
They are just an extension to the two-dimensional case.
Therefore, they need not be discussed further.
B. NON-SEPARABLE FACTORS
The following non-separable canonic factors are now
considered.
(1) [1+ (jaj 1 x 1 ) P (joj 2 T 2 ) q ] a (4.d)
(2) [1 + (ju 1 T 1 + joo 2 T 2 ) P ] a (4.e)
where p, q and a are integers.
(1) First, consider the non-separable canonic
factor of the form:
H(oJ 1 ,aj 2 ) = [1 + (jaj- L T 1 ) P (joJ 2 T 2 ) ]
108
Region I, when
P q 1
1 2 T P T q
z = dB
Region II, when
P q
1 2 T p , q
T l T 2
z = apx + aqy + 20ap log t, + 20 ay log t- dB
The cornerline between the two regions is given by the
equation
or
p q p q
u l w 2 T l T 2
px + qy + 20p log t, + 20q log t 2 =
The cornerline is sketched in Figure 4.6.
The slope of the plane in region II is given by
-4 2
'p + q . Direction numbers of the normal are given
by -ap, -aq and -1.
109
»,= « 2 = I
(*)(*!*
^ x - lo g
I
T
M±)H^
y -lo g
Fig. 4.6. The regions of non-separable case in
log to, , log co^ plane
110
Depending upon the sign of p and q and the magnitude of
t, and x 2 several different orientations of the corner line
are possible.
Example (1)
Discuss this canonic form in the special case with
p = q = 1, a=+l, x, =0.5, and x 2 = 1. The transfer
function gets the following form:
H(o),,a) 2 ) = [1 + (.5jw,) ( ja) 2 ) ]
Take the logarithm of the magnitude square.
i 2 | i+2
z = +20 log|l - .5(jj,co 2 |
In region I, w,t, <<: ^ an< ^ aj ? T 2 <<: "*"
z = dB.
In region II, id,t, >> 1 and w 2 t 2 >> 1
z=x + y+20 log (0.5)
where
x = 20 log co 1 , y = 20 log cu 2 and z = 20 log |H(oj 1 ,u 2 )
The actual frequency response of this factor is shown in
Fig. 4.7 including a contour plot.
Ill
Fig. 4.7a.
The frequency response of
H( joo 1 , joo 2 ) = 1 + ( ja>,) ( jO. 5o3 2 )
z = 20 log|H| = 20 log (1 - 0.5oj,oj 2 )
(dB vs. log go,, log oo^)
112
X-SCRLE=l;OOE+00 UNITS INCH.
T-SCRLE=1..0£iE + 00 units inch.
Fig. 4.7b. The contour plot of the frequency response of
H(jco L ,ja> 2 ) = 1 + ( jco 1 ) (j0.5oj 2 )
z = 20 log|H| = 20 log (1 - O.Su^u^)
in (log co, , log ai 2 ) plane
113
When a = -1 we have a stability problem in this example
The denominator of this transfer function becomes zero
when
or
1 — . 5o)-j 0)2 =
o)-| cop — 2
This is the singularity of the transfer function in (w,,^)
plane- It can be called a singular line of the system. It
is a curve in the (ui,,to 2 ) plane to be compared with a pole
in a one-dimensional transfer function.
Example #2 ;
In this example, we discuss the same canonic form
with different coefficients, i.e.,
p = l, q = -1 , a = +l,
x, = 0.5, x 2 = 1
Then the transfer gets the form of
* 5uj l +1
H( Ul ,« 2 ) = (1 + ~— >
114
Taking the logarithm of its magnitude squared
|H(a), ,u 9 ) = ll + ±\ Z
. 5(jd,
z = +20 log ll +
w 2
. 5a)..
In Region I when - << 1
a) 2
z = dB.
• 5ai,
In Region II, when >> 1
w 9
z = 20x - 20y + 20 log (.5)
where x = log u^, y = log u 2 and z = 20 log |h(u>,/u> 2 ) |. The
boundary between these regions is given by the equation
. 5oi-| = oi~
We also have a stability problem for this canonic factor when
u>2 s 0. This is also a singularity in (u), ,w 2 ) plane in the
form of a straight line. This type of factor is unstable.
(2) Now, consider the non-separable canonic factor
H(o),,a> 2 ) = [1 + (ju),w, + ju^oO ]
We have two regions again, as follows:
Region I, when w,t, + w 2 t 2 < 1
z = dB
115
Region II, when u),x, + w-x- > 1
z = ap log (oa 1 x 1 + aj 2 T 2 ) •
The separation line between the two regions is
given by the equation:
U 1 T 1 + 0) 2 T 2 = 1
or
log (10, x, + o) 2 t 2 ) = . (4.6)
The locus of points of (4.6) is shown in Figure 4.8,
For region II, we have to think of two conditions,
as follows:
a) When w,x, > oo 2 x 2
z = apx + ap20 log x,
b) When oj,t, < ^ 2 x 2
z = aqy + ap20 log x 2 .
Thus, region II is split into two areas with the separation
line given by
W 1 T 1 = W 2 T 2 '
116
Fig. 4.8. Locus points of logC^T^ + w 2 x 2 ) -
117
Frequency, response of this non-separable transfer function
is shown in Fig. 4.9 with p - 1, a = -1, x = x = l.
C. CONCLUSIONS
In this chapter, we discussed the properties of the
separable and non-separable two-dimensional properties of
canonic factors in the analog domain (co,,^ plane). We can
conclude that these canonic factors can be approximated by
planes in the log magnitude versus log co, and log oj^ space.
As a design technique we can obtain any given frequency
specification by proper combinations of planar approximation
of canonic factors. It is obvious that only certain special
cases and symmetries can be obtained using the factors dis-
cussed in this thesis. For example, a wedge-shaped character-
istic in frequency domain could not be synthesized by means
of the factors we have considered.
The other result of this study is the stability problem
which arises in certain transfer functions, that is, the
singularity curve or line. For a separable case, we don't
have a stability problem; but in non-separable cases we may
have stability problems.
118
Fig. 4.9a.
The frequency response of
J 1 J 2 1 + joo, + jog.
z = 20 log|H| = 10 log 1 1+ (oj, +oj 2 )
(dB vs. log to, /log co 2 )
-, 2
119
X-SCRLE=1 . OOE + 00 UNITS INCH
T-SCflLE-1 . OOE+00 UNITS INCH
Fig. 4.9b.
The contour plot of the frequency response of
H(ju>, ,jw ) = . —-.
z = 20 log|H| = 10 log |l + (u^+u^) 2 !
in (log co, , log co-) plane
120
Fig. 4.9c. The frequency response of
1
H(3oa 1 ,:a) 2 ) =
1 + ju).
JO).
z = 20 log|H| = 10 log |l + (u^+a^) |
(dB vs. a) 1 ,oo 2 ) (linear frequency scale)
121
v - TWO-DIMENSIONAL DIGITAL FILTER DESIGN TECHNIQUE
A. INTRODUCTION
The design of two and multi-dimensional recursive digital
filters is difficult due to the absence of the Fundamental
Theorem of Algebra in two or more dimensions [15]. That is,
polynomials in two variables may not, in general, be factored
into lower order polynomials. For this reason, several design
procedures for two-dimensional recursive filters have been
proposed in the literature. We then consider the use of the
double bilinear z-transfer with the semi-infinite plane
approximation in the previous chapters. These are briefly
discussed in the following.
Shank et al suggested three design techniques in their
paper [11] .
1. The use of a series of stability-preserving transforma-
tions to produce a stable two-dimensional recursive discrete
filters from a stable one-dimensional continuous filter.
2. A space domain design technique in which the unit
sample response of the filter being designed is made to
approximate a given desired unit sample response.
3. A stabilization technique in which the denominator
polynomial is replaced by its double planar least-squares
inverse. Planar least-squares can be defined as follows:
We have some given array C. We would like to find
an array P such that C convolved with P approximates the
unit pulse array U. That is,
C * P « U
122
where the symbol * denotes two-dimensional convolution. In
general it will not be possible to make C * P exactly equal
to U. In actuality, C * P will be equal to some other
array, G. If we choose P such that the sum of the squares
of the elements of U-G is minimized, P is called a planar
least square inverse of C. Planar least square inverse of
P is called double planar least square inverse [11] . This
new polynomial is conjuctured to have no zeros on
-2 A , , , ,
u - {z^z^: | ZjJ ^_, |z 2 | >_ 1} and to have a frequency
response which closely approximates that of the original
transfer function.
The first technique is known as a filter design with
rotating axes and these filters are called rotated filters.
This design procedure consists in considering a transfer
function in the analog frequency domain, rotating the axes,
and applying the bilinear transformation to obtain a two-
dimensional recursive filter. In general, the starting
point is a separable filter and stability of the resulting
filter is not guaranteed even if the prototype filter is
stable. This technique was used by Costa and Venetsanopoulus
[23] to design low-pass filters having circular symmetry,
and also they proved that rotated filters are stable if the
angle of rotation is between 270° and 360°. A technique
similar in character but rotated axes in digital frequency
domain ((z,,z 2 ) domain) has been given by Chang and Aggarwal
[24]. Hirona and Aggarwal [25] have accomplished the rota-
tion through the introduction of rational powers of z^ and
123
z 2 , and the interpolation of input and output sequences.
They have proved that, in this technique, the resulting
filter is stable if the prototype filter is stable. Further,
the distortion, a consequence of the bilinear transformation,
is absent in the resulting filter when compared with Shank's
et al. 's filter.
The authors of [11] were unable to prove that the new
denominator polynomial produced by the third technique,
which is an extension of the one-dimensional case, was
stable in the two-dimensional case. Nonetheless, several
hundred filters have been designed using this technique
without a counter example appearing [18]; furthermore, Jury,
Kolavennu and Anderson [26] have shown it to be valid for
certain low order cases. However, Genin and Kamp [27] have
recently shown that counter-examples do, in fact, exist.
Spectral factorization is another technique which can be
used to produce a stable filter from an unstable one. By
spectral factorization we mean a complex map that carries a
stable rational transfer function into another stable transfer
function exhibiting a different frequency response, at the
same time maintaining some desirable characteristic [29] .
Pendergrass, Mitra and Jury [28] used stability-preserving
spectral transformations to modify the frequency response of
quarter-plane filters.
Mitra and Chakrabarti [29] extended this work. Read
and Treitel [30] suggested a two-dimensional Hilbert trans-
form to perform the spectral factorization, while Pister
124
[31] and Ekstrom and Woods [32] made use of the two-
dimensional cepstrum. Ekstrom and Woods also suggested ways
that the factors could be truncated and smoothed to preserve
stability. They also showed that the usual quarter-plane
filter is not the most general two-dimensional recursive
filter. It is instead possible to have a recursive filter
which in the limit has a unit sample response with support
on an asymmetric half-plane, and also they developed
stability conditions for this new class of filters.
B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN PROCEDURE BY
DOUBLE Z-TRANSFORM
1. Double Z-transform
In the one-dimensional recursive digital filter
design, the bilinear z-transform is often applied to an
analog filter transfer function. We proposed that this one-
dimensional design approach can be extended to the two-
dimensional case. Before going into the design procedure,
two-dimensional recursive filter design via the double
bilinear z-transform, it is better to prove that we can
use double z-transform for two-dimensional systems.
Each of the canonic factors, proposed in Chapter
IV, either separable or non-separable, can be converted
into discrete equations by the double z-transform, as
follows:
2 1 - Z I'
125
2 X " z 2
1 + z 2
The single z-transform is equivalent to trapezoidal
integration which can be shown as follows:
X(s) = S Y(s) (5.2)
.<*> = &J§i.
or
n
Y n ~ y n _ x = / x(t) dt
n-1
?(x + x .) (5.4)
2 n n-1
where T is the sampling period.
We write (5.4) in the z transform domain as
(1 - z 1 )Y(z) = |d + z X )X(z)
Hence
X(z) = §( 1 " \[ ) Y(z) (5.5)
1 + z
126
Comparing (5.5) with (5.2), it is seen that
2 1 - z" 1
S = J(i ^- T ) (5.6)
1 + 2
Equation (5.6) is known as a trapezoidal bilinear z -trans form,
We can drive the double z-transform in a similar
way, as follows:
or
Then
,2
W(x,y) = 33^ Z(x,y) (5.7a)
W(s,,s 2 ) = s i s 2 z ( s ]_' s 2^ (5.7b)
/ /d Z(x,y) = / /W(.x,y) dx, dy (5.8)
by a discrete approximation, the left-hand side of (5.8)
becomes (see Fig. 5.1)
/ /d 2 Z(x,y) = [Z(x n ,y n " z ( x n' y n-l )]
" ^V-l'V " Z(x n-l' Y n-l )]
-1 -1 . -1 -1,
(1 - z x - z 2 + z 1 z 2 )Z(.z 1 ,z 2 )
(5.9)
127
and the right-hand side of (5.8) is approximated as follows:
//W(x,y)dxdy = A^£ [w(Xn ,y n ) + W(x n ,y n-1 )
+ "(Vl'V + W(x n-l'*n-l )]
AXAVr, -1 , -1 , "I -1, TT/ V
= — -^- [1 + z, + z 2 + z i z 2 ]W(z 1 ,z 2 )
(5.10)
By equation (5.10) and (5.9) we have
W(z 1 ,z 2 ) 1 - z" 1 - z" 1 + z~ z 2
Z(z,,z ) . -1 -1 -1 -1 AxAy
1 2 1 + z, + z 2 + z, z 2 2
(5.11)
or
-1 i - 1
2 1 " Z l 2 X ~ Z 2
= ^ ( T—*^ ) ^ ( TT7^ ) (5 * 12)
1 + z, J 1 + z 2
By comparing (5.1a) and (5.2b) it is seen that
2 1 " z l
1 " Z 2
Q.E.D.
Thus, this is called a double z-transform. The
double z-transform seems to be a reasonable substitution
128
y(t)
x (O
— n-1 n
► t
(a) One-dimensional case
Z(*,y)
w(x,y)
>W(x,y)
(b) Two-dimensional case
Fig. 5-1- Two-dimensional extension of the
Bilinear Transformation
129
for each canonic factor. However, each factor must be
checked for stability, which is discussed in detail later in
this chapter. If each transformed factor is BIBO stable the
overall transfer function will be BIBO stable. But sometimes the
double z-transform may in certain cases lead to unstable
solutions [33] . Finally, it is noted that the double z-transform
results in a quarter-plane symmetry which is a severe limitation.
2 . Design Procedure
The two-dimensional recursive digital filter design
investigated in this thesis is the double bilinear z-transform.
Proceed as follows:
1. Determine a two-dimensional analog transfer
function,
A(s 1 ,s 9 )
hi«i'-2> = niTTifr (5 - 12>
where A(s,,s 2 ) and B(s,,s 2 ) are mutually prime two-dimensional
polynomials and B(s 1 ,s 2 ) satisfies
B(s 1 ,s 2 ) f
(5.13)
V(s 1 ,s 2 ) z r = {( Sl ,s 2 ): Re (s 1 ) >_0,Re(s 2 ) >_ }
in such a way that, the frequency response of G(s 1 ,s 2 )
meets the given specifications (for example as discussed
in Chapters III and IV) .
130
2. Apply the double bilinear z-transform
2 1 " z l 1
5 1 = A { ^T } (5.14a)
1 + z,
2 1 " z l
5 2 = -(- ij) (5.14b)
1 + z 2
where A is the sampling period, which is assumed the same in both
directions. This yields the two-dimensional digital transfer
function
2 1 -z, 1-2, C(z,,z )
H( f ( TI ) 'f ( =T )} = T ^V Z 2 ] ~ D(z z ) (5 ' 15)
A 1 + z 2 X A 1 + Z2 1 X 2 D( Z;l ,z 2 j
where C(z,,z 2 ) and D(z,,z 2 ) are mutually prime two-dimensional
polynomials in z, and z 2 . (5.15) will have the desired
frequency response with proper sampling time without
aliasing effect. Since
, 1 - z? 1
si = !< *r
1 A 1 + z.
maps the closed unit disc onto the closed right half
plane, one would expect (5.15) to satisfy the following
sufficient stability conditions [11], [20].
D(z 1 ,z 2 ) ?
131
\AW £ u2 - {(.z v z 2 ): | Zl | >_ 1,
z 2 | >_ 1} (5.16)
Unfortunately, however, the conditions (5.13) and (5.16)
as we will show later in this chapter, are not always
satisfied.
In summary, we have two conditions to be satisfied
for stability.
I. B(s lf s 2 ) ? o
\/(s 1 ,s 2 ) E r = {(s 1 ,s 2 ): R e s 1 >_ 0,R e s 2 > 0} (5.17)
II. D(2 1 ,z 2 ) ?
Y(z 1 ,z 2 ) £ u = {(z 1 ,z 2 ): Iz-J > l/|z 2 | >_ 1}
(5.18)
where V stands for "for each", and u is the closed unit bidisc
C. DIGITAL CANONIC FACTORS AND STABILITY CONSIDERATIONS
1 . Introduction
Certain properties of two variable polynomials and
rational functions which will be needed later will now
be discussed. It is well known that a two variable polynomial
is not, in general, factorable into first-order polynomials;
but, a two variable polynomial can be factored into irreduci-
ble factors which are themselves two variable polynomials
but which cannot be further factored. Of course a given
132
polynomial may itself be irreducible. These irreducible
polynomials are unique up to multiplicative constants [33] .
Two polynomials which have no irreducible factors in common
are said to be mutually prime. Consider two variable
rational functions, in general,
P (z, ,z 9
G(z, ,z ) =
l'"2' Q(z 1 ,z 2 )
where P(z^,z 2 ) and Q(z,,z 2 ) are mutually prime. Using
terminology of [34], a point (zj,zl) making Q(z|,z 2 ) =
but P(z',z') ^ is called a pole or a nonessential singu-
larity of the first kind (such a point is analogous to a
pole in the one variable case). A point (z£,z£) making
P(zV,z~) = Q(zV,z 2 ) = is called a nonessential singu-
larity of the second kind (such a point has no analogue in
the one variable case). Clearly, if (zi,zl) is a pole,
G(z,',z') = °°, but if (zV,z") is a nonessential singularity
of the second kind, then G(zV,z") is undefined.
2. Separable Factors
The most common canonic separable factor in Chapter
IV is
H(<jj,,u) ) — "7t — —j, \ 1 1 . a r (5.19a)
1 2 (1 + Joj, x,) (1 + 3oj 2 t 2 )
assuming x = x_ = 1. Substituting
jw-L - S 1 ,
133
ju> 2 = S 2
We get
H(s i' s 2» " (TTT[flTTip < 5 - 19b >
Substituting equations (5.1a and 5.1b) into (5.19b),
performs the double z-transformation. Thus,
8(1 + z" 1 ) (1 + z" 1 )
T(z r z 2 ) = i- £_ (5.20)
(a + z 1 ) (a + z 2 )
where
A 2
(A -2) 2
a =
a + 2
a - 2
A = sampling period
This separable canonic form obviously satisfies the condition
of (5.17). We should be careful about the second condition
of (5.18) .
When z, = z 2 = -1, we have a nonessential singularity
of the second kind if a = 1, but a = 1 never occurs.
Therefore, there is no nonessential singularity of the second
kind. Thus, this separable case is stable by satisfying
134
two essential necessities (5.17) and (5.18). The frequency
response of this separable case is plotted in Fig. 5.2, with
sampling time A = 0.08 which fits that of the analog frequency
response.
The double z-transform for equation (5.19a) is
A 2 (l + z 1 1 ) (l + z^ 1 )
T(Z 1 /Z 2 ) = Ta-2 T;l ) (A-2t 2 ) A+2 T;l ~ aT2~^ ~
( A^2T7 +Z 1 )( A^2TJ +Z 2 >
(5.21)
and this general case for the two-dimensional separable
factor is stable by satisfying the stability necessities
of (5.17) and (5.18)
The second analog canonic factor in Chapter IV is
H(j Ul rJ0. 2 ) (1 + jo) 1 T 1 ) (1 + JW 2 T 2 ) (1 + j w 2 T 3 ) '
When we perform the double bilinear transform for
this, we get
3 -1 -12
A J (1 + z/) (1 + z 9 x r
T(z n ,z„) =
1' 2' A+2t, , A+2T- , A+2t-. ,
(A-2x 1 )(A-2x 2 )(A-2x 3 )[^ + z- i ][ I=2 ^ + z 2 H I=n i+«; ]
(5.22)
The frequency response of (5.22) for t, = 1, t 2 = 2, t^ = 4
is shown in Fig. 5.3 which almost fits the analog case.
135
Fig. 5.2a. The frequency response of
T(z 1 ,z 2 ) =
where 3 =
e(l+z 1 " 1 ) (l+z 2 1 )
(a+z^ ) (a+z~ )
A 2 A+2
~ , a =
(A-2)"
(dB vs. log a), , log u}~)
A-2
= 0.08
136
Fig. 5.2b. The frequency response of
T(z 1 ,z 2 ) =
6(l+z 1 ~ 1 ) (l+z 2 1 )
(a+z 1 " 1 ) (a+z 2 -1 )
3 =
A'
(A-2)
2 '
a = |i| , A = 0.08
(dB vs. a). ,0)-) (linear frequency scale)
137
IT).
0. -,
*4
10.-,
Fig. 5.2c.
The contour plot of the frequency response of
3(l+z 1 " 1 ) (l+z 2 _1 )
T(z x ,z 2 )
6 =
2 '
(a+z 1 ~ ) (a+z 2 )
a = f±§, A= O.OB
where
(A-2)
(linear frequency scale)
138
We can summarize the m xn analog separable factor as
M N
H(ju> ,ju, ) = [ IT , .-*• ] [ n , , / ] (5.23)
After performing double z-transform, (5.23) becomes
M -1 N -1
M a (1 + z, ) N A (1 + z,)
T(z,,z ) = [ n — =i ][ n — | ]
(A-2t.) (-— ii+zT 1 ) 3 L (A-2x.) (_-^J-+z: 1 )
1 A-2x -1 3 A-2t ■ 2.
15.24)
Equations (5.23) and (5.24) satisfy the conditions
of (5.17) and (5.18). As a result, a separable filter
results in separable filters in digital domain and they are
stable.
3. Non-separable Factors
The first non-separable factor in Chapter IV is
4.d with a = -1, p = q = t, = t 2 = 1 as given below
1
H ( joj, t j oj ) - j— - ■
- 1 1 2 1 + 3oj-i](jo~
(5.25a)
or
«<V S 2> ■ TTT^ (5 - 25b:
139
Fig. 5.3a. The frequency response
A 3 (l+z 1 " 1 ) d+z^ 1 )
T(z 1/ z 2 ) =
(A-2) (A-4) (A-S, (f±| .z,- 1 ) (f±| .zf 1 ) <£±f +z 2 -h
where A = 0.0 8
(dB vs. log co, , log u^)
140
Fig. 5.3b. The frequency response of
T( 2l ,z 2 ) =
A 3 (l+z 1 " 1 ) U+z^ 1 )
(A-2) (A-4) (A-8) (|±| +z 1 " 1 ) (|±i +z x X ) (f±f + z 2 1 )
where A = 0.08
(dB vs. u.,ii)J
141
!S.u
Fig. 5.3c. The contour plot of the frequency response of
T(z 1 ,z 2 ) =
A 3 (l + z 1 " 1 ) (l+z^ 1 )
(A-2) (A-4) (A-8) (|±| +z 1 - 1 ) (|±J +a 1 " 1 ) <£±§ -hz^ 1 )
where A = 0.0 8
(linear frequency scale)
142
and by the double 2 -trans form, we get
A 2 (l +Z, 1 ) (1 + Z' 1 )
T(Z 1' Z 2 ) = — 2 =T TT T Tt-ZT (5 ' 26)
(A +4) (A -4) (z L + z l ) + (A +4)z 1 z X
In the analog domain (5.25b) does not satisfy the
condition (5.17). We have a singularity at
Ul u) 2 = 1
For the second condition, we have a nonessential
singularity of the second kind at z, = -1, z_ = -1. This
does not tell us whether it is stable or not, When we
check the frequency response, it is seen that this factor
is unstable.
The second non-separable canonic factor can be
derived from (4.d) with a = q = -1 and T-,=T 2 = las
follows.
H(ju)i ,jui ) = 3 (5.27a)
1 +
2 tu.
or
H < s i< s 2> - — 4r = sr^r < 5 - 27b >
1 + IT
s 2
143
Then, (5.27b) becomes after performing the double z-transform,
as follows:
(l + zT 1 ) (1 -z' 1 )
T(z ,z ) = _ 2 ( 5.28)
2U-« 1 r « 2 I )
This canonic factor does not satisfy either
stability condition. Hence it is unstable.
The other series of canonic factors can be derived
from (4.e) .
For example, with a = -1, t 1 = t- « p » 1, we
get
HtjoiwjwJ = -, , - 1 , ■ (5.29a;
-L ^ 1 + 3 co, +3 u^
or
H(S 1' S 2» = l + s^ + s 2 ( 5 - 29b>
or in the digital domain
A(l + z" 1 ) (1 + z" 1 )
T(z,,z ) = — —, _ (5.30)
(A+4) +A(z x x + z 2 ) + (A-4)z 1 z 2
This canonic factor satisfies the first condition but
it has a nonessential singularity of the second kind, at
z, = z~ = -1 and so we should check the frequency
144
response. We see that there is a singularity at
(u) 1 ,o) 2 ) = (it, it). Thus this non-separable canonic factor
is unstable.
4 . Conclusions
All of the canonic factors discussed in this
chapter are collected in Table I. By comparing the results
shown in this table, we can say that separable stable canonic
factors result in stable separable filters. As seen from
Table I, the separable transfer functions 1, 2 and 3 which
start with stable separable ones in the two-dimensional
continuous frequency domain end up as stable separable trans-
fer functions in the two-dimensional discrete domain. This
is also true for the mxn separable canonic factor 4.
Non-separable canonic factors have very different
characteristics. The unstable non-separable case ends up
with a non-separable case in the discrete domain, which is
also unstable. Although we could not prove this result
generally, we also could not find any counter example to
exist. The transfer functions of 5, and 6 start in the
unstable case, and end up with the unstable case in the
discrete domain. The transfer functions of 7 starts with
the stable case in the continuous frequency domain, ends up
with the unstable case in the discrete domain. Thus we
have to be very careful with the double z-transform since it
does not yield stable results every time. It may yield
unstable results even if it starts from the stable case.
145
UJ
_J
CO
<
>-
CO
<
co
C\J
M
CM
<f)
J)
LU
GQ
<
ea
H
LU
CD
<
CO
I m
w
I
NT
I
4-
CM
LU
co
<
uo
J-flL
If
<<
V
I
i
X
I
00
+
rO
LU
CD
<
CO
r
r ^-^ \
V
♦
*" s
t-r
1
<
~ •
*~<+
«- ^
.**•
C«i
*<
II
rJ
WO
*=*
(/T
+
i ! i
«■
3iaVd3d 3S
146
LU
CD
<
LU
CQ
<
UJ
—I
cn
<
ui
i «
4-
i —
m
+•
-<
fW
J"
«i
1 H
(N
*— «.
\
1 m
—
hi
-.jj
•
1 ^
^
r«
i .
i
fV
—
+
_
N
'-" /
I .4
> *
"7.
»
• -.
I*
ru
+
+
*W
iT)
CD
+
v/T
4»
1^
319VyVd3SNn
147
VI. CONCLUSIONS
A technique for the use of semi-infinite planes to
approximate the log magnitude versus log frequency charac-
teristics of two-dimensional filters is presented. This
technique is applied to quarter plane filters and works
well for separable transfer functions. Other non-
separable canonic (basic) transfer functions are also
studied in terms of their planar approximation.
The following separable and nonseparable canonic factors
have been studied.
(1) (1 + jo) 1 x 1 ) p (l + jw 2 T 2 ) q
(2) (1 + oo 1 x 1 ) P (l + jco 2 T 2 ) q (l + joj 2 T 3 ) q '
(3) [1 + (jco 1 t 1 ) P ] a
q,a
(4) [1 + (joo 1 T 1 )^(l + ja> 2 T 2 )*]
p,a
(5) [1 + (ja) 1 T 1 + jco 2 T 2 )^]
The properties and stability of these canonic factors
are investigated first in the continuous variable domain
and then after the double bilinear z-transform is applied
to extend them to two-dimensional recursive digital filters
These are all summarized in Table I in Chapter V. An
important result is that although the double bilinear
z-transform works well for separable factors, it does not
148
always lead to a stable digital realization starting with a
stable analog transfer function.
The following are suggested for further research:
1. Other types of canonic factors need to be considered.
For example, the general form:
H(j» lf jw 2 ) = [1 + (ja3 1 T 1 ) P +(ja J2 T 2 ) q +(jco 1 T 3 ) r (jw 2 T 4 ) t ] a
2. The work should be extended to non-causal (half-
plane) filters.
3. Transfer functions which have circular symmetry in
the discrete domain need to be investigated.
4. Another approach is to consider half-plane
recursive algorithms, by working backwards from the discrete
domain (z, ,zj to continuous domain (s,,s 2 ).
5. Stability of canonic factors needs to be investigated
in detail .
6. The use of predistortion should also be considered.
149
LIST OF REFERENCES
1. Oppenheim, A.V. and Schafer, R.W., Digital Signal
Processing , Prentice Hall, Englewood Cliffs, New Jersey,
1975.
2. Rabiner, L. and Gold, B. , Theory and Application of
Digital Signal Processing , Prentice Hall, Englewood
Cliffs, New Jersey, 1975.
3. Rabiner, L.R., et al., "Terminology in Digital Signal
Processing," IEEE Trans. Audio Electroacoust . , vol. AU-20,
p. 322-334, December 1972.
4. Stockham, T.G., "High Speed Convolution and Correlation,"
AFIPS Conference Proceedings, 28, pp. 229-233, 1966.
5. Hall, E.L., "A Comparison of Computations for Spatial
Frequency Filtering," Proc. IEEE (Special Issue on Digital
Picture Processing), Vol. 60, pp. 887-891, July 1972.
6. Huang, T.S., "Two-Dimensional Windows," IEEE Trans. Audio
Electroacoust. (Corresp.), Vol. AU-20, pp. 88-89, March
1972.
7. Hu, J.V. and Rabiner, L.R., "Design Techniques for Two-
Dimensional Digital Filters," IEEE Trans. Audio Electro-
acoust. (Special Issue on Digital Filtering), Vol. AU-20,
pp. 249-257, October 1972.
8. Arazi, E., "Two-Dimensional Digital Processing of One-
Dimensional Signal," IEEE Trans. Acoustics, Speech, and
Signal Processing, Vol. ASSP-22, pp. 81-86, April 1974.
9. Kamp, Y. and Thiran, J. P., "Maximally Flat Nonrecursive
Two-dimensional Digital Filters," IEEE Trans. Circuits and
Systems , Vol. CAS-21, pp. 437-449, May 1974.
10. Rabiner, L.R., McClellan, J.H., and Parks, T.W. , "FIR
Digital Filter Design Techniques Using Weighted Chebyshev
Approximation," Proc. IEEE , Vol. 63, pp. 595-609, April
1975.
11. Shanks, J.L., Treitel, S. and Justice, J.H., "Stability
and Synthesis of Two-Dimensional Recursive Filters,"
IEEE Trans. Audio Electroacoust. , Vol AU-20, pp. 115-128,
June 1972.
150
12. Mitra, S.K., Sagar, A.D., and Pendergrass, N.A.,
"Realizations of Two- Dimensional Recursive Digital
Filters," IEEE Trans. Circuits and Systems , Vol. CAS-22,
pp. 177-184, March 1975.
13. Green, B.W., and O'Handley, A.D., "Recent Developments
in Digital Image Processing at the Image Processing
Laboratory at the Jet Propulsion Laboratory," Proc. IEEE ,
Vol. 60, pp. 821-827, July 1972.
14. Treitel, S., Shanks, J.L. , and Frasier, C.W., "Some
Aspects of Far Filtering," Geophysics , Vol. 32,
pp. 789-800, October 1967.
15. Huang, T.S., Schreiber, F.W., and Tretiak, O.J., "Image
Processing," Proc. IEEE , Vol. 59, pp. 1586-1609,
November 1971.
16. Hunt, B.R., "Digital Image Processing," Proc. IEEE ,
Vol. 63, pp. 693-708, April 1975.
17. Wood, L.C. and Treitel, S., "Seismic Signal Processing,"
Proc. IEEE , Vol. 63, pp. 649-661, April 1975.
18. Mersereau, Russell M. and Dudgeon, D.E., "Two-
Dimensional Digital Filtering," Proc. IEEE , Vol. 63,
pp. 610-623, April 1975.
19. Chong, H. and Aggarwal , J.K., "Fundamentals of Two-
Dimensional Recursive Filtering," Digital Signal
Processing , (Edited by Aggarwal, J.K.), pp. 211-260,
Pt. Lobos Press, CA. , 1979.
20. Huang, T.S., "Stability of Two-Dimensional Recursive
Filters," IEEE Trans. Audio Electroacoust . , Vol. AU-20,
pp. 158-163, June 1972.
21. Ansell, H.G., "On Certain Two-Variable Generalizations
of Circuit Theory, with Applications to Network of
Transmission Lines and Lumped Reactances," IEEE Trans.
Circuit Theory , VoL CT-11, pp. 214-223, June 1964.
22. Anderson, B.D.O. and Jury, E.I., "Stability Test for
Two-Dimensional Recursive Filters," IEEE Trans. Audio
Electroacoust . , Vol. AU-21, pp. 366-372, August 1973.
23. Costa, J.M., and Venetsanopoulos , A., "Design of
Circularly Symmetric Two-Dimensional Recursive Filters, 1
IEEE Trans. Acoustics, Speech and Signal Processing ,
Vol. ASSP-22, pp. 432-443, December 1974.
151
24. Chang, H. and Aggarwal, J.K., "Design of Two-Dimensional
Recursive Filters by Interpolation/' IEEE Trans. Circuits
and Systems , Vol. CAS-24, pp. 281-291, June 1977.
25. Hirono, K. and Aggarwal, J.K., "Design of Two-
Dimensional Recursive Digital Filters," IEEE Trans.
Circuits and Systems , Vol. CAS-25, pp. 1066-1076,
December 1978.
26. Jury, I., Kolavennu, V.R., and Anderson, B.D.O.,
"Stabilization of Certain Two-Dimensional Recursive
Digital Filters," Proc. IEEE , Vol. 65, pp. 887-891,
June 1977.
27. Genin, Y.V. and Kamp , Y.G., "Two-Dimensional Stability
and Orthogonal Polynomials on the Hypercircle, " Proc.
IEEE , Vol. 65, pp. 873-881, June 1977.
28. Pendergrass, N.A., Mitra, S.K., and Jury, E.I., "Spectral
Transformations for Two-Dimensional Digital Filters",
IEEE Trans. Circuits and Systems , Vol. CAS-23, pp. 26-
35, January 197 6.
29. Chakrabarti, S. and Mitra, S.K., "Design of Two-Dimensional
Digital Filters via Spectral Transformations", Proc.
IEEE , Vol. 65, pp. 905-914, June 1977.
30. Read, Randol. R. and Treitel , S., "The Stabilization of
Two-Dimensional Recursive Filters via the Discrete
Hilbert Transform" , IEEE Trans. Geoscience Elect .
Vol. GE-11, pp. 153-160. July 1973.
31. Pister, P., "Stability Criterion for Recursive Filters",
IBM Journal of Research and Development , Vol. 18,
pp. 59-71, January 1974.
32. Ekstrom, P.M. and Woods, J.W., "Two-Dimensional Spectral
Factorization with Applications in Recursive Digital
Filtering", IEEE Trans. Acoustics, Speech and Signal
Processing , Vol. ASSP-24, Pp. 115-128, April 1976.
33. Goodman, D., "Some Difficulties with the Double
Bilinear Transformation in Two-Dimensional Recursive
Filter Design", Proc. IEEE , Vol. 66, pp. 796-797,
July 1978.
34. Goodman, D., "Some Stability Properties of Two-Dimensional
Linear Swift-Invariant Digital Filters", IEEE Trans .
Circuits and Systems , Vol. CAS-24, pp. 201-208, April
1977.
152
INITIAL DISTRIBUTION LIST
No. Copies
1. Defense Documentation Center 2
Cameron Station
Alexandria, Virginia 22314
2. Library, Code 014 2 2
Naval Postgraduate School
Monterey, California 93940
3. Department Chairman, Code 62 1
Department of Electrical Engineering
Naval Postgraduate School
Monterey, California 93940
4. Professor S. R. Parker, Code 62Px 10
Department of Electrical Engineering
Naval Postgraduate School
Monterey, California 93940
5. Deniz Kuvvetleri K-ligi 1
Personel Sube Bsk-ligi
ANKARA/TURKEY
6. Deniz Harb Okulu K-ligi 1
Heybeliada
ISTANBUL/TURKEY
7. Coskun Zumrutkaya 3
Yuzbasiler cad. No=20
Degirmendere
KOCAELI/TURKEY
153
I<2 AU^Qo
25787.
Thesis
Z8 Zumrutkaya ] 8381 6
Geometric design tech-
nique for two-dimensional
analog and digital filters.
c.l
l,c A^jdO
- 25 7 8 7^
Thesis
Z8 Zumrutkaya
c.l Geometric design tech-
nique for two-dimensional
analog and digital filters
183816
thesZ8
Geometric design technique for two-dimen
3 2768 001 07029 5
DUDLEY KNOX LIBRARY
TVNews
OpenLibrary