Time Delay Extraction from Frequency Domain Data Using Causal Fourier Continuations for High-Speed Interconnects

We present a new method for time delay estimation using band limited frequency domain data representing the port responses of interconnect structures. The approach is based on the recently developed by the authors spectrally accurate method for causality characterization that employs SVD-based causal Fourier continuations. The time delay extraction is constructed by incorporating a linearly varying phase factor to the system of equations that determines Fourier coefficients. The method is capable of determining time delay using data affected by noise or approximation errors that come from measurements or numerical simulations. It can also be employed when only a limited number of frequency responses is available. The technique can be extended to multi-port and mixed mode networks. Several analytical and simulated examples are used to demonstrate the accuracy and strength of the proposed technique.


I. INTRODUCTION
Identification and extraction of time delay is an important research problem in signal processing and has applications in many fields including radar [25], sonar [33], [8], ultrasonics [10], microwave imaging [26], geophysics [22], seismology [39], [29], wireless communications [36] as well as modeling of passive structures in electronic systems, in particular, transmission line modeling [18], [11], transient simulation of interconnects [24] and co-simulation of passive structures with active devices in a time domain using SPICE. Passive structures in electronic systems have been traditionally analyzed in the frequency domain, while transient simulations are performed in the time domain using suitable models that accurately capture the relevant electromagnetic phenomena. The models are obtained from either direct measurements or electromagnetic simulations. Interconnect models are typically approximated by rational transfer functions using the vector fitting algorithm in various implementations [19], [15], [20], [12], [17], [13], [9], which is the standard macromodeling approach. As clock frequencies increase, the size of passive structures becomes of the same order as the signal wavelength at the operating frequency, which causes the distributed effects Lyudmyla L. Barannyk  such as time delay to play a significant role in the time domain simulations. For this reason, time delay has to be included in macromodeling, in particular, when causality is analyzed. The connection between causality and time delay is in the fact that time delays can pull a non-causal signal into the causal region or vice versa pull a causal signal into the non-causal region, while causality, in turn, can be expressed in terms of the Hilbert transform [32], [38], [31]. Several approaches can be used to extract delays in the frequency domain, for example, using the Hilbert transform [14], [35], [23], the minimum phase all-pass decomposition [28], [27], [24], incorporating an optimal time delay into the vector fitting algorithm [18], [11], employing a modified Lie approximation to develop a passive and compact macromodel [30], using a Gabor transform to develop delayed rational function macromodels for long interconnects [16], [9] or conducting a probabilistic analysis of the cepstrum in the presence of noise [21]. In the time domain, delayed rational functions [7], [6] can be employed to extract delays. In this paper, a novel approach is proposed in which time delay is determined in the frequency domain using a causality argument. Causality is verified using the SVDbased causal Fourier continuation method developed by the authors [4], [2], while the time delay presence is incorporated by a linearly varying phase factor to the system of equations that determines Fourier coefficients. Preliminary results are reported in [5].
The rest of the paper is organized as follows. Section II provides a background on causality for linear time-translation invariant systems and dispersion relations. In Section III, we show main steps in the derivation of causal Fourier continuations using truncated singular value decomposition (SVD) method that was developed to access causality. We also provide error estimates that take into account a possible presence of noise in data. Section IV extends the causality characterization method to develop a technique for time delay extraction. The proposed method is tested in Section V using several analytic and simulated examples. We also analyze the performance of the algorithm when only a limited number of frequency responses is available and when noise/approximation errors are present in data. In Section VI we present our conclusions. The Appendix section is devoted to formulation of error bounds for the causality characterization method based on causal Fourier continuations.

