# Design of Adjustable Square-Shaped 2D IIR Filters

Authors: Radu Matei
Journal: ISRN Signal Processing
 (2013)

---

## Abstract

This paper proposes an analytical design method
for two-dimensional square-shaped IIR filters. The designed
2D filters are adjustable since their bandwidth and orientation
are specified by parameters appearing explicitly in the filter
matrices. The design relies on a zero-phase low-pass 1D
prototype filter. To this filter a frequency transformation is
next applied, which yields a 2D filter with the desired square
shape in the frequency plane. The proposed method combines
the analytical approach with numerical approximations. Since
the prototype transfer function is factorized into partial
functions, the 2D filter also will be described by a factorized transfer
function, which is an advantage in implementation.


---

## Body

## 1. Introduction

Various design methods for 2D filters, both FIR and IIR, have been proposed by many researchers [1]. A frequently used design technique is based on a specified 1D prototype filter, whose transfer function is transformed using various frequency mappings, in order to obtain a 2D filter with a desired frequency response. Some relevant papers approaching 2D filter design using spectral transformations are [2–6]. A class of tunable 2D digital filters is discussed in [7]. The stability problem for 2D filters and stabilization methods are treated in papers like [8–11].Diamond filters have been used as antialiasing filters in the conversion between signals sampled on the rectangular sampling grid and the quincunx sampling grid. Different issues related to design methods for diamond filters were studied in [12–14]. In [15], another design method for diamond-shaped filters was derived, starting from discrete filter transfer functions; two types of filters were obtained, one with complex transfer function and another with zero-phase transfer function. The latter is particularly appropriate for image processing applications as the filters introduce no phase distortions. The technique developed in [15] uses the 2D filter specification in polar coordinates.In this paper an analytical design method is proposed for 2D adjustable zero-phase square-shaped filters, a larger class of filters which may be regarded as a generalization of the common diamond filter. Starting from a 1D prototype filter with factorized transfer function, the corresponding 2D filters are obtained by a particular 1D to 2D frequency mapping. The 2D filter will have a factorized transfer function, which is a useful feature in implementation. This work mainly focuses on presenting the proposed method and describes in detail the design steps. Several design examples are also provided. The typical applications of diamond filters were not approached here, as they are extensively treated in many other works. Image processing applications of the square-shaped filters with arbitrary bandwidth and orientation proposed here will be approached in further work.


## 2. Specification of Square-Shaped Filters

### 2.1. Frequency Plane Specification of Square-Shaped Filters

A particular case of a square-shaped filter is the standard diamond filter; its shape in the frequency plane is shown in Figure1(a). It is a square with a side length of π2, while its axis is tilted by an angle of φ=π/4 radians about the two frequency axes. Next we will consider the orientation angle φ about the ω2 − axis.(a) diamond filter; (b) wide-band oriented filter; (c), (d) wide-band oriented filters with orientations forming an angleφ=π/2; (e) square-shaped filter resulted as product of the above oriented filters; (f) rhomboidal filter.
(a)
(b)
(c)
(d)
(e)
(f)In this work a more general case is approached, that is, a 2D diamond-type filter with a square shape in the frequency plane but with arbitrary side length and axis inclination angle, as shown in Figure1(e). Next we refer to them as square-shaped filters, since they are more general than the diamond filter from Figure 1(a).The square-shaped filter in Figure1(e) is derived as the intersection of two oriented low-pass filters whose axes are perpendicular to each other, for which the shape in the frequency plane is given in Figures 1(c) and 1(d). Correspondingly, the square-shaped filter transfer function HS(z1,z2) results  as a product of two partial transfer functions HD1(z1,z2) and HD2(z1,z2), which we refer to as directional filters:
(1)HS(z1,z2)=HD1(z1,z2)·HD2(z1,z2).The frequency response ofHD2(z1,z2) is ideally identical to the frequency response of HD1(z1,z2), rotated by an angle of φ=π/2. Since this axis rotation implies the frequency variable change, ω1→ω2,ω2→-ω1, the transfer function HD2(z1,z2) can be derived from HD1(z1,z2) as
(2)HD2(z1,z2)=HD1(z2,z1-1).A more general filter belonging to this class is a rhomboidal filter, as shown in Figure1(f). In this case the two oriented low-pass component filters may have different band-widths and their axes are no longer perpendicular to each other.


