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.
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 Figure 1(a). It is a square with a side length of , while its axis is tilted by an angle of radians about the two frequency axes. Next we will consider the orientation angle about the − axis.
(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 Figure 1(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 Figure 1(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 results as a product of two partial transfer functions and , which we refer to as directional filters:
The frequency response of is ideally identical to the frequency response of , rotated by an angle of . Since this axis rotation implies the frequency variable change, , the transfer function can be derived from as
A more general filter belonging to this class is a rhomboidal filter, as shown in Figure 1(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 order is described by the general transfer function in variable : A zero-phase prototype can be obtained from the general filter if the magnitude characteristics 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 filter of order , having the transfer function magnitude as follows: where is the filter cut-off frequency.
We look for a rational expression of the magnitude , 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 : Since a steeper filter would be desirable, let us take a Butterworth low-pass filter of order ; 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 where the constant has the value .
The frequency response of this prototype is plotted in Figure 2(a) and shows a small amplitude ripple in the pass-band. The cut-off frequency is and . 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)
(b)
(c)
Each of the rational expressions () occurring as factors in given by (6), with 4th order polynomials as numerator and denominator, can be generally written as Of course, in a factorized rational approximation like expression (6) there may occur also a second-order term: 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 . In order to return to current frequency values, in the expression of the general factor , we must substitute by . It is more convenient to substitute by , where . Thus the ratio factor from (7) can be rewritten in the following form which includes the parameter , which makes it scalable along the frequency axis: For a value of the parameter different from unity, the characteristic displayed in Figure 2 either stretches (for ) or shrinks (for ). In Figures 2(b) and 2(c) two LP characteristics derived from are plotted, for the indicated values of and , respectively.
Thus the 1D prototype is parametric depending on , and the derived 2D square-shaped filter will have an adjustable pass-band region, specified by the same .
3. Design Method
The main issue in this section is to find the transfer function of the desired 2D square filter using a particular frequency transformation.
From a 1D prototype filter (which varies on one axis only), a 2D oriented filter may be obtained by rotating the axes of the plane by an angle . The rotation is defined by the following linear transformation, where are the original frequency variables and the rotated ones [1]: The spatial orientation is specified by an angle with respect to , defined by the following 1D to 2D frequency mapping: The oriented filter transfer function results in the following: In the complex plane the frequency transformation (11) becomes The 2D oriented filter has the magnitude section along the line , identical with the prototype , and is constant along the perpendicular line: (the filter longitudinal axis).
In our case, since the rational prototype function from (6) contains only even powers of frequency , that is, only powers of , we will derive a convenient expression for , according to mapping (11): Thus using (11) we derived the mapping . Since , , the corresponding expression in the complex variables and is similar, with : Therefore, after applying the 1D to 2D directional mapping (11) by substituting , we reach the following rational expression in the complex frequency variables , which corresponds to each ratio factor of the prototype frequency response, given by (9) Since finally we must reach a transfer function of a discrete 2D filter in the complex variables and , we must determine the discrete counterpart of the general factor .
The currently used method to obtain a discrete filter from an analog prototype is applying the bilinear transform. If the sample interval takes the value , the bilinear transform for and in the complex plane has the form
Even if this method is straightforward, the designed 2D filter, corresponding to the transfer function in , , 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: where is a frequency of the discrete filter and is the corresponding frequency of the analog filter. This error can be corrected by applying a prewarping. Taking in relation (18) we obtain the mappings 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 : Thus the prewarping correction consists in the following substitution, also used in [15]: Even if (21) is useful in general, taking into account that the 2D filter function has even parity in and , that is, the variables appear as , a more efficient pre-warping can be made in this case.
We obtain the following Chebyshev-Padé approximation on the range : Consequently, from (19) and (22), the frequency pre-warping mapping can be written, on both frequency axes: 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 and : The analog-discrete mapping (25) therefore is a corrected form of the bilinear transform including pre-warping in the entire region of interest, namely, . 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 : This mapping is plotted in Figure 3 versus as curve F1 in solid line and is almost linear at least on the range , 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, , the mapping (26) becomes and since we obtain the mapping The next step is to substitute the mappings (25) and (28) into the expression (15) of , and we obtain a rational expression of order 4 in the complex frequency variables and . Thus the discrete approximation of may be expressed in the matrix form
where the variable vectors are , , , and stands for inner product.
The numerator matrix depends on the orientation angle according to the expression where , , and are the following matrices of size with constant elements
The denominator matrix is also a constant matrix: From (16) and (29), after some algebra, the matrices corresponding to the numerator and denominator of result from similar relations: where denotes matrix convolution.
The 2D directional IIR filter will have a real-valued (zero-phase) transfer function in the complex variables and , written in matrix form where the vectors are , and is the directional filter order. If , and are the matrices corresponding to the numerators and denominators of the two factors from (6), the matrices and of the directional filter result as convolution: Finally, the square-shaped filter transfer function is given by an expression similar to (34): where the matrices and result from where the superscript means rotation by 90 degrees; therefore is the clockwise rotated version of . This is due to the fact that the matrix pairs and 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 , and are found; if a second-order factor occurs, and are also found.(b)For a specified orientation angle , the matrix is determined using expression (30) and the partial constant matrices (31).(c)For a desired filter bandwidth , specified by the parameter , the partial matrices and for each of the 4-order factors result from relations (33).(d)The matrices and of the directional filter result as convolution of the matrices determined at the previous step, according to (35).(e)The matrices and 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 is specified by the parameter . As shown before, the matrices , , , and are constant once specified the 1D prototype , in our case given by (6).
The first design example is the diamond-type filter whose frequency response is shown in Figure 4(a). In order to obtain a diamond filter as in Figure 1(a), the orientation angle must be set . Ideally, the bandwidth (cut-off frequency) should be , corresponding to . 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 , corresponding to a bandwidth , 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 , . The square-shaped filter characteristic results as a product of the two directional characteristics.
(a)
(b)
(c)
(d)
In Figures 5 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)
(b)
(c)
(d)
(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 Figure 4(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.