Maximally Flat FIR and IIR Fractional Delay Filters with ... - eurasip
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
MAXIMALLY FLAT FIR AND IIR FRACTIONAL DELAY FILTERS WITH EXPANDED BANDWIDTH Martin Eichler and Arild Lacroix Institute of Applied Physics, J. W. Goethe University Max-von-Laue-Str. 1, D-60438, Frankfurt am Main, Germany phone: + (49) 69-798-47459, fax: + (49) 69-798-47444, email: [email protected] web: www.uni-frankfurt.de
ABSTRACT Fractional Delay (FD) Filters allow shifting of a given signal by a fraction of the sample interval. In the literature, a variety of approaches exist to design such filters, some of which are completely analytic. In this paper, we discuss an alternative analytic formulation for FIR-FD filters, which in one special case coincides with the widely known so-called maximally flat FIR FD filters. However, by variation of the free parameter, significant extension of the bandwith of (approximately) constant group delay is achieved. By deriving the same analytic form by a variable substitution, a way is found to apply this kind of optimization to IIR Allpass FD filters as well. Comparing the optimized FIR and IIR filters, it is shown that the FIR filters reach equivalent bandwidths and are thus attractive with regard to computational load. 1. PROBLEM FORMULATION We consider a linear digital filter given by the transfer function of an ideal FD filter, Hideal (z) = z−α
α ∈ R+ ,
we obtain the {ak } and {bk } (for k = 0, 1, 2, . . . , N):
=
∑
·
(3) (4)
The polynomial T (z) given by the {bk } in (2) now represents a non-causal filter. However, replacing the argument z by z−1 , we get a causal filter H(z) = T (z−1 ) approximating z−α : H(z) =
N
N
k=0
k=0
k −k ≈ z−α . ∑ ak (z−1 − z−1 0 ) = ∑ bk z
(5)
This is possible because for all z on the unit circle the approximation error is complex conjugate to the error when approximating z+α , so the accuracy is the same. Further, the error H(z) − Hideal (z) vanishes for z = z0 . By complete induction, it can be shown that for the {bk } according to (3) and (4) N
α −n n=0 k − n
bk = z0k−α · ∏
(6)
n6=k
due to z = eiω . In case of integer values of α, α = n ∈ N, this Hideal can be realized directly as z−n ; in case of non-integer values of α, a filter approximating (1) has to be developed. As a consequence, the group delay is no longer constant over the full frequency range. A multitude of methods exists for designing such approximating filters by analytic or numeric means ([1, 2]).
in (6) is 1 holds. In the special case z0 = 1, the factor zk−α 0 and the {bk } reduce to the coefficients of the maximally flat design FIR FD filter which is also known as the Lagrange interpolator. In the literature, this filter is derived by setting to zero the first N derivatives of E(z) = H(z)−Hideal (z) with respect to ω at the point ω = ω0 , choosing ω0 = 0 ([1], pp.40). Solving the corresponding equations yields (6) for z0 = 1 directly. If we modify this maximally flat approach for an arbitrary ω0 6= 0, we get N
α −n . n=0 k − n
bk = eiω0 (k−α) · ∏
2. FIR FD FILTER & TAYLOR SERIES EXPANSION In the following, we will consider a Taylor expansion of z+α , which opens an alternative view on the well-known maximally flat design. We consider the series expansion of zα in a polynomial of order N: N
N
k=0
k=0
!
k k α ∑ ak (z − z−1 0 ) = ∑ bk z ≈ z .
Using the binomial theorem and Taylor’s theorem, 1 d k α ak = , z · k! dzk z=z−1 0
©2007 EURASIP
bk
m=k
(1)
d d −∠Hideal (eiω ) = (αω) = α dω dω
T (z) =
= z0
N
having the group delay τ=
1 k−1 · ∏ (α − m), k! m=0 m am (−z0 )k−m . k
−(α−k)
ak
(2)
(7)
n6=k
It is easy to see that equations (6) and (7) are equivalent if z0 = eiω0 . From this, we can conclude that the standard maximally flat approach for ω0 = 0 is equivalent to a Taylor series expansion of z+α about the point z0 = 1, and, in the more general case ω0 6= 0, it is equivalent to a series expansion about −iω0 (in both cases, H(z) = z−α is approximated and z−1 0 =e the error vanishes at z = z0 ). However, together with (5) it now follows that the filter will be complex-valued for real ω0 6= 0. To keep the filter real-valued (for all α ∈ R), z0 must be real and positive, corresponding to purely imaginary ω0 . In section 5, we will see that by varying z0 on the positive real axis an expansion of the filter bandwidth is possible.
1038
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
3. VARIABLE SUBSTITUTION We will now show that the way eq. (6) was obtained is equivalent to applying a variable substitution to the filter resulting for z0 = 1. Let us consider an arbitrary function G1 (z) approximating z−α . We write G1 (z) = z−α − O1 (z) where O1 (z) denotes the approximation error, and assume that the approximation is perfect in z0 = 1 such that O1 (1) = 0. In the argument, we now substitute z by z/z0 and consider z −α Gz0 (z) := z0 · G1 . (8) z0 Using the definition of G1 (z), we can easily see that Gz0 (z) = z0−α ·
z z0
−α
− z−α 0 · O1
z z0
= z−α − Oz0 (z) ≈ z−α , with Oz0 (z) = z−α 0 · O1 (z/z0 ). Therefore, Gz0 (z) approximates z−α perfectly in z = z0 as Oz0 (z0 ) = 0. This means that in Gz0 (z), the point of perfect approximation is shifted from 1 to z0 . Let now G1 (z) be an arbitrary FIR filter polynomial approximating z−α perfectly at z = 1: N
G1 (z) =
∑ hk (1) · z−k ≈ z−α , k=0
where the index of G and the number in parentheses behind hk denote the point of perfect approximation. Inserting this into (8), we directly get the coefficients for Gz0 (z): Gz0 (z) =
z0−α
N
z · ∑ hk (1) · z0 k=0
−k
N
=
∑ hk (z0 ) · z
−k
denoting by ζ the point of perfect approximation. In [3], the coefficients for a purely recursive FD filter are derived under maximally flat conditions fulfilled in z = ζ = 1: N α −N +n N . ak (1) = (−1) k ∏ α −N +n+k k
Using these to build D1 (z) according to (11), a purely recursive filter H(z) = 1/D1 (z) is obtained having a group delay of (α − N)/2. Inserting this into (10), we get an allpass filter with group delay α ([1]). Now, to shift the point of perfect approximation from 1 to z0 , we need to have a closer look at D1 (z), because its phase characteristics approximates that of z(α−N)/2 while its magnitude does not – it is not even normalized ([3] gives a normalization factor which we do not regard here). We now introduce a normalizing function M1 (z) = |D1 (z)| which is fully real-valued all over C. With help of M1 (z), D1 (z) can be written as follows: h i α−N D1 (z) = σ · z 2 − O1 (z) · M1 (z), where σ = sgn(D1 (1)) and the complex-valued error O1 (z) can be obtained from D1 (z) and M1 (z) with O1 (1) = 0, because D1 (z) has real-valued coefficients (namely, (12)) and D1 (1) therefore is real. Thus, D1 (z)/M1 (z) clearly approximates σ z(α−N)/2 . In analogy to (8), we now consider α−N z Dz0 (z) = z0 2 D1 z0 α−N α−N z z 2 = · M1 σ · z 2 − z0 O1 z0 z0 and find that Dz0 (z) approximates the same phase characteristics as D1 (z), with its point of perfect approximation (in terms of phase) at z = ζ = z0 . Using (9) and replacing α by −(α − N)/2, we now directly get the coefficients of Dz0 (z):
,
k=0
with hk (z0 ) = z0k−α · hk (1).
k+ α−N ak (z0 ) = z0 2
(9)
Comparing (9) to (6), we see that multiplying by zk−α in or0 der to shift the point of perfect approximation is exactly what happens in case of the above Taylor series formulation. Most importantly, this modification will also work for any other finite polynomial approximating zβ , β ∈ C. 4. IIR FD ALLPASS FILTERS As for approximating (1) having a constant magnitude is desirable, approaches for designing IIR FD allpass filters also exist. Here, we consider the Thiran allpass filter discussed in [1]. The general form of an allpass filter can be written as z−N D(z−1 ) HAll pass (z) = , D(z)
(12)
n=0
(10)
N α −N +n N · (−1) . k ∏ α −N +n+k k
(13)
n=0
Building the filter (10) from Dz0 (z) with the coefficients (13) will again result in an allpass filter. However, as the filter is recursive, special attention has to be paid to stability. It has been shown analytically in [3] that all zeros of D1 (z) lie inside the unit circle for α > N. However, [4] states that this is the case even for α > N − 1. From this, it follows directly that the modified allpass filter according to (13) will be stable if |z0 | < z0,max =
1 , max |z pole,k |
(14)
where {z pole,k } are the zeros of D1 (z) with its coefficients calculated according to (12) for a given α > N − 1. 5. BANDWIDTH EXPANSION
where D(z) is a polynomial of Nth order, which we write as N
D(z) = Dζ (z) =
∑ ak (ζ ) · z−k , k=0
©2007 EURASIP
(11)
We will now show how the coefficient modifications derived above can be used to optimize the filter bandwidths. We define the effective bandwidth of an FD-Filter of given order N at a given α as the maximum frequency range for
1039
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
Optimum z0 for Maximum Bandwidth 2.4
3
2 1.8 1.6 1.4 1.2
2
0
2.5
α=3 α=2.9 α=2.8 α=2.7 α=2.6 α=2.5 α=2.4 α=2.3 α=2.2 α=2.1 α=2
z (logarithmic scale)
Group Delay [samples]
5th Order: Group Delay Response τ 3.5
standard bandwidth→ optimized bandwidth→
1.5 0
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
−0.5
10th Order: Group Delay Response τ
Group Delay [samples]
5
4.5
standard bandwidth→ optimized bandwidth→
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
Figure 1: Group Delay of 5th (top) and 10th order (bottom) FIR filters. Dashed lines: z0 = 1; solid lines: optimized z0 for ∆τ = 0.01. Circles indicate the effective bandwidths. which the group delay τ does not exceed an error tolerance range [α − ∆τ, α + ∆τ] about the specified group delay α. Now, we define the overall bandwidth of an FD-Filter of given order N as the minimum of all effective bandwidths over the specified optimum range of α. For the FIR filters, this is (N − 1)/2 ≤ α ≤ (N + 1)/2; for the IIR allpass filters, the optimum range is N − 1 < α < N ([4]). Thus, the overall bandwidth represents the frequency range over which the filter can be operated for any α inside the optimum range without exceeding the error tolerance specified by ∆τ. Once the value of ∆τ is defined (based on the requirements of the application), we can optimize z0 such that the effective bandwidth is maximized for each α. From the set of effective bandwidths, the overall bandwidth is obtained. In general, the optimum z0 will be a function of α, depending on the filter order and error tolerance: z0,opt = f (α; N, ∆τ). 5.1 Optimized FIR Filters Two examples for orders 5 and 10 are shown in figure 1; the corresponding curves of z0 are shown in figure 2. Note that already the optimized 5th order filter works almost perfectly in the frequency range [0|π/2]. For integer values of α,
©2007 EURASIP
0 α−N/2
0.5
Figure 2: Optimized z0 for 5th (squares) and 10th order (circles) FIR filters as a function of specified group delay α and error tolerance ∆τ.
α=5.5 α=5.4 α=5.3 α=5.2 α=5.1 α=5 α=4.9 α=4.8 α=4.7 α=4.6 α=4.5
4 0
1 0.9 0.8 0.7 0.6 0.5
6
5.5
N=10, ∆τ=0.05 N=10, ∆τ=0.01
N=5, ∆τ=0.05 N=5, ∆τ=0.01
the group delay is exactly equal to α and the magnitude is H(z) ≡ 1, independent of z0 . z0 is hence arbitrary and can be set to 1. One can observe now that the shape of the z0 curves is different for filters of even and odd orders: Even orders yield discontinuous curves, odd orders yield continuous curves. Further, we observe two symmetries. As can be seen easily from (6), for the coefficients {bk } the relation bk (z0 , α) = bN−k (z−1 0 , N − α) holds for k = 0, 1, 2, . . . , N. As a first consequence of this, both standard and optimized group delay curves are symmetric with respect to the curve for α = N/2 (see figure 1). Second, for the optimized z0 we find that z0,opt (α) = 1/z0,opt (N − α) holds, which causes a point symmetry of log z0,opt with respect to the point (α = N/2, z0 = 1)) (see figure 2). This property can be useful in a practical implementation, as only one branch of the optimized curve needs to be stored or modeled. Figure 7 shows the resulting overall bandwidths for filter orders N = 1 . . . 21. Looking at the magnitude of the FIR filters, we see that they have a low-pass characteristic, with an almost constant value of 1 and then decreasing (see example in figure 3). This general behaviour is kept under optimization of z0 , however, the decrease gets slightly bigger. We determine the maximum magnitude deviation inside the overall bandwidth for standard and z0 -optimized filters by vertically cutting through the set of graphs in figure 7 at the bandwidth frequency and taking the minimum absolute value (see figure 3). Then, for filter orders N = 1 . . . 21, we get the graph shown in figure 4. For ∆τ = 0.01, the deviation is less than ≈ 0.5dB, and for ∆τ = 0.05 it is less than ≈ 2.2dB. 5.2 Optimized IIR Allpass Filters In contrast to the FIR filters, for the IIR allpass, group delay curves and optimum z0 curves have similar shapes for even and odd orders, as can be seen in figure 5. Another big difference is that the z0 curves (figure 6) do not show any symmetry. 5.3 Comparison FIR vs. IIR Allpass Filters Summarizing the overall bandwidths of FIR and IIR filters in figure 7, one can see that the increase in bandwidth gai-
1040
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
5th Order: Magnitude Response
5th Order: Group Delay Response τ
1.05
5.2 α = 2.0, 3.0
Group Delay [samples]
5
Magnitude |H|
1 magnitude deviation inside bandwidth
α = 2.1, 2.9
0.95
α = 2.2, 2.8 standard bandwidth→ optimized bandwidth→
α = 2.3, 2.7 α = 2.4, 2.6
0.9 α = 2.5
4.8 4.6 4.4 4.2 4
0.85 0
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
10 Group Delay [samples]
Magnitude Deviation / dB
α=4.7 α=4.6 α=4.5 α=4.4 α=4.3 α=4.2 α=4.1 α=4.001
0.2
0
∆τ=0.01 (optimized)
∆τ=0.05 (optimized)
9.8 9.6 9.4 9.2 9
−2.5 0
5
10 Order N
15
α=9.9 α=9.8 α=9.7 α=9.6 α=9.5 α=9.4 α=9.3 α=9.2 α=9.1 α=9.001
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
Figure 5: Group Delay of 5th (top) and 10th order (bottom) IIR allpass filters. Dashed lines: z0 = 1; solid lines: optimized z0 for ∆τ = 0.01. Circles indicate the effective bandwidths.
20
Optimum z for Maximum Bandwidth & Maximum Pole Radius
Figure 4: Magnitude deviation in dB of FIR filters of orders 1 through 21, with respect to the ideal magnitude |H(z)| ≡ 1.
0
z
0
1.8
ned by variation of z0 is far bigger in case of the FIR filters. For odd orders, the bandwidths of the optimized FIR filters are almost as high as those of the optimized allpass filters, in some cases even higher. At the same time, the deviation in magnitude exhibited by the FIR filters can be kept quite small for small ∆τ (less than 0.5dB for ∆τ=0.01). If such a magnitude deviation is acceptable under the conditions of the application, we can say that the FIR filters are equivalent to the IIR allpass filters, and may even be preferred due to their bandwidth (especially in case of odd orders) and due to the fact that they only require roughly half of the computational load IIR allpass filters would. 6. PHASE DELAY BANDWIDTHS For additional reference, figure 8 shows the phase delay bandwidths corresponding to the z0 values shown in figure 2 and 6 which were optimized for maximum group delay bandwidth. When optimizing z0 for phase delay, even higher
©2007 EURASIP
1
standard bandwidth→ optimized bandwidth→
α=10.001
0
−1.5 −2
0.4 0.6 0.8 Normalized Frequency / π
10th Order: Group Delay Response τ
∆τ=0.01 (standard) ∆τ=0.05 (standard)
−1
α=4.8
10.2
Bandwidth Evaluation: Magnitude Deviation
−0.5
α=4.9
0
Figure 3: Magnitude of the 5th order FIR filter, corresponding to fig.1 (top). Circles show the frequencies where the respective group delay leaves the error tolerance range ∆τ = 0.01. Dashed lines: z0 = 1; solid lines: optimized z0 . 0.5
standard bandwidth→ optimized bandwidth→
α=5.001
1.6
z0, max for N=5 z0, max for N=10
1.4 Optimized z0 1.2
max{|z
pole
|}
1
N=5, ∆τ=0.05 N=10, ∆τ=0.05 N=5, ∆τ=0.01 N=10, ∆τ=0.01
N=10, ∆τ=0.05 N=10, ∆τ=0.01 N=10, standard filter
0.8 0.6 Maximum Pole Radius
0.4 0.2 0
N=5, ∆τ=0.05 N=5, ∆τ=0.01 N=5, standard filter
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 α − (N−1)
1
Figure 6: Optimized z0 and maximum pole radius for 5th (squares) and 10th order (circles) IIR allpass filters as a function of specified group delay α and error tolerance ∆τ. Dotted curves: maximum pole radius, max{|z pole,k |}, of standard Thiran allpass filter and z0,max , respectively (see eq.(14)).
1041
Bandwidth Evaluation: Normalized Group Delay Bandwidth 1
Bandwidth Evaluation: Normalized Phase Delay Bandwidth 1
0.8
0.8 Bandwidth / π
Bandwidth / π
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
5
10 Order N
15
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
20
5
10 Order N
15
20
Bandwidth Evaluation: Normalized Group Delay Bandwidth 1
Bandwidth Evaluation: Normalized Phase Delay Bandwidth 1
0.8
0.8 Bandwidth / π
Bandwidth / π
0 0
0.6
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
5
10 Order N
15
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
20
5
10 Order N
15
20
Figure 7: Group delay bandwidth comparison of FIR and IIR allpass filters of orders 1 through 21 for ∆τ = 0.01 (top) and ∆τ = 0.05 (bottom).
Figure 8: Phase delay bandwidth of FIR and IIR allpass filters of orders 1 through 21 for ∆τ = 0.01 (top) and ∆τ = 0.05 (bottom), resulting when z0 is optimized for maximum group delay bandwidth.
bandwidths can be expected; however, the deviation in group delay will then increase. Also, as the low-pass characteristic of the FIR filters will persist, the deviation in magnitude next to the border of bandwidth will increase as well.
maximum bandwidths resulting for the FIR filters are equivalent to those resulting for the IIR allpass filters, which suggests that, given a bandwidth required by the application, the FIR filters are advantageous over the IIR allpass filters with regard to computational load.
7. SUMMARY It is shown that the standard maximally flat design approach of FIR FD filters described by the (ideal) transfer function H(z) = z−α is equivalent to a series expansion of z+α about the point 1 and that, rewriting this design for general ω0 , it is equivalent to a Taylor series expansion of z+α about the point −iω0 . Further, it is shown that the resulting filter coz−1 0 =e efficients can also be obtained by a variable substitution in z, mapping z to z/z0 . Applying this principle to the maximally flat IIR allpass filters, a generalized form of their coefficients is obtained which allows the variation of z0 the same way the FIR filter coefficients do. It is shown that, given an error tolerance ∆τ > 0, the bandwidth of both FIR and IIR filters can be increased significantly by varying z0 on the positive real axis, with the filter remaining real-valued. However, the
©2007 EURASIP
REFERENCES [1] T. I. Laakso, V. V¨alim¨aki, M. Karjalainen and U. K. Laine, “Splitting the Unit Delay”, IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30–60, Jan. 1996. [2] J. T. Olkkonen, H. Olkkonen, “Fractional Delay Filter Based on the B-Spline Transform”, IEEE Signal Processing Letters, vol. 14, no. 2, pp. 97–100, Feb. 2007. [3] J.-P. Thiran, “Recursive Digital Filters with Maximally Flat Group Delay”, IEEE Transactions on Circuit Theory, vol. CT-18, no. 6, pp. 659–664, Nov. 1971. [4] V. V¨alim¨aki, PhD Thesis “Discrete-Time Modeling of Acoustic Tubes Using Fractional Delay Filters”, Helsinki University of Technology, Dec. 1995.
1042
MAXIMALLY FLAT FIR AND IIR FRACTIONAL DELAY FILTERS WITH EXPANDED BANDWIDTH Martin Eichler and Arild Lacroix Institute of Applied Physics, J. W. Goethe University Max-von-Laue-Str. 1, D-60438, Frankfurt am Main, Germany phone: + (49) 69-798-47459, fax: + (49) 69-798-47444, email: [email protected] web: www.uni-frankfurt.de
ABSTRACT Fractional Delay (FD) Filters allow shifting of a given signal by a fraction of the sample interval. In the literature, a variety of approaches exist to design such filters, some of which are completely analytic. In this paper, we discuss an alternative analytic formulation for FIR-FD filters, which in one special case coincides with the widely known so-called maximally flat FIR FD filters. However, by variation of the free parameter, significant extension of the bandwith of (approximately) constant group delay is achieved. By deriving the same analytic form by a variable substitution, a way is found to apply this kind of optimization to IIR Allpass FD filters as well. Comparing the optimized FIR and IIR filters, it is shown that the FIR filters reach equivalent bandwidths and are thus attractive with regard to computational load. 1. PROBLEM FORMULATION We consider a linear digital filter given by the transfer function of an ideal FD filter, Hideal (z) = z−α
α ∈ R+ ,
we obtain the {ak } and {bk } (for k = 0, 1, 2, . . . , N):
=
∑
·
(3) (4)
The polynomial T (z) given by the {bk } in (2) now represents a non-causal filter. However, replacing the argument z by z−1 , we get a causal filter H(z) = T (z−1 ) approximating z−α : H(z) =
N
N
k=0
k=0
k −k ≈ z−α . ∑ ak (z−1 − z−1 0 ) = ∑ bk z
(5)
This is possible because for all z on the unit circle the approximation error is complex conjugate to the error when approximating z+α , so the accuracy is the same. Further, the error H(z) − Hideal (z) vanishes for z = z0 . By complete induction, it can be shown that for the {bk } according to (3) and (4) N
α −n n=0 k − n
bk = z0k−α · ∏
(6)
n6=k
due to z = eiω . In case of integer values of α, α = n ∈ N, this Hideal can be realized directly as z−n ; in case of non-integer values of α, a filter approximating (1) has to be developed. As a consequence, the group delay is no longer constant over the full frequency range. A multitude of methods exists for designing such approximating filters by analytic or numeric means ([1, 2]).
in (6) is 1 holds. In the special case z0 = 1, the factor zk−α 0 and the {bk } reduce to the coefficients of the maximally flat design FIR FD filter which is also known as the Lagrange interpolator. In the literature, this filter is derived by setting to zero the first N derivatives of E(z) = H(z)−Hideal (z) with respect to ω at the point ω = ω0 , choosing ω0 = 0 ([1], pp.40). Solving the corresponding equations yields (6) for z0 = 1 directly. If we modify this maximally flat approach for an arbitrary ω0 6= 0, we get N
α −n . n=0 k − n
bk = eiω0 (k−α) · ∏
2. FIR FD FILTER & TAYLOR SERIES EXPANSION In the following, we will consider a Taylor expansion of z+α , which opens an alternative view on the well-known maximally flat design. We consider the series expansion of zα in a polynomial of order N: N
N
k=0
k=0
!
k k α ∑ ak (z − z−1 0 ) = ∑ bk z ≈ z .
Using the binomial theorem and Taylor’s theorem, 1 d k α ak = , z · k! dzk z=z−1 0
©2007 EURASIP
bk
m=k
(1)
d d −∠Hideal (eiω ) = (αω) = α dω dω
T (z) =
= z0
N
having the group delay τ=
1 k−1 · ∏ (α − m), k! m=0 m am (−z0 )k−m . k
−(α−k)
ak
(2)
(7)
n6=k
It is easy to see that equations (6) and (7) are equivalent if z0 = eiω0 . From this, we can conclude that the standard maximally flat approach for ω0 = 0 is equivalent to a Taylor series expansion of z+α about the point z0 = 1, and, in the more general case ω0 6= 0, it is equivalent to a series expansion about −iω0 (in both cases, H(z) = z−α is approximated and z−1 0 =e the error vanishes at z = z0 ). However, together with (5) it now follows that the filter will be complex-valued for real ω0 6= 0. To keep the filter real-valued (for all α ∈ R), z0 must be real and positive, corresponding to purely imaginary ω0 . In section 5, we will see that by varying z0 on the positive real axis an expansion of the filter bandwidth is possible.
1038
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
3. VARIABLE SUBSTITUTION We will now show that the way eq. (6) was obtained is equivalent to applying a variable substitution to the filter resulting for z0 = 1. Let us consider an arbitrary function G1 (z) approximating z−α . We write G1 (z) = z−α − O1 (z) where O1 (z) denotes the approximation error, and assume that the approximation is perfect in z0 = 1 such that O1 (1) = 0. In the argument, we now substitute z by z/z0 and consider z −α Gz0 (z) := z0 · G1 . (8) z0 Using the definition of G1 (z), we can easily see that Gz0 (z) = z0−α ·
z z0
−α
− z−α 0 · O1
z z0
= z−α − Oz0 (z) ≈ z−α , with Oz0 (z) = z−α 0 · O1 (z/z0 ). Therefore, Gz0 (z) approximates z−α perfectly in z = z0 as Oz0 (z0 ) = 0. This means that in Gz0 (z), the point of perfect approximation is shifted from 1 to z0 . Let now G1 (z) be an arbitrary FIR filter polynomial approximating z−α perfectly at z = 1: N
G1 (z) =
∑ hk (1) · z−k ≈ z−α , k=0
where the index of G and the number in parentheses behind hk denote the point of perfect approximation. Inserting this into (8), we directly get the coefficients for Gz0 (z): Gz0 (z) =
z0−α
N
z · ∑ hk (1) · z0 k=0
−k
N
=
∑ hk (z0 ) · z
−k
denoting by ζ the point of perfect approximation. In [3], the coefficients for a purely recursive FD filter are derived under maximally flat conditions fulfilled in z = ζ = 1: N α −N +n N . ak (1) = (−1) k ∏ α −N +n+k k
Using these to build D1 (z) according to (11), a purely recursive filter H(z) = 1/D1 (z) is obtained having a group delay of (α − N)/2. Inserting this into (10), we get an allpass filter with group delay α ([1]). Now, to shift the point of perfect approximation from 1 to z0 , we need to have a closer look at D1 (z), because its phase characteristics approximates that of z(α−N)/2 while its magnitude does not – it is not even normalized ([3] gives a normalization factor which we do not regard here). We now introduce a normalizing function M1 (z) = |D1 (z)| which is fully real-valued all over C. With help of M1 (z), D1 (z) can be written as follows: h i α−N D1 (z) = σ · z 2 − O1 (z) · M1 (z), where σ = sgn(D1 (1)) and the complex-valued error O1 (z) can be obtained from D1 (z) and M1 (z) with O1 (1) = 0, because D1 (z) has real-valued coefficients (namely, (12)) and D1 (1) therefore is real. Thus, D1 (z)/M1 (z) clearly approximates σ z(α−N)/2 . In analogy to (8), we now consider α−N z Dz0 (z) = z0 2 D1 z0 α−N α−N z z 2 = · M1 σ · z 2 − z0 O1 z0 z0 and find that Dz0 (z) approximates the same phase characteristics as D1 (z), with its point of perfect approximation (in terms of phase) at z = ζ = z0 . Using (9) and replacing α by −(α − N)/2, we now directly get the coefficients of Dz0 (z):
,
k=0
with hk (z0 ) = z0k−α · hk (1).
k+ α−N ak (z0 ) = z0 2
(9)
Comparing (9) to (6), we see that multiplying by zk−α in or0 der to shift the point of perfect approximation is exactly what happens in case of the above Taylor series formulation. Most importantly, this modification will also work for any other finite polynomial approximating zβ , β ∈ C. 4. IIR FD ALLPASS FILTERS As for approximating (1) having a constant magnitude is desirable, approaches for designing IIR FD allpass filters also exist. Here, we consider the Thiran allpass filter discussed in [1]. The general form of an allpass filter can be written as z−N D(z−1 ) HAll pass (z) = , D(z)
(12)
n=0
(10)
N α −N +n N · (−1) . k ∏ α −N +n+k k
(13)
n=0
Building the filter (10) from Dz0 (z) with the coefficients (13) will again result in an allpass filter. However, as the filter is recursive, special attention has to be paid to stability. It has been shown analytically in [3] that all zeros of D1 (z) lie inside the unit circle for α > N. However, [4] states that this is the case even for α > N − 1. From this, it follows directly that the modified allpass filter according to (13) will be stable if |z0 | < z0,max =
1 , max |z pole,k |
(14)
where {z pole,k } are the zeros of D1 (z) with its coefficients calculated according to (12) for a given α > N − 1. 5. BANDWIDTH EXPANSION
where D(z) is a polynomial of Nth order, which we write as N
D(z) = Dζ (z) =
∑ ak (ζ ) · z−k , k=0
©2007 EURASIP
(11)
We will now show how the coefficient modifications derived above can be used to optimize the filter bandwidths. We define the effective bandwidth of an FD-Filter of given order N at a given α as the maximum frequency range for
1039
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
Optimum z0 for Maximum Bandwidth 2.4
3
2 1.8 1.6 1.4 1.2
2
0
2.5
α=3 α=2.9 α=2.8 α=2.7 α=2.6 α=2.5 α=2.4 α=2.3 α=2.2 α=2.1 α=2
z (logarithmic scale)
Group Delay [samples]
5th Order: Group Delay Response τ 3.5
standard bandwidth→ optimized bandwidth→
1.5 0
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
−0.5
10th Order: Group Delay Response τ
Group Delay [samples]
5
4.5
standard bandwidth→ optimized bandwidth→
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
Figure 1: Group Delay of 5th (top) and 10th order (bottom) FIR filters. Dashed lines: z0 = 1; solid lines: optimized z0 for ∆τ = 0.01. Circles indicate the effective bandwidths. which the group delay τ does not exceed an error tolerance range [α − ∆τ, α + ∆τ] about the specified group delay α. Now, we define the overall bandwidth of an FD-Filter of given order N as the minimum of all effective bandwidths over the specified optimum range of α. For the FIR filters, this is (N − 1)/2 ≤ α ≤ (N + 1)/2; for the IIR allpass filters, the optimum range is N − 1 < α < N ([4]). Thus, the overall bandwidth represents the frequency range over which the filter can be operated for any α inside the optimum range without exceeding the error tolerance specified by ∆τ. Once the value of ∆τ is defined (based on the requirements of the application), we can optimize z0 such that the effective bandwidth is maximized for each α. From the set of effective bandwidths, the overall bandwidth is obtained. In general, the optimum z0 will be a function of α, depending on the filter order and error tolerance: z0,opt = f (α; N, ∆τ). 5.1 Optimized FIR Filters Two examples for orders 5 and 10 are shown in figure 1; the corresponding curves of z0 are shown in figure 2. Note that already the optimized 5th order filter works almost perfectly in the frequency range [0|π/2]. For integer values of α,
©2007 EURASIP
0 α−N/2
0.5
Figure 2: Optimized z0 for 5th (squares) and 10th order (circles) FIR filters as a function of specified group delay α and error tolerance ∆τ.
α=5.5 α=5.4 α=5.3 α=5.2 α=5.1 α=5 α=4.9 α=4.8 α=4.7 α=4.6 α=4.5
4 0
1 0.9 0.8 0.7 0.6 0.5
6
5.5
N=10, ∆τ=0.05 N=10, ∆τ=0.01
N=5, ∆τ=0.05 N=5, ∆τ=0.01
the group delay is exactly equal to α and the magnitude is H(z) ≡ 1, independent of z0 . z0 is hence arbitrary and can be set to 1. One can observe now that the shape of the z0 curves is different for filters of even and odd orders: Even orders yield discontinuous curves, odd orders yield continuous curves. Further, we observe two symmetries. As can be seen easily from (6), for the coefficients {bk } the relation bk (z0 , α) = bN−k (z−1 0 , N − α) holds for k = 0, 1, 2, . . . , N. As a first consequence of this, both standard and optimized group delay curves are symmetric with respect to the curve for α = N/2 (see figure 1). Second, for the optimized z0 we find that z0,opt (α) = 1/z0,opt (N − α) holds, which causes a point symmetry of log z0,opt with respect to the point (α = N/2, z0 = 1)) (see figure 2). This property can be useful in a practical implementation, as only one branch of the optimized curve needs to be stored or modeled. Figure 7 shows the resulting overall bandwidths for filter orders N = 1 . . . 21. Looking at the magnitude of the FIR filters, we see that they have a low-pass characteristic, with an almost constant value of 1 and then decreasing (see example in figure 3). This general behaviour is kept under optimization of z0 , however, the decrease gets slightly bigger. We determine the maximum magnitude deviation inside the overall bandwidth for standard and z0 -optimized filters by vertically cutting through the set of graphs in figure 7 at the bandwidth frequency and taking the minimum absolute value (see figure 3). Then, for filter orders N = 1 . . . 21, we get the graph shown in figure 4. For ∆τ = 0.01, the deviation is less than ≈ 0.5dB, and for ∆τ = 0.05 it is less than ≈ 2.2dB. 5.2 Optimized IIR Allpass Filters In contrast to the FIR filters, for the IIR allpass, group delay curves and optimum z0 curves have similar shapes for even and odd orders, as can be seen in figure 5. Another big difference is that the z0 curves (figure 6) do not show any symmetry. 5.3 Comparison FIR vs. IIR Allpass Filters Summarizing the overall bandwidths of FIR and IIR filters in figure 7, one can see that the increase in bandwidth gai-
1040
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
5th Order: Magnitude Response
5th Order: Group Delay Response τ
1.05
5.2 α = 2.0, 3.0
Group Delay [samples]
5
Magnitude |H|
1 magnitude deviation inside bandwidth
α = 2.1, 2.9
0.95
α = 2.2, 2.8 standard bandwidth→ optimized bandwidth→
α = 2.3, 2.7 α = 2.4, 2.6
0.9 α = 2.5
4.8 4.6 4.4 4.2 4
0.85 0
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
10 Group Delay [samples]
Magnitude Deviation / dB
α=4.7 α=4.6 α=4.5 α=4.4 α=4.3 α=4.2 α=4.1 α=4.001
0.2
0
∆τ=0.01 (optimized)
∆τ=0.05 (optimized)
9.8 9.6 9.4 9.2 9
−2.5 0
5
10 Order N
15
α=9.9 α=9.8 α=9.7 α=9.6 α=9.5 α=9.4 α=9.3 α=9.2 α=9.1 α=9.001
0.2
0.4 0.6 0.8 Normalized Frequency / π
1
Figure 5: Group Delay of 5th (top) and 10th order (bottom) IIR allpass filters. Dashed lines: z0 = 1; solid lines: optimized z0 for ∆τ = 0.01. Circles indicate the effective bandwidths.
20
Optimum z for Maximum Bandwidth & Maximum Pole Radius
Figure 4: Magnitude deviation in dB of FIR filters of orders 1 through 21, with respect to the ideal magnitude |H(z)| ≡ 1.
0
z
0
1.8
ned by variation of z0 is far bigger in case of the FIR filters. For odd orders, the bandwidths of the optimized FIR filters are almost as high as those of the optimized allpass filters, in some cases even higher. At the same time, the deviation in magnitude exhibited by the FIR filters can be kept quite small for small ∆τ (less than 0.5dB for ∆τ=0.01). If such a magnitude deviation is acceptable under the conditions of the application, we can say that the FIR filters are equivalent to the IIR allpass filters, and may even be preferred due to their bandwidth (especially in case of odd orders) and due to the fact that they only require roughly half of the computational load IIR allpass filters would. 6. PHASE DELAY BANDWIDTHS For additional reference, figure 8 shows the phase delay bandwidths corresponding to the z0 values shown in figure 2 and 6 which were optimized for maximum group delay bandwidth. When optimizing z0 for phase delay, even higher
©2007 EURASIP
1
standard bandwidth→ optimized bandwidth→
α=10.001
0
−1.5 −2
0.4 0.6 0.8 Normalized Frequency / π
10th Order: Group Delay Response τ
∆τ=0.01 (standard) ∆τ=0.05 (standard)
−1
α=4.8
10.2
Bandwidth Evaluation: Magnitude Deviation
−0.5
α=4.9
0
Figure 3: Magnitude of the 5th order FIR filter, corresponding to fig.1 (top). Circles show the frequencies where the respective group delay leaves the error tolerance range ∆τ = 0.01. Dashed lines: z0 = 1; solid lines: optimized z0 . 0.5
standard bandwidth→ optimized bandwidth→
α=5.001
1.6
z0, max for N=5 z0, max for N=10
1.4 Optimized z0 1.2
max{|z
pole
|}
1
N=5, ∆τ=0.05 N=10, ∆τ=0.05 N=5, ∆τ=0.01 N=10, ∆τ=0.01
N=10, ∆τ=0.05 N=10, ∆τ=0.01 N=10, standard filter
0.8 0.6 Maximum Pole Radius
0.4 0.2 0
N=5, ∆τ=0.05 N=5, ∆τ=0.01 N=5, standard filter
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 α − (N−1)
1
Figure 6: Optimized z0 and maximum pole radius for 5th (squares) and 10th order (circles) IIR allpass filters as a function of specified group delay α and error tolerance ∆τ. Dotted curves: maximum pole radius, max{|z pole,k |}, of standard Thiran allpass filter and z0,max , respectively (see eq.(14)).
1041
Bandwidth Evaluation: Normalized Group Delay Bandwidth 1
Bandwidth Evaluation: Normalized Phase Delay Bandwidth 1
0.8
0.8 Bandwidth / π
Bandwidth / π
15th European Signal Processing Conference (EUSIPCO 2007), Poznan, Poland, September 3-7, 2007, copyright by EURASIP
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
5
10 Order N
15
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
20
5
10 Order N
15
20
Bandwidth Evaluation: Normalized Group Delay Bandwidth 1
Bandwidth Evaluation: Normalized Phase Delay Bandwidth 1
0.8
0.8 Bandwidth / π
Bandwidth / π
0 0
0.6
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
5
10 Order N
15
0.6
0.4 IIR Allpass, standard FIR, standard IIR Allpass, optimized FIR, optimized
0.2
0 0
20
5
10 Order N
15
20
Figure 7: Group delay bandwidth comparison of FIR and IIR allpass filters of orders 1 through 21 for ∆τ = 0.01 (top) and ∆τ = 0.05 (bottom).
Figure 8: Phase delay bandwidth of FIR and IIR allpass filters of orders 1 through 21 for ∆τ = 0.01 (top) and ∆τ = 0.05 (bottom), resulting when z0 is optimized for maximum group delay bandwidth.
bandwidths can be expected; however, the deviation in group delay will then increase. Also, as the low-pass characteristic of the FIR filters will persist, the deviation in magnitude next to the border of bandwidth will increase as well.
maximum bandwidths resulting for the FIR filters are equivalent to those resulting for the IIR allpass filters, which suggests that, given a bandwidth required by the application, the FIR filters are advantageous over the IIR allpass filters with regard to computational load.
7. SUMMARY It is shown that the standard maximally flat design approach of FIR FD filters described by the (ideal) transfer function H(z) = z−α is equivalent to a series expansion of z+α about the point 1 and that, rewriting this design for general ω0 , it is equivalent to a Taylor series expansion of z+α about the point −iω0 . Further, it is shown that the resulting filter coz−1 0 =e efficients can also be obtained by a variable substitution in z, mapping z to z/z0 . Applying this principle to the maximally flat IIR allpass filters, a generalized form of their coefficients is obtained which allows the variation of z0 the same way the FIR filter coefficients do. It is shown that, given an error tolerance ∆τ > 0, the bandwidth of both FIR and IIR filters can be increased significantly by varying z0 on the positive real axis, with the filter remaining real-valued. However, the
©2007 EURASIP
REFERENCES [1] T. I. Laakso, V. V¨alim¨aki, M. Karjalainen and U. K. Laine, “Splitting the Unit Delay”, IEEE Signal Processing Magazine, vol. 13, no. 1, pp. 30–60, Jan. 1996. [2] J. T. Olkkonen, H. Olkkonen, “Fractional Delay Filter Based on the B-Spline Transform”, IEEE Signal Processing Letters, vol. 14, no. 2, pp. 97–100, Feb. 2007. [3] J.-P. Thiran, “Recursive Digital Filters with Maximally Flat Group Delay”, IEEE Transactions on Circuit Theory, vol. CT-18, no. 6, pp. 659–664, Nov. 1971. [4] V. V¨alim¨aki, PhD Thesis “Discrete-Time Modeling of Acoustic Tubes Using Fractional Delay Filters”, Helsinki University of Technology, Dec. 1995.
1042