### 2.2. 1D Low-Pass Zero-Phase Prototype Filter

An analog filter of orderN is described by the general transfer function in variable s:
(3)HP(s)=P(s)Q(s)=∑i=0Mpi·si∑j=0Nqj·sj.
A zero-phase prototype can be obtained from the general filter HP(s) if the magnitude characteristics |HP(jω)| are considered. Zero-phase filters are frequently used in image processing as they do not introduce any phase distortions in the filtered image. In order to obtain a maximally flat or low-ripple square-shaped filter, we start from a maximally flat 1D prototype.Let us consider a Butterworth low-pass filterHB(s) of order N, having the transfer function magnitude as follows:
(4)|HB(jω)|=11+(ω/ω0)2N,
where ω0 is the filter cut-off frequency.We look for a rational expression of the magnitude|HB(jω)|, which has to be an approximation as accurate as possible within a specified error on the frequency domain [-π,π]. The most convenient for our purpose is the Chebyshev-Padé expansion, since it yields an efficient approximation of a given function, which is uniform along the entire specified interval. It can be found using a symbolic computation software like MAPLE.Let us first find the approximation for the function (4) in normalized frequency, that is, the function |HBn(jω)|:
(5)|HBn(jω)|=11+ω2N.
Since a steeper filter would be desirable, let us take a Butterworth low-pass filter of order N=12; we obtain on the frequency range [-π,π] the following Chebyshev-Padé rational approximation of order 8 for the magnitude characteristics, as a ratio of polynomials with the same degree 8, which in factorized form looks like
(6)11+ω24≅HP(ω)=ξ·(ω4-5.8589322·ω2+22.393616)(ω4+0.253253·ω2+0.2830378)·(ω4-13.689786·ω2+49.701196)(ω4-1.843889·ω2+0.999034)=ξ·HP1(ω)·HP2(ω),
where the constant has the value ξ=0.0002516.The frequency responseHP(ω) of this prototype is plotted in Figure 2(a) and shows a small amplitude ripple in the pass-band. The cut-off frequency is ω0=1 and p=1. A low-pass filter with a better flatness would obviously have a higher order. The advantage of using this factorized rational approximation is that it is easily scalable on the frequency axis, as will appear clearly in the following sections.(a) Zero-phase low-pass prototype filter characteristic forp=1(ω0=1); (b), (c) scaled LP characteristics for p=0.6(ω0=0.53π) and p=2.7(ω0=0.118π).
(a)
(b)
(c)Each of the rational expressionsHPi(ω) (i=1,2,…) occurring as factors in HP(ω) given by (6), with 4th order polynomials as numerator and denominator, can be generally written as
(7)HPi(ω)=(ω4+b1i·ω2+b0i)(ω4+a1i·ω2+a0i)(i=1,2,…).
Of course, in a factorized rational approximation like expression (6) there may occur also a second-order term:
(8)HPi(ω)=(ω2+c0i)(ω2+d0i)
but this case is not further considered.As specified earlier, the variableω represents here the frequency normalized to the value of the cut-off frequency ω0. In order to return to current frequency values, in the expression of the general factor HPi(ω), we must substitute ω by ω/ω0. It is more convenient to substitute ω by p·ω, where p=1/ω0. Thus the ratio factor HPi(ω) from (7) can be rewritten in the following form which includes the parameter p, which makes it scalable along the frequency axis:
(9)HPi(ω)=p4ω4+b1i·p2ω2+b0ip4ω4+a1i·p2ω2+a0i.
For a value of the parameter p different from unity, the characteristic HP(ω) displayed in Figure 2 either stretches (for p<1) or shrinks (for p>1). In Figures 2(b) and 2(c) two LP characteristics derived from HP(ω) are plotted, for the indicated values of p and ω0, respectively.Thus the 1D prototype is parametric depending onp, and the derived 2D square-shaped filter will have an adjustable pass-band region, specified by the same p.


## 2.1. Frequency Plane Specification of Square-Shaped Filters

