Fractional Delay Digital Filters

1. Introduction

The chapter goal is focused to introduce the concept of fractional delay filters (FDF), as well as a concise description of most of the existing design techniques. For this purpose, several illustrative examples are presented, where each design method is implemented by MATLAB programs.

A fractional delay filter is a filter of digital type having as main function to delay the processed input signal a fractional of the sampling period time. There are several applications where such signal delay value is required, examples of such systems are: timing adjustment in all-digital receivers (symbol synchronization), conversion between arbitrary sampling frequencies, echo cancellation, speech coding and synthesis, musical instruments modelling etc. (Laakson et al., 1996).

In order to achieve the fractional delay filter function, two main frequency-domain specifications must be met by the filter. The filter magnitude frequency response must have an all-pass behaviour in a wide frequency range, as well as its phase frequency response must be linear with a fixed fractional slope through the bandwidth.

Several FIR design methods have been reported during the last two decades. There are two main design approaches: time-domain and frequency-domain design methods. In first one, the fractional delay filter coefficients are easily obtained through classical mathematical interpolation formulas, but there is a small flexibility to meet frequency-domain specifications. On the other hand, the frequency-domain methods are based on frequency optimization process, and a more frequency specification control is available. One important result of frequency-domain methods is a highly efficient implementation structure called Farrow structure, which allows online fractional value update.

The chapter is organized as follows. Next section gives the formal definition of fractional delay filter. In the third section, some design methods are briefly described. Two efficient implementation structures for wideband fractional delay filter, as well as description of recently reported design methods for such structures, are illustrated in fourth section. MATLAB designed examples and concluding remarks are presented in fifth and sixth sections, respectively.

2. Fractional delay filter definition

The continuous-time output signal y_a(t) of a general signal delay system is defined by:

$$y_a\left(t\right)=x\left(t-t_l\right)$$E1

where x(t) is the continuous-time input signal and t_l the obtained time delay value. In a discrete-time system, the input-output relationship of a signal delay system is expressed as:

$$y\left(lT\right)=x\left(nT-DT\right)$$E2

where the delay value is given by DT, y(lT) and x(nT) are the discrete-time versions of output and input signals, respectively, and T is the sampling period time.

A signal delay value equal to a multiple of the sampling period, D as an integer N, can be easily implemented in a discrete-time system by memory elements storing the signal value for a time of NT:

$$y\left(lT\right)=x\left(nT-NT\right)$$E3

In this case, the signal delay value is limited to be only N time the sampling period, t_l=NT. For instance in telephone quality signals, with a sampling frequency of 8 KHz, only delays values multiple of 125μseconds are allowed.

Let us introduce the FDF function using time-domain signals sketched in Fig 1. The FDF output y(lT), squared samples, is obtained a delay time t_l after input x(n_l), with a delay value μ_lT given as a fraction of the sampling period time, 0<μ_l<1. As shown in Fig. 1, the fractional delay value μ_l may be variable; this way, it can be changed at any desired time.

The fundamental design problem of a FDF is to obtain the FDF unit impulse response h_FD(n,μ), in such a way that the obtained output value y(lT) = y_a(DT) be as close as possible to y_a(t_l) for 0<μ_l <1. The simplified block diagram for a FDF is shown in Fig. 2, which output for a no causal FIR FDF filter is given by the discrete-time convolution:

$$y\left(lT\right)={\displaystyle \sum _{k=-N_{FD}/2}^{N_{FD}/2-1}x\left(n_l-k\right)h_{FD}\left(k,{\mu }_l\right)}$$E4

where N_FD is the even length of the FDF. The system function H(z) of the FDF can be expressed as:

$$H\left(z\right)=z^{-D}$$E5

Figure 1.

FDF time-domain behaviour.

Figure 2.

Simplified block diagram for a FDF.

where the delay value is given as: D = D_fix+μ_l, D_fix is a fixed delay value and μ_l is the desired fractional delay value. As a consequence, the ideal frequency response of a FDF H_id(ω,μ_l) is:

$$H_{id}\left(\omega ,{\mu }_l\right)=e^{-j\left(D_{fix}+{\mu }_l\right)\omega }$$E6

Hence the ideal FDF frequency response has an all-band unity magnitude response:

$$\left|H_{id}\left(\omega ,{\mu }_l\right)\right|=1,\left|\omega \right|\le \pi$$E7

and a linear frequency phase response with a constant phase delay given, respectively, by:

$${\theta }_{id}\left(\omega ,{\mu }_l\right)=-D\omega$$E8

$${\tau }_{pid}\left(\omega ,{\mu }_l\right)=D$$E9

The main goal of all existing FDF design methods, based on a frequency design approach, is to obtain the FDF filter coefficients through approximating this ideal frequency performance.

Applying inverse discrete Fourier transform to the ideal FDF frequency response, the ideal FDF filter unit impulse response h_id(n,μ) is obtained as:

$$h_{id}\left(n,\mu \right)=\frac{\sin \left[\pi \left(n-D\right)\right]}{\pi \left(n-D\right)}=\sin c(n-D)$$E10

