Numerical Hilbert Transform Algorithm for Causal Interpolation of Functions Represented by Cubic and Exponential Splines

This paper presents an algorithm for numerical Hilbert transform of functions represented by cubic and exponential splines, which is suitable for causal interpolation of data spanning several frequency decades. It does not suffer from excessive number of data points due to large frequency span, and does not exhibit aliasing, both of which are characteristic for Hilbert transform algorithms based on FFT. The proposed algorithm was used for causal interpolation of characteristic impedance of microstrip line over conductive substrate, and the results were compared to the output of CausalCat.


I. INTRODUCTION
Physical systems are causal, and as a consequence real and imaginary parts of physical complex-valued quantities, such as permittivity = + j , characteristic impedance Z 0 etc., are not independent, but are related by Kramers-Kroning relations [1]. Independent interpolation of real and imaginary parts will violate the Hilbert relations and result in a noncausal response. To preserve causality, only one quantity of the Hilbert transform pair should be interpolated, while the other one should be calculated by the Hilbert transform. In this way the Kramers-Kronig relations are enforced and the causality is preserved.
The importance of Hilbert transform resulted in multitude of algorithms [2]- [6]. Most general-purpose Hilbert transform algorithms are based on the Fast Fourier Transform (FFT), which require fixed frequency step, and are prone to aliasing. Since causal interpolation is usually performed on data spanning several frequency decades, the use of fixed frequency step in FFT-based algorithms results in excessive number of data points. Although FFT is very efficient, excessive number of data points can results in slow computation as shown in [6].
Additional pre-and post-processing of data might be necessary when using FFT-based algorithms to prevent aliasing, such as constructing a function which ensures that input data is periodic, as was done in CausalCat [7]. Construction of additional functions is problem-dependent and has to be addressed on case-by-case basis.
The associate editor coordinating the review of this manuscript and approving it for publication was Giovanni Angiulli .
The algorithm proposed in this paper allows evaluation of Hilbert transform at any frequency. This property allows efficient evaluation of Hilbert transform on logarithmic frequency scale over several decades, where fixed frequency step algorithms would fail due to excessive number of data points. Furthermore, the proposed algorithm does not exhibit aliasing, so no pre-or post-processing is necessary. It can be viewed as an extension of [6] which allows the use of exponential splines to eliminate the undesired inflexion points which occur if only cubic splines are used for interpolation. As in [6], the Hilbert transform of interpolating function is exact, so the quality of results depends solely on how well the interpolating function describes the input data.
The paper is organized as follows. A short review of exponential splines necessary for derivation of the proposed algorithm is given in Sec. II. Derivation of the proposed algorithm is presented in Sec. III, followed by consideration of special cases in Sec. IV. Example of using the proposed algorithm for causal interpolation of characteristic impedance of microstrip line over conductive substrate and comparison with published results is given in Sec. VI. Conclusion and final remarks are in Sec. VII. Some properties of the Hilbert transform used in derivation are given in the Appendix.
can be used to avoid inflexion points which occur when polynomial interpolating functions are used. Coefficients of i-th exponential spline A i , B i , C i and D i and parameter σ i are determined by a fitting algorithm, and frequency points ω i determine the frequency range of each segment. Parameter σ i determines the shape of the interpolating function f i , and can be used to select the between cubic (σ i = 0), exponential (0 < σ i < ∞) or linear interpolating function (σ i → ∞). It has been shown in [8] that there exists a value of parameter σ * i ≥ 0, such that for σ i ≥ σ * i undesired inflexion points do not occur. The problem of determining the values of σ i has been thoroughly discussed in [9], and the following text will use the results of it.
The interpolating function f i (ω) in the interval (ω i , ω i+1 ), equivalent to (1), can be expressed as [9]: where Formulation of interpolating function f i given in (2) is convenient because it depends only on the input data ω i and y i , calculated parameters σ i and first derivatives y i . Parameters σ i and first derivatives y i can be obtained from input data by invoking the algorithm from [9] which is available as FORTRAN and Matlab code.
It is difficult to calculate the Hilbert transform of the interpolating function f i (ω) in the form of (2), but it can be rearranged to a more manageable form by using identities cosh(ω) = (e ω + e −ω )/2 and sinh(ω) = (e ω − e −ω )/2: (13) where coefficients β (i) j , j = 0..3 are given by: β (i) The Hilbert transform of (13) can be calculated, as will be shown in the following section.
The algorithm [9] can output σ i = 0, so it is necessary to determine the form of interpolating function for that case. Letting the value of σ → 0 results in: which is a cubic polynomial. Polynomial (18) can be rearranged to a more convenient form: where constants γ (i) j , j = 0..3 are given by:

III. HILBERT TRANSFORM OF INTERPOLATION FUNCTION
Even or odd function f (ω) known at discrete points can be interpolated as: where f + (ω) is the interpolation of f (ω) for positive values of ω, f + (ω < 0) = 0, and parameter p determines whether the function f has even or odd symmetry: Interpolating function f + (ω) can be constructed from f (ω) known at points ω i ≥ 0, y i = f (ω i ), as a sum of N − 1 piecewise interpolating functions f i (ω) as: where χ i (ω) is the characteristic function defined by: Characteristic function χ i (ω) windows the piecewise inter- Our goal is to determine the Hilbert transform of (24) in order to perform causal interpolation. Hilbert transform H of a function f (ω) is defined by [10, pp. 239]: where P denotes the Cauchy principal value of the integral.
To avoid cumbersome notation, we will adopt the short- where it is unambiguous, which denotes that the Hilbert transform of function f (µ) is evaluated at ω. In the following text we will use the interval (ω i , ω i+1 ) and index i interchangeably so Using the linearity (71) and reflection (73) properties of the Hilbert transform, listed in Appendix together with other properties of interest, F(ω) = Hf (ω) is then: Equation (29) shows that to calculate the Hilbert transform F(ω) = Hf (ω) it is sufficient to calculate F + (ω) = Hf + (ω) and evaluate it for ω and −ω.
Taking the Hilbert transform of (26) and using the linearity property results in: if the sum exists. Special cases when the order of summation and Hilbert transform cannot be exchanged will be considered in Sec. IV. Examining (30) reveals that the problem of calculating the Hilbert transform of f (ω) reduces to calculating the Hilbert transform of interpolating functions multiplied by the characteristic function defining the interval of interpolation.
If the value of parameter σ i is not zero, the i th interpolating function is an exponential spline of the form (13), which is rewritten for clarity: can be calculated term by term by using the linearity property. For the first term, χ i (ω) · const, the Hilbert transform is [10, pp. 243] Hilbert transform of the second term, and using the moment formula (74) on χ i (ω)ω: Using the linearity property The third and fourth term of (13) are of the form e ±σ i (ω i+1 −ω)/h i , and their Hilbert transform can be calculated from H {e µ χ i (µ)} and using the linear scale change property: where the change of variable t = −(µ − ω) was used in the second step. Integral in (34) can be expressed in terms of special functions by noting that: where Ei(ω) is the exponential integral function Using (35) and (36) in (34) the transform H {e µ χ i (µ)} is: The Hilbert transform of T , as shown at the bottom of the next page can be obtained by using the linear scale change property (72) in (37), where + sign is used for T The reason for rearranging (38) is that ExpEi(ω) can be efficiently approximated with rational Chebyshev functions [11], [12]. Numerical code for calculating ExpEi(ω) is available from Netlib repository in the Specfun FORTRAN package.
Finally, the Hilbert transform of exponential spline interpolating function χ i (ω)f i (ω) is: If the value of parameter σ i is zero, the i th interpolating function is a cubic spline of the form (19), rewritten for clarity: can be calculated by using the general moment formula (75) when the sum exists, and the sum and Hilbert transform can be exchanged: Special cases will be treated in Sec. IV. Using Hilbert transform general moment formula (75) for ω n χ i (ω) evaluates to: Using a binomial formula, Hilbert transform of χ i (ω)(ω − ω i ) n can then be expressed as: Using (41) in (42) results in: where P (n) (ω) denotes a polynomial of degree n. The Hilbert transform of χ i (ω)f 0 i (ω) is then: Coefficients of second degree polynomial P (2) (ω − ω i ) can be evaluated by direct summation, and the final result for where coefficients δ

