On the Integration Schemes of Retrieving Impulse Response Functions from Transfer Functions

Abstract

The numerical inverse Laplace transformation (NILM) makes use of numerical integration. Generally, a high-order scheme of numerical integration renders high accuracy. However, surprisingly, this is not true for the NILM to the transfer function. Numerical examples show that the performance of higher-order schemes is no better than that of the trapezoidal scheme. In particular, the solutions from high-order scheme deviate from the exact one markedly over the rear portion of the period of interest. The underlying essence is examined. The deviation can be reduced by decreasing the frequency-sampling interval.

1. Introduction

Some linear dynamic systems are described at first in the frequency domain (FD) via the frequency response function (FRF) 𝐻(𝑗𝜔) or transfer function 𝐻(𝑠). For example, the characteristic of an unbounded media is relatively conveniently described in FD, as are the cases of soil-structure interaction and crack analysis [1–6]. Another salient example is where the transfer function is modified directly in FD to match some special material properties, as is the case of a hysteretic damping model [7–12]. The media property of attenuating wave propagation is also easily described via the FD expression [13].

The unitary impulse response function (UIRF) ℎ(𝑡) is the inverse Fourier transform of the FRF as1ℎ(𝑡)=2𝜋^+∞_−∞𝐻(𝑗𝜔)exp(𝑗𝜔𝑡)𝑑𝜔 (1.1) or the inverse Laplace transform1ℎ(𝑡)=2𝜋𝑗^𝑐+𝑗∞_{𝑐−𝑗∞}𝐻(𝑠)exp(𝑠𝑡)𝑑𝑠𝑐≥𝑐₀. (1.2) Here the real number 𝑐₀ is the convergence abscissa, that is, all the poles lie on the left side of Re(𝑠)=𝑐₀.

Though ℎ(𝑡) and 𝐻(𝑗𝜔) are equivalent for describing a linear time invariant system, sometimes, ℎ(𝑡) is preferred, such as when inspecting the system causality, or computing in time domain [3, 7, 8, 12, 14]. For simple cases, the transform from 𝐻(𝑠) to ℎ(𝑡) can be carried out analytically [15–17]. But from an engineering point of view, a numerical approach is recommended, especially when the closed-form solutions do not exist, for example, in the case of an ideal hysteretic model [8–11]. The numerical inverse Laplace transform (NILT) appears in many engineering problems [18, 19], and lots of algorithms have been constructed. Novel approaches are still developing, for example, Wang’s approach based on the wavelet [20].

The first issue in implementing the NILT is the infinite integral bound of (1.2). This is addressed by choosing a large enough 𝜔_𝑇 to accommodate the essential part of 𝐻(𝑠), and ignoring contribution beyond the [−𝜔_𝑇,𝜔_𝑇], that is,1ℎ(𝑡)≈2𝜋𝑗^{𝑐+𝑗𝜔𝑇}_{𝑐−𝑗𝜔𝑇}𝐻(𝑠)exp(𝑠𝑡)𝑑𝑠𝑐≥𝑐₀. (1.3)

Provided that |𝐻(𝑠)| attenuates very fast as |𝜔|→∞, and 𝜔_𝑇 is large enough, then the approximation (1.3) is acceptable. Otherwise, if |𝐻(𝑠)| attenuates slowly, then, even if |𝐻(𝑠)| is very small for |𝜔|>𝜔_𝑇, we cannot ignore the truncating error arbitrarily, because the contribution from a slowly attenuating function integrating over an infinite bound can be significant [21].

Authors have proposed a semianalytic approach to accommodate the slowly attenuating function [7]. This thanks to that, different from an arbitrary function, the transfer function is the description of a linear system. For most cases in engineering, the systems’ poles are bounded in a circle around the origin. Otherwise, if there were poles at infinity, then the system would vibrate at infinite frequency, which would impose infinite energy. This is not realistic. In light of this characteristic, frequently studied realistic systems possess this property: as |𝜔|→∞, the asymptotic transfer function property becomes simple, and can be approximated with a simple function Ψ(𝑠).