Given a desired factional delay value, the FDF coefficients can be easily obtained with this infinite length delayed sinc function. Due to this infinite length, it is evident that an FIR FDF will be always an approximation to the ideal case.

As an illustrative example, the ideal FDF unit impulse responses for two delay values D= 3.0 (D_fix=3.0 and μ = 0) and D=3.65 (D_fix=3.0 and μ = 0.65) are shown in Fig. 3 and 4, respectively. The unit impulse responses were obtained using MATLAB function sinc. The FDF unit impulse responses are shown as solid lines, and the delayed sinc function as dot line. In the first case, only one three-sample delay is needed, which can be easily implemented with memory components as described above. However, the FDF unit impulse response for the second case has an infinite number of nonzero coefficients (IIR) and it is a no causal sequence, which makes it impractical for implementing in real-time applications.

Figure 3.

Ideal FDF unit impulse response for D=3.0.

Figure 4.

Ideal FDF unit impulse response for D=3.65.

3. FDF Design methods

The existing design methods for FIR FDF use a large range of strategies to approximate as close as possible the ideal FDF unit impulse response h_id(n,μ). It is possible to highlight three main strategies:

• Magnitude frequency response approximation: The FDF unit impulse response is obtained such that its frequency magnitude response is as close as possible to the ideal FDF one, accordingly to some defined error criterion.

• Interpolation design method: The design approach is based on computing FDF coefficients through classical mathematical interpolation methods, such as Lagrange or B-spline. The design is a completely time-domain approach.

• Hybrid analogue-digital model approach: The FDF design is accomplished through the use of an analogue-digital model. The design methods using this strategy are based on a frequency-domain approach.

A concise description of each one of these strategies is presented in the following.

3.1. Magnitude frequency response approximation

The design method goal is to obtain the FDF unit impulse response h_FD(n,μ) based on comparing its magnitude frequency response with the ideal one. The frequency response of the designed FDF with even-length N_FD is given by:

$$H_{FD}\left(\omega ,\mu \right)={\displaystyle \sum _{k=-N_{FD}/2+1}^{N_{FD}/2}{h}_{FD}\left(k,\mu \right)}{e}^{-j\omega k}$$E11

One of the criterions used for the magnitude frequency response comparison is the least squares magnitude error defined as:

$$e_2\left(\omega \right)=\frac{1}{\pi }{\displaystyle \underset{0}{\overset{{\omega }_p}{\int }}{\left|H_{DF}\left(\omega ,\mu \right)-H_{id}\left(\omega ,\mu \right)\right|}^{2}d\omega }$$E12

The error function e₂(ω) is minimized by truncating the ideal unit impulse response to N_FD samples, which can be interpreted as applying a delayed M-length window w(n) to the ideal IIR FDF unit impulse response:

$$h_{FD}\left(n,\mu \right)=h_{id}\left(n-D,\mu \right)w\left(n-D\right)$$E13

where ω(n) is equal to unity in the interval 0≤n≤N_FD-1 and zero otherwise.

The windowing process on the ideal unit impulse response causes not-desired effects on the FDF frequency response, in particular the Gibbs phenomenon for rectangular window (Proakis & Manolakis, 1995).

In general, the performance of a FDF obtained by truncating the sinc function is usually not acceptable in practice. As a design example, the FDF frequency magnitude and phase responses for D=3.65, using a rectangular window with N_FD=50, are shown in Fig 5. We can see that the obtained FDF bandwidth is less than 0.9π and although the IIR sinc function has been truncated up to 50 taps, neither its frequency magnitude nor its phase response are constant.

The windowed unit impulse response h_FD(n,μ) has a low-pass frequency response, in this way it can be modified to approximate only a desired pass-band interval (0,απ) as follows:

$$h_{FD}\left(n,\mu \right)=\{\begin{array}{l}\alpha \sin c\left[\alpha \left(n-D\right)\right]\stackrel{}{}\stackrel{}{}for\stackrel{}{}0\le n\le N_{FD}-1\\ 0\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}\stackrel{}{}otherwise\end{array}$$E14

Figure 5.

FDF frequency response for D=3.65 with rectangular window, N_FD=50.

The magnitude and phase responses of a FDF with N_FD= 8 and α=0.5 are shown in Fig. 6, which were obtained using MATLAB. The phase delay range is from D=3.0 to 3.5 samples with an increment of 0.1. More constant phase delay responses and narrower bandwidth is achieved.

Figure 6.

FDF frequency responses using windowing method for D=3.0 to 3.5 with Ν_FD = 8 and α =0.5.

In principle, window-based design is fast and easy. However, in practical applications it is difficult to meet a desired magnitude and phase specifications by adjusting window parameters. In order to meet a variable fractional delay specification, a real-time coefficient update method is required. This can be achieved storing the window values in memory and computing the values of the sinc function on line, but this would require large memory size for fine fractional delay resolution (Vesma, 1999).

The smallest least squares error can be achieved by defining its response only in a desired frequency band and by leaving the rest as a “don’t care” band. This can be done using a frequency-domain weighting as follows (Laakson et al., 1996):

