Improving accuracy of FIR filters for computing convolution transforms for smooth non-bandlimited signals

We propose to improve the accuracy of FIR filters for computing convolution transforms for smooth non-bandlimited (NBL) signals by designing filters by the identification method with using a pair of bandlimited portions of the chosen NBL input and output signals related with each other by the given transform. A design example of type IV linear phase differentiator is presented, where filter coefficients are calculated from the bandlimited portions of the Cauchy pulse and its derivative. The performance of the designed differentiator is evaluated by comparing the accuracy of computed derivatives for several smooth NBL test signals, such as the Cauchy pulse, the Hilbert transform of Cauchy pulse, the Gaussian function, as well as for a bandlimited sinc-function. Evaluation results demonstrate that the proposed design method generates more accurate differentiators than the Parks-McClellan algorithm, the impulse response truncation method and the identification with using full-band Cauchy pulse and its derivative.


Introduction
Many branches of science and technology, such as material science [1][2][3], relaxation spectroscopy [4,5], geophysics and geophysical prospecting [6,7], etc. face with measuring and processing monotonic or locally monotonic signals, which can be classified as ones belonging to the class of smooth non-bandlimited (NBL) signals. Nevertheless, discrete-time processing of smooth NBL signals is not sufficiently developed yet and has not received much attention in the literature. Currently, most of the research in this area is focused on sampling and reconstruction of NBL signals [8][9][10][11] with a goal to overcome the aliasing effects. In routine practice, however, NBL signals are still often processed by the traditional methods based on bandlimited processing philosophy [12][13][14] without taking into account the non-bandlimitedness of the signals, which results that only a bandlimited portion (BLP) of a NBL signal below the Nyquist frequency is processed, but an out-of-band portion (OBP) above the Nyquist frequency is ignored.
The goal of this study was to investigate how to improve the accuracy of finite impulse response (FIR) filters to compute convolution transforms for smooth NBL signals as accurately as possible with taking into account the non-bandlimitedness of the signals.
The study was motivated by: (i) the significance of convolution transforms and NBL signals for solving relevant technical problems, and (ii) efficiency of FIR filters for computing transforms. The significance of convolution transforms is obvious from the fact that many of classical integral transforms, such as Laplace, Fouriersine, Fourier-cosine, Hankel, Meier, etc. are either in the form of convolution transform or can be put into it by change of variables [15]. Convolution transforms and their inversions with NBL signals are used to solve numerous problems of mathematical physics, they arise in considering many physical systems within the framework of the principle of superposition. Efficiency of FIR filters comes from their operation without interpolation and numerical integration, as well as from the well-developed theory and practice [12,14,16].
The paper is structured as follows. In next Section, we discuss the use of discrete-time filters for computing convolution transforms, describe the basics behind the proposed accuracy improvement of FIR filters, and demonstrate it in an example of designing IV type differentiator. The evaluation results of the designed and some other differentiators are presented in Section 3. Finally, conclusions are given in Section 4.