Assume that |𝜔|>𝜔_𝑇, Ψ(𝑠) matches 𝐻(𝑠) very well. Then (1.2) can be written as1ℎ(𝑡)=2𝜋𝑗^𝑐+𝑗∞_{𝑐−𝑗∞}1[𝐻(𝑠)−Ψ(𝑠)]exp(𝑠𝑡)𝑑𝑠+2𝜋𝑗^𝑐+𝑗∞_{𝑐−𝑗∞}Ψ(𝑠)exp(𝑠𝑡)𝑑𝑠. (1.4) Denote the original function of Ψ(𝑠) as 𝜓(𝑡), then1ℎ(𝑡)≈2𝜋𝑗^{𝑐+𝑗𝜔𝑇}_{𝑐−𝑗𝜔𝑇}[𝐻(𝑠)−Ψ(𝑠)]exp(𝑠𝑡)𝑑𝑠+𝜓(𝑡). (1.5)

Choice of Ψ(𝑠) is case-dependent. However, many systems own an asymptotic form as Ψ(𝑠)=𝜂/𝑠^𝛽, where 𝜂 and 𝛽>0 are constants. The original function is 𝜂𝑡^𝛽−1/Γ(𝛽), where Γ(𝛽) is the Gamma function. It is verified that this approach combining the numerical integration of trapezoidal scheme has effectively inhibited the ripple error [7].

2. Numerical Schemes

To increase accuracy further, another approach is employing high-order integration schemes. However, few authors touched on applying high-order schemes to the Fourier transform and Laplace transform. Only the rectangular scheme and the trapezoidal scheme (TS) are frequently used. One possible reason is that the ripple error might overcast the potential improvement by a high-order scheme. Since the ripple error disappears after compensating the high band contribution by a simple function, the influence of high-order schemes will be examined here.

Denote 𝐺(𝑠)=𝐻(𝑠)−Ψ(𝑠). Assume that the number of sampling intervals is an even 𝑁 (as a consequence, the sampling interval Δ𝜔=2𝜔_𝑇/𝑁). Numerical approximation of the continuous integration in (1.5) is1𝑔(𝑡)=2𝜋𝑗^{𝑐+𝑗𝜔𝑇}_{𝑐−𝑗𝜔𝑇}𝐺(𝑠)exp(𝑠𝑡)𝑑𝑠≈Δ𝜔2𝜋^𝑁_𝑙=0𝐺_𝑙exp(𝑠_𝑙𝑡), (2.1) where 𝑠_𝑙=𝑐+𝑗𝑙Δ𝜔−𝑗𝜔_𝑇. 𝐺_𝑙 in (2.1) depends on the integration scheme. For the trapezoidal scheme (TS)𝐺_𝑙=𝐺(𝑠_𝑙)2𝑙=0or𝑁,𝐺(𝑠_𝑙)0<𝑙<𝑁. (2.2) For the Simpson’s rule (SR)𝐺_𝑙=𝐺(𝑠_𝑙)3𝑙=0or𝑁,4𝐺(𝑠_𝑙)3odd𝑙and0<𝑙<𝑁,2𝐺(𝑠_𝑙)3even𝑙and0<𝑙<𝑁. (2.3) For the Boole's rule(BR) (http://mathworld.wolfram.com/BoolesRule.html):𝐺_𝑙=14𝐺(𝑠_𝑙)45𝑙=0or𝑁,28𝐺(𝑠_𝑙)45𝑙=4𝑘(excluding0and𝑁),64𝐺(𝑠_𝑙)45𝑙=4𝑘+1,24𝐺(𝑠_𝑙)45𝑙=4𝑘+2,64𝐺(𝑠_𝑙)45𝑙=4𝑘+3. (2.4)

To facilitate the fast Fourier transform (FFT), the continuous 𝑡 in (2.1) must be sampled at discrete instants. The FFT prefers a time interval Δ𝑡=2𝜋/(2𝜔_𝑇)=𝜋/𝜔_𝑇. For 𝑡=𝑘Δ𝑡, we have exp(𝑠_𝑁𝑘Δ𝑡)=exp(𝑐𝑘Δ𝑡)exp(−𝑗𝜔_𝑇𝑘Δ𝑡)=exp(𝑠₀𝑘Δ𝑡),exp(𝑠_𝑙𝑘Δ𝑡)=(−1)^𝑘exp(clΔ𝑡+𝑗2𝜋𝑙𝑘𝑁). (2.5) Substituting (2.5) into (2.1) leads to 𝑔(𝑘Δ𝑡)=(−1)^𝑘Δ𝜔exp(𝑐𝑘Δ𝑡)2𝜋^𝑁−1_𝑙=0𝐺_𝑙exp(𝑗2𝜋𝑙𝑘𝑁), (2.6) where 𝐺_𝑙=𝐺_𝑙 for 0<𝑙<𝑁 and 𝐺₀=𝐺₀+𝐺_𝑁. The summation in (2.6) can be sped up via a canonical FFT. Note that (2.6) can automatically avoid the imaginary issue that occurred in the rectangular numerical integration, which was highlighted in [22].