$$e_3\left(\omega \right)=\frac{1}{\pi }{\displaystyle \underset{0}{\overset{{\omega }_p}{\int }}W\left(\omega \right){\left|H_{DF}\left(\omega ,\mu \right)-H_{id}\left(\omega ,\mu \right)\right|}^2}$$E15

where ω_p is the desired pass-band frequency and W(ω) represents the weighting frequency function, which defines the corresponding weight to each band. In this way, the error is defined only in the FDF pass-band, hence the optimization process is applied in this particular frequency range.

In Fig. 7 are shown the FDF frequency responses designed with this method using W(ω)=1, Ν_FD = 8 and α =0.5. We can see a notable improvement in the resulting FDF bandwidth compared with the one obtained using the least square method, Fig. 6.

There is another design method based on the magnitude frequency response approach, which computes the FDF coefficients by minimizing the error function:

$$e_4\left(\omega \right)=\underset{0\le \omega \le {\omega }_p}{\max }\left|H_{FD}\left(\omega ,\mu \right)-H_{id}\left(\omega ,\mu \right)\right|$$E16

The solution to this optimization problem is given by the minimax method proposed by (Oetken, 1979). The obtained FDF has an equiripple pass-band magnitude response. As an illustrative example, the frequency response of an FDF designed through this minimax method is shown in Fig. 8, where N_FD=20 and ω_p=0.9π.

Figure 7.

FDF frequency responses using weighted least square method for D=3.0 to 3.5 with Ν_FD = 8 and α =0.5.

Figure 8.

FDF Frequency responses using minimax method for D=9.0 to 9.5 with Ν_FD = 20 and α =0.9.

3.2. Interpolation design approach

Instead of minimizing an error function, the FDF coefficients are computed from making the error function maximally-flat at ω=0. This means that the derivatives of an error function are equal to zero at this frequency point:

$${\frac{{\partial }^n{e}_c\left(\omega \right)}{\partial {\omega }^n}|}_{\omega =0}=0,\stackrel{}{}n=0,1,2,\mathrm{....}{N}_{FD}-1$$E17

the complex error function is defined as:

$$e_c\left(\omega \right)=H_{FD}\left(\omega ,{\mu }_l\right)-H_{id}\left(\omega ,{\mu }_l\right)$$E18

where H_FD(ω,μ_l) is the designed FDF frequency response, and H_id(ω,μ_l) is the ideal FDF frequency response, given by equation (Eq. 6). The solution of this approximation is the classical Lagrange interpolation formula, where the FDF coefficients are computed with the closed form equation:

$$h_L\left(n\right)={\displaystyle \prod _{\begin{array}{l}k=0\\ k\ne n\end{array}}^{N_{FD}}\frac{D-k}{n-k}}\stackrel{}{}n=0,1,2,\mathrm{....}{N}_{FD}$$E19

where N_FD is the FDF length and the desired delay$$D=\lfloor N_{FD}/2\rfloor +{\mu }_l$$. We can note that the filter length is the unique design parameter for this method.

The FDF frequency responses, designed with Lagrange interpolation, with a length of 10 are shown in Fig. 9. As expected, a flat magnitude response at low frequencies is presented; a narrow bandwidth is also obtained.

Figure 9.

FDF Frequency responses using Lagrange interpolation for D=4.0 to 4.5 with Ν_FD = 10.

The use of this design method has three main advantages (Laakson et al., 1994): 1) the ease to compute the FDF coefficients from one closed form equation, 2) the FDF magnitude frequency response at low frequencies is completely flat, 3) a FDF with polynomial-defined coefficients allows the use of an efficient implementation structure called Farrow structure, which will be described in section 3.3.

On the other hand, there are some disadvantages to be taken into account when a Lagrange interpolation is used in FDF design: 1) the achieved bandwidth is narrow, 2) the design is made in time-domain and then any frequency information of the processed signal is not taken into account; this is a big problem because the time-domain characteristics of the signals are not usually known, and what is known is their frequency band, 3) if the polynomial order is N_FD; then the FDF length will be N_FD, 4) since only one design parameter is used, the design control of FDF specifications in frequency-domain is limited.

The use of Lagrange interpolation for FDF design is proposed in (Ging-Shing & Che-Ho, 1990, 1992), where closed form equations are presented for coefficients computing of the desired FDF filter. A combination of a multirate structure and a Lagrange-designed FDF is described in (Murphy et al., 1994), where an improved bandwidth is achieved.

The interpolation design approach is not limited only to Lagrange interpolation; some design methods using spline and parabolic interpolations were reported in (Vesma, 1995) and (Erup et al., 1993), respectively.

3.3. Hybrid analogue-digital model approach

In this approach, the FDF design methods are based on the hybrid analogue-digital model proposed by (Ramstad, 1984), which is shown in Fig. 10. The fractional delay of the digital signal x(n) is made in the analogue domain through a re-sampling process at the desired time delay t_l. Hence a digital to analogue converter is taken into account in the model, where a reconstruction analog filter h_a(t) is used.

Figure 10.

Hybrid analogue-digital model.

An important result of this modelling is the relationship between the analogue reconstruction filer h_a(t) and the discrete-time FDF unit impulse response h_FD(n,μ), which is given by:

$$h_{FD}\left(n,\mu \right)=h_a\left(\left(n+{\mu }_l\right)T\right)$$E20

where n=-N_FD/2,-N_FD/2+1,…., N_FD/2-1, and T is the signal sampling frequency. The model output is obtained by the convolution expression:

$$y\left(l\right)={\displaystyle \sum _{k=0}^{N_{FD}-1}x\left(n_l-k+N_{FD}/2\right)}{h}_a\left(\left(k+{\mu }_l-N_{FD}/2\right)T\right)$$E21

This means that for a given desired fractional value, the FDF coefficients can be obtained from a designed continuous-time filter.

The design methods using this approach approximate the reconstruction filter h_a(t) in each interval of length T by means of a polynomial-based interpolation as follows:

$$h_a\left(\left(n+{\mu }_l\right)T\right)={\displaystyle \sum _{m=0}^M{c}_m\left(n\right){\mu }_l{}^m}$$E22

for k=-N_FD/2,-N_FD/2+1,…., N_FD/2-1. The c_m(k)’s are the unknown polynomial coefficients and M is the polynomials order.

If equation (Eq. 22) is substituted in equation (Eq. 21), the resulted output signal can be expressed as:

$$y\left(l\right)={\displaystyle \sum _{m=0}^M{v}_m\left(n_l\right){\mu }_l{}^m}$$E23

where:

$$v_m\left(n_l\right)={\displaystyle \sum _{k=0}^{N_{FD}-1}x\left(n_l-k+N_{FD}/2\right)c_m\left(k-N_{FD}/2\right)}$$E24

are the output samples of the M+1 FIR filters with a system function:

$$C_m\left(z\right)={\displaystyle \sum _{k=0}^{N_{FD}-1}{c}_m\left(k-N_{FD}/2\right)z^{-k}}$$E25

The implementation of such polynomial-based approach results in the Farrow structure, (Farrow, 1988), sketched in Fig. 11. This implementation is a highly efficient structure composed of a parallel connection of M+1 fixed filters, having online fractional delay value update capability. This structure allows that the FDF design problem be focused to obtain each one of the fixed branch filters c_m(k) and the FDF structure output is computed from the desired fractional delay given online μ_l.

The coefficients of each branch filter C_m(z) are determined from the polynomial coefficients of the reconstruction filter impulse response h_a(t). Two mainly polynomial-based interpolation filters are used: 1) conventional time-domain design such as Lagrange interpolation, 2) frequency-domain design such as minimax and least mean squares optimization.

Figure 11.

Farrow structure.

As were pointed out previously, Lagrange interpolation has several disadvantages. A better polynomial approximation of the reconstruction filter is using a frequency-domain approach, which is achieved by optimizing the polynomial coefficients of the impulse response h_a(t) directly in the frequency-domain. Some of the design methods are based on the optimization of the discrete-time filter h_FD(n,μ_l)) and others on making the optimization of the reconstruction filter h_a(t). Once that this filter is obtained, the Farrow structure branch filters c_m(k) are related to h_FD(n,m_l) using equations (Eq. 20) and (Eq. 22). One of main advantages of frequency-domain design methods is that they have at least three design parameters: filter length N_FD, interpolation order M, and pass-band frequency ω_p.

There are several methods using the frequency design method (Vesma, 1999). In (Farrow, 1988) a least-mean-squares optimization is proposed in such a way that the squared error between H_FD(ω,μ_l) and the ideal response H_id(ω,μ_l) is minimized for 0≤ω≤ω_p and for 0≤μ_l<1. The design method reported in (Laakson et al., 1995) is based on optimizing c_m(k) to minimize the squared error between h_a(t) and the h_FD(n,μ_l) filters, which is designed through the magnitude frequency response approximation approach, see section 3.1. The design method introduced in (Vesma et al., 1998) is based on approximating the Farrow structure output samples v_m(n_l) as an m^th order differentiator; this is a Taylor series approximation of the input signal. In this sense, C_m(ω) approximates in a minimax or L₂ sense the ideal response of the m^th order differentiator, denoted as D_m(ω), in the desired pass-band frequencies. In (Vesma & Saramaki, 1997) the designed FDF phase delay approximates the ideal phase delay value μ_l in a minimax sense for 0≤ω≤ω_p and for 0≤μ_l<1 with the restriction that the maximum pass-band amplitude deviation from unity be smaller than the worst-case amplitude deviation, occurring when μ=0.5.

4. FDF Implementation structures

As were described in section 3.3, one of the most important results of the analogue-digital model in designing FDF filters is the highly efficient Farrow structure implementation, which was deduced from a piecewise approximation of the reconstruction filter through a polynomial based interpolation. The interpolation process is made as a frequency-domain optimization in most of the existing design methods.

An important design parameter is the FDF bandwidth. A wideband specification, meaning a pass-band frequency of 0.9π or wider, imposes a high polynomial order M as well as high branch filters length N_FD. The resulting number of products in the Farrow structure is given by N_FD(M+1)+M, hence in order to reduce the number of arithmetic operations per output sample in the Farrow structure, a reduction either in the polynomial order or in the FDF length is required.