Discrete-time filters for computing convolution transforms
The idea of computing convolution transforms and their inversions by discrete-time filters is not new. A large amount of examples can be found in literature, where discrete-time (digital) filters have been constructed for computing convolution transforms with the aim to make the computation more effective. The two major areas, where discrete-time filters have been most used for computing integral transforms probably are geophysical electromagnetic prospecting [17][18][19] and material science [20]. In modern signal processing [12][13][14]16], a discrete-time differentiator [21] inverting a convolution transform with the step-function kernel, likely, is the In computing convolution transforms, the shapes of convolved waveforms usually carry information, so attainment of accurate as possible waveforms generally is of primary importance. So, design problem of a filter for computing a convolution transform may be formulated as finding an impulse response that produces convolved waveforms as accurately as possible under the specific processing conditions (filter's length, sampling frequency, etc.). This design differs from that of conventional frequency selective filters [12][13][14]16] intended for removing unwanted frequency parts or extracting useful parts of a signal, where filter coefficients must be found to provide some desired frequency response.
Filters designed in the frequency domain usually are not optimal for computing convolution transforms, in the sense that they do not produce maximum accurate convolved waveforms potentially available for the given signal and the filter's length stated at the specified sampling frequency. A main reason of this is a lack of knowledge how the frequency response of a discretetime filter should deviate from the ideal frequency response to produce a convolved waveforms as accurately as possible. In other words, there are limited possibilities to formulate an optimal design specification in the frequency domain and so to design an optimal filter for computing a convolution transform.
In specific areas, such as geophysical prospecting [17][18][19] and material science [20], it is widespread to design FIR convolution filters from a pair of known input-output functions related with each other by the given transform. This approach has been generalized as an identification (ID) method implementing design of FIR filters in the input-output signal domain [22], although it can be also interpreted as a method based on the learning principle [20,23].
An advantage of ID method is that it effectively eliminates various effects, such as data truncation, rounding-off, etc., from which cannot be avoided by the frequency-domain methods. The identification method has been also used for constructing discrete-time algorithms for non-linear transforms, such as computing the real and imaginary parts [24,25] and the relaxation spectrum [26,27] from the magnitude response, where a non-linear transform is approximated by the convolution transform. A main drawback of the ID method often pointed out in the literature [17][18][19] is dependence of the filter coefficients on the pair of input and output functions chosen for the identification.

Idea behind the filter design
According to the recent studies [28,29] of the aliasing and anti-aliasing effects in processing smooth NBL signals, the total error of a discretely processed smooth NBL signal manifests as where () BLP et is an error of BLP processed in the direct way by the discrete-time algorithm used, whereas () OBP et represents an error of non-bandlimitedness caused by OBP. The non-bandlimitedness error is composed from two components [28,29] where () A et is the aliasing error originating from processing a part or full OBP as a signal transformed back to the Nyquist frequency band, but () AA et is the anti-aliasing error arising from removing a part or full OBP by anti-aliasing filtering (AAF) prior sampling. If AAF is not used, complete OBP is transformed back to the Nyquist frequency band and is processed by the discrete-time algorithm used producing the maximum aliasing error and zero anti-aliasing error. Contrary, in the case of ideal AAF with cut-off at the Nyquist frequency, OBP is completely removed yielding the maximum anti-aliasing error and zero aliasing error. So, the non-bandlimitedness error (2) is restricted by the maximum aliasing and anti-aliasing errors. Depending on how fast the signal spectrum decreases, the nonbandlimitedness error normally decays with increasing sampling frequency. As shown in [28,29], the aliasing error eA(t) weakly depends on the filter specification, but the anti-aliasing error eAA(t) is completely algorithmindependent one. There are very limited possibilities to reduce the non-bandlimitedness error by reducing one of the components (2) because reducing of one component, e.g. eA(t) by AAF, causes an increase of the other component eAA(t), and vice versa [29].
Therefore, the non-bandlimitedness error (2) is nearly algorithm-uncontrollable at the fixed sampling frequency, indicating that the total computation error (1) of a smooth NBL signal can be practically decreased only by two means: (i) by increasing sampling frequency, and (ii) by improving the accuracy of processing BLP. Thus, a key idea behind the proposed filter design is finding maximum accurate algorithms for BLP. To materialize the idea, we propose to design FIR filters in the input-output signal domain by the identification method [20,22] with using pairs of the bandlimited portions of NBL input and output functions instead of commonly used full band functions.

Design example of IV type differentiator
In this Sub-section, we demonstrate designing of FIR filters for computing a convolution transforms according to the proposed idea in an example of constructing IV type linear phase differentiator by using a pair of BLPs of the Cauchy pulse [30] and its derivative (Fig. 1 which have the following expressions [28,29] where Ny = S/2 is the angular Nyquist frequency and S is the angular sampling frequency. according to mean-square error (MSE) criterion [28,29].
In Table in Appendix, filter coefficients are given for 12-point type IV differentiator (denoted as IDBL) designed by the proposed method at sampling frequency S = 15 (Ny = 7.5). Here, the filter coefficients are also shown for a differentiator (denoted as ID) constructed by ID method for full-band Cauchy pulse (3) input with minimizing the error between the differentiator's output and exact full-band derivative (4), as well as for the differentiators designed by the impulse response truncation (IRT) method [14] and the Parks-McClellan (PM) algorithm [12,13,31].

Variation of differentiation errors with sampling frequency
We evaluate the total differentiation errors for several test signals through their MSEs where y(t) is exact analytical derivative of a full-band test signal. MSEs (8) have been calculated for M = 100 over the time interval [0, 10].
Since the non-bandlimitedness error (2) highly depends on sampling frequency, and an increase of sampling frequency is the most effective tool for its reducing, we investigate variation of time-domain error (7) and MSE (8) with sampling frequency.
For the Cauchy pulse (3), we evaluate MSEs also separately for differentiation errors of BLPs and non-bandlimitedness errors representing in this case the maximum aliasing errors [29]. In Fig. 2, differentiation MSEs of the Cauchy pulse with sampling frequency are shown for 12-point differentiators designed by the four above mentioned methods. It can be seen that the four differentiators have practically equal non-bandlimitedness error MSEOBP that exponentially decays with sampling frequency. At the same time, the different differentiators exceedingly differently process BLP producing the error that slightly increases with sampling frequency for IRT differentiator, is practically constant for PM differentiator and decays 3 with sampling frequency for ID and IDBL differentiators. How it was expected, IDBL differentiator has the smallest MSEBLP.
The non-bandlimitedness error is greater than BLP errors at low frequencies resulting in its prevailing contribution in the total errors MSEtotal for all the differentiators (Fig. 2(b)). Due to the exponentially decaying MSEOBP and the slower decaying MSEBLP, BLP errors overcome the non-bandlimitedness error at higher sampling frequencies, which are various for different differentiators. Fig. 2 actually witnesses the rightness of the idea behind the proposed design method that the total processing error for smooth NBL signals can be decreased by two means -increasing sampling frequency and improving the processing accuracy of BLPs.
In Fig. 3, total time-domain errors of derivatives of the Cauchy pulse are shown produced by four 12-point differentiators at two sampling frequencies S = 15 and S = 20 indicated by the arrows in Fig. 2. It can be seen that IDBL and ID differentiators produce significantly more accurate derivatives than IRT and PM differentiators. Contrary to IRT and PM differentiators, which have approximately equal time domain errors at both the sampling frequencies, the errors for IDBL and ID differentiators at sampling frequency S = 20 are significantly smaller than the error at S = 15. If IDBL and ID differentiators have practically equal total time-domain errors at S = 15, the IDBL differentiator has the smaller error than the ID differentiator at sampling frequency S = 20 (Fig. 3(c)).
The observed dissimilar behaviour of time-domain errors can be explained by the fact that the different error components of total error (1) are seen in Fig. 3 for the different differentiators. To return to Fig. 2, one can notice that etotal(t) for IRT and PM differentiators in Fig.  3 actually represents error components eBLP(t) that slightly vary with sampling frequency. Contrary to this, non-bandlimitedness errors eOBP(t) are seen in Fig. 3 for both ID differentiators at sampling frequency S = 15, and total errors etotal(t), that are really composed from both the components of eBLP(t) and eOBP(t), are seen for both ID differentiators only at frequency S = 20. Qualitatively similar behaviour of total MSEs with sampling frequency to that of the Cauchy pulse (see Fig.  2(b)) has been observed for some other smooth NBL test signals, such as the Hilbert transform of the Cauchy pulse (  Variation of the differentiation error with sampling frequency for the Hilbert transform of the Cauchy pulse (see Fig. 4) and the Gaussian function (see Fig. 5) shows that the total error is formed from both of components (1) similarly as for the Cauchy pulse (see Fig. 2) with prevailing of MSEOBP at low sampling frequencies and MSEBLP at the higher ones.
For sinc-function (see Fig. 6) as a bandlimited one, we cannot speak about the non-bandlimitedness error, however, here, also OBP appears, when the Nyquist frequency is smaller than the limit frequency of the frequency band with non-zero spectrum. This OBP creates the aliasing error, which as a quickly decaying error is seen in Fig. 6 at low frequencies.

Magnitude responses
In Fig. 7, errors of magnitude responses () Hj    are compared for the four 12-point differentiators. It can be seen that the deviations of magnitude responses from the ideal magnitude response do not reflect the common differentiation accuracy (see Fig. 2 -Fig. 6). This observation confirms the statement made in the Section 2 that a lack of knowledge about optimal design specification in the frequency domain limits the frequency-domain design methods to construct accurate discrete-time algorithms for computing convolution transforms. Designing FIR filters for computing integral transforms in the specific fields [18][19][20] has shown that smooth magnitude responses are required in the vicinity of zero frequency to compute accurate convolved waveforms. The errors of magnitude responses of ID and IDBL differentiators (Fig. 7 (b)) approve this experience and demonstrate that the more accurate IDBL differentiator (see Fig. 2 -Fig. 6) has the smoother magnitude response than that of the ID differentiator.

Conclusions
Based on the recently findings [28,29] that the error of bandlimited portion (BLP) of a discretely processed smooth non-bandlimited (NBL) signal bellow the Nyquist frequency strongly depends on the discretealgorithm algorithm used, but the error from the nonbandlimitedness caused by out-of-band portion (OBP) above the Nyquist frequency is practically algorithmindependent, we propose to improve the accuracy of FIR filters for computing convolution transforms by designing them by the identification method with using bandlimited portions of NBL signals chosen as input and output ones for the identification. The proposed approach is validated by designing type IV linear phase differentiator with using BLPs of the Cauchy pulse and its derivative. The performance of the designed differentiator is evaluated by estimating the time-domain errors of computed derivatives through the mean-square errors (MSEs) at different sampling frequencies for several smooth NBL test signals, such as the Cauchy pulse, the Hilbert transform of Cauchy pulse, the Gaussian function, as well as for a bandlimited test signal -sinc-function. Evaluation results demonstrate that the proposed design method generates more accurate differentiators than other design methods, such as the Parks-McClellan algorithm, the impulse response truncation method and the identification with using the full-band Cauchy pulse and its derivative.