Mathematical Problems in Engineering

Volume 2010, Article ID 143582, 9 pages

http://dx.doi.org/10.1155/2010/143582

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

^{1}College of Science, China Agricultural University, P.B. 74, East Campus, Beijing 100083, China^{2}The First Affiliated Hospital of China People's Liberation, Army General Hospital, Beijing 100037, China

Received 7 February 2010; Accepted 28 December 2010

Academic Editor: Moran Wang

Copyright © 2010 Kui Fu Chen and Yan Feng Li. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 as or the inverse Laplace transform Here the real number is the convergence abscissa, that is, all the poles lie on the left side of .

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,

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 as Denote the original function of as , then

Choice of is case-dependent. However, many systems own an asymptotic form as , where and are constants. The original function is , 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 ). Numerical approximation of the continuous integration in (1.5) is where . in (2.1) depends on the integration scheme. For the trapezoidal scheme (TS) For the Simpson’s rule (SR) For the Boole's rule(BR) (http://mathworld.wolfram.com/BoolesRule.html):

To facilitate the fast Fourier transform (FFT), the continuous in (2.1) must be sampled at discrete instants. The FFT prefers a time interval . For , we have Substituting (2.5) into (2.1) leads to where for 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 , and its transfer function is . 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 , 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 .

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 , and are completely unjustifiable. From the view of numerical integration, the integrand 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.

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 does not necessarily coincide with this assumption.

The significant error can also be appreciated along the following theme. In light of (2.1), 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,

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 . 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 . 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 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 , , , , are examined. The truncating frequency keeps a constant of , 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 , especially approaching . This is taken for granted, since the ripples, if significant, have been expelled rightward.

However, the feature around , 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 and are almost comparable at , although the former is better than the latter over most part of . This lead to the guess that decreasing cannot improve the accuracy at 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, computed with fives cases of and two cases of are illustrated in Figure 5.

This figure shows that, firstly, the error with is comparable to that . 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 , 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 .

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

#### References

