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

Authors: Kui Fu Chen; Yan Feng Li
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 1. Introduction

Some linear dynamic systems are described at first in the frequency domain (FD) via the frequency response function (FRF)H(jω) or transfer function H(s). 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)h(t) is the inverse Fourier transform of the FRF as(1.1)h(t)=12π∫-∞+∞H(jω)exp(jωt)dω
or the inverse Laplace transform(1.2)h(t)=12πj∫c-j∞c+j∞H(s)exp(st)dsc≥c0.
Here the real number c0 is the convergence abscissa, that is, all the poles lie on the left side of Re(s)=c0.Thoughh(t) and H(jω) are equivalent for describing a linear time invariant system, sometimes, h(t) 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 H(s) to h(t) 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 ωT to accommodate the essential part of H(s), and ignoring contribution beyond the [-ωT,ωT], that is,(1.3)h(t)≈12πj∫c-jωTc+jωTH(s)exp(st)dsc≥c0.Provided that|H(s)| attenuates very fast as |ω|→∞, and ωT is large enough, then the approximation (1.3) is acceptable. Otherwise, if |H(s)| attenuates slowly, then, even if |H(s)| is very small for |ω|>ωT, 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 Ψ(s).Assume that|ω|>ωT, Ψ(s) matches H(s) very well. Then (1.2) can be written as(1.4)h(t)=12πj∫c-j∞c+j∞[H(s)-Ψ(s)]exp(st)ds+12πj∫c-j∞c+j∞Ψ(s)exp(st)ds.
Denote the original function of Ψ(s)  as ψ(t), then(1.5)h(t)≈12πj∫c-jωTc+jωT[H(s)-Ψ(s)]exp(st)ds+ψ(t).Choice ofΨ(s) is case-dependent. However, many systems own an asymptotic form as Ψ(s)=η/sβ, where η and β>0 are constants. The original function is ηtβ-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.DenoteG(s)=H(s)-Ψ(s). Assume that the number of sampling intervals is an even N (as a consequence, the sampling interval Δω=2ωT/N). Numerical approximation of the continuous integration in (1.5) is(2.1)g(t)=12πj∫c-jωTc+jωTG(s)exp(st)ds≈Δω2π∑l=0NGlexp(slt),
where sl=c+jlΔω-jωT. Gl in (2.1) depends on the integration scheme. For the trapezoidal scheme (TS)(2.2)Gl={G(sl)2l=0orN,G(sl)0<l<N.
For the Simpson’s rule (SR)(2.3)Gl={G(sl)3l=0orN,4G(sl)3oddland0<l<N,2G(sl)3evenland0<l<N.
For the Boole's rule(BR) (http://mathworld.wolfram.com/BoolesRule.html):(2.4)Gl={14G(sl)45l=0orN,28G(sl)45l=4k(excluding0andN),64G(sl)45l=4k+1,24G(sl)45l=4k+2,64G(sl)45l=4k+3.To facilitate the fast Fourier transform (FFT), the continuoust in (2.1) must be sampled at discrete instants. The FFT prefers a time interval Δt=2π/(2ωT)=π/ωT. For t=kΔt, we have (2.5)exp(sNkΔt)=exp(ckΔt)exp(-jωTkΔt)=exp(s0kΔt),exp(slkΔt)=(-1)kexp(clΔt+j2πlkN).
Substituting (2.5) into (2.1) leads to (2.6)g(kΔt)=(-1)kΔωexp(ckΔt)2π∑l=0N-1G̃lexp(j2πlkN),
where G̃l=Gl for 0<l<N and G̃0=G0+GN. 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 ish(t)=exp(-αt), and its transfer function is H(s)=1/(s+α). 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].Ψ(s) is chosen as 1/s, and the corresponding original function ψ(t) 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 Figure1. The results of SR and BR deviate from the exact one after t>0.5, and are completely unjustifiable. From the view of numerical integration, the integrand H(s)exp(st) fluctuates too much, which does not satisfy the flatness premise of a numerical integration.The parameterc is very important in computation. The computational condition for Figure 2 is the same as those in Figure 1, except that c 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 c, the flatter H(s) is, as indicated in Figure 3.Figure 2
Results with a largerc (other parameters cf. Figure 1).The flatness of original functions depending onc.
(a)
real part(b)
imaginary partEquation (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 G(s)exp(st)does not necessarily coincide with this assumption.The significant error can also be appreciated along the following theme. In light of (2.1), G(s)exp(st) is multiplied by a periodic series κi. For SR (Similar analysis can be applied to the Boole's rule scheme), the inverse Fourier Transform of κi is δ(t)-δ(t+π/Δω)), where δ(t) is the Dirac delta function. The multiplication in FD changes into the convolution in TD, that is,(3.1)h̃(t)≈∫-∞∞h(t-τ)[δ(τ)-13δ(τ+π/Δω)]dτ=h(t)-13h(t+π/Δω).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 T as TΔω=2π. Thus, regardless of the accuracy, the FFT yields a computational UIRF up to a length of T. 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 TΔω=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 ωT keeps a constant of 256π, as a consequence, the corresponding N 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 t=1. This is taken for granted, since the ripples, if significant, have been expelled rightward.Comparison of the influence of the sampling interval and integration schemes. TS: trapezoidal scheme; SR: Simpson’s rule; BR: Boole's rule.
(a)
Error trajectories over[0,1](b)
Error trajectories zoom aroundt=0However, the feature aroundt=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 t=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 t=0 any longer. In other words, improving the accuracy around this discontinuity necessitates increasing ωT, which means weakening the contribution from side bands. To validate this guess, h(0) computed with fives cases of ωT(256π,512π,1024π,2048π,and4096π) and two cases of Δω(π/2andπ/4) are illustrated in Figure 5.Dependence of errors att=0 on the frequency truncation bound.
(a)
Δω=π/2(b)
Δω=π/4This 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 ωT increases. This reflects that the influence of truncating the integral limits is dominating in the error of h(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 h(0).Surprisingly, the result from the TS is even better than the results of higher-order integrations over the front portion in Figure1, 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 Δω.


## 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 ish(t)=exp(-αt), and its transfer function is H(s)=1/(s+α). 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].Ψ(s) is chosen as 1/s, and the corresponding original function ψ(t) 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 Figure1. The results of SR and BR deviate from the exact one after t>0.5, and are completely unjustifiable. From the view of numerical integration, the integrand H(s)exp(st) fluctuates too much, which does not satisfy the flatness premise of a numerical integration.The parameterc is very important in computation. The computational condition for Figure 2 is the same as those in Figure 1, except that c 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 c, the flatter H(s) is, as indicated in Figure 3.Figure 2
Results with a largerc (other parameters cf. Figure 1).The flatness of original functions depending onc.
(a)
real part(b)
imaginary partEquation (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 G(s)exp(st)does not necessarily coincide with this assumption.The significant error can also be appreciated along the following theme. In light of (2.1), G(s)exp(st) is multiplied by a periodic series κi. For SR (Similar analysis can be applied to the Boole's rule scheme), the inverse Fourier Transform of κi is δ(t)-δ(t+π/Δω)), where δ(t) is the Dirac delta function. The multiplication in FD changes into the convolution in TD, that is,(3.1)h̃(t)≈∫-∞∞h(t-τ)[δ(τ)-13δ(τ+π/Δω)]dτ=h(t)-13h(t+π/Δω).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 T as TΔω=2π. Thus, regardless of the accuracy, the FFT yields a computational UIRF up to a length of T. 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 TΔω=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 ωT keeps a constant of 256π, as a consequence, the corresponding N 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 t=1. This is taken for granted, since the ripples, if significant, have been expelled rightward.Comparison of the influence of the sampling interval and integration schemes. TS: trapezoidal scheme; SR: Simpson’s rule; BR: Boole's rule.
(a)
Error trajectories over[0,1](b)
Error trajectories zoom aroundt=0However, the feature aroundt=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 t=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 t=0 any longer. In other words, improving the accuracy around this discontinuity necessitates increasing ωT, which means weakening the contribution from side bands. To validate this guess, h(0) computed with fives cases of ωT(256π,512π,1024π,2048π,and4096π) and two cases of Δω(π/2andπ/4) are illustrated in Figure 5.Dependence of errors att=0 on the frequency truncation bound.
(a)
Δω=π/2(b)
Δω=π/4This 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 ωT increases. This reflects that the influence of truncating the integral limits is dominating in the error of h(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 h(0).Surprisingly, the result from the TS is even better than the results of higher-order integrations over the front portion in Figure1, 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.



---

*Source: 143582-2011-02-07.xml*