2.1

The implementation principle of digital fractional delay

It is not appropriate to treat the delay of a uniformly sampled band-limited digital signal as a continuous-time signal. To address this issue, the continuous-time signal must be converted into a discrete-time signal through sampling at time instant t = nT_s and delaying by a positive real number D, consisting of an integer and fractional part d [1][2][3]. The resulting delayed discrete signal can be expressed as Eq. (1) and Eq. (2),

Where x[n], y[n] donate the input and output signal (delayed version), L is a linear operator. The integer delay part corresponds to one of the previous outputs, while the fractional delay is challenging due to its potential location between two samples [1][4][13]. Considering the delay as a linear time-invariant operation, the potential location can be achieved by bandlimited interpolation, firstly treating the delay system in z transfer domain, it can be obtained formally as,

However, this approach is only effective for integer delays, and noninteger delays must be approximated in some way. One traditional method is to construct a series expansion for Z^–D, while a more common and effective approach is to formulate the design objective in the frequency domain by substituting Z = e^jω [4]. Therefore, the frequency response of an ideal fractional delay filter is given as follow,

By applying the inverse discrete-time Fourier transforms. The corresponding impulse response is shown below,

From Eq. (6), the ideal output after digital fractional delay in time domain is achieved by spreading it over all the discrete-time signal values, weighted by appropriate values of the sin c[n – D] function. However, the impulse response of the fractional delay (FD) filter h[n] = sin c[n – D] is not only non-causal but also infinite, which is not realistic, the computer does not have an infinite number of buffers to store filter coefficients of infinite length and the time cannot be negative [4]. To account for the engineering value, the most straightforward method is truncating the impulse response of ideal sin c[n – D] function to a finite length [4][7]. Here is equivalent to applying a truncation with length L = N+1, as shown in Eq. (7),

M is the optimal integer time index of the first nonzero value of the impulse response, which can be calculated from Eq. (8),

If utilize a 21-order filter with a fractional delay of 0.25 samples, the value of M can be calculated as -10. After inputting x[n] = {1, 2,3, 4,5}, the output can be generated as the convolution of x[n] and h[n], its impulse response is shown in Figure 1. The result is reasonable. For example, y[0] should be equivalent to x[-0.25]. Since the input signal x[n] increases with n (0 ≤ n ≤ 4), x[-0.25] should be slightly smaller than x[0]=1, y[0]=0.8281 is a reasonable value. Figure 2 shows the magnitude response of the FD filter, which exhibits ripples, especially at the edge, due to the truncation with two discontinuous points applied in the time domain, thereby causes smearing or known as Gibbs phenomenon in frequency domain [1]. The following content proposes solutions to improve digital fractional delay based on truncated sinc function method and focuses on mitigating the Gibbs phenomenon caused by discontinuity.

Figure 1.

Impulse response of factional delay FIR filters in Truncated Sinc Function Method.

Figure 2.

Magnitude response of the fractional delay FIR filters in Truncated Sinc Function Method.

3.1

Windowed sinc functions

The use of window functions for time-domain weighting has been shown to be effective in reducing the phenomenon of amplitude discontinuities at the boundaries of finite sequence [4][9]. A bell-shaped window function can be used to emphasize the middle values of the impulse response and reduce peak magnitude errors. However, this comes at the expense of a wider transition band. The windowed impulse response of the FD FIR filter can be expressed as Eq. (9).

Various window functions are available in signal processing, such as Hanning window, Hamming window, Kaiser window, Chebyshev window, Blackman window, among others. While Hanning and Hamming windows are simple and easy to use, they cannot explicitly control stop-band attenuation and main-lobe width. Kaiser and Chebyshev window can help control them using a parameter. Chebyshev window can achieve lower stopband attenuation with the same window length compared to Kaiser window [14]. This results in smaller peak sidelobe ripple, but at the cost of a wider main-lobe width. A higher main lobe and lower side lobe can bring the frequency response of the FD filter closer to its ideal impulse response.

3.2

Smooth transition functions

Alternatively, applying smooth transition functions is possible to achieve the same outcome by modifying the desired magnitude response. To eliminate the Gibbs phenomenon while maintaining optimal approximation, a smooth transition function is utilized to establish a continuous between passband and stopband, which can effectively eliminate the discontinuity associated with the rectangular window. By choosing a suitable smooth transition function, the impulse response of the FD filter can asymptotically decay as n increases, resulting in reduced approximation error without windowing explicitly [5]. In smooth transition function, different responses are formulated for different bands, with the transition band being defined by a Pth-order spline [15], as outlined in Table 1.

