Description
Digital Signal Processing
1. Consider a medical CT (computerized tomography) scan which results in a set of twodimensional lateral image ìslicesî taken along a longitudinal axis. With reference to
(x; y; z) coordinates, x; y are the lateral coordinates, and z the longitudinal coordinate.
For simplicity assume the spacing between the slices has unit length.
We will model the resulting 3-D image as a signal g (~r; n) with coordinate system
~r 2 R
2
(the continuous lateral coordinates), and n 2 Z (the discrete longitudinal
coordinate). In what follows, do not separate ~r into x; y coordinates, or join it with n
to form a 3-D vector: keep ~r,n separate. Let the corresponding frequency coordinates
be denoted
~kxy; kz
, where ~kxy = (kx; ky) is kept separate from kz. Do not explicitly
write x; y; kx; ky in any of the following,
(a) Specify the domains of ~kxy and kz, respectively.
(b) Write the Fourier and inverse Fourier transforms.
(c) State Parsevalís theorem.
(d) Write the formulas for convolution in the spatial and wavenumber domains.
(e) Prove that convolution in the spatial domain leads to multiplication in the wavenumber domain.
(f) Let h (~r; n) denote the impulse response of an LTI system applied to such a 3-D
object. Write the condition for stability as an explicit mathematical expression
(i.e., do not simply write h 2 L
1
).
2. McClellan Transform Design Method for Multidimensional FIR Filters
The problem presents a technique for designing multidimensional FIR Ölters based on
1-D ìprototypesî. The approach described here can be used for 2-D, 3-D (or higher!)
Ölters, but we will stay with a 2-D example. The idea is to start with a 1-D FIR Ölter
“prototype”, H (z), and make a substitution z F (z1; z2) to obtain a 2-D Ölter:
G (z1; z2) = H (F (z1; z2))
The McClellan transform proposes the form of F () that we chose. Usually, the target
Ölter G should be zero-phase with real coe¢ cients, so we chose H to be zero-phase
with real coe¢ cients. The function F is also zero-phase with real coe¢ cients, and
it is usually designed to map “DC” to “DC” and high frequencies (say the points
(!1; !2) = (; )) to ! = . Thus, for example, if H is lowpass, so is G, if
H is bandpass, so is G. In addition, F is often chosen to achieve certain desirable
symmetries. For example, if F is nearly isotropic, so is G.
If we start with a zero-phase FIR Ölter with real coe¢ cients, it can be expressed (in
the frequency domain) as a polynomial1
in cos !. SpeciÖcally, if h [] has support
fN n Ng, h [n] = h [n], then we can write:
H (!) = h0 + 2h1 cos ! + + 2hN cos N!
= h0 + 2X
N
m=1
hmTm (cos !)
where Tm () is the Chebyshev polynomial of degree m. Note that H (!) = H (!) is
real. With the mapping cos ! ! (z + z
1
) =2 we get:
H (z) = h0 + 2X
N
m=1
hmTm
z + z
1
2
(1)
Next, we make a substitution:
cos ! = F (!1; !2)
where F (!1; !2) is a low degree polynomial in cos !1; cos !2;sin !1;sin !2. This results
in F (z1; z2) being a 2-D FIR Ölter (i.e., involving relatively low order positive and
negative powers of z1; z2). For example, cos !1 !
z1+z
1
1
2
, and sin !1 !
z1z
1
1
2j
, and
similarly for cos !2;sin !2. Note that symmetries in F will ensure its Önal form will
involve only real coe¢ cients, however. This gives us the following formulas for the 2-D
Ölter G in the frequency and transform domains:
G (!1; !2) = h0 + 2X
N
m=1
hmTm (F (!1; !2))
G (z1; z2) = h0 + 2X
N
m=1
hmTm (F (z1; z2))
The function F (!1; !2) is real and must satisfy:
1 F (!1; !2) 1 for (!1; !2) 2 [; ]
2
Then, under the mapping cos ! = F (!1; !2), for each (!1; !2) 2 [; ]
2
, there is a
unique ! 2 [0; ] to which it corresponds. Let us denote that as !(!1;!2)
. Then the
value of the frequency response G (!1; !2) at (!1; !2) is the value of the frequency
response H at !(!1;!2)
.
So here is some help with MATLAB. If x is a symbolic variable, T=chebyshevT(1:N,x)
will create a vector of symbolic Chebyshev polynomials Tm (x) for 1 m N.
Also, I have emailed you a MATLAB function sym2z2d that will map a symbolic
2D FIR transfer function H (z1; z2) into the impulse response matrix h [n1; n2]; the
script testsym2z2d gives an example of its use.
Let us take a speciÖc example. Take:
F (!1; !2) = 1
2
(1 + cos !1 + cos !2 + cos !1 cos !2)
1Since cos (m!) = Tm (cos !) where Tm is the Chebyshev polynomial of degree m, the given formula is
indeed a polynomi
(a) Check the following:
1. F has 8-fold symmetry, in that F (!1; !2) = F (!2; !1) (i.e., replacing
either !1 or !2 by its negative, or swapping them, does not change F). This
is what leads to ìapproximatelyîisotropic behavior.
2. Either !1 = or !2 = corresponds to ! = .
3. F (!1; 0) = cos !1, which means G (!1; 0) = H (!1), i.e., the ìsliceîof G along
the !1 axis is the same as the 1-D Ölter H, and similarly G (0; !2) = H (!2).
(b) Compute F (!1; !2) on a 128×128 grid over [; ]
2
, and graph a contour plot of
it. Check that its values are bounded between 1 and 1.
(c) Find F (z1; z2) as a 2-D polynomial.
(d) If the 1-D Ölter has support fN n Ng, what is the support of the 2-D Ölter
designed using this method?
(e) Write a MATLAB function that will do the following: given a zero-phase real
1-D FIR Ölter h [n] with support fN n Ng, it will use THIS SPECIFIC
F (!1; !2) (donít try to write a general purpose McClellan transform function!)
and:
1. Compute and graph the magnitude response of H (!) over [; ] (128 points
equally spaced is Öne). Use a linear (not decibel) scale.
2. Compute and graph the magnitude response of G (!1; !2) over [; ]
2
. (a
128 128 grid). Also on a linear scale.
3. Use the symbolic toolbox to Önd G (z1; z2), and compute g [n1; n2] (represented as a matrix). These two quantities should be returned by the function.
(f) Apply your function to the following. Let h be a Hamming window of length 7,
but centered so as to have support f3 n 3g. This will serve as a lowpass
prototype Ölter.
(g) Now apply your function to the following bandpass Ölter: Önd the impulse response h0 [n] of an ideal bandpass Ölter with H0 (!) = 1 if
3
2
3
, and 0 otherwise
(it will have the property h0 [n] = h0 [n]); I expect you to use the exact formula
(you can look it up); truncate it with your Hamming window function to obtain
the length 7 prototype Ölter hBP