3. Numerical Examples

3.1. Influence of the Integration Schemes

We have examined several systems, including the exponentially attenuating model, viscous damping model, hysteretic damping model, and infinite elastic layer model. The ripple error can be inhibited significantly. Nonetheless, the performance of high-order integration schemes is counterintuitive, that is, high-order integration schemes yield rather coarse results. To stress this point, only the results from the exponentially attenuating model, the simplest case, are presented here.

The original function is ℎ(𝑡)=exp(−𝛼𝑡), and its transfer function is 𝐻(𝑠)=1/(𝑠+𝛼). This is the URIF of the first order linear differential equation, and appears ubiquitously in electrical circuit analysis. This function also occurs in the one-dimension scaled boundary finite-element analysis of unbounded medium [23].

Ψ(𝑠) is chosen as 1/𝑠, and the corresponding original function 𝜓(𝑡) is 1. Other computational parameters are listed in Figure 1. Inspecting this figure, we can see the result from the TS is free of ripples, and the computational trajectory follows the exact solution very well over the whole period [0,1].

Figure 1 

Comparison of the results from different numerical integration schemes.

Contradicting to the general impression, high-order integrations do not fatefully lead to good results. This is shown in Figure 1. The results of SR and BR deviate from the exact one after 𝑡>0.5, and are completely unjustifiable. From the view of numerical integration, the integrand 𝐻(𝑠)exp(𝑠𝑡) fluctuates too much, which does not satisfy the flatness premise of a numerical integration.

The parameter 𝑐 is very important in computation. The computational condition for Figure 2 is the same as those in Figure 1, except that 𝑐 is 9.22, twice that for Figure 1. We can see that the errors of the SR and BR are greater than those in Figure 1 by about 10 times. Still, there are significant ripples pertaining to the right end of curves from the trapezoidal scheme. It should be pointed that, the greater is 𝑐, the flatter 𝐻(𝑠) is, as indicated in Figure 3.

Figure 2 

Results with a larger 𝑐 (other parameters cf. Figure 1).

(a) real part

(b) imaginary part

Figure 3 

The flatness of original functions depending on 𝑐.

Equation (2.1) relates with the numerical inverse Fourier Transform. Insofar as authors know, few literature utilizes the high-order integrations to approximate inverse Fourier Transform, thus the error knowledge is rare. The frequently used numerical integrations are based on approximating the integrand with low-order polynomials. The belief that high-order integration schemes can improve accuracy is that the integrand between the sampling points is flat enough, so that a low-order polynomial can accommodate it. However, the integral 𝐺(𝑠)exp(𝑠𝑡)does not necessarily coincide with this assumption.

