- Email: [email protected]
The Hilbert transform of cubic splines
Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Research paper
The Hilbert transform of cubic splines Mina B. Abd-el-Malek∗, Samer S. Hanna Department of Engineering Mathematics and Physics, Faculty of Engineering, Alexandria University, Alexandria 21544, Egypt
a r t i c l e
i n f o
Article history: Received 12 March 2018 Revised 24 December 2018 Accepted 26 August 2019 Available online 30 August 2019 Keywords: Hilbert transform Splines Piecewise cubic polynomial Hilbert–Huang transform
a b s t r a c t In this work, we propose an efficient method to calculate the Hilbert transform of cubic splines. We start by developing an analytical derivation for the problem. Using the structure of the obtained formula, we suggest an efficient algorithm to evaluate it. Our method rewrites the cubic spline and reorders calculations to increase the temporal locality of related tasks. Then, it creates lookup tables to reuse calculations instead of repeating them. Results show that this method accurately calculates Hilbert transform when applied to functions with a known analytical Hilbert transform. We also show that it significantly reduces the execution time compared to the direct substitution. Our suggested algorithm can be used for the Hilbert–Huang Transform, whose calculation is based on cubic splines and Hilbert transform. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The Hilbert transform (HT) has countless applications in the fields of signal processing, applied mathematics, and physical sciences. It has been used to solve several problems in these fields [1–3]. The significance of some of these problems has led to an interest in efficient methods for the evaluation of the Hilbert transform. Most recently, a transform based on Hilbert transform and cubic splines was introduced. The Hilbert–Huang Transform (HHT) is a time-frequency transform suitable for non-linear and non-stationary signals [4]. HHT uses a method called Empirical Mode Decomposition (EMD) to decompose a given time series into signals with magnitudes and frequencies of close scales called Intrinsic Mode Functions (IMFs). Then, it applies the Hilbert transform to these IMFs, in order to derive a timefrequency representation of the signal. The steps of EMD relies on the use of cubic splines. The Computationally Adaptive Empirical Mode Decomposition (CAEMD) [5] was developed as an enhancement to EMD. When using CAEMD, some of the IMFs are obtained in the form of Piecewise Cubic Polynomials (PCPs). In the second step of HHT, we need to evaluate the Hilbert transform of these PCPs. In this work, we derive an efficient method to calculate the Hilbert transform of a given piecewise cubic polynomial. The Hilbert transform of a PCP is of interest for methods like the Hilbert–Huang transform. To obtain the Hilbert transform of a PCP IMF using any of the existing methods in the literature, the PCPs needs to be evaluated at all points before applying the Hilbert transform. This due to the fact that none of these methods has addressed the problem of obtaining the Hilbert transform of a function given in the form of a PCP. Other than the HHT, piecewise cubic polynomials are useful to create a mathematical model for data [6]. They have applications in many fields like material science [7], robotics [8], and computer graphics [9]. This further motivates the development of a method to calculate the Hilbert transform of piecewise cubic polynomials. ∗
Corresponding author. E-mail addresses: [email protected] (M.B. Abd-el-Malek), [email protected] (S.S. Hanna).
https://doi.org/10.1016/j.cnsns.2019.104983 1007-5704/© 2019 Elsevier B.V. All rights reserved.
2
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
In this work, we present an algorithm that calculates the Hilbert transform of a piecewise cubic spline directly. This algorithm uses the properties of the Hilbert transform and the nature of cubic splines to simplify the calculations. It uses the known Hilbert transform of a cubic polynomial to efficiently calculates the HT of a PCP. Although this method was originally developed to calculate the HT of PCPs obtained from CAEMD, it can be used to calculate the Hilbert transform of any set of points that can be fitted or approximated using a PCP. The rest of this paper is organized as follows. In Section 2, we derive an analytical formula for the Hilbert transform of a PCP. The problem, the obtained formula, and the assumptions of this work are discussed in Section 3. Section 4 starts analyzes the case of a PCP with uniform knots, while Section 5 addresses a more general case. In Section 6, we propose reformulation the PCP to further improve computational efficiency. The proposed algorithm is then presented in Section 7. In Section 8, the suggested method is validated against functions with known Hilbert transform and applied to realistic signals. The execution times of the direct and proposed methods are compared in Section 9. A survey of existing methods for the evaluation of the Hilbert transform are presented in Section 10, while Section 11 concludes this work. 2. Analytical derivation The Hilbert transform of a function f (t ) ∈ L p (R ) (1 ≤ p < ∞) is defined by the following principle value integral
F (x ) = H { f (t )} =
1
π
∞
PV −∞
f (t ) dt x−t
(1)
In this work, we only consider the case where f is a real-valued function belonging to the set of square integrable functions defined on the real line i.e, f ∈ L2 (R ). The Hilbert transform is well defined for this type of functions and defines a function in the same class [10, p. 311]. We start by considering a piecewise cubic polynomial p(t) having knots ki with i = 0, . . . , n + 1 defined as follows
⎧ 3 2 ⎪ ⎨a0 (t − k0 )3 + b0 (t − k0 )2 + c0 (t − k0 ) + d0 k0 ≤ t < k1 a1 (t − k1 ) + b1 (t − k1 ) + c1 (t − k1 ) + d1 k1 ≤ t < k2 p(t ) = , ⎪ ⎩. . . an (t − kn )3 + bn (t − kn )2 + cn (t − kn ) + dn kn ≤ t < kn+1
(2)
where the set of knots ki ∈ R is bounded and strictly increasing (−∞ < k0 < . . . < kn+1 < ∞) and ai , bi , ci , and di are the coefficients of the polynomial defined in the interval X = [ki , ki+1 ). We want to calculate
P (x ) = H { p(t )}.
(3)
To do this we start by defining pi (t) to be the ith piece polynomial of p(t) such that
pi (t ) = ai (t − ki )3 + bi (t − ki )2 + ci (t − ki ) + di ,
(4)
and the shifted rectangular pulse function srect(t) to be equal to
srect(t ) =
1 0≤t <1 0 Otherwise
(5)
Now, we can rewrite p(t) in terms of srect as
p(t ) =
n
pi (t − ki ) srect
i=0
t − ki wi
(6)
where the width of the ith piece wi is defined as wi = ki+1 − ki . By solving the Hilbert transform integration, we find that
H {(at 3 + bt 2 + ct + d ) srect(t/w )} =
x (ax3 + bx2 + cx + d ) log π x−w 2
3 1
w −w x +a − wx2 − 2 3
w2 + b −wx − 2
− cw.
(7)
Using the time shift property of Hilbert transform, which states if H { f (t )} = F (x ), then H { f (t + a )} = F (x + a ), and by substituting in Eq. (7), we get
H { p( x ) } =
n 1
π
Qi ( x − ki )
0 < x ≤ wi
(8)
i=1
where Qi (x) is the Hilbert transform of the ith polynomial which is given by
Qi (x ) = (ai x3 + bi x2 + ci x + di ) log
x x − w i
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
+ ai
−w2i x wi − wi x2 − 2 3
3
+ bi −wi x −
wi 2
2
− ci w i .
3
(9)
Looking at Eq. (9), we can see that there is a singularity at x = 0 and x = wi . But, since the PCP is continuous, i.e. pi (wi ) = pi+1 (0 ) which is equivalent to ai w3i + bi w2i + ci wi + di = di+1 , the limit at the knots in Eq. (8) takes finite values. 3. Problem formulation and motivation We start by introducing the problem and some conventions and assumptions we use throughout this work. Then we discuss the drawback of using the analytic formula given by Eq. (8) to evaluate the numerical Hilbert transform. 3.1. Problem formulation In this work, we assume the signal for which we want to calculate the Hilbert transform is provided in the form of a PCP as in Eq. (2). Many methods to fit or approximate a set of sampled points using PCPs exist in the literature [5,11–14] and are out of the scope of this work. In this work, we assume that we have a piecewise cubic polynomial p(t) consisting of n pieces, having knots ki , where i = 0, . . . , n + 1. Knots ki are assumed to be integers for all i (ki ∈ Z+ ). We want to evaluate p(t) at m + 1 points t = 0, 1, . . . , m within the range of definition of p(t), i.e, k0 = 0 and kn+1 = m. The assumption of integer knots and evaluation points satisfies the requirements for calculating the Hilbert transform of IMFs in HHT. Yet, this work is still valid as long as all the knots and evaluation points can be expressed as integer multiples of any real constant δ ∈ R. 3.2. Motivation While, we can evaluate the numerical Hilbert transform of a cubic spline by directly evaluating Eq. (9), this incurs a high computational complexity. By counting the number of operations needed for the direct evaluation, we get that this requires
(10 add + 9 mul + 1 log +1 div ) × m × n. From this, we see that the order of computationally expensive operations like the logarithm and the division is m × n. As we will discuss later, a big portion of these calculations is redundant and can be avoided by reusing some of the values. We start by discussing the special case where the knots are uniformly spaced, then we discuss the more general case. 4. Uniform knots Let us assume that the breaks of p(t), ki where i = 0, . . . , n + 1, are equidistant, i.e., wi = w for all i. This assumption enables us to rewrite Eq. (9) as follows
x −w2 x w3 2 + − wx − x − w 2 3
x 2 w + bi x2 log + −wx − 2 x−w x
x
+ ci log + w + di log ,
Qi (x ) = ai x3 log
x−w
x−w
(10)
from which we can see that when evaluating the piece given by Qi (x − ki ), we will be substituting in Qi (y) for y = −ki , . . . , m − ki . From Eq. (10),the dependence of Qi (x) on the piece index i is due to the coefficients ai , bi , ci and di only. Hence, we can rewrite Qi (x) as follows
Qi (x ) = ai (F1 (x )) + bi (F2 (x ))) + ci (F3 (x )) + di (F4 (x )) where
x −w2 x w3 2 F1 (x ) = x log − wx − , + x−w 2 3 3
x w2 F2 (x ) = x log + −wx − 2 , x−w 2
F3 (x ) = log
x x − w + w,
(11)
(12)
(13) (14)
4
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
and
F4 (x ) = log
x . x−w
(15)
As we can see, F1 (x), F2 (x), F3 (x), and F4 (x) do not depend on the piece index i and depend only on the value of x. This enables us to precompute all these functions for x = −kn , . . . , m. Then, when evaluating Qi (x), we will read the values of these functions for y = −ki , . . . , m − ki and multiply them by the appropriate values of ai , bi , ci , and di . Using this property, instead of substituting in Eq. (9) m × n times, we only need to precompute F1 (x), F2 (x), F3 (x), and F4 (x) using 2m points and then substitute in Eq. (11). The number of operations needed in such case is as follows: For precomputing, we need
(7 add + 6 mul + 1 log +1 div ) × 2 × m and for calculating P(x)
(3 add + 4 mul ) × m × n. Hence, we have reduced the order of the log and div operations and most of the mul and add from m × n to m. Thus, reducing the computational complexity of evaluating the Hilbert transform of a cubic polynomial. 5. General knots Now, we discuss the general case, where the difference between successive breaks is not necessarily constant for all pieces. We can still use the same technique we used in the previous section, if we divide the n pieces into sets, where the pieces in each set have the same width. For a set of pieces having the same width whose indexes are i1 , i2 , . . . , in , such that wi1 = wi2 = . . . = win where i1 < i2 < . . . < in , we need to evaluate the four function, F1 (x), F2 (x), F3 (x), and F4 (x), for y = −kin , . . . , m − ki1 . This way, we would be avoiding the recalculation of the four functions for m − ki1 + kin points. By precomputing the four functions and reusing the values, we are guaranteed to reduce the number of points needed for calculations or at least retain the same number of operations as the worst case scenario. To show that, we consider the worst case scenario where we only have two pieces having the same width and they are the first and the last piece in the PCP. In that case, we would have ki1 = 0 and kin = m. In that case, we will need to calculate 2m points for these two pieces. This is the same number of points we would have calculated if we had used the direct substitution method. So, in that case, despite we did not gain any benefits from our proposed method, we did not perform any extra calculations compared to the direct substitution method. In all other cases, the overlap between the pieces will reduce the number of points that need to be calculated to a number less than 2m. The more pieces having the same width exist in our PCP, the more reduction in calculations will occur. Now, we try to calculate the number of operations needed by this technique. Using this method, the number of operations depends on the number of distinct values of w, we will refer to as nw , and the average of the number of coefficients that needs to be computed for a length wi , which we will refer to as Yav given by
n w
Yav =
i=1
m − kimin + kimax nw
(16)
where imin and imin are the index of the first and last breaks having a width of wi , and we have m ≤ Yav ≤ 2m. The number of operation to precompute the four functions using this method is
(7 add + 6 mul + 1 log +1 div ) × Yav × nw and to calculate P(x) it would require an additional
(3 add + 4 mul ) × m × n. As we discussed before this method decreases the required calculations since Yav × nw ≤ n × m. Hence, a reduction in operations is expected. In the special case where all breaks are equidistant, we would have Yav = 2 × m and nw = 1, which is the case discussed in Section 4. 6. Inserting additional knots The previously discussed method relies on having multiple pieces with the same width. When we have a piece with a unique width, that is, no other piece in PCP has the same width, our proposed method works similar to direct substitution. One way to avoid unique width pieces is to divide them into pieces with already existing widths. This will decrease the number of distinct widths of the pieces, nw . The method for breaking a piece into two is discussed in Appendix A. Assuming we have a single piece of width wl , for this single width, we will have to precompute m coefficients for all of the four functions and then perform m evaluations of Eq. (11). The cost of this piece would be
(7 add + 6 mul + 1log + 1 div ) × m + (3 add + 4 mul ) × m.
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
5
Fig. 1. Plotting of Eq. (18) for n = 10 0 0, i = 250, and n = 10 0 0. This equation shows the relation between the probability of having |ki − kl | > m/2 and the probability of having probability equal w (p)
If we are capable of rewriting this piece into two pieces of width wi and wj , such that there are at least one piece of width wi and one piece of width wj already existing in p(t), we will avoid the calculations for the unique piece. The effect of splitting this piece will cause an increase in the precomputation of the ith piece by |kl − ki | and a similar increase for the precomputation of the jth piece by |kl − k j |. The two new pieces will require (3 add + 4 mul) × m × 2 for P(x). Hence, this division would have saved us
m−
|kl − k j | + |kl − ki |
× (7 add + 6 mul + 1 log +1 div )
and replaced them with an additional (3 add + 4 mul) × m. This method saves calculations if (|kl − k j | + |kl − ki | ) is substantially less than m, else if it is equal or bigger than m this method will incur additional computations. Now, we evaluate the probability distribution of having (|kl − k j | + |kl − ki | > m ). Assuming the width of the pieces of the piecewise polynomial are independent and p is the probability of element i to have width wi , the probability that the distance between two pieces having the same width wi equals r is
prob(|kl − ki | = r such that
1
wi = wl ) =
c 1 c
( p(1 − p)2r−1 + p2 (1 − p)2r−2 )
1≤r≤i
p( 1 − p )
i
i+r−1
(17)
where c is a normalization constant. This relation was calculated on the assumption that i is within the first half of the PCP, 1 ≤ i ≤ m/2. The case where i is in the second half can similarly be developed. The probability that the distance between two pieces having the same width wi bigger than r is given by
prob(|kl − ki | ≥ r
such that
wi = wl ) =
1 p(1 − p)i−1 c
( 1 − p )r − ( 1 − p )n p
+ ( 1 − p )n
(18)
Plotting this equation for r = m/2 in Fig. 1, we can see that for n = 10 0 0 and a probability of occurrence of width w, p, equals 10−3 , the probability of having |ki − kl | > m/2 equals 0.2. For the w with a bigger value of p, the probability of having |ki − kl | > m/2 is close to zero. Hence, the probability of having |kl − ki | ≥ m/2 for pieces with wi = wl is low for most likely to occur widths. So, the probability of finding two rare widths, wi and wj , such that their sum equals wl is not likely to happen. Even in that case, the probability that the distance between our break kl and then new breaks ki and ki exceeds m is less likely to happen. Hence, the condition that makes modifying the breaks not feasible is not likely to occur.
7. Algorithm In this section, we combine the previously discussed techniques and present an algorithm to efficiently calculate the Hilbert transform of a cubic spline. In the proposed method, we start by dividing the pieces having a unique width into the sum of two widths already existing in the polynomial as demonstrated in Appendix A. Afterwards, for each width, we precompute the functions F1 (x), F2 (x), F3 (x) and F4 (x) as we explained in Section 5. Then, we use these values to calculate the contribution of each piece to the Hilbert transform of the piecewise cubic polynomial. The algorithm is shown in Table 1. From lines 1 to 14, we divide the pieces having a unique value of w. From lines 15 to 24, we calculate the contribution of each piece to the Hilbert transform.
6
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983 Table 1 Hilbert transform of PCP.
Table 2 Analytical Hilbert transform. Function
Hilbert transform
1 (t 2 +1 )
x) (x2 +1 )
1 (t 4 +1 )
2 √x(x +1 ) 2(x4 +1 )
sin (t ) (t 2 +1 )
e−1 −cos(x ) x2 +1
Fig. 2. Hilbert transform of f (t ) = (t 21+1) .
8. Validation To test the validity of our approach, we use our proposed method to evaluate the Hilbert transform of some functions with a known Hilbert transform. These functions are shown in Table 2. We also compare our method with the Hilbert transform calculated using FFT which was introduced in [15]. To facilitate further comparisons, we selected some of the example functions to be the same as the ones used in [16–18]. To use our approach, we approximated the given set of points of each of the functions using a spline. The method we used for this approximation is the one introduced in [5]. Then, we calculated the Hilbert transform of the piece-wise cubic spline according to our method. Figs. 2–4 show that our method is capable of correctly calculating the Hilbert transform. This is demonstrated by the closeness of the obtained Hilbert transform (HT) to the analytic one. An evaluation of the
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
7
Fig. 3. Hilbert transform of f (t ) = (t 41+1) .
(t ) . Fig. 4. Hilbert transform of f (t ) = (sin t 2 +1 )
Table 3 % RRMSE of FFT Hilbert transform and proposed method Function
FFT HT
HT of cubic splines
1
19.51%
1.10%
1
13.52%
0.11%
3.49%
1.89%
(t 2 +1 ) (t 4 +1 ) sin(t ) (t 2 +1 )
Relative Root Mean Square Error (RRMSE) is shown in Table 3. From this table, we can conclude that our method is more accurate than the FFT based method for the given functions. To further validate our work, we compared it with the method in [15] using more realistic signals, which don’t have a known closed form of the Hilbert transform. In Fig. 5, we show the Hilbert transform of a sine wave added to Gaussian noise with a mean of 0.2 using our method and FFT. For comparison, we also plotted the known Hilbert transform of a sine wave without the added noise. In Fig. 6, we compared the FFT and the proposed Hilbert transform algorithms on a speech signal. These results further validate our proposed algorithm and show its applicability to realistic signals. 9. Performance evaluation Now, we evaluate the gain obtained from using our proposed method. To do that, we use points obtained from applying the CAEMD method on 64 realizations of white Gaussian noise of length 214 . We use the piecewise cubic polynomial of the 10th Intrinsic Mode Function (IMF) to test our method. We compare the execution time of the direct method for calculating the Hilbert transform and compare it to the optimized method. To ensure that the accuracy of our measurements, for each IMF, we timed the execution of each method eight times and then selected the median. Fig. 7 shows the average and the standard deviation of the execution time of both methods. We can see that our optimized method is more than four times faster than the original method on the average. These calculations were carried on a computer using an Intel Core i7-3612QM CPU and equipped with 8GB of RAM. This CPU has a 2.10 GHz clock with 6MB of cache. The evaluated implementation of the algorithm is single threaded although the algorithm can be parallelized by performing the calculations for each set of values of w in a separate thread. To understand the reason behind this improvement, we take a closer look at the distribution of piece widths of one of the PCPs. This distribution is shown in Fig. 8. From this figure, we can see that the majority of pieces share the same width
8
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
Fig. 5. Hilbert transform of sin(t) added to white Gaussian noise of mean 0.2.
Fig. 6. Hilbert transform of a speech signal.
Fig. 7. Comparison of direct substitution method and the proposed method in terms of execution time. The X presents the average and the length of the horizontal bar represents two standard deviations.
Fig. 8. Histogram of the distribution of w for one of the IMFs.
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
9
except for a small number of pieces which have a unique value of widths. Our method is capable of reusing the majority of the calculations for pieces having the same width and it avoids repeating expensive operations like the logarithm and division. For the few unique widths, our method is capable of rewriting them as the sum of two existing widths to avoid making calculations for a unique piece. Thus, our method is capable of reformulating the calculations in a manner that significantly reduces execution time. 10. Related work Several methods exist for the evaluation of the Hilbert transform. Existing work in the literature can be divided into three types. The first type uses the relation between the Hilbert transform and other transforms for which there exists efficient numerical algorithms. In [19, p.203], Henrici suggested a method to compute the Hilbert transform based on its relation with the Fourier transform. This method is widely used in mathematical software packages according to the equations suggested in [15]. The second type works by approximating the Hilbert transform integration. Several methods to approximate Cauchy integrals exist in the literature, some of them work over bounded intervals [20–27], while other work over unbounded intervals [28,29]. All these methods to approximate the Cauchy integral can be used for the calculation of Hilbert transform. The third type works by approximating any given function or fitting a set of given samples using functions with a known Hilbert transform. The work in [30] falls within this category, where a method based on the expansion in rational eigenfunctions is used to implement Hilbert transform on the real line. Other methods calculate the Hilbert transform of discrete samples by fitting these samples to a function with known Hilbert transform. From this category, Zhou et al. [16] presented a method that approximates given samples using the Haar transform. Then the Hilbert transform is calculated using using the analytical transform of the Haar basis. Chen et al. presented a method to calculate the Hilbert transform based on the B-spline approximation of the samples as part of an alternative method for EMD in [31]. Continuing on the work in [31], Micchelli et al. provided a fast algorithm for the evaluation of the Hilbert transform of the B-spline approximation of the sample data [17]. In [18], an algorithm to calculate the Hilbert transform of a real-valued function by interpolating using a piecewise-linear function is proposed. Methods to solve special cases of Hilbert transform were also developed. In [32], a numerical procedure to calculate the Hilbert transform for oscillatory functions using a convergence accelerator technique was developed. The Hilbert transform defined over the interval [−1,1] is discussed in [33]. For signal processing applications, methods for the design of Finite Impulse Response (FIR) filters to calculate the Hilbert transform were proposed in [34,35]. Our proposed method falls within the third category. Any set of points can be approximated using a PCP. By using the analytic Hilbert transform of one of the pieces of the PCP, we obtain the Hilbert transform of any given PCP. From this relation, we derived an efficient algorithm to evaluate the Hilbert transform using this analytic expression. As for our assumptions of integer knots and evaluation points. These assumptions are less restrictive than the ones made in [18], where the authors assume that the function is approximated on a uniform grid, while we are assuming that the knots are integer multiples of a real number, but not necessarily uniform. In [17], the authors do not impose any conditions on the evaluation points. However, in [17], the grid of the B-splines knots has to be different from the sample points to avoid singularities in the evaluation. Our proposed method does not have this problem and can have the same knots and the evaluation points. 11. Conclusion In this paper, we introduced a method to calculate the Hilbert transform of a cubic spline. This method optimizes the calculation of the analytical expression for the Hilbert transform of splines. It starts by breaking the pieces having unique widths. Then, it creates lookup tables for calculating the Hilbert transform of pieces having the same width to avoid recalculations. This method is valid as long as the points to be evaluated and the knots of the spline are multiples of some real value δ . Evaluations have shown that the proposed technique gives accurate results with a significant performance gain compared to direct substitution. The proposed work is applicable in the calculation of the Hilbert–Huang transform. Acknowledgments The authors would like to thank the reviewers for their valuable feedback, which helped improve the quality of this work. Appendix A. Dividing a polynomial to pieces The polynomial
a p,1 (x − x p )3 + b p,1 (x − x p )2 + c p,1 (x − x p ) + d p,1
x p ≤ x < x2
(A.1)
10
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
is to be rewritten into two branches separated at q as
aip,1 (x − x p )3 + bip,1 (x − x p )2 + cip,1 (x − x p ) + dip,1 aip,1 (x − x p )3 + bip,1 (x − x p )2 + cip,1 (x − x p ) + dip,1
x p ≤ x < xq xq ≤ x < x2
(A.2)
The first branch would have the same coefficients. The coefficients of the second branch can be obtained by substituting for (x − p) by ((x − q ) + q − p) in the original polynomial and expanding in terms of (x − q ). The values of aip , bip , cip and dip can easily be deduced using some algebraic simplifications as follows
aip = a p bip = 3a p r + b p cip = 3a p r 2 + 2b p r + c p dip = a p r 3 + b p r 2 + c p r + d p
(A.3)
where r = q − p. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]
Hahn SL. Hilbert transforms in signal processing, 2. Artech House Boston; 1996. Pandey JN. The Hilbert transform of Schwartz distributions and applications, 27. John Wiley & Sons; 2011. King FW. Hilbert transforms, 2. Cambridge University Press Cambridge; 2009. Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc R Soc A 1998;454(1971):903–95. doi:10.1098/rspa.1998.0193. Abd-el-Malek MB, Hanna SS. Using filter bank property to simplify the calculations of Emprical Mode Decompostion. Commun. Nonlinear Sci. Numer. Simul. 2018;62:429–44. De Boor C, Rice JR. Least squares cubic spline approximation, II-variable knots; 1968. Teichert GH, Gunda NH, Rudraraju S, Natarajan AR, Puchala B, Garikipati K, et al. A comparison of Redlich-Kister polynomial and cubic spline representations of the chemical potential in phase field computations. Comput Mater Sci 2017;128:127–39. Bruno Siciliano OK. Handbook of robotics. Springer; 2016. Mortenson ME. Mathematics for computer graphics applications. Industrial Press Inc.; 1999. ISBN 978-0-8311-3111-1. Google-Books-ID: YmQy799flPkC Butzer PL, Nessel RJ. Fourier analysis and approximation. Basel: Birkhuser Basel; 1971. doi:10.1007/978- 3- 0348- 7448- 9. ISBN 978-3-0348-7450-2 9783-0348-7448-9 De Boor C, De Boor C, Mathmaticien E-U, De Boor C, De Boor C. A practical guide to splines, 27. Springer-Verlag New York; 1978. Akima H. A new method of interpolation and smooth curve fitting based on local procedures. J ACM 1970;17(4):589–602. doi:10.1145/321607.321609. Prenter PM, Westwater ER. Three adaptive discrete least squares cubic spline procedures for the compression of data. Comput Vis Graph Image Process 1986;33(3):327–45. doi:10.1016/0734-189X(86)90181-7. Can E. Piecewise cubic approximation for data. Am J Appl Math 2013;1(2):24. doi:10.11648/j.ajam.20130102.11. Marple SL Jr. Computing the discrete-time analytic signal via FFT. Signal Process IEEE Trans 1999;47(9):2600–3. Zhou C, Yang L, Liu Y, Yang Z. A novel method for computing the Hilbert transform with Haar multiresolution approximation. J Comput Appl Math 2009;223(2):585–97. doi:10.1016/j.cam.2008.02.006. Micchelli CA, Xu Y, Yu B. On computing with the Hilbert spline transform. Adv Comput Math 2013;38(3):623–46. doi:10.1007/s10444- 011- 9252- x. Bilato R, Maj O, Brambilla M. An algorithm for fast Hilbert transform of real functions. Adv Comput Math 2014;40(5-6):1159–68. doi:10.1007/ s10444- 014- 9345- 4. Henrici P. Applied and computational complex analysis. Discrete Fourier analysis, Cauchy integrals, construction of conformal maps, univalent functions, vol. 3. John Wiley & Sons; 1993. Dagnino C, Santi E. Spline product quadrature rules for cauchy singular integrals. J Comput Appl Math 1990;33(2):133–40. doi:10.1016/0377-0427(90) 90363-5. Monegato G. The numerical evaluation of one-dimensional Cauchy principal value integrals. Computing 1982;29(4):337–54. doi:10.1007/BF02246760. Diethelm K. Modified compound quadrature rules for strongly singular integrals. Computing 1994;52(4):337–54. doi:10.1007/BF02276881. Diethelm K. Peano kernels and bounds for the error constants of Gaussian and related quadrature rules for Cauchy principal value integrals. Numer Math 1996;73(1):53–63. doi:10.10 07/s0 02110 050183. Elliott D, Paget DF. Gauss type quadrature rules for Cauchy principal value integrals. Math Comput 1979;33(145):301–9. Diethelm K. Uniform convergence of optimal order quadrature rules for Cauchy principal value integrals. J Comput Appl Math 1994;56(3):321–9. Criscuolo G. A new algorithm for Cauchy principal value and Hadamard finite-part integrals. J Comput Appl Math 1997;78(2):255–75. Orsi AP. Spline approximation for Cauchy principal value integrals. J Comput Appl Math 1990;30(2):191–201. Damelin S, Diethelm K. Interpolatory product quadratures for Cauchy principal value integrals with Freud weights. Numer Math 1999;83(1):87–105. doi:10.10 07/s0 02110 050440. Damelin SB, Diethelm K. Boundedness and uniform numerical approximation of the weighted Hilbert transform on the real line. Numer Funct Anal Optim 2001;22(1-2):13–54. doi:10.1081/NFA-100103786. Weideman JAC. Computing the Hilbert transform on the real line. Math Comput 1995;64(210):745–62. Chen Q, Huang N, Riemenschneider S, Xu Y. A B-spline approach for empirical mode decompositions. Adv Comput Math 2006;24(1-4):171–95. doi:10. 1007/s10444- 004- 7614- 3. King FW, Smethells GJ, Helleloid GT, Pelzl PJ. Numerical evaluation of Hilbert transforms for oscillatory functions: a convergence accelerator approach. Comput Phys Commun 2002;145(2):256–66. doi:10.1016/S0010-4655(02)00155-8. Hasegawa T. Uniform approximations to finite Hilbert transform and its derivative. J Comput Appl Math 2004;163(1):127–38. Schussler H, Weith J. On the design of recursive Hilbert-transformers. In: Acoustics, speech, and signal processing, IEEE international conference on ICASSP’87., 12. IEEE; 1987. p. 876–9. Rabiner LR, Schafer RW. On the behavior of minimax FIR digital Hilbert transformers. Bell Labs Tech J 1974;53(2):363–90.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Research paper
The Hilbert transform of cubic splines Mina B. Abd-el-Malek∗, Samer S. Hanna Department of Engineering Mathematics and Physics, Faculty of Engineering, Alexandria University, Alexandria 21544, Egypt
a r t i c l e
i n f o
Article history: Received 12 March 2018 Revised 24 December 2018 Accepted 26 August 2019 Available online 30 August 2019 Keywords: Hilbert transform Splines Piecewise cubic polynomial Hilbert–Huang transform
a b s t r a c t In this work, we propose an efficient method to calculate the Hilbert transform of cubic splines. We start by developing an analytical derivation for the problem. Using the structure of the obtained formula, we suggest an efficient algorithm to evaluate it. Our method rewrites the cubic spline and reorders calculations to increase the temporal locality of related tasks. Then, it creates lookup tables to reuse calculations instead of repeating them. Results show that this method accurately calculates Hilbert transform when applied to functions with a known analytical Hilbert transform. We also show that it significantly reduces the execution time compared to the direct substitution. Our suggested algorithm can be used for the Hilbert–Huang Transform, whose calculation is based on cubic splines and Hilbert transform. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The Hilbert transform (HT) has countless applications in the fields of signal processing, applied mathematics, and physical sciences. It has been used to solve several problems in these fields [1–3]. The significance of some of these problems has led to an interest in efficient methods for the evaluation of the Hilbert transform. Most recently, a transform based on Hilbert transform and cubic splines was introduced. The Hilbert–Huang Transform (HHT) is a time-frequency transform suitable for non-linear and non-stationary signals [4]. HHT uses a method called Empirical Mode Decomposition (EMD) to decompose a given time series into signals with magnitudes and frequencies of close scales called Intrinsic Mode Functions (IMFs). Then, it applies the Hilbert transform to these IMFs, in order to derive a timefrequency representation of the signal. The steps of EMD relies on the use of cubic splines. The Computationally Adaptive Empirical Mode Decomposition (CAEMD) [5] was developed as an enhancement to EMD. When using CAEMD, some of the IMFs are obtained in the form of Piecewise Cubic Polynomials (PCPs). In the second step of HHT, we need to evaluate the Hilbert transform of these PCPs. In this work, we derive an efficient method to calculate the Hilbert transform of a given piecewise cubic polynomial. The Hilbert transform of a PCP is of interest for methods like the Hilbert–Huang transform. To obtain the Hilbert transform of a PCP IMF using any of the existing methods in the literature, the PCPs needs to be evaluated at all points before applying the Hilbert transform. This due to the fact that none of these methods has addressed the problem of obtaining the Hilbert transform of a function given in the form of a PCP. Other than the HHT, piecewise cubic polynomials are useful to create a mathematical model for data [6]. They have applications in many fields like material science [7], robotics [8], and computer graphics [9]. This further motivates the development of a method to calculate the Hilbert transform of piecewise cubic polynomials. ∗
Corresponding author. E-mail addresses: [email protected] (M.B. Abd-el-Malek), [email protected] (S.S. Hanna).
https://doi.org/10.1016/j.cnsns.2019.104983 1007-5704/© 2019 Elsevier B.V. All rights reserved.
2
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
In this work, we present an algorithm that calculates the Hilbert transform of a piecewise cubic spline directly. This algorithm uses the properties of the Hilbert transform and the nature of cubic splines to simplify the calculations. It uses the known Hilbert transform of a cubic polynomial to efficiently calculates the HT of a PCP. Although this method was originally developed to calculate the HT of PCPs obtained from CAEMD, it can be used to calculate the Hilbert transform of any set of points that can be fitted or approximated using a PCP. The rest of this paper is organized as follows. In Section 2, we derive an analytical formula for the Hilbert transform of a PCP. The problem, the obtained formula, and the assumptions of this work are discussed in Section 3. Section 4 starts analyzes the case of a PCP with uniform knots, while Section 5 addresses a more general case. In Section 6, we propose reformulation the PCP to further improve computational efficiency. The proposed algorithm is then presented in Section 7. In Section 8, the suggested method is validated against functions with known Hilbert transform and applied to realistic signals. The execution times of the direct and proposed methods are compared in Section 9. A survey of existing methods for the evaluation of the Hilbert transform are presented in Section 10, while Section 11 concludes this work. 2. Analytical derivation The Hilbert transform of a function f (t ) ∈ L p (R ) (1 ≤ p < ∞) is defined by the following principle value integral
F (x ) = H { f (t )} =
1
π
∞
PV −∞
f (t ) dt x−t
(1)
In this work, we only consider the case where f is a real-valued function belonging to the set of square integrable functions defined on the real line i.e, f ∈ L2 (R ). The Hilbert transform is well defined for this type of functions and defines a function in the same class [10, p. 311]. We start by considering a piecewise cubic polynomial p(t) having knots ki with i = 0, . . . , n + 1 defined as follows
⎧ 3 2 ⎪ ⎨a0 (t − k0 )3 + b0 (t − k0 )2 + c0 (t − k0 ) + d0 k0 ≤ t < k1 a1 (t − k1 ) + b1 (t − k1 ) + c1 (t − k1 ) + d1 k1 ≤ t < k2 p(t ) = , ⎪ ⎩. . . an (t − kn )3 + bn (t − kn )2 + cn (t − kn ) + dn kn ≤ t < kn+1
(2)
where the set of knots ki ∈ R is bounded and strictly increasing (−∞ < k0 < . . . < kn+1 < ∞) and ai , bi , ci , and di are the coefficients of the polynomial defined in the interval X = [ki , ki+1 ). We want to calculate
P (x ) = H { p(t )}.
(3)
To do this we start by defining pi (t) to be the ith piece polynomial of p(t) such that
pi (t ) = ai (t − ki )3 + bi (t − ki )2 + ci (t − ki ) + di ,
(4)
and the shifted rectangular pulse function srect(t) to be equal to
srect(t ) =
1 0≤t <1 0 Otherwise
(5)
Now, we can rewrite p(t) in terms of srect as
p(t ) =
n
pi (t − ki ) srect
i=0
t − ki wi
(6)
where the width of the ith piece wi is defined as wi = ki+1 − ki . By solving the Hilbert transform integration, we find that
H {(at 3 + bt 2 + ct + d ) srect(t/w )} =
x (ax3 + bx2 + cx + d ) log π x−w 2
3 1
w −w x +a − wx2 − 2 3
w2 + b −wx − 2
− cw.
(7)
Using the time shift property of Hilbert transform, which states if H { f (t )} = F (x ), then H { f (t + a )} = F (x + a ), and by substituting in Eq. (7), we get
H { p( x ) } =
n 1
π
Qi ( x − ki )
0 < x ≤ wi
(8)
i=1
where Qi (x) is the Hilbert transform of the ith polynomial which is given by
Qi (x ) = (ai x3 + bi x2 + ci x + di ) log
x x − w i
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
+ ai
−w2i x wi − wi x2 − 2 3
3
+ bi −wi x −
wi 2
2
− ci w i .
3
(9)
Looking at Eq. (9), we can see that there is a singularity at x = 0 and x = wi . But, since the PCP is continuous, i.e. pi (wi ) = pi+1 (0 ) which is equivalent to ai w3i + bi w2i + ci wi + di = di+1 , the limit at the knots in Eq. (8) takes finite values. 3. Problem formulation and motivation We start by introducing the problem and some conventions and assumptions we use throughout this work. Then we discuss the drawback of using the analytic formula given by Eq. (8) to evaluate the numerical Hilbert transform. 3.1. Problem formulation In this work, we assume the signal for which we want to calculate the Hilbert transform is provided in the form of a PCP as in Eq. (2). Many methods to fit or approximate a set of sampled points using PCPs exist in the literature [5,11–14] and are out of the scope of this work. In this work, we assume that we have a piecewise cubic polynomial p(t) consisting of n pieces, having knots ki , where i = 0, . . . , n + 1. Knots ki are assumed to be integers for all i (ki ∈ Z+ ). We want to evaluate p(t) at m + 1 points t = 0, 1, . . . , m within the range of definition of p(t), i.e, k0 = 0 and kn+1 = m. The assumption of integer knots and evaluation points satisfies the requirements for calculating the Hilbert transform of IMFs in HHT. Yet, this work is still valid as long as all the knots and evaluation points can be expressed as integer multiples of any real constant δ ∈ R. 3.2. Motivation While, we can evaluate the numerical Hilbert transform of a cubic spline by directly evaluating Eq. (9), this incurs a high computational complexity. By counting the number of operations needed for the direct evaluation, we get that this requires
(10 add + 9 mul + 1 log +1 div ) × m × n. From this, we see that the order of computationally expensive operations like the logarithm and the division is m × n. As we will discuss later, a big portion of these calculations is redundant and can be avoided by reusing some of the values. We start by discussing the special case where the knots are uniformly spaced, then we discuss the more general case. 4. Uniform knots Let us assume that the breaks of p(t), ki where i = 0, . . . , n + 1, are equidistant, i.e., wi = w for all i. This assumption enables us to rewrite Eq. (9) as follows
x −w2 x w3 2 + − wx − x − w 2 3
x 2 w + bi x2 log + −wx − 2 x−w x
x
+ ci log + w + di log ,
Qi (x ) = ai x3 log
x−w
x−w
(10)
from which we can see that when evaluating the piece given by Qi (x − ki ), we will be substituting in Qi (y) for y = −ki , . . . , m − ki . From Eq. (10),the dependence of Qi (x) on the piece index i is due to the coefficients ai , bi , ci and di only. Hence, we can rewrite Qi (x) as follows
Qi (x ) = ai (F1 (x )) + bi (F2 (x ))) + ci (F3 (x )) + di (F4 (x )) where
x −w2 x w3 2 F1 (x ) = x log − wx − , + x−w 2 3 3
x w2 F2 (x ) = x log + −wx − 2 , x−w 2
F3 (x ) = log
x x − w + w,
(11)
(12)
(13) (14)
4
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
and
F4 (x ) = log
x . x−w
(15)
As we can see, F1 (x), F2 (x), F3 (x), and F4 (x) do not depend on the piece index i and depend only on the value of x. This enables us to precompute all these functions for x = −kn , . . . , m. Then, when evaluating Qi (x), we will read the values of these functions for y = −ki , . . . , m − ki and multiply them by the appropriate values of ai , bi , ci , and di . Using this property, instead of substituting in Eq. (9) m × n times, we only need to precompute F1 (x), F2 (x), F3 (x), and F4 (x) using 2m points and then substitute in Eq. (11). The number of operations needed in such case is as follows: For precomputing, we need
(7 add + 6 mul + 1 log +1 div ) × 2 × m and for calculating P(x)
(3 add + 4 mul ) × m × n. Hence, we have reduced the order of the log and div operations and most of the mul and add from m × n to m. Thus, reducing the computational complexity of evaluating the Hilbert transform of a cubic polynomial. 5. General knots Now, we discuss the general case, where the difference between successive breaks is not necessarily constant for all pieces. We can still use the same technique we used in the previous section, if we divide the n pieces into sets, where the pieces in each set have the same width. For a set of pieces having the same width whose indexes are i1 , i2 , . . . , in , such that wi1 = wi2 = . . . = win where i1 < i2 < . . . < in , we need to evaluate the four function, F1 (x), F2 (x), F3 (x), and F4 (x), for y = −kin , . . . , m − ki1 . This way, we would be avoiding the recalculation of the four functions for m − ki1 + kin points. By precomputing the four functions and reusing the values, we are guaranteed to reduce the number of points needed for calculations or at least retain the same number of operations as the worst case scenario. To show that, we consider the worst case scenario where we only have two pieces having the same width and they are the first and the last piece in the PCP. In that case, we would have ki1 = 0 and kin = m. In that case, we will need to calculate 2m points for these two pieces. This is the same number of points we would have calculated if we had used the direct substitution method. So, in that case, despite we did not gain any benefits from our proposed method, we did not perform any extra calculations compared to the direct substitution method. In all other cases, the overlap between the pieces will reduce the number of points that need to be calculated to a number less than 2m. The more pieces having the same width exist in our PCP, the more reduction in calculations will occur. Now, we try to calculate the number of operations needed by this technique. Using this method, the number of operations depends on the number of distinct values of w, we will refer to as nw , and the average of the number of coefficients that needs to be computed for a length wi , which we will refer to as Yav given by
n w
Yav =
i=1
m − kimin + kimax nw
(16)
where imin and imin are the index of the first and last breaks having a width of wi , and we have m ≤ Yav ≤ 2m. The number of operation to precompute the four functions using this method is
(7 add + 6 mul + 1 log +1 div ) × Yav × nw and to calculate P(x) it would require an additional
(3 add + 4 mul ) × m × n. As we discussed before this method decreases the required calculations since Yav × nw ≤ n × m. Hence, a reduction in operations is expected. In the special case where all breaks are equidistant, we would have Yav = 2 × m and nw = 1, which is the case discussed in Section 4. 6. Inserting additional knots The previously discussed method relies on having multiple pieces with the same width. When we have a piece with a unique width, that is, no other piece in PCP has the same width, our proposed method works similar to direct substitution. One way to avoid unique width pieces is to divide them into pieces with already existing widths. This will decrease the number of distinct widths of the pieces, nw . The method for breaking a piece into two is discussed in Appendix A. Assuming we have a single piece of width wl , for this single width, we will have to precompute m coefficients for all of the four functions and then perform m evaluations of Eq. (11). The cost of this piece would be
(7 add + 6 mul + 1log + 1 div ) × m + (3 add + 4 mul ) × m.
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
5
Fig. 1. Plotting of Eq. (18) for n = 10 0 0, i = 250, and n = 10 0 0. This equation shows the relation between the probability of having |ki − kl | > m/2 and the probability of having probability equal w (p)
If we are capable of rewriting this piece into two pieces of width wi and wj , such that there are at least one piece of width wi and one piece of width wj already existing in p(t), we will avoid the calculations for the unique piece. The effect of splitting this piece will cause an increase in the precomputation of the ith piece by |kl − ki | and a similar increase for the precomputation of the jth piece by |kl − k j |. The two new pieces will require (3 add + 4 mul) × m × 2 for P(x). Hence, this division would have saved us
m−
|kl − k j | + |kl − ki |
× (7 add + 6 mul + 1 log +1 div )
and replaced them with an additional (3 add + 4 mul) × m. This method saves calculations if (|kl − k j | + |kl − ki | ) is substantially less than m, else if it is equal or bigger than m this method will incur additional computations. Now, we evaluate the probability distribution of having (|kl − k j | + |kl − ki | > m ). Assuming the width of the pieces of the piecewise polynomial are independent and p is the probability of element i to have width wi , the probability that the distance between two pieces having the same width wi equals r is
prob(|kl − ki | = r such that
1
wi = wl ) =
c 1 c
( p(1 − p)2r−1 + p2 (1 − p)2r−2 )
1≤r≤i
p( 1 − p )
i
i+r−1
(17)
where c is a normalization constant. This relation was calculated on the assumption that i is within the first half of the PCP, 1 ≤ i ≤ m/2. The case where i is in the second half can similarly be developed. The probability that the distance between two pieces having the same width wi bigger than r is given by
prob(|kl − ki | ≥ r
such that
wi = wl ) =
1 p(1 − p)i−1 c
( 1 − p )r − ( 1 − p )n p
+ ( 1 − p )n
(18)
Plotting this equation for r = m/2 in Fig. 1, we can see that for n = 10 0 0 and a probability of occurrence of width w, p, equals 10−3 , the probability of having |ki − kl | > m/2 equals 0.2. For the w with a bigger value of p, the probability of having |ki − kl | > m/2 is close to zero. Hence, the probability of having |kl − ki | ≥ m/2 for pieces with wi = wl is low for most likely to occur widths. So, the probability of finding two rare widths, wi and wj , such that their sum equals wl is not likely to happen. Even in that case, the probability that the distance between our break kl and then new breaks ki and ki exceeds m is less likely to happen. Hence, the condition that makes modifying the breaks not feasible is not likely to occur.
7. Algorithm In this section, we combine the previously discussed techniques and present an algorithm to efficiently calculate the Hilbert transform of a cubic spline. In the proposed method, we start by dividing the pieces having a unique width into the sum of two widths already existing in the polynomial as demonstrated in Appendix A. Afterwards, for each width, we precompute the functions F1 (x), F2 (x), F3 (x) and F4 (x) as we explained in Section 5. Then, we use these values to calculate the contribution of each piece to the Hilbert transform of the piecewise cubic polynomial. The algorithm is shown in Table 1. From lines 1 to 14, we divide the pieces having a unique value of w. From lines 15 to 24, we calculate the contribution of each piece to the Hilbert transform.
6
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983 Table 1 Hilbert transform of PCP.
Table 2 Analytical Hilbert transform. Function
Hilbert transform
1 (t 2 +1 )
x) (x2 +1 )
1 (t 4 +1 )
2 √x(x +1 ) 2(x4 +1 )
sin (t ) (t 2 +1 )
e−1 −cos(x ) x2 +1
Fig. 2. Hilbert transform of f (t ) = (t 21+1) .
8. Validation To test the validity of our approach, we use our proposed method to evaluate the Hilbert transform of some functions with a known Hilbert transform. These functions are shown in Table 2. We also compare our method with the Hilbert transform calculated using FFT which was introduced in [15]. To facilitate further comparisons, we selected some of the example functions to be the same as the ones used in [16–18]. To use our approach, we approximated the given set of points of each of the functions using a spline. The method we used for this approximation is the one introduced in [5]. Then, we calculated the Hilbert transform of the piece-wise cubic spline according to our method. Figs. 2–4 show that our method is capable of correctly calculating the Hilbert transform. This is demonstrated by the closeness of the obtained Hilbert transform (HT) to the analytic one. An evaluation of the
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
7
Fig. 3. Hilbert transform of f (t ) = (t 41+1) .
(t ) . Fig. 4. Hilbert transform of f (t ) = (sin t 2 +1 )
Table 3 % RRMSE of FFT Hilbert transform and proposed method Function
FFT HT
HT of cubic splines
1
19.51%
1.10%
1
13.52%
0.11%
3.49%
1.89%
(t 2 +1 ) (t 4 +1 ) sin(t ) (t 2 +1 )
Relative Root Mean Square Error (RRMSE) is shown in Table 3. From this table, we can conclude that our method is more accurate than the FFT based method for the given functions. To further validate our work, we compared it with the method in [15] using more realistic signals, which don’t have a known closed form of the Hilbert transform. In Fig. 5, we show the Hilbert transform of a sine wave added to Gaussian noise with a mean of 0.2 using our method and FFT. For comparison, we also plotted the known Hilbert transform of a sine wave without the added noise. In Fig. 6, we compared the FFT and the proposed Hilbert transform algorithms on a speech signal. These results further validate our proposed algorithm and show its applicability to realistic signals. 9. Performance evaluation Now, we evaluate the gain obtained from using our proposed method. To do that, we use points obtained from applying the CAEMD method on 64 realizations of white Gaussian noise of length 214 . We use the piecewise cubic polynomial of the 10th Intrinsic Mode Function (IMF) to test our method. We compare the execution time of the direct method for calculating the Hilbert transform and compare it to the optimized method. To ensure that the accuracy of our measurements, for each IMF, we timed the execution of each method eight times and then selected the median. Fig. 7 shows the average and the standard deviation of the execution time of both methods. We can see that our optimized method is more than four times faster than the original method on the average. These calculations were carried on a computer using an Intel Core i7-3612QM CPU and equipped with 8GB of RAM. This CPU has a 2.10 GHz clock with 6MB of cache. The evaluated implementation of the algorithm is single threaded although the algorithm can be parallelized by performing the calculations for each set of values of w in a separate thread. To understand the reason behind this improvement, we take a closer look at the distribution of piece widths of one of the PCPs. This distribution is shown in Fig. 8. From this figure, we can see that the majority of pieces share the same width
8
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
Fig. 5. Hilbert transform of sin(t) added to white Gaussian noise of mean 0.2.
Fig. 6. Hilbert transform of a speech signal.
Fig. 7. Comparison of direct substitution method and the proposed method in terms of execution time. The X presents the average and the length of the horizontal bar represents two standard deviations.
Fig. 8. Histogram of the distribution of w for one of the IMFs.
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
9
except for a small number of pieces which have a unique value of widths. Our method is capable of reusing the majority of the calculations for pieces having the same width and it avoids repeating expensive operations like the logarithm and division. For the few unique widths, our method is capable of rewriting them as the sum of two existing widths to avoid making calculations for a unique piece. Thus, our method is capable of reformulating the calculations in a manner that significantly reduces execution time. 10. Related work Several methods exist for the evaluation of the Hilbert transform. Existing work in the literature can be divided into three types. The first type uses the relation between the Hilbert transform and other transforms for which there exists efficient numerical algorithms. In [19, p.203], Henrici suggested a method to compute the Hilbert transform based on its relation with the Fourier transform. This method is widely used in mathematical software packages according to the equations suggested in [15]. The second type works by approximating the Hilbert transform integration. Several methods to approximate Cauchy integrals exist in the literature, some of them work over bounded intervals [20–27], while other work over unbounded intervals [28,29]. All these methods to approximate the Cauchy integral can be used for the calculation of Hilbert transform. The third type works by approximating any given function or fitting a set of given samples using functions with a known Hilbert transform. The work in [30] falls within this category, where a method based on the expansion in rational eigenfunctions is used to implement Hilbert transform on the real line. Other methods calculate the Hilbert transform of discrete samples by fitting these samples to a function with known Hilbert transform. From this category, Zhou et al. [16] presented a method that approximates given samples using the Haar transform. Then the Hilbert transform is calculated using using the analytical transform of the Haar basis. Chen et al. presented a method to calculate the Hilbert transform based on the B-spline approximation of the samples as part of an alternative method for EMD in [31]. Continuing on the work in [31], Micchelli et al. provided a fast algorithm for the evaluation of the Hilbert transform of the B-spline approximation of the sample data [17]. In [18], an algorithm to calculate the Hilbert transform of a real-valued function by interpolating using a piecewise-linear function is proposed. Methods to solve special cases of Hilbert transform were also developed. In [32], a numerical procedure to calculate the Hilbert transform for oscillatory functions using a convergence accelerator technique was developed. The Hilbert transform defined over the interval [−1,1] is discussed in [33]. For signal processing applications, methods for the design of Finite Impulse Response (FIR) filters to calculate the Hilbert transform were proposed in [34,35]. Our proposed method falls within the third category. Any set of points can be approximated using a PCP. By using the analytic Hilbert transform of one of the pieces of the PCP, we obtain the Hilbert transform of any given PCP. From this relation, we derived an efficient algorithm to evaluate the Hilbert transform using this analytic expression. As for our assumptions of integer knots and evaluation points. These assumptions are less restrictive than the ones made in [18], where the authors assume that the function is approximated on a uniform grid, while we are assuming that the knots are integer multiples of a real number, but not necessarily uniform. In [17], the authors do not impose any conditions on the evaluation points. However, in [17], the grid of the B-splines knots has to be different from the sample points to avoid singularities in the evaluation. Our proposed method does not have this problem and can have the same knots and the evaluation points. 11. Conclusion In this paper, we introduced a method to calculate the Hilbert transform of a cubic spline. This method optimizes the calculation of the analytical expression for the Hilbert transform of splines. It starts by breaking the pieces having unique widths. Then, it creates lookup tables for calculating the Hilbert transform of pieces having the same width to avoid recalculations. This method is valid as long as the points to be evaluated and the knots of the spline are multiples of some real value δ . Evaluations have shown that the proposed technique gives accurate results with a significant performance gain compared to direct substitution. The proposed work is applicable in the calculation of the Hilbert–Huang transform. Acknowledgments The authors would like to thank the reviewers for their valuable feedback, which helped improve the quality of this work. Appendix A. Dividing a polynomial to pieces The polynomial
a p,1 (x − x p )3 + b p,1 (x − x p )2 + c p,1 (x − x p ) + d p,1
x p ≤ x < x2
(A.1)
10
M.B. Abd-el-Malek and S.S. Hanna / Commun Nonlinear Sci Numer Simulat 80 (2020) 104983
is to be rewritten into two branches separated at q as
aip,1 (x − x p )3 + bip,1 (x − x p )2 + cip,1 (x − x p ) + dip,1 aip,1 (x − x p )3 + bip,1 (x − x p )2 + cip,1 (x − x p ) + dip,1
x p ≤ x < xq xq ≤ x < x2
(A.2)
The first branch would have the same coefficients. The coefficients of the second branch can be obtained by substituting for (x − p) by ((x − q ) + q − p) in the original polynomial and expanding in terms of (x − q ). The values of aip , bip , cip and dip can easily be deduced using some algebraic simplifications as follows
aip = a p bip = 3a p r + b p cip = 3a p r 2 + 2b p r + c p dip = a p r 3 + b p r 2 + c p r + d p
(A.3)
where r = q − p. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35]
Hahn SL. Hilbert transforms in signal processing, 2. Artech House Boston; 1996. Pandey JN. The Hilbert transform of Schwartz distributions and applications, 27. John Wiley & Sons; 2011. King FW. Hilbert transforms, 2. Cambridge University Press Cambridge; 2009. Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc R Soc A 1998;454(1971):903–95. doi:10.1098/rspa.1998.0193. Abd-el-Malek MB, Hanna SS. Using filter bank property to simplify the calculations of Emprical Mode Decompostion. Commun. Nonlinear Sci. Numer. Simul. 2018;62:429–44. De Boor C, Rice JR. Least squares cubic spline approximation, II-variable knots; 1968. Teichert GH, Gunda NH, Rudraraju S, Natarajan AR, Puchala B, Garikipati K, et al. A comparison of Redlich-Kister polynomial and cubic spline representations of the chemical potential in phase field computations. Comput Mater Sci 2017;128:127–39. Bruno Siciliano OK. Handbook of robotics. Springer; 2016. Mortenson ME. Mathematics for computer graphics applications. Industrial Press Inc.; 1999. ISBN 978-0-8311-3111-1. Google-Books-ID: YmQy799flPkC Butzer PL, Nessel RJ. Fourier analysis and approximation. Basel: Birkhuser Basel; 1971. doi:10.1007/978- 3- 0348- 7448- 9. ISBN 978-3-0348-7450-2 9783-0348-7448-9 De Boor C, De Boor C, Mathmaticien E-U, De Boor C, De Boor C. A practical guide to splines, 27. Springer-Verlag New York; 1978. Akima H. A new method of interpolation and smooth curve fitting based on local procedures. J ACM 1970;17(4):589–602. doi:10.1145/321607.321609. Prenter PM, Westwater ER. Three adaptive discrete least squares cubic spline procedures for the compression of data. Comput Vis Graph Image Process 1986;33(3):327–45. doi:10.1016/0734-189X(86)90181-7. Can E. Piecewise cubic approximation for data. Am J Appl Math 2013;1(2):24. doi:10.11648/j.ajam.20130102.11. Marple SL Jr. Computing the discrete-time analytic signal via FFT. Signal Process IEEE Trans 1999;47(9):2600–3. Zhou C, Yang L, Liu Y, Yang Z. A novel method for computing the Hilbert transform with Haar multiresolution approximation. J Comput Appl Math 2009;223(2):585–97. doi:10.1016/j.cam.2008.02.006. Micchelli CA, Xu Y, Yu B. On computing with the Hilbert spline transform. Adv Comput Math 2013;38(3):623–46. doi:10.1007/s10444- 011- 9252- x. Bilato R, Maj O, Brambilla M. An algorithm for fast Hilbert transform of real functions. Adv Comput Math 2014;40(5-6):1159–68. doi:10.1007/ s10444- 014- 9345- 4. Henrici P. Applied and computational complex analysis. Discrete Fourier analysis, Cauchy integrals, construction of conformal maps, univalent functions, vol. 3. John Wiley & Sons; 1993. Dagnino C, Santi E. Spline product quadrature rules for cauchy singular integrals. J Comput Appl Math 1990;33(2):133–40. doi:10.1016/0377-0427(90) 90363-5. Monegato G. The numerical evaluation of one-dimensional Cauchy principal value integrals. Computing 1982;29(4):337–54. doi:10.1007/BF02246760. Diethelm K. Modified compound quadrature rules for strongly singular integrals. Computing 1994;52(4):337–54. doi:10.1007/BF02276881. Diethelm K. Peano kernels and bounds for the error constants of Gaussian and related quadrature rules for Cauchy principal value integrals. Numer Math 1996;73(1):53–63. doi:10.10 07/s0 02110 050183. Elliott D, Paget DF. Gauss type quadrature rules for Cauchy principal value integrals. Math Comput 1979;33(145):301–9. Diethelm K. Uniform convergence of optimal order quadrature rules for Cauchy principal value integrals. J Comput Appl Math 1994;56(3):321–9. Criscuolo G. A new algorithm for Cauchy principal value and Hadamard finite-part integrals. J Comput Appl Math 1997;78(2):255–75. Orsi AP. Spline approximation for Cauchy principal value integrals. J Comput Appl Math 1990;30(2):191–201. Damelin S, Diethelm K. Interpolatory product quadratures for Cauchy principal value integrals with Freud weights. Numer Math 1999;83(1):87–105. doi:10.10 07/s0 02110 050440. Damelin SB, Diethelm K. Boundedness and uniform numerical approximation of the weighted Hilbert transform on the real line. Numer Funct Anal Optim 2001;22(1-2):13–54. doi:10.1081/NFA-100103786. Weideman JAC. Computing the Hilbert transform on the real line. Math Comput 1995;64(210):745–62. Chen Q, Huang N, Riemenschneider S, Xu Y. A B-spline approach for empirical mode decompositions. Adv Comput Math 2006;24(1-4):171–95. doi:10. 1007/s10444- 004- 7614- 3. King FW, Smethells GJ, Helleloid GT, Pelzl PJ. Numerical evaluation of Hilbert transforms for oscillatory functions: a convergence accelerator approach. Comput Phys Commun 2002;145(2):256–66. doi:10.1016/S0010-4655(02)00155-8. Hasegawa T. Uniform approximations to finite Hilbert transform and its derivative. J Comput Appl Math 2004;163(1):127–38. Schussler H, Weith J. On the design of recursive Hilbert-transformers. In: Acoustics, speech, and signal processing, IEEE international conference on ICASSP’87., 12. IEEE; 1987. p. 876–9. Rabiner LR, Schafer RW. On the behavior of minimax FIR digital Hilbert transformers. Bell Labs Tech J 1974;53(2):363–90.