Table 1.

Response in different bands for Smooth Transition Function design.

If a first-order spline is employed as the transition function, its magnitude response |H₁(ω)| can be expressed as Eq. (10).

By applying inverse Fourier transform to Eq. (10), the impulse response can be obtained as follows.

If a second-order spline is utilized as the transition function, its frequency response can be expressed as the convolution of an ideal no-transition-band frequency response and a triangle pulse with a width of (ω_s – ω_p). This can be considered as the convolution of two rectangular pulses each with a width of [4][16]. To achieve this in the time domain, one can multiply Eq. (11) by the square of the inverse transform of the half width-pulses [2][3]. Therefore, the Pth order spline can be easily expended, as demonstrated in Eq. (12),

After delaying D samples, the impulse response of FD FIR filter can be derived as follow:

The spline power, donated by P, can be utilized to regulate the character of the approximation.

3.3

Complex general ls design

Theoretically, FIR filters can be designed by minimizing the least weighted squared error. The fractional delay filter with the smallest LS (Least Square) error within a specific approximate band can define the response only in this frequency band while disregarding any errors generated in other frequency bands as “don’t care” band [4][5][10]. The general LS method, in comparison to previous methods, can effectively minimize the waste of resources outside the approximate bandwidth. The error formulation of General Least Square design method is provided in Eq. (14).

Where E represents the Least Square Error and is defined within the frequency band of [0,απ]. The non-negative weighting function in the frequency domain is donated as W(ω). For simplicity, it is assumed to be equal to 1 in this paper, signifying that all frequencies are of equal importance within the approximate bandwidth. To obtain the filter coefficients that can minimize the LS error, the partial derivatives with respect to filter coefficient can be set to zero. After derivation and simplification, filter coefficients h can be derived as follows [4][12]:

Where the parameters of P and p1 can be formulated as Eq. (16) and Eq. (17) [4]:

L donates for the filter length.

The technique under consideration can be extended to frequency discretization designs. If appropriate trapezoidal integration method applied formulas (16) and (17), along with the division of the frequency range [0,απ] into n segments for computation, it will enable effective frequency-domain discretization. However, a potential challenging of this discretization is that it may lead to prolonged computation time.

3.4

Maximally flat design: lagrange interpolation

In contrast to the approach of minimizing the weighted Least Square error, an alternative method for achieving optimal approximation is maximizing the flattening of the error function at a specific frequency, typically at ω₀ = 0 [4][5][8].

This can be obtained by setting the derivatives of error function at frequency ω = ω₀ to zero [4][8],

Where E(e^jω) is the complex error function, has a desired response H_id (e^jω) = e^–jωD. Then it can be expressed in the form of impulse response (ω₀ = 0),

or using matrix,

h donates for the vector of FD FIR filter coefficients, V is an LxL Vandermonde matrix and v is the vector of delay.

The classical Lagrange interpolation formula is a suitable approach to tackle the resolution of Eq. (21). By applying this method, a straightforward solution for an FD FIR filter coefficient can be obtained, this is demonstrated in the following equation,

Where D stands for the fractional delay, N is the order of the filter.

4.1

Windowed sinc functions

The implementation of the FD FIR filter with Chebyshev, Hamming, Hanning and Blackman windows is demonstrated, and their impulse response y[n]=x[n-0.25] is shown in Table 2, the magnitude response is illustrated in Figure 3. It is observed that the impulse response is reasonable and windowed filter significantly reduces ripples in the magnitude response compared to truncated sinc function method. This is attributed to the reduction of discontinuities at the boundaries, thereby mitigating the Gibbs phenomenon. While hamming and hanning windows are straightforward to use and yields comparable amplitude responses, they do not allow for explicit control over the stopband attenuation and mainlobe width. On the other hand, the implementation of Blackman window produces almost no ripple but a narrow passband, approximately 0.75 (normalized), while Chebyshev windows exhibits a wider passabnd (around 0.83) but more ripples and slight Gibbs phenomenon at the edge.

Figure 3.

The magnitude response of FD FIR filters with Chebyshev, hamming, hanning and blackman windowed method.

Table 2.