II. CAUSALITY OF LINEAR TIME-INVARIANT SYSTEMS
Consider a linear and time-invariant physical system with the impulse response h(t) subject to a time-dependent input f (t), to which it responds by an output x(t). Denote by the Fourier transform of h(t), which is also called the transfer function.
The system is causal if the output cannot precede the input, i.e. if f (t) = 0 for t < T , the same must be true for x(t). This primitive causality condition in the time domain implies h(t) = 0, t < 0. Hence, domain of integration in (1) can be reduced to [0, ∞).
Assume H(w) ∈ L 2 (R). Then starting from Cauchy's theorem and using contour integration, one can show [31] that for any point w on the real axis, H(w) can be written 1 as where denotes Cauchy's principal value. Separating the real and imaginary parts of (2), we get Expressions (3) and (4)   Evaluation of the Hilbert transform requires integration on (−∞, ∞), which can be reduced to [0, ∞) by spectrum symmetry of H(w) if h(t) is real valued. In practice, only a limited number of discrete values of H(w) is available on [w min , w max ]. Thus, the domain of integration has to be truncated. This usually causes serious boundary artifacts due to the lack of out-of-band frequency responses. To reduce or even completely remove boundary artifacts, the authors recently developed periodic polynomial [1], [3] and causal Fourier continuation [2], [4], respectively, based methods for causality characterization. The approach was motivated by the example H(w) = e −iaw , a > 0, that is not square integrable but still satisfies the dispersion relations. The causality characterization method based on causal Fourier continuations allows one to construct highly accurate approximations of a given transfer function on the original frequency interval [w min , w max ] with the uniform error that decreases as the number of Fourier coefficients increases. The technique is applicable to both baseband and bandpass cases and capable of detecting very small localized causality violations. The method can also be extended to multidimensional cases.
In the next section, for completeness of presentation, we show main steps in the derivation of the causal Fourier continuation method that can be used to access causality of a given transfer function whose values are available at a discrete set of frequencies. We also provide upper bounds of reconstruction error between the given function and its causal Fourier continuation. We use these error estimates to understand how to extract time delay when data with different resolutions are available and when data are affected by noise or other approximation errors.

III. CAUSAL FOURIER CONTINUATIONS
Consider a transfer function H(w) = Re H +i Im H, whose N discrete values are available on [w min , w max ], w min ≥ 0. The idea of a causal Fourier continuation is to construct an accurate Fourier series approximation of H(x) by allowing the Fourier series to be periodic and causal in an extended domain. The result is the Fourier continuation of H that we denote by C(H), and it is defined by for even number 2M of terms, whereas for odd number 2M +1 of terms, the index k varies from −M to M . Throughout this paper, we assume that the number M of Fourier coefficients is even, for simplicity. When M is odd, analogous results can be formulated. Here b is the period of approximation. For SVDbased periodic continuations b is normally chosen as twice the length of the domain on which function H is given though the value b = 2 is not necessarily optimal. The optimal value b depends on a function being approximated. In practice, several values b ∈ (1, 4) may be tried to get a better reconstruction of H(x) with a Fourier series.
. It can be shown that H{φ k (x)} = i sgn(k)φ k (x), which implies that functions {φ k (x)} are the eigenfunctions of the Hilbert transform H with associated . For a causal periodic continuation, according to (5), we need Im C(H)(x) to be the Hilbert transform of − Re C(H)(x). It can be shown [4] that this implies α k = 0 for k ≤ 0 in (6). Hence, a causal Fourier continuation has the form Evaluating H(x) at points x j , j = 1, . . . ,Ñ , x j ∈ [−0.5, 0.5], produces a complex valued system Since Re H(x) and Im H(x) are even and odd functions of x, respectively, the Fourier coefficients are real. Here¯denotes the complex conjugate. To ensure that numerically computed Fourier coefficients α k are real, instead of solving complex-valued system (8), one can separate the real and imaginary parts of C(H)(x j ) to obtain real-valued system We show in [4] that real formulation (9) provides slightly more accurate results than complex.
To access the quality of approximation of H(x) with its causal Fourier continuation C(H)(x), we introduce reconstruction errors E R (x) and E I (x), on the original interval [−0.5, 0.5].
The error analysis performed in [4] (see also Appendix) shows that the error between H(x) and its causal Fourier continuation C(H + ε) under the presence of a noise ε, has the following upper bound: Here is the error due to approximation of H with a causal Fourier series and it decays as O(M −k+1 ), where k is the smoothness order of the transfer function H(x).
is the error due to the truncation of singular values and it is typically small and close to the cut-off value ξ. As (14) indicates, ǫ T depends on b and the function H being approximated.
is the error due to the presence of a noise or approximation errors in the given data and it shows a level of causality violation. In practice the size of ǫ n is close to the size of noise in data. FunctionĤ M and constants Λ 1 , Λ 2 and K are defined in Appendix. These constants depend only on the continuation parameters N , M , b and ξ as well as location of discrete points x j , and not on the function H. The error bound (12) shows that the reconstruction errors E R and E I decrease as M increases due to the causal Fourier series approximation error with the error bound ǫ F until either the level ǫ of a noise or level ǫ T due to truncation of singular values is reached. If only round-off errors are present in data, the errors will level off at ǫ T . If reconstruction errors level off at some value ǫ > ǫ T as the resolution increases, the data are declared non-causal with the error approximately at the order of ǫ. More information about the error analysis for the causality characterization methods based on causal Fourier continuations can be found in [4].