- J. L. Humar, A. Bagchi, and H. Xia, “Frequency domain analysis of soil-structure interaction,”
*Computers & Structures*, vol. 66, no. 2-3, pp. 337–351, 1998. View at Google Scholar · View at Scopus - D. Polyzos, K. G. Tsepoura, and D. E. Beskos, “Transient dynamic analysis of 3-D gradient elastic solids by BEM,”
*Computers & Structures*, vol. 83, no. 10-11, pp. 783–792, 2005. View at Publisher · View at Google Scholar · View at Scopus - J. Yan, C. Zhang, and F. Jin, “A coupling procedure of FE and SBFE for soil-structure interaction in the time domain,”
*International Journal for Numerical Methods in Engineering*, vol. 59, no. 11, pp. 1453–1471, 2004. View at Publisher · View at Google Scholar · View at Scopus - N. A. Dumont and R. De Oliveira, “From frequency-dependent mass and stiffness matrices to the dynamic response of elastic systems,”
*International Journal of Solids and Structures*, vol. 38, no. 10-13, pp. 1813–1830, 2001. View at Publisher · View at Google Scholar · View at Scopus - M. C. Genes and S. Kocak, “Dynamic soil-structure interaction analysis of layered unbounded media via a coupled finite element/boundary element/scaled boundary finite element model,”
*International Journal for Numerical Methods in Engineering*, vol. 62, no. 6, pp. 798–823, 2005. View at Publisher · View at Google Scholar · View at Scopus - S. W. Liu and J. H. Huang, “Transient dynamic responses of a cracked solid subjected to in-plane loadings,”
*International Journal of Solids and Structures*, vol. 40, no. 18, pp. 4925–4940, 2003. View at Publisher · View at Google Scholar · View at Scopus - K. F. Chen and S. W. Zhang, “A semi-analytic approach to retrieve the unitary impulse response function,”
*Earthquake Engineering & Structural Dynamics*, vol. 36, no. 3, pp. 429–437, 2007. View at Publisher · View at Google Scholar · View at Scopus - D. K. Kim and C. B. Yun, “Earthquake response analysis in the time domain for 2D soil-structure systems using analytical frequency-dependent infinite elements,”
*International Journal for Numerical Methods in Engineering*, vol. 58, no. 12, pp. 1837–1855, 2003. View at Publisher · View at Google Scholar · View at Scopus - J. A. Inaudi and J. M. Kelly, “Linear hysteretic damping and the Hilbert transform,”
*Journal of Engineering Mechanics*, vol. 121, no. 5, pp. 626–632, 1995. View at Publisher · View at Google Scholar · View at Scopus - N. Makris and J. Zhang, “Time-domain viscoelastic analysis of earth structures,”
*Earthquake Engineering & Structural Dynamics*, vol. 29, no. 6, pp. 745–768, 2000. View at Publisher · View at Google Scholar · View at Scopus - H. C. Tsai and T. C. Lee, “Dynamic analysis of linear and bilinear oscillators with rate-independent damping,”
*Computers & Structures*, vol. 80, no. 2, pp. 155–164, 2002. View at Publisher · View at Google Scholar · View at Scopus - K. F. Chen and S. W. Zhang, “On the impulse response precursor of an ideal linear hysteretic damper,”
*Journal of Sound and Vibration*, vol. 312, no. 4-5, pp. 576–583, 2008. View at Publisher · View at Google Scholar · View at Scopus - N. V. Sushilov and R. S. C. Cobbold, “Wave propagation in media whose attenuation is proportional to frequency,”
*Wave Motion*, vol. 38, no. 3, pp. 207–219, 2003. View at Publisher · View at Google Scholar · View at Scopus - M. Schanz, “A boundary element formulation in time domain for viscoelastic solids,”
*Communications in Numerical Methods in Engineering*, vol. 15, no. 11, pp. 799–809, 1999. View at Google Scholar · View at Scopus - A. S. Omar and A. H. Kamel, “Frequency- and time-domain expressions for transfer functions and impulse responses related to the waveguide propagation,”
*IEE Proceedings: Microwaves, Antennas and Propagation*, vol. 151, no. 1, pp. 21–25, 2004. View at Publisher · View at Google Scholar · View at Scopus - J. Suk and E. J. Rothwell, “Transient analysis of TE-plane wave reflection from a layered medium,”
*Journal of Electromagnetic Waves and Applications*, vol. 16, no. 2, pp. 281–297, 2002. View at Google Scholar · View at Scopus - B. Gatmiri and K. van Nguyen, “Time 2D fundamental solution for saturated porous media with incompressible fluid,”
*Communications in Numerical Methods in Engineering*, vol. 21, no. 3, pp. 119–132, 2005. View at Publisher · View at Google Scholar · View at Scopus - A. Elzein, “A three-dimensional boundary element/laplace transform solution of uncoupled transient thermo-elasticity in non-homogeneour rock media,”
*Communications in Numerical Methods in Engineering*, vol. 17, no. 9, pp. 639–646, 2001. View at Publisher · View at Google Scholar · View at Scopus - K. E. Barrett and A. C. Gotts, “Finite element analysis of a compressible dynamic viscoelastic sphere using FFT,”
*Computers & Structures*, vol. 80, no. 20-21, pp. 1615–1625, 2002. View at Publisher · View at Google Scholar · View at Scopus - Y. Y. Wang, K. Y. Lam, and G. R. Liu, “Strip element method for the transient analysis of symmetric laminated plates,”
*International Journal of Solids and Structures*, vol. 38, no. 2, pp. 241–259, 2001. View at Publisher · View at Google Scholar · View at Scopus - A. Maffucci and G. Miano, “An accurate time-domain model of transmission lines with frequency-dependent parameters,”
*International Journal of Circuit Theory and Applications*, vol. 28, no. 3, pp. 263–280, 2000. View at Publisher · View at Google Scholar · View at Scopus - F. V. Filho, A. M. Claret, and F. S. Barbosa, “Frequency and time domain dynamic structural analysis: convergence and causality,”
*Computers & Structures*, vol. 80, no. 18-19, pp. 1503–1509, 2002. View at Publisher · View at Google Scholar · View at Scopus - C. Song and J. P. Wolf, “Scaled boundary finite-element method—a primer: solution procedures,”
*Computers & Structures*, vol. 78, no. 1, pp. 211–225, 2000. View at Publisher · View at Google Scholar · View at Scopus