IV. SPECIAL CASES FOR ω = ω I
Inspecting the Hilbert transform of exponential (40), (31)-(38), and cubic (45) spline interpolating functions reveals that they contain terms log |ω i − ω| and ExpEi(±σ i (ω − ω i )/h i ) which are singular at ω i . The singularity at ω = ω i arises because the Hilbert transform was performed over individual interpolation intervals (ω i , ω i+1 ), and taking the principal value of integral in (28) could not remove it. Consequently, for ω = ω i the sum in (30) does not exist, and the order of Hilbert transform and summation cannot be exchanged. The value of F(ω i ) can be calculated by observing that only F i−1 and F i contain singular terms, and that calculating the Hilbert transform over two interpolation intervals (ω i−1 , ω i+1 ) will remove the singularity at ω i by taking the principal value of the integral. Using this insight, the sum (30) can be rearranged to: The Hilbert transform over two intervals is where the Cauchy principal value of integral excludes the singular point ω i by splitting the integration interval (ω i−1 , ω i+1 ) to intervals (ω i−1 , ω i − ) and (ω i − , ω i+1 ), and taking the limit → 0. By introducing F i (ω) and F i (ω), defined as equation (50) can be expressed as: Functions F i (ω) and F i (ω) can be thought of as the Hilbert transform of interpolating function f i (ω) where the interpolation interval is reduced by from the lower or upper side. Therefore, F i (ω) can be directly obtained from F i (ω) by replacing ω i with ω i + , while F i (ω) is obtained by replacing ω i+1 with ω i+1 − . Taking the Cauchy principal value of integrals in (51) and (52) is not needed, since the evaluation points ω i will always be outside of the integration interval and they will always be well-defined. Values of F i (ω i ) and F i (ω i+1 ) are needed to evaluate (53), and they are derived for cases of cubic and exponential splines in the following text.
Function F 0 i (ω i ) can be derived from (45) by replacing the lower integration limit ω i with ω i + and substituting ω = ω i .
and can be rearranged as: where the non-singular term F 0 i (ω i ) is: Similarly, F i (ω i+1 ) can be obtained by replacing the upper integration limit ω i+1 with ω i+1 − and substituting ω = ω i+1 , which results in: The Hilbert transform of exponential spline (38) contains both logarithmic and exponential integral terms which are singular at ω i and ω i+1 . Logarithmic terms can be split into singular and non-singular terms, as was previously done for cubic spline case, while the exponential integral terms can be expanded to power series at singular points. Series expansion where γ E ≈ 0.5772 is the Euler-Mascheroni constant.
Function F i (ω i ) can be obtained from (40) by replacing the lower integration limit ω i with ω i + , substituting ω = ω i , and using the series expansion (59) around y = ω i in (38) for terms containing ω − ω i . After collecting singular and non-singular terms the result is: where the non-singular term F i (ω i ) is In a similar manner, function F i (ω i+1 ) can be obtained from (40) by replacing the upper integration limit ω i+1 with ω i+1 − , substituting ω = ω i+1 , and using the series expansion (59) around y = ω i+1 in (38) for terms containing ω − ω i+1 . After collecting singular and non-singular terms the result is: where the non-singular term F i (ω i+i ) is C. EVALUATION AT POINT ω i Using the previous results, the Hilbert transform over two return β, γ , δ 13: end procedure intervals (53) can be expressed as: Since the interpolating function is assumed to be continuous, it follows that f i−1 (ω i ) = f i (ω i ), and the singular term log | | vanishes. Finally, the value of the Hilbert transform over two intervals (ω i−1 , ω i+1 ) at point ω i is: Result (65) is valid for any combination of cubic and exponential splines, because both of (56) and (61) contain the term −π −1 f i (ω i ) log | |, and both of (58) and (63) contain the term π −1 f i (ω i ) log | |, which cancel out in (64).