IV. TIME DELAY ESTIMATION
The above approach for causality assessment can be transformed into a delay estimation algorithm by observing the following. Suppose that h(t) is non-zero only from time T 0 ≥ 0, and we would like to identify the time delay T 0 . Consider the Fourier transform H(w) of h(t): where we used the substitution x = aw, a = 0.5 wmax . Introducing τ = t a , we can write where G(x) is the Fourier transform of a causal function with no time delay. This implies that when 0 ≤ T ≤ T 0 , the transfer function H(x) e ix T a is causal, but when T ≥ T 0 , the transfer function H(x) e ix T a has a non-causal component. Therefore, T 0 = T 0 /a is the time delay for H(x), and the delay T 0 for the original function H(w) is recovered by multiplyingT 0 by a. Since one can add any integral multiple of 2π to xT /a, it is enough to restrict our investigations to the interval Then for each potential time delay 0 ≤ T a ≤ T max , we solve the following modified system or its equivalent real-valued formulation. For T < T 0 , the reconstruction errors E R , E I should be small and approximately of the same order. As T increases and becomes greater than some critical transition time close to the time delay T 0 , the reconstruction errors should start increase. The goal is to approximate T 0 . The difficulty is that the reconstruction errors grow gradually as T ≥ T 0 , so transition is not sharp. Moreover, the order of reconstruction errors for T < T 0 depends on the resolution of data and threshold ξ used in the truncated SVD method, which, in turn, affects a transition time. In addition, a noise in data, if present, also affects when reconstruction errors start growing. A similar approach was used in [23] to estimate the time delay for square integrable transfer functions. In this contribution, we extend the approach to more general transfer functions. In addition, we use a different causality measure than in [23] and take into account different resolutions of given data and a possible presence of noise. The approach can be extended to multi-port and mixed mode networks by applying it to each element of the transfer matrix.

V. NUMERICAL EXAMPLES
In this section, we apply a proposed technique to several analytic and simulated examples when the time delay is either known exactly or can be estimated using other techniques. We also consider the effect of noise presence on the accuracy of timed delay estimation.

