- Email: [email protected]
A Design Technique for 2-D Linear Phase Frequency Sampling Filters with Fourfold Symmetry
A Design Technique for 2-D Linear Phase Frequency Sampling Filters with Fourfold Symmetry Peter A. Stubberud University of Nevada, Las Vegas
Abstract In this chapter, system functions are developed for two dimensional (2-D) frequency sampling filters that have real impulse responses and linear phase and for 2-D frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. Under certain conditions, these frequency sampling filters can implement narrowband 2-D linear phase filters and narrowband 2-D linear phase filters with fourfold symmetry much more efficiently than direct convolution implementations. Also, a technique for determining optimal frequency sampling filter coefficients is developed for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. This technique approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response.
I. Introduction Some two dimensional (2-D) signal processing systems, including image processing systems, require linear phase or zero phase filters. A 2-D linear phase or zero phase filter implemented by direct convolution uses the filter's impulse response as coefficients. If a 2-D linear phase filter has a region of support, R N, where CONTROL AND DYNAMIC SYSTEMS, VOL. 77 Copyright 9 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.
117
118
PETER A. STUBBERUD
R N = {(nl,n2)'0_< n I _
and N 1, N 2, n 1 and n 2 are members of the set of integers ( N 1, N 2, n 1, n 2 ~ / ) , then the filter's impulse response has the form h(nl,n2) = h(N 1 -1-nl,N 2 -l-n2) [ 1]. If the filter is implemented using direct convolution, approximately N 1 N 2 / 2 multiplies are required to compute each output sample. To reduce the
number of multiplies required by a direct convolution implementation of a linear phase filter, fourfold symmetry conditions are often imposed on linear phase filters. If a 2-D fourfold symmetric linear phase filter has support over the region, R N, then the filter's impulse response has the form h ( n I ,n2) = h ( n l , N 2 - 1 - n2) = h ( N 1 - 1 - n 1, n 2 ) = h(N 1 -1-n
1,N 2 -1-n
2)
[2], and a direct convolution implementation of the filter requires approximately N 1 N 2 / 4 multiplies per output sample. Unlike direct convolution implementations which use the filter's impulse response as coefficients in the filter's implementation, frequency sampling filters use frequency samples, which are specific values from the filter's frequency response, as coefficients in the filter's implementation. The frequency sampling filters discussed in this chapter interpolate a frequency response through a set of (NIN2) frequency samples. When using a frequency sampling filter to implement a frequency selective filter, the frequency samples that lie in the filter's stopbands can be set to zero and do not have to be implemented. Thus, only the non-zero frequency samples which lie in the illter's passbands and transition bands are required in the filter's implementation. As a filter's passbands or transition bands narrow or its stopband specifications become more stringent, the filter's region of support increases which implies that the values of N 1 and N 2 increase. Because the number of multiplies per output sample of a direct convolution filter is approximately proportional to (N1N2), even a small increase in the filter's specifications can substantially increase a direct convolution filter's computational requirements. However, as a frequency sampling filter's passbands and transition bands narrow and its stopbands increase, the number of nonzero frequency samples decrease; and thus, the number of coefficients required in the filter's imple-
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
119
mentation decreases. Thus, frequency sampling filters have the potential to implement narrowband frequency selective filters much more efficiently than direct convolution filters. In this chapter, the 2-D Type 1-1 frequency sampling filter system function described in reference [3] is further developed for frequency sampling filters that have real impulse responses and linear phase and for frequency sampiing filters that have real impulse responses, linear phase and fourfold symmetry. This chapter also introduces three new types, Type 1-2, Type 2-1 and Type 2-2, of frequency sampling filter system functions. These system functions are also developed for frequency sampling filters that have real impulse responses and linear phase and for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. The resulting system functions for frequency sampling filters that have a real impulse responses and linear phase interpolate a frequency response through a set of approximately
N1N2/2independent samples, and the resulting system functions for frequency sampling filters that have a real impulse responses, linear phase and fourfold symmetry interpolate a frequency response through a set of approximately
N1N2/4independent samples.
As a result these system functions are computationally more efficient for frequency sampling filters that have real impulse responses, linear phase and for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry than the system function described in reference [3]. Although frequency sampling filters interpolate a frequency response through a set of frequency samples, the frequency response may not be well behaved between samples. References [3; 4] describe two of the design methods currently used to control the interpolation errors between frequency samples and approximate a desired frequency response. In this chapter, a design technique is developed that controls the in-
terpolation errors and approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response. This results in a constrained optimization problem which can be solved by using the Lagrange multiplier optimization method.
120
PETER A. STUBBERUD
II. 2-D Frequency Sampling Filters A. Type 1-1 Frequency Sampling Filters Consider a filter that has an impulse response h(n) where n - (nl,n2), a region of support R N where
RN = { ( n l , n 2 ) ' O < n 1 < N l - l , 0 _ < n 2 < N 2 - 1 } and a frequency response, Hro(to), where
n l,n 2 ~ 1,
0)= (091,0)2) ~ RK~ = {(0)1,0)2)'0 < 0)1 < 2g,0 <_0)2 < 21r} C~ 0)2 e the set of real numbers (9t). Let Hk(k) where k - ( k l , k 2 ) ~ R K - { ( k l , k 2 ) ' O _ < k 1 <_N l - l , 0 < k
2<_N 2 - 1 }
kl,k 2 ~ I
represent a set of (NlN 2) frequency samples taken from the filter's frequency response, H (to) for to e R For a Type l-1 frequency sampling filter, the (o .(_2. set, Hk(k) for k ~ R K, of frequency samples are selected such that 2rt 2rt U k (k) = Ho~ (0)1,0)2)10~1=.~l_lk 1,0~2=.~_2k2 where Hk(k ) can be written as
H k (k) - [ H k (k)leJO(k) where k ~ R K, j = ~ - l and 0(k) is the phase of Hk(k ). Thus, for a Type 1-1 frequency sampling filter, the set, Hk(k) for k ~ R K, of frequency samples corresponds to the Discrete Fourier Transform (DFT) coefficients. Also, the impulse response, h(n) for n e R N, which interpolates a frequency response through the set of frequency samples, Hk(k) for k ~ R K, can be determined from the inverse discrete Fourier transform (IDFT),
~,~'~H k (k)e j(2rt/NI )nlkl e j(2rt/N 2 )n2k 2 , h(n) = ~ 1 NIN2 k~RK
(1)
and a set, Hk(k) for k ~ R K, of frequency samples can be determined from the filter's impulse response by the discrete Fourier transform (DFT),
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS Hk(k) -
121
Z h ( n ) e - J ( 2 g / N 1 )nlkl e-J(2rr'/N2 )n2k2 " neR N
If we let H(z 1,z2) represent the z transform of h(n), then the system function,
H(z 1,z2), of the filter that interpolates a frequency response through the set, Hk(k) for k e R K, of frequency samples is N 1-1N2-1 H(Zl, Z2 ) = Z Z h(nl'n2 )z? nl z2 n2 "
(2)
n 1=0 n2 =0 By substituting Equation (1) into Equation (2), interchanging the order of the summations and performing the summation over the n 1 and n 2 indices,
H(z 1,z2) can be written as 1 - Zl N1 1 - z2 N2 H(Zl,Z2) = ~ x
N1
N2
N 1-1N2-1 (3)
Equation (3), originally developed in [3], describes the general form of a Type 1-1 frequency sampling filter, and has the form of an interpolation formula. The complex function, H(Zl,Z2), interpolates a polynomial through the frequency samples in the set, Hk(k) for k ~ R K, so that n ( z l , z2)] Zl =eJ(2rc / N1)kl ,z2 =eJ(2rc / N2 )k2 = n k (kl' k2 )" Thus, as desired, the frequency sampling filter's frequency response passes through the (N1N2) frequency samples in the set Hk(k) for k ~ R K. The frequency sampling filter in Equation (3) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A filter with a real impulse response has a DFT of the form
IHk(kl,k2)l-]Hk(N1-
kl,N2 -
k2)l
0(k 1, k2 ) = - 0 ( N 1 - k 1, N 2 - k 2 ). It can be shown[ 1] that the phase of a linear phase FIR filter with a region of support R N is
122
PETER A. STUBBERUD
arg[Hco ((01 (02)] = -(01
(N,-1/ (N2-1/
'
- (02
2
Therefore, the phase of the frequency samples is
O(kl'k2) =- N 1
2
--~2 k2
2
9
(4)
(N2-1 /" 2
Substituting these constraints into Equation (3), the filter's system function can be written as
Hz(Zl,Z2) =
1- Z?NI 1- Z2N2 I
Hk(O'O)
~k 1 ~l ~-l~k'21Hk~kl,0~lcos(~/tlZl' t
l--'II--2COS kl/zll'N1 +ZI2]I--z;I +2;
~2
["Trk2 - z;l ) (-1~2 ~l.~(0,~l cos~w)(1
k2=l(1-zll) N 1-1
M2
kl=l
k2=l
+ZZ
1 - 2 c o s -~22k2 z 2 + z 2 2
(_l)kl +k2 21Hk
cos~?
+ -~2 ~k2
z-l z" I" - cos NI
N2
+
(5)
where M1 - (N l- 1)/2, M 2 - (N 2-1 )/2, and N1 and N 2 are odd. The frequency sampling filter in Equation (5) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 1-1 frequency sampling filter with fourfold symmetry has a DFT of the form
IHk (kl' k2)l - I H k
(kl, N2 - k2)l - I H k (N1 - kl, k2 )1 - I H k (N1 - kl' N2 - k2)l"
Substituting these constraints into Equation (5), the filter's system function can be written as
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
Hz(Zl,Z2)= 1- ZlN1 1- z2N2 I
.,,,,
123
Hk(O'O)
.,,,2 (l_z~l)(l_z;,)
o,icos(nkl)(l_ z l)
.,,,,,
+Z
kl=l[1-2c~
+Z12](l-z21)
~k2 )(1-Z21) M2 (-1)k22[Hk(O,k2)lcos(-~2 = I ( 1 - z l l ) 1 - 2 c o s (t,2u2 ~k2
M1 M2
+ZE
z21
z-2
(_l)kl +k2 4]Hk(kl,k2)lx
k1=1 k2=l Irk2 ~l~osr~~(~_z~l)(~_z; k. 2 J
COSC~)
)
1
+z;21
(6)
where N 1 and N 2 are odd.
B. Type 2-2 Frequency Sampling Filters The Type 1-1 frequency sampling filters described in Equations (5) and (6) interpolate a frequency response through a set of frequency samples for ~ R Q starting at col - 092 = 0. Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (N1N2) frequency samples starting at (.01 -- ~ N 1 and (.02 - ~ N 2 instead of col -- (-02 -- 0. This type of frequency sampling filter is called a Type 2-2 frequency sampling filter. To develop Type 2-2 frequency sampling filters, the system function,
H(z1,z2),and the set, Hk(k)for k ~ RK,of frequency samples are related such that
124
PETERA. STUBBERUD
Hk (k) = IHk (k) Ie jO(k)
= H(Zl,Z2)lZl =e j(2rc/N1 )(kl +l/2),z2=eJ(2rc/N2 )(k2 +1/2) This relationship modifies the DFT into the form
"~, h(n)e-J(2r~/N1 )(kl +1/2)nl e-J(2n/N2 )(k2 +1/2)n2
Hk(k)
(7)
n~_RN where the modified IDFT is Z Hk(k)eJ(2x/N1 )(kl +1/2)n I eJ(2rr,/N2)(k 2 +1/2)n2 . (8)
h(n) = ~
NIN2 keRK Using a development similar to the one used for deriving the system function for Type 1-1 frequency sampling filters, Equation (8) can be substituted into the filter's z transform so that the system function, H(Zl,Z2), of a Type 2-2 frequency sampling filter can be written as 1+ Zl NI 1+ Z2 N2
H(Zl,Z2)=~
~ • N2
NI
NI-I N2-1
Hk(kl'k2)
k~-o k2-0 1-
zi-
l-
(9)
2)zl
The frequency sampling filter in Equation (9) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 2-2 frequency sampling filter with a real impulse response has frequency samples of the form
IHk(k,,k2)l-Ink(N1 - 1 -
kl,N 2 - 1 - k2) I
O(kl,k2) =-O(N 1 - 1 - kl,N 2 - 1 - k2). Using a development similar to the one used for deriving the system function for linear phase Type 1-1 frequency sampling filters, the system function of a linear phase Type 2-2 frequency sampling filter can be written as
H(zl,z2) =
1+
Z? ul
1+
Z2N2 f~,
l ]~,
z)
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
+Z
125
(,+ z~')(,- 2c~(k'LN, + ,, 2,]z l + z 21
k1=0 NI-1 M2-1
+Z
Z (-1)k'+k22[Hk(kl'k2)l x
kl=0 k2=0
c~
,~ l~ltz-1I + I) - ~-~-2 (k2 + -~1_],
1 +z2)
(I-2cos[-~-i (kl + l)]zll + Zl 2)(I-2cos[-~-2 (k2 + l)]z21 + z22) c~
+ 1 ) + ~ 2 (k2 +1)](1 +zllz21)
(I- 2cosE--~-i(kl + I 'lz-I
- 2cos[-~-2(k2 1'Iz-I
(lo)
where N 1 and N 2 are odd. The frequency sampling filter in Equation (10) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 2-2 frequency sampling filter with fourfold symmetry has a DFT of the form
k2)l =[Hk(N1-1-kl,k2)l
IHk(kl, k2)l- IHk(kl,N2
-
= In k (N 1 -
1 -
1 - k 1, N 2 - 1 - k2)1.
Substituting these constraints into Equation (10) yields,
H(Zl,Z2) = ml_l
+Z
k1=0
1 + Zl N1 1+ Z2N2 f_ ml
N2
( l + z l l ) ( l + z 2 l)
(_l)klNlHk(kl,M2)]cosI_~ll(kl
k'LNI
__ + 1 ) _ ~N2]( 12 +zll)
126
PETERA.STUBBERUD
M2_1 (-1)k22[Hk(ml,k2)lcosI--~2(k2 + 89
+2;
+Z21)
M1 M2
+Z
Z (-1)kl+k24lSk(kl'k2)[x
k1=0 k2 =0 sin[-~ll (kl + 89 sin[-~22 (k2 + 89
+ Zl 1)(1 + z21 )
(1-2COS[-~1(kl +l)]zll +Zl2)(1- 2cosI-~-2(k2 + 89
2)
(11)
where N 1 and N 2 are odd.
C. Type 1-2 Frequency Sampling Filters Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (NIN2) frequency samples starting at 091 = 0 and c02 - g/N2. This type of frequency sampling filter is called a Type 1-2 frequency sampling filter. To develop Type 1-2 frequency sampling filters, the system function, H(Zl,Z2), and the set, Hk(k ) for k ~ R K, of frequency samples are related such that Hk(k ) = IHk(k)leJO(k)
: H(zl,z )lz,
=eJ(2n/N1)kl =eJ(2r:/N2)(k2+1/2)
This relationship modifies the DFT into the form
Hk(k) = Z h(n)e-J(Zx/ul )klnle-J(Zx/N2)(k2+1/2)n2 n~_RN where the modified IDFT is h(n)= 1-----~ ZHk(k)eJ(2x/N1)klnleJ(2x/N2)(k2+l/2)n2.
(12)
N1 N2 k~RK Using a development similar to the one used for deriving the system function
DESIGN FOR 2-D FREQUENCY SAMPLINGFILTERS
12/
for Type 1-1 frequency sampling filters, Equation (12) can be substituted into the the filter's z transform so that the system function,
H(z 1,z2), of a Type 1-2
frequency sampling filter can be written as
H(Zl,Z2) =
1- Zl N1 l + z 2 N2
N1 N 1-1 N 2 -1 kl=0 k2=0
N2
•
(1-eJ(2rc/N1)klzll)(1-eJ(2rc/N2)(k2+l/2)z21)
The frequency sampling filter in Equation (13) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 1-2 frequency sampling filter with a real impulse response has frequency samples of the form
[Hk(kl,k2) ]= [Hk(N 1 - kl,N 2 - 1 - k2) ] O(kl,k2) = -O(N 1 - kl,N 2 - 1 - k2). Using a development similar to the one used for deriving the system function for linear phase Type 1-1 frequency sampling filters, the system function of a linear phase Type 1-2 frequency sampling filter can be written as
1- z7 N1 1 + z2 N2 f
+
(1- 1)(l_ cos[
Hk(O'M2 )
+ )1z;1+
M1 N 2-1
+Z Z
k1=1 k2 =0
(_l)kl +k2 2]Hk(kl,k2)]•
)
128
PETER A. STUBBERUD
where N 1 and N2 are odd. The frequency sampling filter in Equation ( 1 4 ) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 1-2 frequency sampling filter with fourfold symmetry has a DIT of the form I ~ k ( k l , k : ! ) l = I ~ k ( k l- ,1~- k2z ) I = IHk (
.
~ -1kl k2 I)
= ( H ~ ( N- Ik l , N 2 - 1 - k 2 ) I .
Substituting these constraints into Equation ( 1 4 ) yields,
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS M1
129
M 2-1
Z (-1)kl +k2 41Hk(kl'k2)l•
+Z k 1=1
k2 =0
c o s ( k k, )sin[~2-2 (k2+ 89
z~-1)(1+ z21)
[1-2 cos(-~l k1) z l l + Zl 2 ](1-2 c o s I ~2 (k2+ 89
(15) + z2 2 )
where N 1 and N2 are odd. D. Type 2-1 Frequency Sampling Filters Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (N1N2) frequency samples starting at 091 - ~N 1 and co2 - 0. This type of frequency sampling filter is called a Type 2-1 frequency sampling filter. To develop Type 2-1 frequency sampling filters, the system function, H(Zl,Z2), and the set, Hk(k) for k ~ RK, of frequency samples are related such that H k (k) = IHk (k)[ e jO(k)
= n(Zl,Z 2 )lzl=eJ(2rc/N1)(kl +l/2),z2=eJ(2rt/N2)k 2 " This relationship modifies the DFT into the form Hk(k) = Z h(n)e-J(2rt/N1 )(kl +1/2)n 1e-J(27r,/N2)k2n2 n~R N
where the modified IDFT is h(n) = ----~1 ZHk(k)eJ(2rr'/N1)(kl+l/2)nle j(ert/Na)kan2 .
(16)
N1N2 keRK Using a development similar to the one used for deriving the system function for Type 1-2 frequency sampling filters, Equation (1 6) can be substituted into the the filter's z transform so that the system function, H(z 1,z2), of a Type 2-1 frequency sampling filter can be written as
130
PETERA. STUBBERUD
H(Zl,Z2) =
1+ Zl N1 1- z2 N2 N1 N2
N 1-1N2-1
Z
Z
:o
:o
Hk(kl'k2)
(1
+l ' 'Zl II(1
z;l )
(17)
The frequency sampling filter in Equation (17) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 2-1 frequency sampling filter with a real impulse response has frequency samples of the form IHk(kl,k2)l = [ H k ( U l
- 1
- kl,N 2 -
k2) I
0(k 1, k2 ) = - 0 ( N 1 - 1 - k 1, N 2 - k2 ). Using a development similar to the one used for deriving the system function for linear phase Type 1-2 frequency sampling filters, the system function of a linear phase Type 2-1 frequency sampling filter can be written as
H(Zl,Z2) = 1+ Zl ul 1-Z2 N2 f _ MI-1
_
__
(-1)kl2]Hk(kl,O)lsinI-~ll(kl +1)]( 1 + z ?
1)
~_~cosI~(~l +l)z-,+Zl ~- (l-z;,) [.NI 1
N1-1 M2 + Z Z (-1)k' +~2 21nk(kl,k2)I• k1=1 k2 =0
sin[-~l (kl ( 1 - 2 cosI--~-1 (kl+ 89
_f (1 - 2 c o s
+ 89
)
+ Zl 2 ) I 1 - 2 c o s ( ~ 2 k2/z21+ z2 2 ]
sin[~ll(kl + 89 k21(zll -z21.! k + z +z 1 2cos I2~l-~l-~-~l';-1S)~----(~2k21z21+z221
where N 1 and N 2 are odd.
(18)
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
131
The frequency sampling filter in Equation (18) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 2-1 frequency sampling filter with fourfold symmetry has a DFr of the form
IHk(kl' k2 )l = [Hk(kl' N2 - k2 )l =[Hk(N1- 1 - kl,k2)l =[Hk (N 1 - 1 - kl, N 2 - k2)[. Substituting these constraints into Equation (18) yields,
H(Zl,Z2)= 1+ Z?N1 1-Z2 N2 f _
__
M,-, (-1' kl 2[Hk(kl,O)[sin[-~l(kl+ l ) ] ( l + z l l )
+Z
kl=O (1-2cos[2~(kl+l)lzll LN1
+ z12 ) (1 - z21 )
+Z MI-1 M 2
+ Z Z (-1~' +~ 41"~(~1,~1 • k 1 =1 k 2 = 0
sin[~ll (kl+ 1)] c o s ( ~ 2 k2)(1+ Zl 1)(1- z21)
(19)
where N 1 and N 2 are odd.
E. ComputationalAdvantage of Frequency Sampling Filters Each of the frequency sampling filter system functions described by Equations (3), (5), (6), (9)-(11), (13)-(15) and (17)-(19) describe an architec-
132
P E T E R A. STUB B E R U D
ture that has a nonrecursive comb filter created by the terms, Ill+ z,N1 \] and (l + z2N2 ) , in series with a number of parallel resonators 1. When ~' "most! of a frequency sampling filter's frequency samples, the Hk(k)'s, are exactly zero, most of the filter's resonators are not implemented. Therefore, in the case of a narrowband filter where most of the frequency samples are in the stopband and are set equal to zero, only a small number of the filter's resonators need to be implemented. In such instances, frequency sampling implementations usually require fewer arithmetic operations than the direct convolution implementations. Each of the resonators in the frequency sampling filters described by Equations (3), (5), (6), (9)-(11), (13)-(15) and (17)-(19) can be represented by a separable system function 2. If separable system function, H(Zl,Z2), can be written as H(Zl,Z2) - HI(Zl)H2(z2), then it can be shown that that system is stable iff each of the 1-D system functions, HI(Zl) and H2(z2), are stable[2]. Because each of the resonator system functions are separable, it is straight forward to observe that for each of the frequency sampling filters the zeroes createdby the terms, (I _+zlNI ) and (l _+z2N2 ), exactly cancel the resonators ' poles that lie on the 4-D unit sphere. When implementing these filters with finite word length systems, exact pole zero cancellation generally does not occur, and uncanceled poles on the unit sphere will cause the filter to be unstable. To prevent this instability, each z.-1 in each of the frequency sampling filter system functions is replaced by rzi -~1",r < 1", i - 1,2. This moves all the poles inside the 4-D unit sphere. To keep the frequency response of these modified filters as close as possible to the frequency response of the original filters, r should be chosen near to the value 1. Substituting rz.-1 for z- 1 where r < 1 guarantees the stability of a frequent l cy sampling filter at the cost of increasing its computational requirements. If implemented with z i- 1 replaced by rz i 1, a linear phase frequency sampling filter requires at most 8 multiplies per resonator, and a linear phase frequency sampling filter with fourfold symmetry requires at most 5 multiplies per resonator If only K of the frequency samples are non-zero, then a linear phase fre1. The term resonator is used in this paper to denote a system which has its poles on or near the 4-D unit sphere or has its poles on or near the 2-D unit circle. 2. A system function is said to be separable iff it can be expressed as the product of one dimensional (l-D) system functions.
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
133
quency sampling structure requires approximately 8K multiplies per output sample, and a linear phase frequency sampling filter with fourfold symmetry requires approximately 5K multiplies per output sample. If a linear phase FIR filter with region of support RN is implemented using direct convolution, it requires approximately N1N2/2 multiplies per output sample; and if a linear phase FIR filter with fourfold symmetry and region of support RN is implemented using direct convolution, it requires approximately N1N2/4 multiplies per output sample. Because multiplications typically require more time to compute and are more complex to implement than adds, a frequency sampling structure becomes computationally more efficient than a direct convolution implementation when it uses fewer multiplies. Therefore, if we wish to implement a linear phase filter with a region of support R N, a frequency sampling filter can implement the filter more efficiently (in the sense of fewer multiplies) than a direct convolution implementation when 8K < N1N2/2 or
K < N1N2/16. If we wish to implement a linear phase filter with fourfold symmetry and region of support RN, a frequency sampling filter can implement the filter more efficiently than a direct convolution implementation when 5K < N1N2/4 or K < N1N2/20.
III. The Design of 2-D Frequency Sampling Filters It was shown in section II that a frequency sampling filter approximates a desired frequency response by interpolating a frequency response through a set of frequency samples. In this section, a design technique that controls interpolation errors between samples and approximates a desired frequency response is developed for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. This design technique approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response. To maintain the computational efficiency of the filter, this error criterion is minimized while constraining the frequency samples in the stopband to be identically zero. This design technique is developed as a constrained optimization problem which is solved by the method of Lagrange multipliers. The Lagrange multipliers optimization method results in a set of linear equations which are solved to determine the filter's
134
PETER A. STUBBERUD
coefficients. The design method applies to each of the four types of 2-D frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. The frequency response of a frequency sampling filter which has an impulse response, h(n) for n ~ R N, is Hco ( tO) = H ( eJ(-Ol , eJCO2 ) =
h(n )e-J~ nl e-JOJ2 n2
n~R N If the filter has a real impulse response, fourfold symmetry and linear phase then h(n 1, n 2) = h(n 1,N 2 - 1 - n 2) = h ( N 1 - 1 - n 1, n 2 ) = h(N 1 -1-n
1,N 2 - 1 - n
2),
and it can then be shown that Hco(a)) can be written as H(o (o)1, (.o2) = e-J[~176 where Hr((_o 1,(.o2) is a real function given by M1 M2 Z Z 2sgn(nl)2sgn(n2)h(Ml-nl'M2-n2)cOs((-Olnl)cOs(OJ2n2) n 1=0 n2 =0 M 1 - (N 1-1)/2, M 2 - (N 2-1)/2, N 1 is odd, N 2 is odd and sgn(o) is the signum function. Hr(O~) is the amplitude of Hoj(m) and is often referred to as the zero phase frequency response or the amplitude response of the filter. If we let Hr(t~
x=[h(MI,M2)
2h(MI-I,M2)...
4h(1,0)
4h(0,0)] r
where the superscript T denotes transpose, and let s(col ,co2 )=[Sr (oJ1 ,co2 )]=[cos(al (r)col)cos(a2 (r)co2)] where r = 0, 1 .. (MI+I)(M2+I) - 1 s (co) represents the rth element in the column vector, s(m), '
"
'
~
r
a I ( r ) = r m o d ( M 1 +1)
a2 (r)=int( M-~+I )
DESIGNFOR 2-D FREQUENCYSAMPLINGFILTERS
135
and int(x) truncates x, then
Hr(tO) = xTs(to) = s T (to)x. (20) In Equations (1), (8), (12) and (16), the impulse response, h(n), of each of the frequency sampling filters is expressed as a function of its DFT. By constraining the filters to have linear phase, the IDFT's can be written as Z 2sgn(kl+k2)(_l)(kl+k2) N1N2 [H~(kl'k2)[ • k~Rll cos[2~ kl(nl + 89 27r )] LN1 ~22 k2(n2 + 89 Z k6R12
2sgn(kl+M2_k2)(_l)(kl+k2)
N1N2
Type 1-1
Igk(kl,k=~l•
sinF2Zc 2zr (k2 +l)(n2 + 89)] LN 1 kl(n 1 + 89 "~2 h(n) =
~, k6R21
2sgn(Ml-kl+k2)(-1)(kl+k2)
N1N2
IH~(kl,k=~l•
sinF~ (~, +~)(n,+~)+ ~ Z
k6R22
2sgn(M1 -k 1+ M2 -k 2 ) (_ 1)(kl +k2 + 1)
N1N2
cos[~ (~1+~)(nl+~)+ ~
Type 1- 2
)]
Type 2-1
IHk(kl,k2~l•
)1
Type 2- 2
where R l l = { ( k l , k 2 ) ' 0 < k I <_N1-1,1_< k2
IA IA
~
,=...
IA
~
~
IA
o.
il
I
!
i
•
,-
+
7"
~,..-
I
I
~.
~-
~I--
,
+
i
, t'4
m
I
i
~
i
x
+ I
~
~
~
*
+
+
i
~
b~
I,,./
, t'4 Ill
I
I
II
-I~I--
+
I
I
m
x
'~
b~
I
~1-"
m
x
,~
b~
b~
t'rl
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
2sgn(kl )2sgn(k2)(_l)(kl +k2)
N1N2
137
x
cosFIN 1 l(n + 2" cos "~2 k2 (n2 + 89
y el-1
2 sgn(kl)2sgn(M2-k 2) (_ 1)(kl +k2 ) x
N1N2
COSI2' klLUl (nl+l]]sin[ ( k 2 z ' LNe' . ]27r
f(n) =
2 sgn(M1-kl )2sgn(k2)(_ 1)(kl +k2 )
N1N2
Type 1-2 + 1)(n2 +1)]
•
sin[2n'(klLN 1, + 89 nl +l)]c~ n 2 z[_N2, j
+1)]
Type 2-1
2sgn(Ml -kl !2sgn(M2--k2_!(-1)(kl +k2! • N1N2 l
2re sin[ "~1"1(kl + ~)(nl + 2")lsln[ "~2-2(k2 + 1)(n2 + 1)] 9
2zr
1
1
Type 2- 2
9
for r - 0 , 1, (MI+I)(M2+I) - 1; k 1 = int[r/(M2+l)]; and k2 = rmod(M2+l), then h(n) = xTf(n)= fT (n)X where the appropriate expression for f(n) is used depending upon the type of frequency sampling filter being designed. If we also define F-[f(M1,M2)
2f(M 1 - 1,M2) -.. 4f(1,0) 4f(0,0)] T
then x - FX
(21)
Substituting (21) into (20),
Hr(t~) - xTFTs(to) - s T (to)FX. (22) An approximation of a linear phase filter with a desired passband amplitude, Hal(to), can be obtained by minimizing the mean square error between the desired amplitude response, Hd(tO), and the FIR filter's amplitude response, H(~). If we let
138
PETER A. STUBBERUD
~pb
= {(COl,CO2)" COl E passband, Co2 ~ passband},
and let J b represent the mean square error over the region of passband frequencies, then
Jpb = m(~pb~1
SftoE~pbIHr (to) -
Hd (to)l 2 dto
(23)
where m(o) is the 2-D Lebesgue measure. Substituting Equation (22) into Equation (23), Jpb can be written as
1pb-"----~x T F T Jpb (X) = m(~"2
SS s(~
( ~ ) d ~ FX
~E~pb _2xTF T
SS Hd(~)s(~)dm + ~-~ pb
S H (m)dr
(24)
mE~ pb
If we let W(m) = s(m)s r ( ~ )
then w(m) = [Wrc(O)] = [cos(a 1(r)co 1)cos(a 2 (r)co 2) cos(a I (c)Co1)cos(a 2 (c)Co2)] where Wrc(m) represent the element in the rth row and the cth column of the matrix, W(to), andr, c - 0 , 1..... (MI+I)(M2+I) - 1. Defining Q(m) as the matrix,
the element in the rth row and the cth column of the matrix, Q(~), can be written as
Qrc ( C~ 0)2 ) - Qlc (col)Q2c ( o92 ) where
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
ai(r)=ai(c)=O
Ogi i Qrc ( Ogi) =
Ogi
139
sin 2ai(r)O9i
2 4ai(r) sin[ai(r) + ai(c)]o9i sin[ai(r ) - ai(c)]O9i 2[ai(r)+ai(c)] + 2[ai(r ) - ai(c)]
ai(r)=ai(c)•O ai(r) r ai(c)
and r, c = O, 1..... (M 1+ 1)(M2+ 1) - 1. If we define gg
QP
= JJm ~ p b W(m) dm
where Qp is calculated by evaluating
and
Q2 (O92)= [Q2c (o92)] at the appropriate limits and performing an element by element product of the resulting matrices, and also define the terms,
R(m) = f l Hd (m)s(m)do~
~(~)=
~IH2(~)doJ
Rp=
II Hd(m)s(m)dm
~'p=
IIH2(o~)dm
O~pb then Equation (24) can be written as
Jpb(X) =
1.......~-[xTFTQpF x _ 2xTFTRp + )'p] m(~pb)
(25)
For filters which approximate constant values in their passbands, the expression for R(m) can be simplified. If Hd(O~) is equal to a constant in the passband, then without loss of generality, we can let Hd(m) -- 1, and Equation (25) becomes
Jpb(X) =
I~[xTFTQpFX m(~pb)
- 2xTFTRp]
+1
(26)
Because the elements in the vector, s(~), are separable, the matrix, R(~), can be written as
140
PETER A. STUBBERUD
where
Coi Ric(f~ - sin[ai(r)ooi]/ai(r )
ai(r) = O ai(r ) r 0
This expression for R(m) eliminates the need for integrating Hd(m)s(m) when calculating Rp. The matrix, Rp, in Equation (26) can be calculated by evaluating Rl(O)l)= [Rlc(COl)] and R2(fo2)=[R2c(OO2)] atthe appropriate limits and performing an element by element product of the resulting matrices 9 An approximation of a desired amplitude response can also be obtained by imposing constraints on the amplitude response and its derivatives at particular frequencies 9 For Example, nr(m)lm=to0 = K000
o(m+n)Hr(~) I o300m~o)~
-..
Hr(m)lm=m t =
tg(m+ n) Hr (OJ) = gmn 0 ""
= Kmn I tO=O~1
~=~0
c92PHr(~) I OcOp OOoff
02 P Hr (O~) = K ppO
KO01
= Kpp l
"'"
~=~0
~=~1
where the K's are constants 9 Using Equation (22), these passband constraints can be written as
o(m+n)Hr(~ ) = o(m+n)s T (O~) FX.
oo,Too, If we let cT=Is(0~0 ) S(t~l )
o3(2P)s(0~)
o32s(t0) 00910(.02
ta=o~0
and KT=[K000
K001 ...
KI10
...
1
Kppl],
the passband constraints can be written in the matrix form, C F X - K. (27) If we let Jsb represent the mean square error over the stopband frequencies, then
D E S I G N F O R 2-D F R E Q U E N C Y S A M P L I N G FILTERS
141
1
Jsb - m(ns b) IIcoef~sb [nr(cO) - Hd(~)l 2 dm
(28)
where
~sb Because
= {(CO1,0)2)" CO1E stopband, 092 e stopband}.
Hd(m) equals zero in the stopband, Equation (28) can be written as Jsb = m(f2sb )
(29)
~-f~sb
Substituting Equation (22) into Equation (29), Jsb(X) =
1
xTF Tff,,,
s(to)s T(tO) dto
FX
(30)
Recall that earlier, we let W(to) -s(~)sT(to) and Q(m) = f ~ w(m) dto. Thus, if we define -
Q s IIco~f~sb
W(m)
dm
where Qs is calculated by evaluating Q l(O~l) and Q2(co2) at the appropriate limits and performing an element by element product of the resulting matrices, Equation (30) can be written as Jsb(X) =
~m----xTFTQ s FX. m(asb)
(31)
To generate an efficient realization of a linear phase FIR filter with a frequency sampling structure, the frequency samples in the stopband are constrained equal to zero. Recall that the vector, X, contains the frequency sampies. In Equations (25), (27) and (31), constraining frequency samples to zero is equivalent to removing elements in the vector, X, that contain the corresponding frequency samples and also removing the corresponding columns from the matrix, F. If we define these reduced terms as X and F C
C~
then the
design problem can now be stated as follows. Minimize the error function
142
PETER A. STUBBERUD
J(X c) = t~pb (Xc) + (1 - a)Jsb (X c) =
~
1-o~
m(~sb)
XcT FcT Qs Fc x m(~pb)
c
[Xc c Q cXc
where 0 < ~ < 1 subject to the constraint, CF X - K. The scalar term, a, in C
C
Equation (32) allows the designer to weight the relative importance between the mean square error over the passbands and the mean square error over the stopbands. If the number of non-redundant constraints in Equation (27) is less than the number of elements in the vector X c, then the remaining degrees of freedom are used to minimize J(Xc). If the number of non-redundant constraints is equal to the number of elements in the vector X c, then the design becomes a linear algebra problem and the frequency samples are unique for the constraints. And if the number of non-redundant constraints is greater than the number of elements in the vector X c, then no solution will exist which satisfies all the constraints. Assuming that the number of non-redundant constraints in Equation (27) is less than the number of elements in the vector X c, the design problem can be solved by defining a Lagrange multiplier vector as ~T =1.[zOo0
zoo1
"'" '~'110 "'"
Zppl]j
and minimizing the augmented cost function,
Ja(Xc,~) = OClpb(Xc) + (1 - a)Jsb(Xc) + ~T(CFcXc - K).
(33)
Substituting the appropriate expressions for %b and Jsb into Equation (33),
Ja(Xc ~ ) = ~ [ X ' m(~pb)
T T c Fc QpFcXc - 2 x T F T R p + Yp]
1-a XcT FcT Qs FcXc + ~T (CFcXc _ K) m(~2sb ) The necessary conditions for an optimal solution are ,)J a ( X c , ~.)
O~c
=0
(34)
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
143
and
tgJa ( X c , ~,)
Ok
(35)
"-0.
Because F c TQpFc and Fc TQsFc are symmetric matrices, Equation (34) becomes 2 ( 1 - a) FTQs Fc X + FTcT2~ = 0. (36)
m(~sb)
m(f2pb)
c
Equation (35) implies
OJa(Xc,~)
ovk
= CFcX c - K = 0.
(37)
Equations (36) and (37) can be written in the matrix form,
I
2a
m(~pb )
0 lrjL-2j=,r r
FTQpFc + 2 ( l - a ) F t Q s Fc ,, F T c T
m(f2sb )
~
Xc
I m(f2Pb )
K
(38)
A solution to Equation (38) will exist when C has full rank which occurs when C contains no redundant or trivial constraints.
IV. Example A. Example of a Type 1-1 Frequency Sampling Filter In this example, a Type 1-1 linear phase fourfold symmetric frequency sampling filter that approximates the amplitude response, Hd (~ - {10
~~2Pb OiE~sb
where
f2pb - {(091,092) .0 < 091 -< 0.1r~,0 _<092 _<0.1~} ~sb - {(091,092)" 0-13rt <
091 -< rtU0.13rt < 092 < rt}
is designed. The filter will be designed to have an impulse response that has a region of support R N where
R N = {(nl,n2)" 0 < n 1 < 100,0 < n 2 < 100}. To approximate the passband response, the mean square error over the passband will be minimized subject to the constraints
144
PETER A. STUBBERUD Hr(to)lm-(O, - ,0) = 1
t?(m+ n) Hr (to)
=0
where n, m = 0, 1, 2, 3, 4 excluding m = n = 0. To approximate the stopband response, the scalar term, o~, which weights the relative importance between the mean square error over the passbands and the mean square error over the stopbands, will be set to 0.1. Substituting c~ = 0.1 into Equation (38) gives 0.2
m(~pb ) where
FTQpFc +
m(~sb )
s
~
FTcT
Xc ]L-~..j=
0.2
FTRp
m(~pb)
l.....
(39) .....
m(~pb)=(O.l~) 2 and m(~sb)=~2-(O.13~;) 2 .
To solve Equation (39), the matrix, C, must have full rank, and the nonzero frequency samples must be determined. Because the odd ordered partial derivatives of a linear phase filter evaluated at to - (0,0) are zero, constraining them is redundant, and they must be omitted. Thus, for the matrix, C, to have full rank, only the constraints Hr(to)]m=(O,O ) = 1
tg(m +n ) Hr ( to )
=0
where n, m - 0, 2, 4 excluding m = n - 0 are considered. The nonzero frequency samples include all of the frequency samples not in the stopbands. In this example, the stopband frequency samples are H(k) for k ~ {(kl,k2) 97 < k I _< 5 0 U 7 < k 2 < 50} which implies that the nonzero frequency samples are H(k) for k ~ {(kl,k2)" 0 < k 1 _< 6,0 _
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
145
arg[Ha~(fOl, co2 )] = -50091 - 50o92 .
V. Summary and Conclusions In this chapter, four types of f r e q u e n c y s a m p l i n g filters were developed. F o r each type of filter, s y s t e m functions were d e v e l o p e d for filters with arbitrary impulse responses and arbitrary frequency responses; filters with real impulse responses and linear phase; and filters with real i m p u l s e responses, linear phase and fourfold s y m m e t r y . Also, in this chapter, a technique for designing frequency s a m p l i n g filters with real i m p u l s e responses, linear phase and fourfold s y m m e t r y was developed. This design t e c h n i q u e a p p r o x i m a t e s a desired frequency r e s p o n s e by m i n i m i z i n g a w e i g h t e d m e a n square error over the stopbands and passbands subject to constraints on the filter's a m p l i t u d e response. This a p p r o x i m a t i o n technique results in a set of linear equations which are solved to d e t e r m i n e the filter's coefficients. This design m e t h o d applies to each o f the four types o f f r e q u e n c y s a m p l i n g filters d e v e l o p e d in this chapter.
Table 1. Nonzero Frequency Samples for the Example. H(0,0) - 1.0000000
H(2,3) - H(3,2) - 0 . 9 9 6 7 5 4 4
H(0,1) - H(1,0) - 0.9999942
H(2,4) - H(4,2) - 0.9926905
H(0,2) - H(2,0) - 0 . 9 9 9 8 7 4 4
H(2,5) - H(5,2) - 0.9701703
H(0,3) - H(3,0) - 0.9940741
H(2,6) - H(6,2) - 0.5011298
H(0,4) - H(4,0) - 0.9953709
H(3,3) -- 1.0131842
H(0,5) = H(5,0) - 0.9704475
H(3,4) - H(4,3) - 0.9991182
n ( 0 , 6 ) - H(6,0) -- 0 . 5 2 6 6 0 8 0
H(3,5) -- H(5,3) -- 0 . 9 7 7 7 6 2 4
H(1,1) - 1.0005610
H(3,6) - H(6,3) = 0 . 4 5 3 0 7 1 0
H(1,2) - H(2,1) - 1.0064317
H(4,4) -- 0.9915667
H(1,3) - H(3,1) - 1.0057280
H(4,5) - H(5,4) -- 0 . 9 7 1 9 1 5 4
H(1,4) = H(4,1) - 1.0035502
H(4,6) - H(6,4) -- 0 . 4 9 4 1 6 3 6
H(1,5) - H(5,1) -- 0 . 9 7 9 1 7 3 4
H(5,5) -- 0.9312007
H(1,6) - H(6,1) - 0.5152585
H(5,6) - H(6,5) -- 0.3930366
H(2,2) = 0.9951068
H(6,6) - 0.0884227
146
PETER A. STUBBERUD
1 0.8
~,
~, 0.6
0.4 0.2 O= /Z"
-0.2tr ~ O~( . 0.21r " ~ tra~"/#~ 0.6x ~/e)
0.2zr tr
-n:
0.6zt"
-0.6n: OX(,~&.ls~ga'p\e')
Figure 1. Magnitude Response of the Type 1-1 Frequency Sampling Filter in the Example. 0 -50
~
-100 -150 -200 0 ~0 . 4~0 ' 2' ~ ~ r_~r 0.61r ~
e)
~ ' <0.8'1r~
8~
O.
. . > ~ i ~ 0.2zr o
~ o3x
0.6~ ~ 0.4/1: ~$~6.1s~9\e')
Figure 2. Magnitude Response of the Type 1-1 Frequency Sampling Filter in the Example.
zc
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
~"
_'g
147
0
o ~o -0.5
-1 0
0-047r~"~-..~
o~ ~
o.6~ ~
~'~'~~ O.06zr
0.08zr
'~'4,~,,,. 0 . 0 8 ~ r ~ ~ 0.04zr '"~ o"l x " ~ " i ~ O O.02zrcox(~6.1 s~px'e')
Figure 3. Magnitude Response of the Passband for the Type 1-1 Frequency Sampling Filter in the Example.
VI.
References
1. L.R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing, Prentice Hall, New Jersey (1975). 2. D.E. Dudgeon and R. M. Mersereau, Multidimensional Digital Signal Processing, Prentice Hall, New Jersey (1984). 3. J.V. Hu and L. R. Rabiner, "Design Techniques for Two-Dimensional Filters," IEEE Transactions on Audio and Electroacoustics, Vol. AU-20, No. 4, pp. 249-257 (1972). 4. J.S. Lim, Two-Dimensional Signal and Image Processing, Prentice Hall, New Jersey (1990).
Abstract In this chapter, system functions are developed for two dimensional (2-D) frequency sampling filters that have real impulse responses and linear phase and for 2-D frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. Under certain conditions, these frequency sampling filters can implement narrowband 2-D linear phase filters and narrowband 2-D linear phase filters with fourfold symmetry much more efficiently than direct convolution implementations. Also, a technique for determining optimal frequency sampling filter coefficients is developed for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. This technique approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response.
I. Introduction Some two dimensional (2-D) signal processing systems, including image processing systems, require linear phase or zero phase filters. A 2-D linear phase or zero phase filter implemented by direct convolution uses the filter's impulse response as coefficients. If a 2-D linear phase filter has a region of support, R N, where CONTROL AND DYNAMIC SYSTEMS, VOL. 77 Copyright 9 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.
117
118
PETER A. STUBBERUD
R N = {(nl,n2)'0_< n I _
and N 1, N 2, n 1 and n 2 are members of the set of integers ( N 1, N 2, n 1, n 2 ~ / ) , then the filter's impulse response has the form h(nl,n2) = h(N 1 -1-nl,N 2 -l-n2) [ 1]. If the filter is implemented using direct convolution, approximately N 1 N 2 / 2 multiplies are required to compute each output sample. To reduce the
number of multiplies required by a direct convolution implementation of a linear phase filter, fourfold symmetry conditions are often imposed on linear phase filters. If a 2-D fourfold symmetric linear phase filter has support over the region, R N, then the filter's impulse response has the form h ( n I ,n2) = h ( n l , N 2 - 1 - n2) = h ( N 1 - 1 - n 1, n 2 ) = h(N 1 -1-n
1,N 2 -1-n
2)
[2], and a direct convolution implementation of the filter requires approximately N 1 N 2 / 4 multiplies per output sample. Unlike direct convolution implementations which use the filter's impulse response as coefficients in the filter's implementation, frequency sampling filters use frequency samples, which are specific values from the filter's frequency response, as coefficients in the filter's implementation. The frequency sampling filters discussed in this chapter interpolate a frequency response through a set of (NIN2) frequency samples. When using a frequency sampling filter to implement a frequency selective filter, the frequency samples that lie in the filter's stopbands can be set to zero and do not have to be implemented. Thus, only the non-zero frequency samples which lie in the illter's passbands and transition bands are required in the filter's implementation. As a filter's passbands or transition bands narrow or its stopband specifications become more stringent, the filter's region of support increases which implies that the values of N 1 and N 2 increase. Because the number of multiplies per output sample of a direct convolution filter is approximately proportional to (N1N2), even a small increase in the filter's specifications can substantially increase a direct convolution filter's computational requirements. However, as a frequency sampling filter's passbands and transition bands narrow and its stopbands increase, the number of nonzero frequency samples decrease; and thus, the number of coefficients required in the filter's imple-
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
119
mentation decreases. Thus, frequency sampling filters have the potential to implement narrowband frequency selective filters much more efficiently than direct convolution filters. In this chapter, the 2-D Type 1-1 frequency sampling filter system function described in reference [3] is further developed for frequency sampling filters that have real impulse responses and linear phase and for frequency sampiing filters that have real impulse responses, linear phase and fourfold symmetry. This chapter also introduces three new types, Type 1-2, Type 2-1 and Type 2-2, of frequency sampling filter system functions. These system functions are also developed for frequency sampling filters that have real impulse responses and linear phase and for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. The resulting system functions for frequency sampling filters that have a real impulse responses and linear phase interpolate a frequency response through a set of approximately
N1N2/2independent samples, and the resulting system functions for frequency sampling filters that have a real impulse responses, linear phase and fourfold symmetry interpolate a frequency response through a set of approximately
N1N2/4independent samples.
As a result these system functions are computationally more efficient for frequency sampling filters that have real impulse responses, linear phase and for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry than the system function described in reference [3]. Although frequency sampling filters interpolate a frequency response through a set of frequency samples, the frequency response may not be well behaved between samples. References [3; 4] describe two of the design methods currently used to control the interpolation errors between frequency samples and approximate a desired frequency response. In this chapter, a design technique is developed that controls the in-
terpolation errors and approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response. This results in a constrained optimization problem which can be solved by using the Lagrange multiplier optimization method.
120
PETER A. STUBBERUD
II. 2-D Frequency Sampling Filters A. Type 1-1 Frequency Sampling Filters Consider a filter that has an impulse response h(n) where n - (nl,n2), a region of support R N where
RN = { ( n l , n 2 ) ' O < n 1 < N l - l , 0 _ < n 2 < N 2 - 1 } and a frequency response, Hro(to), where
n l,n 2 ~ 1,
0)= (091,0)2) ~ RK~ = {(0)1,0)2)'0 < 0)1 < 2g,0 <_0)2 < 21r} C~ 0)2 e the set of real numbers (9t). Let Hk(k) where k - ( k l , k 2 ) ~ R K - { ( k l , k 2 ) ' O _ < k 1 <_N l - l , 0 < k
2<_N 2 - 1 }
kl,k 2 ~ I
represent a set of (NlN 2) frequency samples taken from the filter's frequency response, H (to) for to e R For a Type l-1 frequency sampling filter, the (o .(_2. set, Hk(k) for k ~ R K, of frequency samples are selected such that 2rt 2rt U k (k) = Ho~ (0)1,0)2)10~1=.~l_lk 1,0~2=.~_2k2 where Hk(k ) can be written as
H k (k) - [ H k (k)leJO(k) where k ~ R K, j = ~ - l and 0(k) is the phase of Hk(k ). Thus, for a Type 1-1 frequency sampling filter, the set, Hk(k) for k ~ R K, of frequency samples corresponds to the Discrete Fourier Transform (DFT) coefficients. Also, the impulse response, h(n) for n e R N, which interpolates a frequency response through the set of frequency samples, Hk(k) for k ~ R K, can be determined from the inverse discrete Fourier transform (IDFT),
~,~'~H k (k)e j(2rt/NI )nlkl e j(2rt/N 2 )n2k 2 , h(n) = ~ 1 NIN2 k~RK
(1)
and a set, Hk(k) for k ~ R K, of frequency samples can be determined from the filter's impulse response by the discrete Fourier transform (DFT),
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS Hk(k) -
121
Z h ( n ) e - J ( 2 g / N 1 )nlkl e-J(2rr'/N2 )n2k2 " neR N
If we let H(z 1,z2) represent the z transform of h(n), then the system function,
H(z 1,z2), of the filter that interpolates a frequency response through the set, Hk(k) for k e R K, of frequency samples is N 1-1N2-1 H(Zl, Z2 ) = Z Z h(nl'n2 )z? nl z2 n2 "
(2)
n 1=0 n2 =0 By substituting Equation (1) into Equation (2), interchanging the order of the summations and performing the summation over the n 1 and n 2 indices,
H(z 1,z2) can be written as 1 - Zl N1 1 - z2 N2 H(Zl,Z2) = ~ x
N1
N2
N 1-1N2-1 (3)
Equation (3), originally developed in [3], describes the general form of a Type 1-1 frequency sampling filter, and has the form of an interpolation formula. The complex function, H(Zl,Z2), interpolates a polynomial through the frequency samples in the set, Hk(k) for k ~ R K, so that n ( z l , z2)] Zl =eJ(2rc / N1)kl ,z2 =eJ(2rc / N2 )k2 = n k (kl' k2 )" Thus, as desired, the frequency sampling filter's frequency response passes through the (N1N2) frequency samples in the set Hk(k) for k ~ R K. The frequency sampling filter in Equation (3) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A filter with a real impulse response has a DFT of the form
IHk(kl,k2)l-]Hk(N1-
kl,N2 -
k2)l
0(k 1, k2 ) = - 0 ( N 1 - k 1, N 2 - k 2 ). It can be shown[ 1] that the phase of a linear phase FIR filter with a region of support R N is
122
PETER A. STUBBERUD
arg[Hco ((01 (02)] = -(01
(N,-1/ (N2-1/
'
- (02
2
Therefore, the phase of the frequency samples is
O(kl'k2) =- N 1
2
--~2 k2
2
9
(4)
(N2-1 /" 2
Substituting these constraints into Equation (3), the filter's system function can be written as
Hz(Zl,Z2) =
1- Z?NI 1- Z2N2 I
Hk(O'O)
~k 1 ~l ~-l~k'21Hk~kl,0~lcos(~/tlZl' t
l--'II--2COS kl/zll'N1 +ZI2]I--z;I +2;
~2
["Trk2 - z;l ) (-1~2 ~l.~(0,~l cos~w)(1
k2=l(1-zll) N 1-1
M2
kl=l
k2=l
+ZZ
1 - 2 c o s -~22k2 z 2 + z 2 2
(_l)kl +k2 21Hk
cos~?
+ -~2 ~k2
z-l z" I" - cos NI
N2
+
(5)
where M1 - (N l- 1)/2, M 2 - (N 2-1 )/2, and N1 and N 2 are odd. The frequency sampling filter in Equation (5) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 1-1 frequency sampling filter with fourfold symmetry has a DFT of the form
IHk (kl' k2)l - I H k
(kl, N2 - k2)l - I H k (N1 - kl, k2 )1 - I H k (N1 - kl' N2 - k2)l"
Substituting these constraints into Equation (5), the filter's system function can be written as
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
Hz(Zl,Z2)= 1- ZlN1 1- z2N2 I
.,,,,
123
Hk(O'O)
.,,,2 (l_z~l)(l_z;,)
o,icos(nkl)(l_ z l)
.,,,,,
+Z
kl=l[1-2c~
+Z12](l-z21)
~k2 )(1-Z21) M2 (-1)k22[Hk(O,k2)lcos(-~2 = I ( 1 - z l l ) 1 - 2 c o s (t,2u2 ~k2
M1 M2
+ZE
z21
z-2
(_l)kl +k2 4]Hk(kl,k2)lx
k1=1 k2=l Irk2 ~l~osr~~(~_z~l)(~_z; k. 2 J
COSC~)
)
1
+z;21
(6)
where N 1 and N 2 are odd.
B. Type 2-2 Frequency Sampling Filters The Type 1-1 frequency sampling filters described in Equations (5) and (6) interpolate a frequency response through a set of frequency samples for ~ R Q starting at col - 092 = 0. Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (N1N2) frequency samples starting at (.01 -- ~ N 1 and (.02 - ~ N 2 instead of col -- (-02 -- 0. This type of frequency sampling filter is called a Type 2-2 frequency sampling filter. To develop Type 2-2 frequency sampling filters, the system function,
H(z1,z2),and the set, Hk(k)for k ~ RK,of frequency samples are related such that
124
PETERA. STUBBERUD
Hk (k) = IHk (k) Ie jO(k)
= H(Zl,Z2)lZl =e j(2rc/N1 )(kl +l/2),z2=eJ(2rc/N2 )(k2 +1/2) This relationship modifies the DFT into the form
"~, h(n)e-J(2r~/N1 )(kl +1/2)nl e-J(2n/N2 )(k2 +1/2)n2
Hk(k)
(7)
n~_RN where the modified IDFT is Z Hk(k)eJ(2x/N1 )(kl +1/2)n I eJ(2rr,/N2)(k 2 +1/2)n2 . (8)
h(n) = ~
NIN2 keRK Using a development similar to the one used for deriving the system function for Type 1-1 frequency sampling filters, Equation (8) can be substituted into the filter's z transform so that the system function, H(Zl,Z2), of a Type 2-2 frequency sampling filter can be written as 1+ Zl NI 1+ Z2 N2
H(Zl,Z2)=~
~ • N2
NI
NI-I N2-1
Hk(kl'k2)
k~-o k2-0 1-
zi-
l-
(9)
2)zl
The frequency sampling filter in Equation (9) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 2-2 frequency sampling filter with a real impulse response has frequency samples of the form
IHk(k,,k2)l-Ink(N1 - 1 -
kl,N 2 - 1 - k2) I
O(kl,k2) =-O(N 1 - 1 - kl,N 2 - 1 - k2). Using a development similar to the one used for deriving the system function for linear phase Type 1-1 frequency sampling filters, the system function of a linear phase Type 2-2 frequency sampling filter can be written as
H(zl,z2) =
1+
Z? ul
1+
Z2N2 f~,
l ]~,
z)
DESIGN FOR 2-D FREQUENCYSAMPLINGFILTERS
+Z
125
(,+ z~')(,- 2c~(k'LN, + ,, 2,]z l + z 21
k1=0 NI-1 M2-1
+Z
Z (-1)k'+k22[Hk(kl'k2)l x
kl=0 k2=0
c~
,~ l~ltz-1I + I) - ~-~-2 (k2 + -~1_],
1 +z2)
(I-2cos[-~-i (kl + l)]zll + Zl 2)(I-2cos[-~-2 (k2 + l)]z21 + z22) c~
+ 1 ) + ~ 2 (k2 +1)](1 +zllz21)
(I- 2cosE--~-i(kl + I 'lz-I
- 2cos[-~-2(k2 1'Iz-I
(lo)
where N 1 and N 2 are odd. The frequency sampling filter in Equation (10) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 2-2 frequency sampling filter with fourfold symmetry has a DFT of the form
k2)l =[Hk(N1-1-kl,k2)l
IHk(kl, k2)l- IHk(kl,N2
-
= In k (N 1 -
1 -
1 - k 1, N 2 - 1 - k2)1.
Substituting these constraints into Equation (10) yields,
H(Zl,Z2) = ml_l
+Z
k1=0
1 + Zl N1 1+ Z2N2 f_ ml
N2
( l + z l l ) ( l + z 2 l)
(_l)klNlHk(kl,M2)]cosI_~ll(kl
k'LNI
__ + 1 ) _ ~N2]( 12 +zll)
126
PETERA.STUBBERUD
M2_1 (-1)k22[Hk(ml,k2)lcosI--~2(k2 + 89
+2;
+Z21)
M1 M2
+Z
Z (-1)kl+k24lSk(kl'k2)[x
k1=0 k2 =0 sin[-~ll (kl + 89 sin[-~22 (k2 + 89
+ Zl 1)(1 + z21 )
(1-2COS[-~1(kl +l)]zll +Zl2)(1- 2cosI-~-2(k2 + 89
2)
(11)
where N 1 and N 2 are odd.
C. Type 1-2 Frequency Sampling Filters Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (NIN2) frequency samples starting at 091 = 0 and c02 - g/N2. This type of frequency sampling filter is called a Type 1-2 frequency sampling filter. To develop Type 1-2 frequency sampling filters, the system function, H(Zl,Z2), and the set, Hk(k ) for k ~ R K, of frequency samples are related such that Hk(k ) = IHk(k)leJO(k)
: H(zl,z )lz,
=eJ(2n/N1)kl =eJ(2r:/N2)(k2+1/2)
This relationship modifies the DFT into the form
Hk(k) = Z h(n)e-J(Zx/ul )klnle-J(Zx/N2)(k2+1/2)n2 n~_RN where the modified IDFT is h(n)= 1-----~ ZHk(k)eJ(2x/N1)klnleJ(2x/N2)(k2+l/2)n2.
(12)
N1 N2 k~RK Using a development similar to the one used for deriving the system function
DESIGN FOR 2-D FREQUENCY SAMPLINGFILTERS
12/
for Type 1-1 frequency sampling filters, Equation (12) can be substituted into the the filter's z transform so that the system function,
H(z 1,z2), of a Type 1-2
frequency sampling filter can be written as
H(Zl,Z2) =
1- Zl N1 l + z 2 N2
N1 N 1-1 N 2 -1 kl=0 k2=0
N2
•
(1-eJ(2rc/N1)klzll)(1-eJ(2rc/N2)(k2+l/2)z21)
The frequency sampling filter in Equation (13) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 1-2 frequency sampling filter with a real impulse response has frequency samples of the form
[Hk(kl,k2) ]= [Hk(N 1 - kl,N 2 - 1 - k2) ] O(kl,k2) = -O(N 1 - kl,N 2 - 1 - k2). Using a development similar to the one used for deriving the system function for linear phase Type 1-1 frequency sampling filters, the system function of a linear phase Type 1-2 frequency sampling filter can be written as
1- z7 N1 1 + z2 N2 f
+
(1- 1)(l_ cos[
Hk(O'M2 )
+ )1z;1+
M1 N 2-1
+Z Z
k1=1 k2 =0
(_l)kl +k2 2]Hk(kl,k2)]•
)
128
PETER A. STUBBERUD
where N 1 and N2 are odd. The frequency sampling filter in Equation ( 1 4 ) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 1-2 frequency sampling filter with fourfold symmetry has a DIT of the form I ~ k ( k l , k : ! ) l = I ~ k ( k l- ,1~- k2z ) I = IHk (
.
~ -1kl k2 I)
= ( H ~ ( N- Ik l , N 2 - 1 - k 2 ) I .
Substituting these constraints into Equation ( 1 4 ) yields,
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS M1
129
M 2-1
Z (-1)kl +k2 41Hk(kl'k2)l•
+Z k 1=1
k2 =0
c o s ( k k, )sin[~2-2 (k2+ 89
z~-1)(1+ z21)
[1-2 cos(-~l k1) z l l + Zl 2 ](1-2 c o s I ~2 (k2+ 89
(15) + z2 2 )
where N 1 and N2 are odd. D. Type 2-1 Frequency Sampling Filters Another type of frequency sampling filter can be designed by interpolating a frequency response through a set of (N1N2) frequency samples starting at 091 - ~N 1 and co2 - 0. This type of frequency sampling filter is called a Type 2-1 frequency sampling filter. To develop Type 2-1 frequency sampling filters, the system function, H(Zl,Z2), and the set, Hk(k) for k ~ RK, of frequency samples are related such that H k (k) = IHk (k)[ e jO(k)
= n(Zl,Z 2 )lzl=eJ(2rc/N1)(kl +l/2),z2=eJ(2rt/N2)k 2 " This relationship modifies the DFT into the form Hk(k) = Z h(n)e-J(2rt/N1 )(kl +1/2)n 1e-J(27r,/N2)k2n2 n~R N
where the modified IDFT is h(n) = ----~1 ZHk(k)eJ(2rr'/N1)(kl+l/2)nle j(ert/Na)kan2 .
(16)
N1N2 keRK Using a development similar to the one used for deriving the system function for Type 1-2 frequency sampling filters, Equation (1 6) can be substituted into the the filter's z transform so that the system function, H(z 1,z2), of a Type 2-1 frequency sampling filter can be written as
130
PETERA. STUBBERUD
H(Zl,Z2) =
1+ Zl N1 1- z2 N2 N1 N2
N 1-1N2-1
Z
Z
:o
:o
Hk(kl'k2)
(1
+l ' 'Zl II(1
z;l )
(17)
The frequency sampling filter in Equation (17) can be expressed in a computationally more efficient form by constraining the filter to have a real impulse response and linear phase. A Type 2-1 frequency sampling filter with a real impulse response has frequency samples of the form IHk(kl,k2)l = [ H k ( U l
- 1
- kl,N 2 -
k2) I
0(k 1, k2 ) = - 0 ( N 1 - 1 - k 1, N 2 - k2 ). Using a development similar to the one used for deriving the system function for linear phase Type 1-2 frequency sampling filters, the system function of a linear phase Type 2-1 frequency sampling filter can be written as
H(Zl,Z2) = 1+ Zl ul 1-Z2 N2 f _ MI-1
_
__
(-1)kl2]Hk(kl,O)lsinI-~ll(kl +1)]( 1 + z ?
1)
~_~cosI~(~l +l)z-,+Zl ~- (l-z;,) [.NI 1
N1-1 M2 + Z Z (-1)k' +~2 21nk(kl,k2)I• k1=1 k2 =0
sin[-~l (kl ( 1 - 2 cosI--~-1 (kl+ 89
_f (1 - 2 c o s
+ 89
)
+ Zl 2 ) I 1 - 2 c o s ( ~ 2 k2/z21+ z2 2 ]
sin[~ll(kl + 89 k21(zll -z21.! k + z +z 1 2cos I2~l-~l-~-~l';-1S)~----(~2k21z21+z221
where N 1 and N 2 are odd.
(18)
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
131
The frequency sampling filter in Equation (18) can be expressed in a computationally more efficient form if we also constrain the filter to have fourfold symmetry. A Type 2-1 frequency sampling filter with fourfold symmetry has a DFr of the form
IHk(kl' k2 )l = [Hk(kl' N2 - k2 )l =[Hk(N1- 1 - kl,k2)l =[Hk (N 1 - 1 - kl, N 2 - k2)[. Substituting these constraints into Equation (18) yields,
H(Zl,Z2)= 1+ Z?N1 1-Z2 N2 f _
__
M,-, (-1' kl 2[Hk(kl,O)[sin[-~l(kl+ l ) ] ( l + z l l )
+Z
kl=O (1-2cos[2~(kl+l)lzll LN1
+ z12 ) (1 - z21 )
+Z MI-1 M 2
+ Z Z (-1~' +~ 41"~(~1,~1 • k 1 =1 k 2 = 0
sin[~ll (kl+ 1)] c o s ( ~ 2 k2)(1+ Zl 1)(1- z21)
(19)
where N 1 and N 2 are odd.
E. ComputationalAdvantage of Frequency Sampling Filters Each of the frequency sampling filter system functions described by Equations (3), (5), (6), (9)-(11), (13)-(15) and (17)-(19) describe an architec-
132
P E T E R A. STUB B E R U D
ture that has a nonrecursive comb filter created by the terms, Ill+ z,N1 \] and (l + z2N2 ) , in series with a number of parallel resonators 1. When ~' "most! of a frequency sampling filter's frequency samples, the Hk(k)'s, are exactly zero, most of the filter's resonators are not implemented. Therefore, in the case of a narrowband filter where most of the frequency samples are in the stopband and are set equal to zero, only a small number of the filter's resonators need to be implemented. In such instances, frequency sampling implementations usually require fewer arithmetic operations than the direct convolution implementations. Each of the resonators in the frequency sampling filters described by Equations (3), (5), (6), (9)-(11), (13)-(15) and (17)-(19) can be represented by a separable system function 2. If separable system function, H(Zl,Z2), can be written as H(Zl,Z2) - HI(Zl)H2(z2), then it can be shown that that system is stable iff each of the 1-D system functions, HI(Zl) and H2(z2), are stable[2]. Because each of the resonator system functions are separable, it is straight forward to observe that for each of the frequency sampling filters the zeroes createdby the terms, (I _+zlNI ) and (l _+z2N2 ), exactly cancel the resonators ' poles that lie on the 4-D unit sphere. When implementing these filters with finite word length systems, exact pole zero cancellation generally does not occur, and uncanceled poles on the unit sphere will cause the filter to be unstable. To prevent this instability, each z.-1 in each of the frequency sampling filter system functions is replaced by rzi -~1",r < 1", i - 1,2. This moves all the poles inside the 4-D unit sphere. To keep the frequency response of these modified filters as close as possible to the frequency response of the original filters, r should be chosen near to the value 1. Substituting rz.-1 for z- 1 where r < 1 guarantees the stability of a frequent l cy sampling filter at the cost of increasing its computational requirements. If implemented with z i- 1 replaced by rz i 1, a linear phase frequency sampling filter requires at most 8 multiplies per resonator, and a linear phase frequency sampling filter with fourfold symmetry requires at most 5 multiplies per resonator If only K of the frequency samples are non-zero, then a linear phase fre1. The term resonator is used in this paper to denote a system which has its poles on or near the 4-D unit sphere or has its poles on or near the 2-D unit circle. 2. A system function is said to be separable iff it can be expressed as the product of one dimensional (l-D) system functions.
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
133
quency sampling structure requires approximately 8K multiplies per output sample, and a linear phase frequency sampling filter with fourfold symmetry requires approximately 5K multiplies per output sample. If a linear phase FIR filter with region of support RN is implemented using direct convolution, it requires approximately N1N2/2 multiplies per output sample; and if a linear phase FIR filter with fourfold symmetry and region of support RN is implemented using direct convolution, it requires approximately N1N2/4 multiplies per output sample. Because multiplications typically require more time to compute and are more complex to implement than adds, a frequency sampling structure becomes computationally more efficient than a direct convolution implementation when it uses fewer multiplies. Therefore, if we wish to implement a linear phase filter with a region of support R N, a frequency sampling filter can implement the filter more efficiently (in the sense of fewer multiplies) than a direct convolution implementation when 8K < N1N2/2 or
K < N1N2/16. If we wish to implement a linear phase filter with fourfold symmetry and region of support RN, a frequency sampling filter can implement the filter more efficiently than a direct convolution implementation when 5K < N1N2/4 or K < N1N2/20.
III. The Design of 2-D Frequency Sampling Filters It was shown in section II that a frequency sampling filter approximates a desired frequency response by interpolating a frequency response through a set of frequency samples. In this section, a design technique that controls interpolation errors between samples and approximates a desired frequency response is developed for frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. This design technique approximates a desired frequency response by minimizing a weighted mean square error over the passbands and stopbands subject to constraints on the filter's amplitude response. To maintain the computational efficiency of the filter, this error criterion is minimized while constraining the frequency samples in the stopband to be identically zero. This design technique is developed as a constrained optimization problem which is solved by the method of Lagrange multipliers. The Lagrange multipliers optimization method results in a set of linear equations which are solved to determine the filter's
134
PETER A. STUBBERUD
coefficients. The design method applies to each of the four types of 2-D frequency sampling filters that have real impulse responses, linear phase and fourfold symmetry. The frequency response of a frequency sampling filter which has an impulse response, h(n) for n ~ R N, is Hco ( tO) = H ( eJ(-Ol , eJCO2 ) =
h(n )e-J~ nl e-JOJ2 n2
n~R N If the filter has a real impulse response, fourfold symmetry and linear phase then h(n 1, n 2) = h(n 1,N 2 - 1 - n 2) = h ( N 1 - 1 - n 1, n 2 ) = h(N 1 -1-n
1,N 2 - 1 - n
2),
and it can then be shown that Hco(a)) can be written as H(o (o)1, (.o2) = e-J[~176 where Hr((_o 1,(.o2) is a real function given by M1 M2 Z Z 2sgn(nl)2sgn(n2)h(Ml-nl'M2-n2)cOs((-Olnl)cOs(OJ2n2) n 1=0 n2 =0 M 1 - (N 1-1)/2, M 2 - (N 2-1)/2, N 1 is odd, N 2 is odd and sgn(o) is the signum function. Hr(O~) is the amplitude of Hoj(m) and is often referred to as the zero phase frequency response or the amplitude response of the filter. If we let Hr(t~
x=[h(MI,M2)
2h(MI-I,M2)...
4h(1,0)
4h(0,0)] r
where the superscript T denotes transpose, and let s(col ,co2 )=[Sr (oJ1 ,co2 )]=[cos(al (r)col)cos(a2 (r)co2)] where r = 0, 1 .. (MI+I)(M2+I) - 1 s (co) represents the rth element in the column vector, s(m), '
"
'
~
r
a I ( r ) = r m o d ( M 1 +1)
a2 (r)=int( M-~+I )
DESIGNFOR 2-D FREQUENCYSAMPLINGFILTERS
135
and int(x) truncates x, then
Hr(tO) = xTs(to) = s T (to)x. (20) In Equations (1), (8), (12) and (16), the impulse response, h(n), of each of the frequency sampling filters is expressed as a function of its DFT. By constraining the filters to have linear phase, the IDFT's can be written as Z 2sgn(kl+k2)(_l)(kl+k2) N1N2 [H~(kl'k2)[ • k~Rll cos[2~ kl(nl + 89 27r )] LN1 ~22 k2(n2 + 89 Z k6R12
2sgn(kl+M2_k2)(_l)(kl+k2)
N1N2
Type 1-1
Igk(kl,k=~l•
sinF2Zc 2zr (k2 +l)(n2 + 89)] LN 1 kl(n 1 + 89 "~2 h(n) =
~, k6R21
2sgn(Ml-kl+k2)(-1)(kl+k2)
N1N2
IH~(kl,k=~l•
sinF~ (~, +~)(n,+~)+ ~ Z
k6R22
2sgn(M1 -k 1+ M2 -k 2 ) (_ 1)(kl +k2 + 1)
N1N2
cos[~ (~1+~)(nl+~)+ ~
Type 1- 2
)]
Type 2-1
IHk(kl,k2~l•
)1
Type 2- 2
where R l l = { ( k l , k 2 ) ' 0 < k I <_N1-1,1_< k2
IA IA
~
,=...
IA
~
~
IA
o.
il
I
!
i
•
,-
+
7"
~,..-
I
I
~.
~-
~I--
,
+
i
, t'4
m
I
i
~
i
x
+ I
~
~
~
*
+
+
i
~
b~
I,,./
, t'4 Ill
I
I
II
-I~I--
+
I
I
m
x
'~
b~
I
~1-"
m
x
,~
b~
b~
t'rl
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
2sgn(kl )2sgn(k2)(_l)(kl +k2)
N1N2
137
x
cosFIN 1 l(n + 2" cos "~2 k2 (n2 + 89
y el-1
2 sgn(kl)2sgn(M2-k 2) (_ 1)(kl +k2 ) x
N1N2
COSI2' klLUl (nl+l]]sin[ ( k 2 z ' LNe' . ]27r
f(n) =
2 sgn(M1-kl )2sgn(k2)(_ 1)(kl +k2 )
N1N2
Type 1-2 + 1)(n2 +1)]
•
sin[2n'(klLN 1, + 89 nl +l)]c~ n 2 z[_N2, j
+1)]
Type 2-1
2sgn(Ml -kl !2sgn(M2--k2_!(-1)(kl +k2! • N1N2 l
2re sin[ "~1"1(kl + ~)(nl + 2")lsln[ "~2-2(k2 + 1)(n2 + 1)] 9
2zr
1
1
Type 2- 2
9
for r - 0 , 1, (MI+I)(M2+I) - 1; k 1 = int[r/(M2+l)]; and k2 = rmod(M2+l), then h(n) = xTf(n)= fT (n)X where the appropriate expression for f(n) is used depending upon the type of frequency sampling filter being designed. If we also define F-[f(M1,M2)
2f(M 1 - 1,M2) -.. 4f(1,0) 4f(0,0)] T
then x - FX
(21)
Substituting (21) into (20),
Hr(t~) - xTFTs(to) - s T (to)FX. (22) An approximation of a linear phase filter with a desired passband amplitude, Hal(to), can be obtained by minimizing the mean square error between the desired amplitude response, Hd(tO), and the FIR filter's amplitude response, H(~). If we let
138
PETER A. STUBBERUD
~pb
= {(COl,CO2)" COl E passband, Co2 ~ passband},
and let J b represent the mean square error over the region of passband frequencies, then
Jpb = m(~pb~1
SftoE~pbIHr (to) -
Hd (to)l 2 dto
(23)
where m(o) is the 2-D Lebesgue measure. Substituting Equation (22) into Equation (23), Jpb can be written as
1pb-"----~x T F T Jpb (X) = m(~"2
SS s(~
( ~ ) d ~ FX
~E~pb _2xTF T
SS Hd(~)s(~)dm + ~-~ pb
S H (m)dr
(24)
mE~ pb
If we let W(m) = s(m)s r ( ~ )
then w(m) = [Wrc(O)] = [cos(a 1(r)co 1)cos(a 2 (r)co 2) cos(a I (c)Co1)cos(a 2 (c)Co2)] where Wrc(m) represent the element in the rth row and the cth column of the matrix, W(to), andr, c - 0 , 1..... (MI+I)(M2+I) - 1. Defining Q(m) as the matrix,
the element in the rth row and the cth column of the matrix, Q(~), can be written as
Qrc ( C~ 0)2 ) - Qlc (col)Q2c ( o92 ) where
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
ai(r)=ai(c)=O
Ogi i Qrc ( Ogi) =
Ogi
139
sin 2ai(r)O9i
2 4ai(r) sin[ai(r) + ai(c)]o9i sin[ai(r ) - ai(c)]O9i 2[ai(r)+ai(c)] + 2[ai(r ) - ai(c)]
ai(r)=ai(c)•O ai(r) r ai(c)
and r, c = O, 1..... (M 1+ 1)(M2+ 1) - 1. If we define gg
QP
= JJm ~ p b W(m) dm
where Qp is calculated by evaluating
and
Q2 (O92)= [Q2c (o92)] at the appropriate limits and performing an element by element product of the resulting matrices, and also define the terms,
R(m) = f l Hd (m)s(m)do~
~(~)=
~IH2(~)doJ
Rp=
II Hd(m)s(m)dm
~'p=
IIH2(o~)dm
O~pb then Equation (24) can be written as
Jpb(X) =
1.......~-[xTFTQpF x _ 2xTFTRp + )'p] m(~pb)
(25)
For filters which approximate constant values in their passbands, the expression for R(m) can be simplified. If Hd(O~) is equal to a constant in the passband, then without loss of generality, we can let Hd(m) -- 1, and Equation (25) becomes
Jpb(X) =
I~[xTFTQpFX m(~pb)
- 2xTFTRp]
+1
(26)
Because the elements in the vector, s(~), are separable, the matrix, R(~), can be written as
140
PETER A. STUBBERUD
where
Coi Ric(f~ - sin[ai(r)ooi]/ai(r )
ai(r) = O ai(r ) r 0
This expression for R(m) eliminates the need for integrating Hd(m)s(m) when calculating Rp. The matrix, Rp, in Equation (26) can be calculated by evaluating Rl(O)l)= [Rlc(COl)] and R2(fo2)=[R2c(OO2)] atthe appropriate limits and performing an element by element product of the resulting matrices 9 An approximation of a desired amplitude response can also be obtained by imposing constraints on the amplitude response and its derivatives at particular frequencies 9 For Example, nr(m)lm=to0 = K000
o(m+n)Hr(~) I o300m~o)~
-..
Hr(m)lm=m t =
tg(m+ n) Hr (OJ) = gmn 0 ""
= Kmn I tO=O~1
~=~0
c92PHr(~) I OcOp OOoff
02 P Hr (O~) = K ppO
KO01
= Kpp l
"'"
~=~0
~=~1
where the K's are constants 9 Using Equation (22), these passband constraints can be written as
o(m+n)Hr(~ ) = o(m+n)s T (O~) FX.
oo,Too, If we let cT=Is(0~0 ) S(t~l )
o3(2P)s(0~)
o32s(t0) 00910(.02
ta=o~0
and KT=[K000
K001 ...
KI10
...
1
Kppl],
the passband constraints can be written in the matrix form, C F X - K. (27) If we let Jsb represent the mean square error over the stopband frequencies, then
D E S I G N F O R 2-D F R E Q U E N C Y S A M P L I N G FILTERS
141
1
Jsb - m(ns b) IIcoef~sb [nr(cO) - Hd(~)l 2 dm
(28)
where
~sb Because
= {(CO1,0)2)" CO1E stopband, 092 e stopband}.
Hd(m) equals zero in the stopband, Equation (28) can be written as Jsb = m(f2sb )
(29)
~-f~sb
Substituting Equation (22) into Equation (29), Jsb(X) =
1
xTF Tff,,,
s(to)s T(tO) dto
FX
(30)
Recall that earlier, we let W(to) -s(~)sT(to) and Q(m) = f ~ w(m) dto. Thus, if we define -
Q s IIco~f~sb
W(m)
dm
where Qs is calculated by evaluating Q l(O~l) and Q2(co2) at the appropriate limits and performing an element by element product of the resulting matrices, Equation (30) can be written as Jsb(X) =
~m----xTFTQ s FX. m(asb)
(31)
To generate an efficient realization of a linear phase FIR filter with a frequency sampling structure, the frequency samples in the stopband are constrained equal to zero. Recall that the vector, X, contains the frequency sampies. In Equations (25), (27) and (31), constraining frequency samples to zero is equivalent to removing elements in the vector, X, that contain the corresponding frequency samples and also removing the corresponding columns from the matrix, F. If we define these reduced terms as X and F C
C~
then the
design problem can now be stated as follows. Minimize the error function
142
PETER A. STUBBERUD
J(X c) = t~pb (Xc) + (1 - a)Jsb (X c) =
~
1-o~
m(~sb)
XcT FcT Qs Fc x m(~pb)
c
[Xc c Q cXc
where 0 < ~ < 1 subject to the constraint, CF X - K. The scalar term, a, in C
C
Equation (32) allows the designer to weight the relative importance between the mean square error over the passbands and the mean square error over the stopbands. If the number of non-redundant constraints in Equation (27) is less than the number of elements in the vector X c, then the remaining degrees of freedom are used to minimize J(Xc). If the number of non-redundant constraints is equal to the number of elements in the vector X c, then the design becomes a linear algebra problem and the frequency samples are unique for the constraints. And if the number of non-redundant constraints is greater than the number of elements in the vector X c, then no solution will exist which satisfies all the constraints. Assuming that the number of non-redundant constraints in Equation (27) is less than the number of elements in the vector X c, the design problem can be solved by defining a Lagrange multiplier vector as ~T =1.[zOo0
zoo1
"'" '~'110 "'"
Zppl]j
and minimizing the augmented cost function,
Ja(Xc,~) = OClpb(Xc) + (1 - a)Jsb(Xc) + ~T(CFcXc - K).
(33)
Substituting the appropriate expressions for %b and Jsb into Equation (33),
Ja(Xc ~ ) = ~ [ X ' m(~pb)
T T c Fc QpFcXc - 2 x T F T R p + Yp]
1-a XcT FcT Qs FcXc + ~T (CFcXc _ K) m(~2sb ) The necessary conditions for an optimal solution are ,)J a ( X c , ~.)
O~c
=0
(34)
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
143
and
tgJa ( X c , ~,)
Ok
(35)
"-0.
Because F c TQpFc and Fc TQsFc are symmetric matrices, Equation (34) becomes 2 ( 1 - a) FTQs Fc X + FTcT2~ = 0. (36)
m(~sb)
m(f2pb)
c
Equation (35) implies
OJa(Xc,~)
ovk
= CFcX c - K = 0.
(37)
Equations (36) and (37) can be written in the matrix form,
I
2a
m(~pb )
0 lrjL-2j=,r r
FTQpFc + 2 ( l - a ) F t Q s Fc ,, F T c T
m(f2sb )
~
Xc
I m(f2Pb )
K
(38)
A solution to Equation (38) will exist when C has full rank which occurs when C contains no redundant or trivial constraints.
IV. Example A. Example of a Type 1-1 Frequency Sampling Filter In this example, a Type 1-1 linear phase fourfold symmetric frequency sampling filter that approximates the amplitude response, Hd (~ - {10
~~2Pb OiE~sb
where
f2pb - {(091,092) .0 < 091 -< 0.1r~,0 _<092 _<0.1~} ~sb - {(091,092)" 0-13rt <
091 -< rtU0.13rt < 092 < rt}
is designed. The filter will be designed to have an impulse response that has a region of support R N where
R N = {(nl,n2)" 0 < n 1 < 100,0 < n 2 < 100}. To approximate the passband response, the mean square error over the passband will be minimized subject to the constraints
144
PETER A. STUBBERUD Hr(to)lm-(O, - ,0) = 1
t?(m+ n) Hr (to)
=0
where n, m = 0, 1, 2, 3, 4 excluding m = n = 0. To approximate the stopband response, the scalar term, o~, which weights the relative importance between the mean square error over the passbands and the mean square error over the stopbands, will be set to 0.1. Substituting c~ = 0.1 into Equation (38) gives 0.2
m(~pb ) where
FTQpFc +
m(~sb )
s
~
FTcT
Xc ]L-~..j=
0.2
FTRp
m(~pb)
l.....
(39) .....
m(~pb)=(O.l~) 2 and m(~sb)=~2-(O.13~;) 2 .
To solve Equation (39), the matrix, C, must have full rank, and the nonzero frequency samples must be determined. Because the odd ordered partial derivatives of a linear phase filter evaluated at to - (0,0) are zero, constraining them is redundant, and they must be omitted. Thus, for the matrix, C, to have full rank, only the constraints Hr(to)]m=(O,O ) = 1
tg(m +n ) Hr ( to )
=0
where n, m - 0, 2, 4 excluding m = n - 0 are considered. The nonzero frequency samples include all of the frequency samples not in the stopbands. In this example, the stopband frequency samples are H(k) for k ~ {(kl,k2) 97 < k I _< 5 0 U 7 < k 2 < 50} which implies that the nonzero frequency samples are H(k) for k ~ {(kl,k2)" 0 < k 1 _< 6,0 _
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
145
arg[Ha~(fOl, co2 )] = -50091 - 50o92 .
V. Summary and Conclusions In this chapter, four types of f r e q u e n c y s a m p l i n g filters were developed. F o r each type of filter, s y s t e m functions were d e v e l o p e d for filters with arbitrary impulse responses and arbitrary frequency responses; filters with real impulse responses and linear phase; and filters with real i m p u l s e responses, linear phase and fourfold s y m m e t r y . Also, in this chapter, a technique for designing frequency s a m p l i n g filters with real i m p u l s e responses, linear phase and fourfold s y m m e t r y was developed. This design t e c h n i q u e a p p r o x i m a t e s a desired frequency r e s p o n s e by m i n i m i z i n g a w e i g h t e d m e a n square error over the stopbands and passbands subject to constraints on the filter's a m p l i t u d e response. This a p p r o x i m a t i o n technique results in a set of linear equations which are solved to d e t e r m i n e the filter's coefficients. This design m e t h o d applies to each o f the four types o f f r e q u e n c y s a m p l i n g filters d e v e l o p e d in this chapter.
Table 1. Nonzero Frequency Samples for the Example. H(0,0) - 1.0000000
H(2,3) - H(3,2) - 0 . 9 9 6 7 5 4 4
H(0,1) - H(1,0) - 0.9999942
H(2,4) - H(4,2) - 0.9926905
H(0,2) - H(2,0) - 0 . 9 9 9 8 7 4 4
H(2,5) - H(5,2) - 0.9701703
H(0,3) - H(3,0) - 0.9940741
H(2,6) - H(6,2) - 0.5011298
H(0,4) - H(4,0) - 0.9953709
H(3,3) -- 1.0131842
H(0,5) = H(5,0) - 0.9704475
H(3,4) - H(4,3) - 0.9991182
n ( 0 , 6 ) - H(6,0) -- 0 . 5 2 6 6 0 8 0
H(3,5) -- H(5,3) -- 0 . 9 7 7 7 6 2 4
H(1,1) - 1.0005610
H(3,6) - H(6,3) = 0 . 4 5 3 0 7 1 0
H(1,2) - H(2,1) - 1.0064317
H(4,4) -- 0.9915667
H(1,3) - H(3,1) - 1.0057280
H(4,5) - H(5,4) -- 0 . 9 7 1 9 1 5 4
H(1,4) = H(4,1) - 1.0035502
H(4,6) - H(6,4) -- 0 . 4 9 4 1 6 3 6
H(1,5) - H(5,1) -- 0 . 9 7 9 1 7 3 4
H(5,5) -- 0.9312007
H(1,6) - H(6,1) - 0.5152585
H(5,6) - H(6,5) -- 0.3930366
H(2,2) = 0.9951068
H(6,6) - 0.0884227
146
PETER A. STUBBERUD
1 0.8
~,
~, 0.6
0.4 0.2 O= /Z"
-0.2tr ~ O~( . 0.21r " ~ tra~"/#~ 0.6x ~/e)
0.2zr tr
-n:
0.6zt"
-0.6n: OX(,~&.ls~ga'p\e')
Figure 1. Magnitude Response of the Type 1-1 Frequency Sampling Filter in the Example. 0 -50
~
-100 -150 -200 0 ~0 . 4~0 ' 2' ~ ~ r_~r 0.61r ~
e)
~ ' <0.8'1r~
8~
O.
. . > ~ i ~ 0.2zr o
~ o3x
0.6~ ~ 0.4/1: ~$~6.1s~9\e')
Figure 2. Magnitude Response of the Type 1-1 Frequency Sampling Filter in the Example.
zc
DESIGN FOR 2-D FREQUENCY SAMPLING FILTERS
~"
_'g
147
0
o ~o -0.5
-1 0
0-047r~"~-..~
o~ ~
o.6~ ~
~'~'~~ O.06zr
0.08zr
'~'4,~,,,. 0 . 0 8 ~ r ~ ~ 0.04zr '"~ o"l x " ~ " i ~ O O.02zrcox(~6.1 s~px'e')
Figure 3. Magnitude Response of the Passband for the Type 1-1 Frequency Sampling Filter in the Example.
VI.
References
1. L.R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing, Prentice Hall, New Jersey (1975). 2. D.E. Dudgeon and R. M. Mersereau, Multidimensional Digital Signal Processing, Prentice Hall, New Jersey (1984). 3. J.V. Hu and L. R. Rabiner, "Design Techniques for Two-Dimensional Filters," IEEE Transactions on Audio and Electroacoustics, Vol. AU-20, No. 4, pp. 249-257 (1972). 4. J.S. Lim, Two-Dimensional Signal and Image Processing, Prentice Hall, New Jersey (1990).