A particular case of a square-shaped filter is the standard diamond filter; its shape in the frequency plane is shown in Figure1(a). It is a square with a side length of π2, while its axis is tilted by an angle of φ=π/4 radians about the two frequency axes. Next we will consider the orientation angle φ about the ω2 − axis.(a) diamond filter; (b) wide-band oriented filter; (c), (d) wide-band oriented filters with orientations forming an angleφ=π/2; (e) square-shaped filter resulted as product of the above oriented filters; (f) rhomboidal filter.
(a)
(b)
(c)
(d)
(e)
(f)In this work a more general case is approached, that is, a 2D diamond-type filter with a square shape in the frequency plane but with arbitrary side length and axis inclination angle, as shown in Figure1(e). Next we refer to them as square-shaped filters, since they are more general than the diamond filter from Figure 1(a).The square-shaped filter in Figure1(e) is derived as the intersection of two oriented low-pass filters whose axes are perpendicular to each other, for which the shape in the frequency plane is given in Figures 1(c) and 1(d). Correspondingly, the square-shaped filter transfer function HS(z1,z2) results  as a product of two partial transfer functions HD1(z1,z2) and HD2(z1,z2), which we refer to as directional filters:
(1)HS(z1,z2)=HD1(z1,z2)·HD2(z1,z2).The frequency response ofHD2(z1,z2) is ideally identical to the frequency response of HD1(z1,z2), rotated by an angle of φ=π/2. Since this axis rotation implies the frequency variable change, ω1→ω2,ω2→-ω1, the transfer function HD2(z1,z2) can be derived from HD1(z1,z2) as
(2)HD2(z1,z2)=HD1(z2,z1-1).A more general filter belonging to this class is a rhomboidal filter, as shown in Figure1(f). In this case the two oriented low-pass component filters may have different band-widths and their axes are no longer perpendicular to each other.


## 2.2. 1D Low-Pass Zero-Phase Prototype Filter

An analog filter of orderN is described by the general transfer function in variable s:
(3)HP(s)=P(s)Q(s)=∑i=0Mpi·si∑j=0Nqj·sj.
A zero-phase prototype can be obtained from the general filter HP(s) if the magnitude characteristics |HP(jω)| are considered. Zero-phase filters are frequently used in image processing as they do not introduce any phase distortions in the filtered image. In order to obtain a maximally flat or low-ripple square-shaped filter, we start from a maximally flat 1D prototype.Let us consider a Butterworth low-pass filterHB(s) of order N, having the transfer function magnitude as follows:
(4)|HB(jω)|=11+(ω/ω0)2N,
where ω0 is the filter cut-off frequency.We look for a rational expression of the magnitude|HB(jω)|, which has to be an approximation as accurate as possible within a specified error on the frequency domain [-π,π]. The most convenient for our purpose is the Chebyshev-Padé expansion, since it yields an efficient approximation of a given function, which is uniform along the entire specified interval. It can be found using a symbolic computation software like MAPLE.Let us first find the approximation for the function (4) in normalized frequency, that is, the function |HBn(jω)|:
(5)|HBn(jω)|=11+ω2N.
Since a steeper filter would be desirable, let us take a Butterworth low-pass filter of order N=12; we obtain on the frequency range [-π,π] the following Chebyshev-Padé rational approximation of order 8 for the magnitude characteristics, as a ratio of polynomials with the same degree 8, which in factorized form looks like
(6)11+ω24≅HP(ω)=ξ·(ω4-5.8589322·ω2+22.393616)(ω4+0.253253·ω2+0.2830378)·(ω4-13.689786·ω2+49.701196)(ω4-1.843889·ω2+0.999034)=ξ·HP1(ω)·HP2(ω),
where the constant has the value ξ=0.0002516.The frequency responseHP(ω) of this prototype is plotted in Figure 2(a) and shows a small amplitude ripple in the pass-band. The cut-off frequency is ω0=1 and p=1. A low-pass filter with a better flatness would obviously have a higher order. The advantage of using this factorized rational approximation is that it is easily scalable on the frequency axis, as will appear clearly in the following sections.(a) Zero-phase low-pass prototype filter characteristic forp=1(ω0=1); (b), (c) scaled LP characteristics for p=0.6(ω0=0.53π) and p=2.7(ω0=0.118π).
(a)
(b)
(c)Each of the rational expressionsHPi(ω) (i=1,2,…) occurring as factors in HP(ω) given by (6), with 4th order polynomials as numerator and denominator, can be generally written as
(7)HPi(ω)=(ω4+b1i·ω2+b0i)(ω4+a1i·ω2+a0i)(i=1,2,…).
Of course, in a factorized rational approximation like expression (6) there may occur also a second-order term:
(8)HPi(ω)=(ω2+c0i)(ω2+d0i)
but this case is not further considered.As specified earlier, the variableω represents here the frequency normalized to the value of the cut-off frequency ω0. In order to return to current frequency values, in the expression of the general factor HPi(ω), we must substitute ω by ω/ω0. It is more convenient to substitute ω by p·ω, where p=1/ω0. Thus the ratio factor HPi(ω) from (7) can be rewritten in the following form which includes the parameter p, which makes it scalable along the frequency axis:
(9)HPi(ω)=p4ω4+b1i·p2ω2+b0ip4ω4+a1i·p2ω2+a0i.
For a value of the parameter p different from unity, the characteristic HP(ω) displayed in Figure 2 either stretches (for p<1) or shrinks (for p>1). In Figures 2(b) and 2(c) two LP characteristics derived from HP(ω) are plotted, for the indicated values of p and ω0, respectively.Thus the 1D prototype is parametric depending onp, and the derived 2D square-shaped filter will have an adjustable pass-band region, specified by the same p.