The significant error can also be appreciated along the following theme. In light of (2.1), 𝐺(𝑠)exp(𝑠𝑡) is multiplied by a periodic series 𝜅_𝑖. For SR (Similar analysis can be applied to the Boole's rule scheme), the inverse Fourier Transform of 𝜅_𝑖 is 𝛿(𝑡)−𝛿(𝑡+𝜋/Δ𝜔)), where 𝛿(𝑡) is the Dirac delta function. The multiplication in FD changes into the convolution in TD, that is,ℎ(𝑡)≈^∞_−∞1ℎ(𝑡−𝜏)[𝛿(𝜏)−13𝛿(𝜏+𝜋/Δ𝜔)]d𝜏=ℎ(𝑡)−3ℎ(𝑡+𝜋/Δ𝜔). (3.1)

Equation (3.1) explains the significant difference between the computational result and the original. Underlying the SR is that, the rear portion is folded to the front, with a weight 1/3. However this error will be zoomed by the exponential term in (2.6).

3.2. Influence of the Frequency Interval

Another approach to improve the computational result is decreasing the sampling interval Δ𝜔. For FFT, Δ𝜔 correlates with the retrieved time interval 𝑇 as 𝑇Δ𝜔=2𝜋. Thus, regardless of the accuracy, the FFT yields a computational UIRF up to a length of 𝑇. However, for most case, only a portion of FFT output, especially the front portion is of interest. For this case, the pruning FFT can decrease the computation burden. In addition, the constraint 𝑇Δ𝜔=2𝜋 can be uncoupled completely via the chirp-Z transform (CZT). With the CZT, the starting instant and terminating instant can be arbitrary. But the cost is that, the computation burden of one CZT is about six times of a compatible length FFT.

The influence of Δ𝜔 on the accuracy is depicted in Figure 4. Three cases of Δ𝜔, 2𝜋, 𝜋, 𝜋/2, are examined. The truncating frequency 𝜔_𝑇 keeps a constant of 256𝜋, as a consequence, the corresponding 𝑁 are 256, 512, 1024, respectively. Only the first 256 points are retained for computing the error trajectories shown in Figure 4(a). Decreasing the sampling interval Δ𝜔 can dramatically increase the accuracy over most part of [0,1], especially approaching 𝑡=1. This is taken for granted, since the ripples, if significant, have been expelled rightward.

(a) Error trajectories over [0,1]

(b) Error trajectories zoom around 𝑡=0

Figure 4 

Comparison of the influence of the sampling interval and integration schemes. TS: trapezoidal scheme; SR: Simpson’s rule; BR: Boole's rule.

However, the feature around 𝑡=0, which is zoomed in Figure 4(b), is interesting. First, the result is more erroneous than at other instants, which is probably due to the discontinuity around this instant. Another point is that, the results from Δ𝜔=𝜋/2 and Δ𝜔=𝜋 are almost comparable at 𝑡=0, although the former is better than the latter over most part of [0,1]. This lead to the guess that decreasing Δ𝜔 cannot improve the accuracy at 𝑡=0 any longer. In other words, improving the accuracy around this discontinuity necessitates increasing 𝜔_𝑇, which means weakening the contribution from side bands. To validate this guess, ℎ(0) computed with fives cases of 𝜔_𝑇(256𝜋,512𝜋,1024𝜋,2048𝜋,and4096𝜋) and two cases of Δ𝜔(𝜋/2and𝜋/4) are illustrated in Figure 5.

(a) Δ𝜔=𝜋/2

(b) Δ𝜔=𝜋/4

Figure 5 

Dependence of errors at 𝑡=0 on the frequency truncation bound.

This figure shows that, firstly, the error with Δ𝜔=𝜋/2 is comparable to that Δ𝜔=𝜋/4. This indicates again that the influence Δ𝜔 on accuracy is secondary if Δ𝜔 is greater than some threshold. Secondly, the error decreases as 𝜔_𝑇 increases. This reflects that the influence of truncating the integral limits is dominating in the error of ℎ(0), at least, with the computational parameters adopted in Figure 5. Thirdly, this figure also shows that three numerical schemes make no much difference to the error of ℎ(0).

Surprisingly, the result from the TS is even better than the results of higher-order integrations over the front portion in Figure 1, the result deviates significantly from the exact on the rear half, but not the front half. However, in Figure 4, the sampling interval Δ𝜔=𝜋/2 is so small that the significant erroneous occurring over the rear half is expelled from the period of interest. But, the accuracy of high-order schemes still is inferior to that of the trapezoidal scheme with an identical Δ𝜔.

4. Summary

In general, numerically retrieving the unitary impulse response function from the transfer function is contaminated by two kinds of errors, the truncating error due to ignoring the high-frequency band contribution, and the discrete error due to the numerical integration. The truncating error can be mitigated by approximating the transfer function over the high frequency band as a simple function.

Regard to the second error, to be counterintuitive is that, high-order integration schemes, rather than being accurate, produces results deviating from the exact dramatically. Another approach to improve performance is decreasing the sampling interval, and then discarding the computed UIRF over the rear portion of the retrieved period.

Acknowledgments

The authors thank Mr. Stephen Crowsen for his help with editing this paper. The project is jointly sponsored by the National Science Foundation of China (nos. 30801299 and 30971692), the Excellent Scholar Cultivating Fund of Beijing Municipality, Scientific and Technological New Star Foundation of Beijing  (no. 2009A40), China Postdoctoral Science Foundation (no. 20090450261), and Chinese Universities Scientific Fund   (no. 2009-1-30).