Some design approaches for efficient implementation structures have been proposed to reduce the number of arithmetic operations in a wideband FDF. A modified Farrow structure, reported in (Vesma & Samaraki, 1996), is an extension of the polynomial based interpolation method. In (Johansson & Lowerborg, 2003), a frequency optimization technique is used a modified Farrow structure achieving a lower arithmetic complexity with different branch filters lengths. In (Yli-Kaakinen & Saramaki, 2006a, 2006a, 2007), multiplierless techniques were proposed for minimizing the number of arithmetic operations in the branch filters of the modified Farrow structure. A combination of a two-rate factor multirate structure and a time-domain designed FDF (Lagrange) was reported in (Murphy et al., 1994). The same approach is reported in (Hermanowicz, 2004), where symmetric Farrow structure branch filters are computed in time-domain with a symbolic approach. A combination of the two-rate factor multirate structure with a frequency-domain optimization process was firstly proposed in (Jovanovic-Docelek & Diaz-Carmona, 2002). In subsequence methods (Hermanowicz & Johansson, 2005) and (Johansson & Hermanowicz &, 2006), different optimization processes were applied to the same multirate structure. In (Hermanowicz & Johansson, 2005), a two stage FDF jointly optimized technique is applied. In (Johansson & Hermanowicz, 2006) a complexity reduction is achieved by using an approximately linear phase IIR filter instead of a linear phase FIR in the interpolation process.

Most of the recently reported FDF design methods are based on the modified Farrow structure as well as on the multirate Farrow structure. Such implementation structures are briefly described in the following.

4.1. Modified Farrow structure

The modified Farrow structure is obtained by approximating the reconstruction filter with the interpolation variable 2μ_l -1 instead of μ_l in equation (Eq. 22):