## 3. Design Method

The main issue in this section is to find the transfer function of the desired 2D square filterHS(z1,z2) using a particular frequency transformation.From a 1D prototype filterHP(s)=HP(jω) (which varies on one axis only), a 2D oriented filter may be obtained by rotating the axes of the plane (ω1,ω2) by an angle φ. The rotation is defined by the following linear transformation, where ω1,ω2 are the original frequency variables and ω-1,ω-2 the rotated ones [1]:
(10)[ω1ω2]=[cosφsinφ-sinφcosφ]·[ω-1ω-2].
The spatial orientation is specified by an angle φ with respect to ω1-axis, defined by the following 1D to 2D frequency mapping:
(11)ω⟶ω1cosφ+ω2sinφ.
The oriented filter transfer function results in the following:
(12)Hφ(ω1,ω2)=HP(ω1cosφ+ω2sinφ).
In the complex plane (s1,s2) the frequency transformation (11) becomes
(13)s⟶s1cosφ+s2sinφ.
The 2D oriented filter Hφ(ω1,ω2) has the magnitude section along the line ω1cosφ+ω2sinφ=0, identical with the prototype HP(ω), and is constant along the perpendicular line: ω1sinφ-ω2cosφ=0 (the filter longitudinal axis).In our case, since the rational prototype functionHP(ω) from (6) contains only even powers of frequency ω, that is, only powers of ω2, we will derive a convenient expression for Ωφ=(ω1cosφ+ω2sinφ)2, according to mapping (11):
(14)Ωφ=0.5sin2φ·(ω1+ω2)2+0.5(1+cos2φ-sin2φ)·ω12+0.5(1-cos2φ-sin2φ)·ω22.
Thus using (11) we derived the mapping ω2→Ωφ. Since s12=-ω12, s22=-ω22, the corresponding expression Sφ in the complex variables s1 and s2 is similar, with Sφ=-Ωφ:
(15)Sφ=(s1cosφ+s2sinφ)2=0.5sin2φ·(s1+s2)2+0.5(1+cos2φ-sin2φ)·s12+0.5(1-cos2φ-sin2φ)·s22.
Therefore, after applying the 1D to 2D directional mapping (11) by substituting ω2→Ωφ=-Sφ, we reach the following rational expression in the complex frequency variables s1,s2, which corresponds to each ratio factor HPi(ω) of the prototype frequency response, given by (9)
(16)FDi(s1,s2)=p4Sφ2-b1i·p2Sφ+b0ip4Sφ2-a1i·p2Sφ+a0i.
Since finally we must reach a transfer function of a discrete 2D filter in the complex variables z1 and z2, we must determine the discrete counterpart HDi(z1,z2) of the general factor FDi(s1,s2).The currently used method to obtain a discrete filter from an analog prototype is applying the bilinear transform. If the sample interval takes the valueT=1, the bilinear transform for s1 and s2 in the complex plane (s1,s2) has the form
(17)s1=2(z1-1)(z1+1),s2=2(z2-1)(z2+1).Even if this method is straightforward, the designed 2D filter, corresponding to the transfer function inz1, z2, will inherently present substantial linearity distortions towards the limits of the frequency plane as compared to the ideal frequency response. This fact is mainly due to the so-called frequency warping effect of the bilinear transform, expressed by the continuous to discrete frequency mapping:
(18)ω=2T·arctg(ωaT2),
where ω is a frequency of the discrete filter and ωa is the corresponding frequency of the analog filter. This error can be corrected by applying a prewarping. Taking T=1 in relation (18) we obtain the mappings
(19)ω1⟶2·arctg(ω12),ω2⟶2·arctg(ω22).
In order to include the nonlinear mappings (19) into the frequency transformation, a rational approximation would be needed. Using again the efficient Chebyshev-Padé method, we get the following accurate approximation on the range [-π,π]:
(20)arctg(ω2)≅0.4751·ω(1+0.05·ω2).
Thus the prewarping correction consists in the following substitution, also used in [15]:
(21)ω⟶0.95·ω(1+0.05·ω2).
Even if (21) is useful in general, taking into account that the 2D filter function FDi(s1,s2) has even parity in s1 and s2, that is, the variables appear as s12,s22, a more efficient pre-warping can be made in this case.We obtain the following Chebyshev-Padé approximation on the rangeω∈[-π,π]:
(22)(arctg(ω2))2≅(0.017251+0.218196·ω2)(1+0.116048·ω2).
Consequently, from (19) and (22), the frequency pre-warping mapping can be written, on both frequency axes:
(23)ω1,22⟶(0.069+0.872783·ω1,22)(1+0.116048·ω1,22),(24)s1,22⟶-(0.069-0.872783·s1,22)(1-0.116048·s1,22).
Let us now apply the bilinear transform on both axes; that is, we substitute relations (17) into mapping (24); we get the following simple mapping, valid on both axes ω1 and ω2:
(25)s1,22⟶6.3868·(z1,22-2.080651·z1,2+1)(z1,22+5.465361·z1,2+1).
The analog-discrete mapping (25) therefore is a corrected form of the bilinear transform including pre-warping in the entire region of interest, namely, ω1,2∈[-π,π]. It is noteworthy that this pre-warping does not increase the order. The mapping (25) can be written on each axis in terms of frequency ω:
(26)ω2⟶6.3868·(1.0403-cosω)(2.7326+cosω).
This mapping is plotted in Figure 3 versus ω2 as curve F1 in solid line and is almost linear at least on the range ω∈[0,0.7π], as compared to the mapping (21)—curve F2, plotted in dash line. Therefore the proposed correction compensates the distortions introduced by the bilinear transform. As shown further in the design examples, using the simple frequency mapping (25) which includes pre-warping, we obtain 2D square filters with good linearity along their axes. For the sum of frequency variables, ω1+ω2, the mapping (26) becomes
(27)(ω1+ω2)2⟶6.3868·(1.0403-cos(ω1+ω2))(2.7326+cos(ω1+ω2)),
and since cos(ω1+ω2)=0.5·(z1z2+z1-1z2-1) we obtain the mapping
(28)(s1+s2)2⟶6.3868·(z12z22-2.0806·z1z2+1)(z12z22+5.4653·z1z2+1).
The next step is to substitute the mappings (25) and (28) into the expression (15) of Sφ, and we obtain a rational expression of order 4 in the complex frequency variables z1 and z2. Thus the discrete approximation Fφ(z1,z2) of Sφ(s1,s2) may be expressed in the matrix form
(29)Sφ(s1,s2)⟶Fφ(z1,z2)=k·Mφ(z1,z2)N(z1,z2)=k·(z1×Mφ×z2T)(z1×N×z2T),Figure 3
Plots of mapping (26)—F1 in solid line, and mapping (21)—F2 in dash line.where the variable vectors arez1=[1z1z12z13z14], z2=[1z2z22z23z24], k=3.1934, and × stands for inner product.The numerator matrixMφ depends on the orientation angle φ according to the expression
(30)Mφ=M0+sin2φ·M1+cos2φ·M2,
where M0, M1, and M2 are the following matrices of size 5×5 with constant elements(31)M0=[-0.25-0.42307-0.2500-0.423071.47642-2.73526-1.366300-0.25-2.7352615.036-2.73526-0.250-1.36630-2.735261.47642-0.4230700-0.25-0.42307-0.25],M1=[0.125-0.26010.12500-0.2601-4.94993.473471.6263700.1253.47347-7.518033.473470.12501.626373.47347-4.9499-0.2601000.125-0.26010.125],M2=[0-0.943220000.943220-4.211690004.2116904.21169000-4.2116900.94322000-0.943220].The denominator matrixN is also a constant 5×5 matrix:
(32)N=[0.1250.683150.125000.683154.41674.41670.6831500.1254.416720.65464.41670.12500.683154.41674.41670.68315000.1250.683150.125].
From (16) and (29), after some algebra, the matrices Bi,Ai(i=1,2,…) corresponding to the numerator and denominator of FDi(s1,s2) result from similar relations:
(33)Bi=p4k2·Mφ*Mφ-b1i·p2k·Mφ*N+b0i·N*N,Ai=p4k2·Mφ*Mφ-a1i·p2k·Mφ*N+a0i·N*N,
where * denotes matrix convolution.The 2D directional IIR filter will have areal-valued (zero-phase) transfer function HD(z1,z2) in the complex variables z1=ejω1 and z2=ejω2, written in matrix form
(34)HD(z1,z2)=BD(z1,z2)AD(z1,z2)=z1×BD×z2Tz1×AD×z2T,
where the vectors z1,z2 are z1=[1z1⋯z1N-1z1N], z2=[1z2⋯z2N-1z2N] and N is the directional filter order. If B1,B2,A1, and A2 are the matrices corresponding to the numerators and denominators of the two factors from (6), the matrices BD and AD of the directional filter result  as convolution:
(35)BD=B1*B2,AD=A1*A2.
Finally, the square-shaped filter transfer function is given by an expression similar to (34):
(36)HS(z1,z2)=BS(z1,z2)AS(z1,z2)=z1×BS×z2Tz1×AS×z2T,
where the matrices BS and AS result from
(37)BS=BD*BDR,AS=AD*ADR,
where the superscript R means rotation by 90 degrees; therefore ADR is the clockwise rotated version of AD. This is due to the fact that the matrix pairs (AD,BD) and (ADR,BDR) correspond to two directional filters which form a right angle in the frequency plane.To summarize, the design steps for the adjustable square-shaped filter are the following.(a)
Once adopted an 1D low-pass filter prototype of the form given by (4), its 4th order factors are determined; therefore the coefficients a0i,a1i,b0i, and b1i are found; if a second-order factor occurs, c0i and d0i are also found.(b)
For a specified orientation angleφ, the matrix Mφ is determined using expression (30) and the partial constant matrices (31).(c)
For a desired filter bandwidthω0, specified by the parameter p=1/ω0, the partial matrices Bi and Ai for each of the 4-order factors result from relations (33).(d)
The matricesBD and AD of the directional filter result as convolution of the matrices determined at the previous step, according to  (35).(e)
The matricesBS and AS of the square-shaped filter result as convolution of the matrices determined at the previous step, according to  (37).