Impulse response of FD FIR filters with different windowed sinc design.

In summary, the window-based design approach is an efficient and expeditious approach that aids in mitigating the Gibbs phenomenon, albeit at the expense of a wider transition band. This method is particularly suitable for real-time coefficient updates if multiple samples of the windowed sinc function are stored in the memory for interpolation filter coefficients computation [5][11]. However, controlling the magnitude error through window parameter adjustments, particularly for filters with shorter lengths, is a challenging task when employing window-based methods.

4.2

Smooth transition functions

The simulation result shows that the impulse response is y[n] = {0.6006,1.862, 2.698,3.67,5.202}, which is reasonable. The magnitude response of the FD FIR filter with Pth-order spline in transition band exhibits significantly lower ripples compared to the implementation of Chebyshev and hanning windows, as depicted in Figure 4. This is due to the smooth transition function, which continuously connects passband and stopband, removing the discontinuity present in the rectangular window. However, the smooth transition function method has a narrower passband, approximately 0.7.

Figure 4.

The impulse response (left) and magnitude response (right) of FD FIR filters with smooth transition function method.

In summary, the smooth transition function method is employed to reduce Gibbs phenomenon by smoothing the transition band. The use of Pth-order spline offers the simplicity of window methods and explicit control of transition band, promising least squared error optimality [16]. Moreover, this method is more suitable for online coefficient calculations. It is important to note that the application of optimal P, p1 and transition bandwidth significantly affect the performance of filter. However, the frequency response should be explicitly designed in the whole Nyquist band, and the impulse response is still infinitely long, albeit decaying fast and truncation always leads to some error.

4.3

Complex general LS design

Figure 5 demonstrates the impulse response and the magnitude response of FD FIR filters with Complex General LS design and Complex General LS design by trapezoidal rule method. The impulse response of Complex LS design is y[n] = {0.7619,1.672, 2.907,3.457,5.405} and Complex General LS design by trapezoidal rule method is y[n] = {0.7622,1.672, 2.907,3.457,5.405}, which are almost same. The results are satisfied. It can be observed that the magnitude response of these two design exhibits minimal ripple and almost no Gibbs phenomenon. Additionally, the peak magnitude error has been significantly reduced. This can be attributed to the emphasis placed by this method on the response in particular part of the frequency band within the error measurement, which at the expense of an increased error outside the approximation band. Furthermore, this method facilitates a wider passband, approximately 0.8. The design with trapezoidal rule works slightly worse at the edge, that possibly resulted from the error generated by frequency discretization.

Figure 5.

The impulse response (left) and the magnitude response (right) of FD FIR filters with Complex General LS design and Complex General LS design by trapezoidal rule method.

In summary, the Complex General LS design method offers a closed-form solution by minimizing the least weighted squared error in the frequency domain, leveraging matrix properties in the design process. This method enables efficient computation of filter coefficients, making it suitable for designing various types of 2D FIR filters without imposing symmetry constraints on their frequency response [6]. However, this method incurs higher computational cost and prolonged computation time due to the increased desired number of additions and multiplications during design process, as matrix inversion is typically proportional to (N + 1)³ [12]. It also suffers from numerical errors in the approximation process, particularly in narrowband approximation.

It worth noting that our design employs a simple implementation with unity weight function W(ω) = 1, which only roughly characterizes the performance of the signal response. In practice, signals are usually random and characterized by their average power spectrum of input signal x[n] instead of W(ω) = 1 to obtain a smaller LS error.

4.4

Maximally flat design: lagrange interpolation

In contrast to the previous results, the magnitude response of the FD filter with Langrage Interpolation method exhibits smoothness and minimal ripples, as shown in Figure 6. However, the bandwidth of the magnitude response is relatively narrow bandwidth, approximately 0.6, as the flatness constraint is only enforced at ω₀ = 0. The impulse response is reasonable which is y[n] = {0.7105,1.737, 2.834,3.524,5.363 }.

Figure 6.

The impulse response (left) and magnitude response (right) of FD FIR filter with Lagrange interpolation method.

In conclusion, the design method offers simple formulas for computing coefficients, avoiding numerical problems encountered in the Complex General LS method. Additionally, this method yields excellent performance at low frequencies, resulting in a smooth magnitude response and small approximation error. However, the response at high frequencies is compromised. Overall, this design method is suitable for lower-order FD FIR filter design at low frequencies.