V. ALGORITHM FOR EVALUATION
The Hilbert transform of a function interpolated with a mix of cubic and exponential splines can be expressed as: Evaluation of (66) can be divided into setup and evaluation steps. In the setup setup, shown in Alg. 1, coefficients which do not depend on frequency are evaluated, and needs to be performed only once for a given set of interpolating splines.
Algorithm 2 evaluates the Hilbert transform of function interpolated by a mix of cubic and exponential splines at a given frequency. The algorithm evaluates the contribution of individual interpolating functions at ±ω and accumulates Algorithm 2 Hilbert Transform of Interpolation Function 1: procedure Hilbert(p, ω) microstrip transmission line over conductive substrate [13]. The characteristic impedance data used in this example, taken from CausalCat [14], is calculated by power-voltage definition in [15] from full-wave simulation of 5 µm wide microstrip over a 1 µm oxide stacked on top of 100 µm conductive silicon substrate. Theory of causal interpolation is explained in detail in [14]- [16], and a short overview is given here for completeness. Characteristic impedance is causal and minimum-phase function with magnitude-phase relation: Using the property of inverse Hilbert transform (76) results in: where λ is a constant. Magnitude of characteristic impedance is then Causal interpolation of characteristic impedance can be performed by interpolating the phase arg(Z 0 (ω)) and calculating the magnitude (69). Additionally, using (69) can reveal whether the voltage integration path is correctly chosen [13].
Causal interpolation of characteristic impedance has been implemented in NIST software CausalCat [7], which uses FFT based Hilbert transform, and in [6], which uses analytic solution for piecewise polynomial interpolating function. CausalCat computes the causal characteristic impedance magnitude by interpolating the input data phase, sampling the interpolated function with a fixed frequency step and performing the Hilbert transform by FFT. Using the results from [6] for n = 3 is equivalent to forcing the proposed algorithm to use cubic splines for all interpolating functions by setting σ i = 0.
Hilbert transform is non-local, meaning that the value of Hilbert transform Hf (ω) at a given frequency depends on the entire function f . In practice, measured quantity will be known up to some frequency ω max , which will require to make assumptions about the frequency behavior above ω max . In the case of causal characteristic impedance interpolation, the error at frequency ω due to unknown behavior above ω max is bounded by [16]: where Z 0 is the true characteristic impedance andZ 0 is the calculated characteristic impedance. The error due to unknown high frequency behavior can be made arbitrarily small by restricting to frequencies where the error is acceptable.
In the microstrip example, phase of characteristic impedance approaches ±π/4 as ω → 0, tends to zero as frequency increases and might exhibit complicated behavior at very high frequencies [14]. Since the frequency range of measurement equipment is not infinite, an assumption about   the phase at frequencies above the measurement range must be made. By assuming that the phase is zero when ω > ω max an error bounded by (70) is introduced. To avoid aliasing, CausalCat corrects the phase at ω = 0 by summing the input data with a function f = ±(π/4)/(1 + a|ω|/8), so ϕ(0) = ϕ(ω max ) = 0, and by subtracting Hf from the result. Such correction is not needed in the proposed algorithm because it does not suffer from aliasing.
Interpolation of arg(Z 0 (ω)) with cubic and mix of cubic and exponential splines is shown in Fig. 1. In both cases 23 interpolating functions were used, where for the case of mixed splines 11 interpolating functions were exponential splines. From zoomed inset in Fig. 1 it can be seen that cubic splines exhibit oscillations due to undesired inflexion points. The use of exponential splines eliminates these oscillations.
Phase of characteristic impedance is an odd function of frequency, so symmetry parameter p should be set to p = 1 when calculating its Hilbert transform. Causal interpolation of characteristic impedance magnitude was performed by calculating the interpolated phase Hilbert transform and using (69), and the results are shown in Fig. 2. It can be seen from zoomed inset that using only cubic splines results in oscillations, while using a mix of cubic and exponential splines does not.
Relative error of calculated |Z 0 | for the case of cubic (dashed) and mix of cubic and exponential (solid) splines with CausalCat results as reference is shown in Fig. 3. It can be seen that both cubic and mixed splines interpolation are in good agreement with referent results. The use of exponential splines results in lower overall relative error compared to using only cubic splines. Lower relative error is obtained by eliminating undesired inflexion points introduced by interpolating with cubic splines.

VII. CONCLUSION
Algorithm for numerical Hilbert transform of functions represented by cubic and exponential splines is presented in this paper. Allowing the interpolating function to contain both cubic and exponential splines can eliminate undesired inflexion points which can occur if only cubic splines are used. The algorithm calculates exact Hilbert transform of interpolating function, so the accuracy of results depends solely on the quality of interpolation. It is suitable for causal interpolation of data spanning several decades without requiring excessive number of data points, and does not suffer from aliasing. Proposed algorithm was used for causal interpolation of characteristic impedance spanning five frequency decades, and the results are compared to the output of CausalCat.