## 4. Design Examples

Following the previous steps, the design of a 2D square-shaped filter using this method is straightforward. The bandwidthω0 is specified by the parameter p=1/ω0. As shown before, the 5×5 matrices M0, M1, M2, and N are constant once specified the 1D prototype HP(ω), in our case given by (6).The first design example is the diamond-type filter whose frequency responseHS1(ω1,ω2) is shown in Figure 4(a). In order to obtain a diamond filter as in Figure 1(a), the orientation angle must be set φ=π/4. Ideally, the bandwidth (cut-off frequency) should be ω0=π/2, corresponding to p=0.45. However, due to the inherent errors and the finite slope of the filter transition region, using these values would lead to a filter with large distortions towards the margins of the frequency plane. Using instead the value p=0.6, corresponding to a bandwidth ω0=0.53π, the filter will have a slightly narrower bandwidth but also smaller distortions. Also the contour plot shows a satisfactory linearity on the four sides of the square-shaped section. In Figures 4(c) and 4(d), the frequency responses and contour plots are displayed for the two directional component filters HD1(z1,z2), HD2(z1,z2). The square-shaped filter characteristic results as a product of the two directional characteristics.(a) Frequency response of a square-shaped filter with orientationφ=π/4 and p=0.6 (bandwidth ω0=1.66≅0.53π); (b) contour plot of this filter; (c), (d) frequency responses and contour plots of the two partial oriented filters.
(a)
(b)
(c)
(d)In Figures5 and 6, the frequency responses and contour plots are shown, for other square-shaped filters with different specifications, that is, bandwidth and orientation. It can be noticed that the filter characteristics have satisfactory steepness and relatively low ripple distortions in the pass-band region. Of course, better characteristics can be obtained using a higher order 1D prototype filter, that is, a more accurate approximation of the Butterworth LP filter. This however would increase the order of the 2D square-shaped filter.(a) Frequency response of a square-shaped filter with orientationφ=3π/8 and p=0.54 (bandwidth ω0=1.85≅0.589π); (b) contour plot of this filter; (c), (d) frequency response and contour plot of one oriented filter.
(a)
(b)
(c)
(d)(a) Frequency response of a square-shaped filter with orientationφ=3π/8 and p=1.6 (bandwidth ω0=0.625≅0.2π); (b) contour plot of this filter; (c), (d) frequency response and contour plot of one oriented filter; (e), (f) frequency response and contour plot for a square-shaped filter with φ=π/3 and p=1 (bandwidth ω0=1≅0.318π).
(a)
(b)
(c)
(d)
(e)
(f)A more general filter of this kind can also be derived, namely, a rhomboidal-shaped filter, in which the two component filters are not in a right angle to each other. Furthermore, a directional filter like the ones in Figure4(c) or Figure 6(c) can be used independently, especially if it has    a narrow bandwidth and can be used to extract lines with different spatial orientations from an image by tuning it to specified directions.Even if this analytical design method may not yield optimal filters in terms of complexity, as is generally the case with techniques based on numerical optimization, the advantage of having parametric, closed form expressions for the filter matrices is the possibility of tuning or adjusting the 2D filter for a new set of specifications—in our case bandwidth and orientation—without the need to resume the design process every time.Stability of the proposed filters is also an important issue. The stability problem for 2D filters is much more difficult than for 1D filters. There exist various stability criteria [7, 8] and also stabilization methods which can be applied for some unstable filters [9, 10]. A full stability analysis of the proposed adjustable 2D filters will be approached in further work.The main goal of the paper was to present in detail this analytical design method and some relevant design examples for these filters. Since they can select adjustable square or rhomboidal regions in the frequency plane, they may have interesting applications in image processing. Their applications in image filtering will be investigated in future work on this topic.


## 5. Conclusions

An analytical design method for recursive zero-phase square-shaped filters was proposed. The design is based on analog zero-phase low-pass prototypes with specified parameters. As 1D prototype, a Butterworth filter was used, and, using an efficient rational approximation, a zero-phase filter small pass-band ripple was derived. The resulted 2D filter is parametric or adjustable in the sense that the specified parameters—bandwidth and orientation angle—appear explicitly in the final filter matrices. The method includes the bilinear transform applied along the two axes and uses a frequency pre-warping to compensate for distortions. This leads to a particular frequency mapping, which is applied to the 1D prototype in order to obtain the 2D filter. This design technique is relatively simple, efficient and versatile, in the sense that, by changing the specifications, the new 2D filter matrices result directly, without the need to resume the entire design process all the way from the start.



---

*Source: 796830-2013-09-11.xml*