A. Four-Pole Example
Consider a transfer function with four poles and time delay T 0 , defined by where r 1 = 1 + 2i, p 1 = 1 + 3i, r 2 = 2 3 + 1 2 i, p 2 = 1 2 + 5i, and T 0 = 0.25. Since the poles ofH(w) are located in the upper half w-plane at ±3 + i and ±5 + 1 2 i , this function is causal as a sum of four causal transforms, and has no time delay. Therefore, the function H(w) is a causal function delayed with offset T 0 . H is sampled on [0, w max ] at N frequency points varying from 50 to 1500 with w max = 6. The real and imag-  Fig. 2. Even though given H(x) and its causal Fourier continuation approximation look indistinguishable, the actual reconstruction errors E R and E I in both real and imaginary parts, that are defined in (10), (11), are on the order of 10 −6 and they decreases as M increases (with M = N ). For example, with M = 800, the errors are on the order of 10 −13 . Since both errors E R and E I are of the same order, it is enough to analyze one of the errors, for example, E R . The results using E I are similar.
To estimate the time delay, we analyze the evolution of the ||E R ||∞, shown in Fig. 3, for various values M . Since the error due to a causal Fourier series approximation decreases with M (see error bound (13)), the reconstruction error between the given transfer function H and its causal Fourier continuation C(H) also decreases as M increases until it either reaches the level ξ of filtering of singular values or a level ǫ of noise/causality violations (see error bounds (14) and (15)  errors are dominated by a causal Fourier series approximation error and then for T greater than some transition time -by causality violations since this value T provides large enough negative time delay and shifts a causal function into a noncausal area. A transition value T = T c , we call it a critical time, from a plateau region to a growth region, is different for each M and it decreases as the resolution or number of Fourier coefficients increases if the error is dominated by the causal Fourier series approximation error. The critical times T c approach the time delay T 0 as M increases. The goal is to estimate T 0 using the error curves shown in Fig. 3. Analyzing graphs of the error curves for M > 800, we observe some non-monotonic behavior at T close to T 0 . This behavior is due to the filtering of the singular values below the threshold ξ = 10 −13 that we used in our experiments. By increasing the value of ξ, the non-monotonic behavior will be present at smaller values of M . This suggests that portions of error curves close to threshold ξ are affected by filtering and may be inaccurate and difficult to use for time delay estimation as we find in our experiments. To estimate critical times T c of transition from the plateau region to the growth region, we approximate the growing region by a quadratic function on the loglog scale. Specifically, we assume that where coefficients a 0 , a 1 , and a 2 are determined in the least squares sense. The resulting quadratic function f (ln ||E R || ∞ ) is then evaluated at the value of ||E R || ∞ at T = 0 that is assumed to be the "most causal" time.  Table I. The results indicate that the approximations become more accurate as M increases. The error with M = 900 is less than 1 %. At the same time, the error with M = 1500 is about 3%, which is due to the fact that the results in this case are more affected by the filtering of singular values. In the cases when M is high and the resulting error is not flat for T < T 0 , instead of evaluating a fitted quadratic curve at the value of ||E R || ∞ at T = 0 we evaluate it at ξ, the threshold of filtering singular values, to avoid using results affected by filtering. In practice the number N of samples of the transfer function H(w) is usually limited, which sets the bound for number M = N of Fourier coefficients, so it may not always be possible to use large enough M to obtain critical time T c close enough to the actual time delay T 0 . A good method should be capable of producing an accurate approximation of T 0 even with a small number of data points. We achieve this by employing another approach for time delay estimation. Using the obtained fitted quadratic error curves, we extrapolate them to the value ξ of filtering of singular values, which is typically chosen to be close to the machine precision. This corresponds to finding time T at which the error reaches the value ξ. This choice is natural since the errors below ξ are most likely affected by filtering and may not be accurate enough to use. The results of such extrapolation are shown in Fig.  5 for M = 200, 400, 600 and 800. An intersection of the extrapolated curve corresponding to M = 200 is at a value T = 0.45451, which is a bit far from the exact T 0 = 0.25. At the same time, intersections of extrapolated curves with higher values of M are much closer -see Table II for details. Results    shown in this table indicate that as M increases, extrapolated quadratic curve intersect the horizontal line the value ξ at times closer to T 0 . Obtained approximations of T 0 can be averaged producing T 0 ≈ 0.24805. The approach with extrapolation provides a faster convergence and good approximations of T 0 even for small values of M , i.e. less data points are needed to approximate T 0 .
We also consider the effect of noise on the time delay estimation. To study this, we impose a sine perturbation a sin(10πx) of various amplitudes a that we add to Re H, while keeping Im H unchanged. We choose N = 800 and vary a from 10 −10 to 10 −3 . The reconstruction error E R with no perturbation for early times T < T 0 is of the order of 10 −13 , as shown in Fig. 6, that corresponds to the level of filtering of singular values. When the perturbation is added, the reconstruction errors for T < T 0 are higher and approximately of the order of a. Once some critical transition time greater than T 0 is passed, reconstruction errors start growing and they grow at the same rate and coincide almost perfectly with each other. This observation suggests that the proposed approach can also be used in the cases when data have a noise, which is typical in real-life applications, when data have either measurement or simulation errors. For noise with a smaller amplitude, the region close to T 0 will be less affected by noise and a bigger growing region will be available for fitting, so we expect better accuracy of time delay estimation in such cases. When more noise in data is present, less growing region will be available for fitting and extrapolation of fitted quadratic error curves may be less accurate. We demonstrate this by considering two cases: with a = 10 −5 (noisier case) and a = 10 −8 (less noisy case). The error curves with a higher amplitude a = 10 −5 are presented in Fig. 7. It is clear that the error does not become smaller than 10 −5 as M ≥ 300 gets larger because of the noise. We use available growing regions and extrapolate fitted error curves to find their intersection with the horizontal line with value ξ. This gives us time T when the error reaches the value ξ for each considered M . The results of such extrapolation for M = 200, 400, 600 and 800 are shown in Fig. 8. Clearly, extrapolated error curves reach value ξ at times around T 0 but not close enough to T 0 and without established convergence but rather in a spread-out manner around T 0 . Approximations of T 0 for values of M that we investigated are shown in Table III. Averaging these approximations we obtain T       Next we show results when a smaller noise of amplitude a = 10 −8 is added. The evolution of ||E R || ∞ as T increases is shown in Fig. 10. We can see that the plateau error region in this case is at about 10 −9 level, so the error growth region is bigger than in the previous case, which should make fitting and extrapolation more accurate. Indeed, extrapolated quadratic curves intersect the horizontal line with value ξ in a more localized region about T 0 as shown in Fig. 11, while averaging of obtained approximation to T 0 produces T (aver) 0 = 0.25436, which is more accurate than in the case with a higher amplitude a = 10 −5 .

B. Transmission Line Example
We consider a uniform transmission line segment with the following per-unit-length parameters: L = 7.574 nH/inch, C = 2.61166 pF/inch, R = 16.278 mΩ/inch, G = 5.58 µS/inch and length L = 5 inches. The frequency is sampled on the interval (0, 5.0] GHz. The scattering matrix of the structure is computed using Matlab function rlgc2s. We consider the elementH(w) = S 11 (w) and impose the time delay T 0 = 1.25 ns by multiplyingH(w) by exp(−iwT 0 ) to get the delayed transfer function H(w) = exp(−iwT 0 )H(w). The real and imaginary parts of H(w) are given in Fig. 12. The intersections with ||E R || ∞ at T = 0 or finding times when these fitted error curves reach the value ξ of the error for M ≥ 600, we get a sequence of critical transition times T c , that we show in Fig. 14. Clearly, critical times T c converge to T 0 and provide a good approximation of T 0 for M ≥ 500.
Using an alternative approach when we extrapolate the fitted quadratic error curves to find their intersections with the error threshold ξ, we find approximations of T 0 . Some of these curves for M = 200, 400, 600 and 800 are depicted in Fig.  15. Approximations of T 0 using extrapolation procedure for various values of M ranging from M = 200 to 1500 are given in Table V. A good approximation of T 0 is obtained even with M = 300. As before, approximations of T 0 become better as M increases, but for very large values of M ≥ 1000 when the reconstruction error falls below the filtering threshold ξ and

C. Dawson's Integral Example
We consider here another analytic example [23] modeled by the transfer function  Fig. 16. The evolution of ||E R || ∞ for various M is shown in Fig. 17, where it is clear that critical transition times approach T 0 . Constructing fitted quadratic error curves and extrapolating them to find their intersections with the horizontal line corresponding to the error value ξ produces a set of approximation of T 0 , shown in Table VI. Averaging obtained approximations of T 0 for M ≥ 200, once some convergence is established, gives T (aver) 0 = 0.12528. It is interesting to note behavior of the relative error E rel R in this example. The evolution of its ∞ norm is shown in Fig. 18. behavior. Even though the behavior of the relative error E rel R can be used to determine the time delay in this example, we did not find the same pronounced behavior in other examples we considered. At the same time, extrapolating fitted quadratic curves of ∞ norms of the absolute error E R was a robust approach in all examples we considered.

D. Stripline Example
We simulated an asymmetric stripline modeled in [37] with length L = 8 in, width W = 14 mils, distances from the trace to reference planes H 1 = 10 mils, H 2 = 20 mils, substrate dielectric Megtron6-1035, Laminate with a dielectric constant ǫ r = 3.45 using a Cadence software tool with FEM full-wave field solver. The stripline was simulated on [0, w max ] with w max = 2 GHz. We analyzed element H(w) = S 11 (w) of the transfer matrix. The real and imaginary parts of H are shown in Fig. 19. The evolution of ||E R || ∞ for various M is depicted in Fig. 20. Even for high values of M , the error in causality does not go to machine precision and instead levels off around 10 −6 . This indicates that our finite element simulation results are accurate only within 10 −7 − 10 −6 . For causality characterization, this implies that data have noise/approximation errors with amplitude around 10 −7 −10 −6 . Graphs of ||E R || ∞ suggest that for M ≤ 2000, the error is dominated by Fourier series approximation error, while for higher of M , the error is dominated by the noise/approximation errors from the finite element method.
In this example, the time delay was estimated using a closed form microwave theory approximation as T 0 = 8 × 0.0254/(c 0 / √ ǫ r ) = 1.25809 ns, where c 0 = 3 × 10 8 m/s is the speed of light, 0.0254 is a conversion factor to convert from inches to meters. The error curves were fitted to quadratic curves as explained above. Because of relatively high errors in data, the fitted regions are not long enough. Besides, there is more nonlinear behavior of the error curves for high values of T > T 0 . All this makes it difficult to estimate the time delay as shown in Fig. 21. As can be seen, extrapolated quadratic curves do not focus at T 0 but instead spread out around T 0 similar to the four-pole example with an imposed noise considered in Section V-A. This problem can be corrected   by decreasing the fitting range and going more away from transition regions. The results are shown in Fig. 22. The approximations of T 0 are given in Table VII

VI. CONCLUSIONS
We proposed a new method for time delay extraction from tabulated frequency responses. The approach uses the spectrally accurate causality enforcement technique constructed using SVD-based causal Fourier continuations, that was recently developed by the authors. The time delay is incorporated to the causality characterization approach by introducing a linear varying phase factor to the system of equations that defines Fourier coefficients. Varying time until a threshold time, that depends on the maximum frequency at which the transfer function is available, results in the reconstruction error between the given data and their causal Fourier continuations to go from an almost constant small value to a rapidly growing function at some critical transition time. The critical transition times depend on resolution and approach the time delay as resolution increases. Several sets of frequency responses with increasing resolution can be used to establish convergence and get an approximation of the time delay. Alternatively, when only a limited number of samples is available, a growing portion of the error curve can be extrapolated to find an approximation of the time delay. The method is applicable to data that have noise or other approximation errors. A few sets of frequency responses can be used to improve the accuracy of time delay approximation by averaging the obtained results. The technique can be extended for multi-port and mixed mode networks. The performance of the method is demonstrated using several analytic and simulated examples, including data with noise, for which time delay is known exactly or can be estimated using other approaches.
The following result is true [4]. Theorem : Consider a rescaled transfer function H(x) defined by symmetry on Ω = [−0.5, −a] ∪ [a, 0.5], where a = 0.5 wmin wmax , whose values are available at points x j ∈ Ω, j = 1, . . . , N . Then the error in approximation of H(x), that is known with some error ε, by its causal Fourier continuation C(H)(x) defined in (7)  and holds for all functions of the form (19). Here V kj φ k (x) are each an up to M term causal Fourier series with coefficients given by the jth column of V ; K denotes the number of singular values that are discarded, i.e. the number of j for which σ j < ξ, where ξ is the cut-off tolerance.
It can be seen that constants Λ 1 , Λ 2 and K depend only on the continuation parameters N , M , b and ξ as well as location of discrete points x j , and not on the function H.
For brevity, we can write the above error estimate as ||H − C(H + ε)|| L2(Ω) ≤ ǫ F + ǫ n + ǫ T . Here is the error due to a causal Fourier series approximation and it decays at least as fast as O(M −k+1 ), k is the smoothness order of H(x), which can be estimated numerically using reconstruction errors with different resolution (see [4]).
that is the error due to the truncation of singular values. It is typically small and close to the cut-off value ξ.
is the error due to noise ǫ in data. Numerical experiments reveal that ǫ n has the order of noise ǫ in data.