$$h_a\left(\left(n+{\mu }_l\right)T\right)={\displaystyle \sum _{m=0}^M{c}_m^{\text{'}}\left(k\right){\left(2{\mu }_l-1\right)}^m}$$E26

for k=-N_FD/2,-N_FD/2+1,…., N_FD/2-1. The first four basis polynomials are shown in Fig. 12. The symmetry property h_a(-t)= h_a(t) is achieved by:

$$c_m^{\text{'}}\left(n\right)={\left(-1\right)}^m{c}_m^{\text{'}}\left(-n-1\right)$$E27

for m= 0, 1, 2,…,M and n=0, 1,….,N_FD/2. Using this condition, the number of unknowns is reduced to half.

The reconstruction filter h_a(t) can be now approximated as follows:

$$h_a\left(t\right)={\displaystyle \sum _{n=0}^{N_{FD}/2}{\displaystyle \sum _{m=0}^M{c}_m^{\text{'}}\left(n\right)g\left(n,m,t\right)}}$$E28

where c_m(n) are the unknown coefficients and g(n,m,t)’s are basis functions reported in (Vesma & Samaraki, 1996).

Figure 12.

Basis polynomials for modified Farrow structure for 0≤ m ≤ 3.

The modified Farrow structure has the following properties: 1) polynomial coefficients c_m(n) are symmetrical, according to equation (Eq. 27); 2) The factional value μ_l is substituted by 2μ_l -1, the resulting implementation of the modified Farrow structure is shown in Fig. 13; 3) the number of products per output sample is reduced from N_FD(M+1)+M to N_FD(M+1)/2+M.

The frequency design method in (Vesma et al., 1998) is based on the following properties of the branch digital filters C_m(z):

• The FIR filter C_m(z), 0≤m≤M, in the original Farrow structure is the m^th order Taylor approximation to the continuous-time interpolated input signal.

• In the modified Farrow structure, the FIR filters C’_m(z) are linear phase type II filters when m is even and type IV when m is odd.

Each filter C_m(z) approximates in magnitude the function K_mw^m, where K_m is a constant. The ideal frequency response of an m^th order differentiator is (jω)^m, hence the ideal response of each C_m(z) filter in the Farrow structure is an m^th order differentiator.

In same way, it is possible to approximate the input signal through Taylor series in a modified Farrow structure for each C’_m(z), (Vesma et al., 1998). The m^th order differential approximation to the continuous-time interpolated input signal is done through the branch filter C’_m(z), with a frequency response given as:

Figure 13.

Modified Farrow structure.

$$C{\text{'}}_m\left(\omega \right)=e^{j\omega \left(N_{FD}-1\right)/2}\frac{{\left(-j\omega \right)}^m}{2^{m}m!}$$E29

The input design parameters are: the filter length N_FD, the polynomial order M, and the desired pass-band frequency ω_p.

The N_FD coefficients of the M+1 C’_m(z) FIR filters are computed in such a way that the following error function is minimized in a least square sense through the frequency range [0,ω_p]:

$$e_m\left(\omega \right)=\left|{\displaystyle \sum _{n=o}^{N_{FD}/2-1}\left[c_m\left(N_{FD}/2-n\right)\right]\gamma \left(m,n,\omega \right)-D\left(m,\omega \right)}\right|$$E30

where:

$$\begin{array}{l}D\left(m,\omega \right)=\frac{{\left(-\omega \right)}^m}{2^{m}m!},\underset{}{}\gamma \left(m,n,\omega \right)=2\cos \left[\left(n+1/2\right)\omega \right],\underset{}{}m\underset{}{}even\\ \underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\underset{}{}\gamma \left(m,n,\omega \right)=2\sin \left[\left(n+1/2\right)\omega \right],\underset{}{}m\underset{}{}odd\end{array}$$E31

Hence the objective function is given as:

$$E_1={\displaystyle \underset{0}{\overset{{\omega }_p}{\int }}\left[{\displaystyle \sum _{n=0}^{N_{FD}/2-1}{c}_m\left(N_{FD}/2-1\right)}\gamma \left(m,n,\omega \right)-D\left(m,\omega \right)\right]}d\omega$$E32

From this equation it can be observed that the design of a wide bandwidth FDF requires an extensive computing workload. For high fractional delay resolution FDF, high precise differentiator approximations are required; this imply high branch filters length, N_FD, and high polynomial order, M. Hence a FDF structure with high number of arithmetic operations per output sample is obtained.

4.2. Multirate Farrow structure

A two-rate-factor structure in (Murphy et al., 1994), is proposed for designing FDF in time-domain. The input signal bandwidth is reduced by increasing to a double sampling frequency value. In this way Lagrange interpolation is used in the filter coefficients computing, resulting in a wideband FDF.

The multirate structure, shown in Fig. 14, is composed of three stages. The first one is an upsampler and a half-band image suppressor H_HB(z) for incrementing twice the input sampling frequency. Second stage is the FDF H_DF(z), which is designed in time-domain through Lagrange interpolation. Since the signal processing frequency of H_DF(z) is twice the input sampling frequency, such filter can be designed to meet only half of the required bandwidth. Last stage deals with a downsampler for decreasing the sampling frequency to its original value. Notice that the fractional delay is doubled because the sampling frequency is twice. Such multirate structure can be implemented as the single-sampling-frequency structure shown in Fig. 15, where H₀(z) and H₁(z) are the first and second polyphase components of the half-band filter H_HB(z), respectively. In the same way H_FD0(z) and H_FD1(z) are the polyphase components of the FDF H_FD(z) (Murphy et al, 1994).

The resulting implementation structure for H_DF(z) designed as a modified Farrow structure and after some structure reductions (Jovanovic-Dolecek & Diaz-Carmona, 2002) is shown in Fig. 16. The filters C_m,0(z) and C_m,1(z) are the first and second polyphase components of the branch filter C_m(z), respectively.

Figure 14.

FDF Multirate structure.

Figure 15.

Single-sampling-frequency structure.

Figure 16.

Equivalent single-sampling-frequency structure.

The use of the obtained structure in combination with a frequency optimization method for computing the branch filters C_m(z) coefficients was exploited in (Jovanovic-Dolecek & Diaz-Carmona, 2002). The approach is a least mean square approximation of each one of the m^th differentiator of input signal, which is applied through the half of the desired pass-band. The resulting objective function, obtained this way from equation (Eq. 32), is:

$$E_2={\displaystyle \underset{0}{\overset{{}_{\frac{{\omega }_p}{2}}}{\int }}\left[{\displaystyle \sum _{n=0}^{N_{FD}/2-1}{c}_m\left(N_{FD}/2-1\right)}\gamma \left(m,n,\omega \right)-D\left(m,\omega \right)\right]}d\omega$$E33

The decrease in the optimization frequency range allows an abrupt reduction in the coefficient computation time for wideband FDF, and this less severe condition allows a resulting structure with smaller length of filters C_m(z).

The half-band H_HB(z) filter plays a key role in the bandwidth and fractional delay resolution of the FDF filter. The higher stop-band attenuation of filter H_HB(z), the higher resulting fractional delay resolution. Similarly, the narrower transition band of H_HB(z) provides the wider resulting bandwidth.

In (Ramirez-Conejo, 2010) and (Ramirez-Conejo et al., 2010a), the branch filters coefficients c_m(n) are obtained approximating each m^th differentiator with the use of another frequency optimization method. The magnitude and phase frequency response errors are defined, for 0≤w≤w_p and 0≤μ_l≤1, respectively as:

$$e_{mag}\left(\omega \right)=\left|H_{FD}\left(\omega \right)\right|-1,\underset{}{}$$E34

$$e_{pha}\left(\omega \right)=-\frac{\varphi \left(\omega \right)}{\omega }-\left(D_{fix}+{\mu }_l\right),\underset{}{}$$E35

where H_FD(ω) and ϕ(ω) are, respectively, the frequency and phase responses of the FDF filter to be designed. In the same way, this method can also be extended for designing FDF with complex specifications, where the complex error used is given by equation (Eq. 18).

The coefficients computing of the resulting FDF structure, shown in Fig. 16, is done through frequency optimization for global magnitude approximation to the ideal frequency response in a minimax sense. The objective function is defined as:

$${\Delta }_m=\underset{0\le {\mu }_l\le 1}{\max }\left[\underset{0\le \omega \le {\omega }_p}{\max }\left|e_m\left(\omega \right)\right|\right]$$E36

The objective function is minimized until a magnitude error specification δ_m is met. In order to meet both magnitude and phase errors, the global phase delay error is constrained to meet the phase delay restriction:

$${\Delta }_p=\underset{0\le {\mu }_l\le 1}{\max }\left[\underset{0\le \omega \le {\omega }_p}{\max }\left|e_p\left(\omega \right)\right|\right]\le {\delta }_p$$E37

where δ_p is the FDF phase delay error specification. The minimax optimization can de performed using the function fminmax available in the MATLAB Optimization Toolbox.

As is well known, the initial solution plays a key role in a minimax optimization process, (Johansson & Lowenborg, 2003), the proposed initial solution is the individual branch filters approximations to the m^th differentiator in a least mean squares sense, accordingly to (Jovanovic-Delecek & Diaz-Carmona, 2002):

$$E_m={\displaystyle \underset{0}{\overset{\raisebox{1ex}{${\omega }_p$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\int }}{\left[e_m\left(\omega \right)\right]}^{2}d\omega }$$E38

The initial half-band filter H_HB(z) to the frequency optimization process can be designed as a Doph-Chebyshev window or as an equirriple filter. The final hafband coefficients are obtained as a result of the optimization.

The fact of using the proposed optimization process allows the design of a wideband FDF structure with small arithmetic complexity. Examples of such designing are presented in section 5.

An implementation of this FDF design method is reported in (Ramirez-Conejo et al., 2010b), where the resulting structure, as one shown in Fig. 16, is implemented in a reconfigurable hardware platform.

5. FDF Design examples

The results obtained with FDF design methods described in (Diaz-Carmona et al., 2010) and (Ramirez-Conejo et al., 2010) are shown through three design examples, that were implemented in MATLAB.

Example 1:

The design example is based on the method described in (Diaz-Carmona et al., 2010). The desired FDF bandwidth is 0.9π, and a fractional delay resolution of 1/10000.

A half-band filter H_HB(z) with 241 coefficients was used, which was designed with a Dolph-Chebyshev window, with a stop-band attenuation of 140 dBs. The design parameters are: M=12 and N_FD=10 with a resulting structure arithmetic of 202 products per output sample.

The frequency optimization is applied up to only ω_p=0.45π, causing a notably computing workload reduction, compared with an optimization on the whole desired bandwidth (Vesma et al., 1998). As a matter of comparison, the MATLAB computing time in a PC running at 2GHz for the optimization on half of the desired pass-band is 1.94 seconds and 110 seconds for the optimization on the whole pass-band. The first seven differentiator approximations for both cases are shown in Fig. 17 and Fig. 18.

The frequency responses of the resulted FDF from μ=0.008 to 0.01 samples for the half pass-band and for the whole pass-band optimization process, are shown in Fig. 19 and Fig. 20, respectively.

The use of the optimization process (Vesma et al., 1998) with design parameters of M=12 and N_FD=104 results in a total number of 688 products per output sample. Accordingly to the described example in (Zhao & Yu, 2006), using a weighted least squares design method, an implementation structure with N_FD=67 and M=7 is required to meet ω_p=0.9π, which results in arithmetic complexity of 543 products per output sample.

Figure 17.

Frequency responses of the first seven ideal differentiators (dotted line) and the obtained approximations (solid line) in 0≤ω≤0.45π with N_FD=10 and M=12.

Figure 18.

First seven differentiator ideal frequency responses (dotted line) and obtained approximations (solid line) in 0≤ω≤0.9π with N_FD=104 and M=12.

Figure 19.

FDF frequency responses using half band frequency optimization method for μ_l=0.0080 to 0.0100 with Ν_FD = 10 and M=12.

Figure 20.

FDF frequency responses, using all-bandwidth frequency optimization method for μ_l=0.0080 to 0.0100 with N_FD=104 and M=12.

In order to compare the frequency-domain approximation achieved by the described method with existing design methods results, the frequency-domain absolute error e(ω,μ), the maximum absolute error e_max, and the root mean square error e_rms are defined, like in (Zhao & Yu, 2006), by:

$$e\left(\omega ,\mu \right)=\left|H_{FD}\left(\omega ,\mu \right)-H_{id}\left(\omega ,\mu \right)\right|$$E39

$$e_{\max }=\max \left\{e\left(\omega ,\mu \right)\right\},\underset{}{}0\le \omega \le {\omega }_p,\underset{}{}0\le \mu \le 1$$E40

$$e_{rms}={\left[{\displaystyle \underset{0}{\overset{{\omega }_p}{\int }}{\displaystyle \underset{0}{\overset{1}{\int }}{e}^2\left(\omega ,\mu \right)d\mu d\omega }}\right]}^{1/2}$$E41

The maximum absolute magnitude error and the root mean square error obtained are shown in Table 1, reported in (Diaz-Carmona et al., 2010), as well as the results reported by some design methods.

Method	emax(dBs)	erms
(Tarczynski et al., 1997)	-100.0088	2.9107x10^-6
(Wu-Sheng, & Tian-Bo, 1999)	-100.7215	2.7706x10^-6
(Tian-Bo, 2001)	-99.9208	4.931x10^-4
(Zhao & Yu, 2006)	-99.3669	2.8119x10^-6
(Vesma et al., 1998)	-93.69	4.81x10^-4
(Diaz-Carmona et al., 2010)	-86.17	2.78x10^-4

Table 1.

Magnitude frequency response error comparison.

Example 2:

The FDF is designed using the explained minimax optimization approach applied on the single-sampling-frequency structure, Fig. 16, according to (Ramirez et al., 2010a). The FDF specifications are: ω_p = 0.9π, δ_m = 0.01 and δ_p =0.001, the same ones as in the design example of (Yli-Kaakinen & Saramaki, 2006a). The given criterion is met with N_FD = 7 and M = 4 and a half-band filter length of 55. The overall structure requires Prod = 32 multipliers, Add = 47 adders, resulting in a Δ_m = 0.0094448 and Δ_p = 0.00096649. The magnitude and phase delay responses obtained for μ_l = 0 to 0.5 with 0.1 delay increment are depicted in Fig. 21. The results obtained, and compared with those reported by other design methods, are shown in Table 2. The design described requires less multipliers and adders than (Vesma & Saramaki, 1997), (Johansson & Lowenborg, 2003), the same number of multipliers and nine less adders than (Yli-Kaakinen & Saramaki, 2006a), one more multiplier and three less adders than (Yli-Kaakinen & Saramaki, 2006b), and two more multipliers than (Yli-Kaakinen, & Saramaki, 2007).

Method	Arithmetic complexity
NFD	M	Prod	Add	Δm	Δp
(Vesma & Saramaki, 1997)	26	4	69	91	0.006571	0.0006571
(Johansson, & Lowenborg, 2003)	28	5	57	72	0.005608	0.0005608
(Yli-Kaakinen & Saramaki, 2006a)	28	4	32	56	0.009069	0.0009069
(Yli-Kaakinen & Saramaki, 2006b)	28	4	31	50	0.009742	0.0009742
(Yli-Kaakinen & Saramaki, 2007)	28	4	30	-	0.009501	0.0009501
(Ramirez-Conejo et al.,2010)	7	4	32	47	0.0094448	0.0009664

Table 2.

Arithmetic complexity results for example 2.- Not reported

Figure 21.

FDF frequency responses, using minimax optimization approach in example 2.

Figure 22.

FDF frequency response errors, using minimax optimization approach in example 2.

Example 3:

This example shows that the same minimax optimization approach can be extended for approximating a global complex error. For this purpose, the filter design example described in (Johansson & Lowenborg 2003) is used, which specifications are ω_p = 0.9π, and maximum global complex error of δc = 0.0042. Such specifications are met with N_FD = 7 and M = 4 and a half-band filter length of 69. The overall structure requires Prod = 35 multipliers with a resulting maximum complex error Δ_c = 0.0036195. The results obtained are compared in Table 3 with the reported ones in existing methods. The described method requires less multipliers than (Johansson & Lowenborg 2003), (Hermanowicz, 2004) and case A of (Hermanowicz & Johansson, 2005). Reported multipliers of (Johansson & Hermanowicz, 2006) and case B of (Hermanowicz & Johansson, 2005) are less than the obtained with the presented design method. It should be pointed out that in (Johansson & Hermanowicz, 2006) an IIR half-band filter is used and in case B of (Hermanowicz & Johansson, 2005) and (Johansson & Hermanowicz, 2006) a switching technique between two multirate structures must be implemented. The resulted complex error magnitude is shown in Fig. 23 for fractional delay values from D =17.5 to 18.0 with 0.1 increment, magnitude response of the designed FDF is shown in Fig. 24 and errors of magnitude and phase frequency responses are presented in Fig 25.

Method	Arithmetic complexity
NFD	M	Prod
(Johansson & Lowenborg 2003)	39	6	73
(Johansson & Lowenborg 2003)^a	31	5	50
(Hermanowicz, 2004)	11	6	60(54)
(Hermanowicz & Johansson, 2005)	7	5	36
(Hermanowicz & Johansson, 2005)^b	7	3	26
(Johansson & Hermanowicz, 2006)	-	6	32
(Johansson & Hermanowicz, 2006)^b	-	3	22
(Ramirez-Conejo et al., 2010)	7	4	35

Table 3.

Arithmetic complexity results for example 3.a. Minimax design with subfilters jointly optimized.

Figure 23.

FDF frequency complex error, using minimax optimization approach in example 3.

Figure 24.

FDF frequency response using minimax optimization approach in example 3.

Figure 25.

FDF frequency response errors using minimax optimization approach in example 3.

6. Conclusion

The concept of fractional delay filter is introduced, as well as a general description of most of the existing design methods for FIR fractional delay filters is presented. Accordingly to the explained concepts and to the results of recently reported design methods, one of the most challenging approaches for designing fractional delay filters is the use of frequency-domain optimization methods. The use of MATLAB as a design and simulation platform is a very useful tool to achieve a fractional delay filter that meets best the required frequency specifications dictated by a particular application.

Acknowledgments

Authors would like to thank to the Technological Institute of Celaya at Guanajuato State, Mexico for the facilities in the project development, and CONACYT for the support.