CRITERIA FOR DETERMINING THE COMPUTATIONAL E R R O R OF NUMERICALLY CALCULATED FOURIER INTEGRAL TRANSFORMS C H A R L E S J. M E S S A , W I L L I A M L. L U Y B E N , A N D GARY W. P O E H L E I N Department of Chemical Engineering, Lehigh University, Bethlehem, Pa. 18015 Commonly used criteria for determining the upper frequency limit of numerically calculated Fourier integral transforms of continuous time records are incorrect applications of the “sampled-data theorem.” This theorem assumes no knowledge of the sampled function between data points. The upper frequency limit depends only on how well the approximating function used in the numerical integration represents the real continuous function. This i s illustrated with several common functions. Complete equations are presented for trapezoidal anld parabolic numerical integration which permit multiple sampling intervals.
HE use of mathematical models of chemical processes in Tprocess control and optimization studies has led to the development of techniques for measuring the transient behavior of processes to known disturbances-e.g., pulse tests. Such information can be used to check mathematical models and to determine constants in semiempirical transfer functions. This paper discusses some of the problems associated with the analgsis of these transient responses and their traiisformation into the frequency domain. The numerical evalurition of Fourier integral transforms (FIT’S) and transfer functions (TF’s) from continuous records of a process disturbance and response always involves the question of how often niust the data recordings be sampled to achieve reasonable accuracy in the FIT’S in the upper frequency range-i.e., how many data points must be used? Sampling criteria which are commonly applied are listed in Table I. All criteria can be rispresented as some fraction of the sampled data theorem. This theorem is rigorous for sampleddata systems where one is dealing with discrete data points and nothing is known about the function between sampling intervals. T e consider the calculation of FIT’S of continuously recorded functions, the usual situation in pulse tests of chemical processes. The criteria given in Table I yield poor estimates of the number of data points rt:quired. The upper frequency limit depends only on how well the approximating function describes the continuous function; only two data points are required to calculate esactly the F I T of a square pulse a t all frequencies if a first-order approximating function is used. The same result is obtained for a triangular pulse. The response of a first-order system to an impulse-driving function is also tested s ith first- and second-order approximating functions. Development of Error Criiterion
The analysis to determine the error criterion is straightforward, using the following approach. The FIT of a function, f ( t ) , of finite duration, T,, is given by :
In general, pulse functions are not describable analytically over their entire time duration. Therefore, the above integral cannot be evaluated analytically. However, wit’hout loss of rigor me can divide the time interval 0 to T p into a for the purpose of integration. number of subintervals, This subdivision permits a valuable mathematical synthesis-namely, j ( t ) can be represented by simple functions in every subinterval Atk. Unfortunately, for most experimentally generated pulses, the formulation of an exact’ analytical espression for f ( t ) in every subinterval is prohibitive. Hence, one considers the use of a simple analytical expression (withvariable parameters) to approximate j ( t ) in every interval. Traditionally, a polynomial is used. This polynomial has variable parameters which are calculated from the values of f ( t ) in each interval: ffOk
+
fflhTk
FIT actual - FITcalculated
(1)
The error normally increases as the value of the frequency increases and the number of data points decreases. The upper frequency limit is defined as that frequency a t which a specified error is exceeded for a given sample spacing, At (the same as specifying the number of data points used).
+
ff2kTk2
+ .. . +
ffnkTkn
constants t-
(3 1
8k
k
Ati i=l
sample intervals Once the parameters are calculated for some interval, the polynomial is used in the integrating formula to calculate the value of the F I T in that interval. This process is continued for all intervals with the F I T given by: x-1
=
( 21
j(t)e+U‘dt
0
Sumerical calculation of the F I T is an approximation and has an error associated with it, defined as: 6
TP
J
FIT( j ( t ) ) =
BX+l
F I T SZ
Yk (l)e-j”‘ dt k=O
(4)
ok
The computational error (Equation 1) which results from this approximation is: E
= k=O
r+’
{ j ( t ) - pk(t)~e+-jwldt
8k
(5 1
As pointed out by Filon (1928), the error obviously depends only on how vvell pk(f) approximates f ( t ) in each VOL.
8
NO.
4
NOVEMBER
1 9 6 9
745
Table 1.
1. 2. 3. 4. 5.
Source Draper, McKay, and Lees (1953) Dreifke 11961) Hougen (1964) Lees and Dougherty (1967) Sampled datatheorem a
w
=
factor which determines this limit is only a function of how well the first-order approximating function represents f ( t ) in every interval and cannot be expressed by a simple quantitative formulation. This limit is different for every $oh ( t ) and f ( t ) combination.
Upper Frequency limit Criterion
frequency. A1
=
Criteria“ w ( A t ) 5 */2 w(At) 5 x/2 u ( A t ) 5 x/2
5 a/10
W(At)
w(At)
Some Ideal Cases
2x
data-sampling interval.
interval. Thus, no unique limiting value of the product In fact, if $oh (t) is identically equal to f ( t ) in every interval, no error results a t any frequency. To estimate the error, both an approximating function and the pulse function f ( t ) must be specified. Even after specifying these functions, a simple quantitative formulation of an upper frequency limit seems improbable. However, a qualitative criterion to determine the correct number of data points can be stated as follows: “For a given approximating function (pk ( t ) the tabulated values of the pulse function S ( t ) should be chosen so that in the interval A h , ‘ p k ( t ) represents the function f ( t ) sufficiently well.” A more quantitative and workable criterion might be that $ok ( t ) should representf(t) as accurately as the experimental data permit. Once an approximating function is chosen, the constant parameters CY%, can be determined from the sampled values of f ( t ) for the interval A t k . Equation 5 can then be integrated analytically over the interval. Summing the results of the integration as shown in Equation 4-i.e., from t = 0 to t = ?“,-results in a simplified algorithm. For a first-order polynomial approximating function this algorithm is given by Equation 6 if the sampling intervals-Le., At,-are constant. w ( A t ) can be postulated.
The accuracy of Equation 6 depends on both o and At, which could be expressed as the w ( A t ) limit. However, the Table II.
I. Time domain representation 11. Analytical expression of transfer function 111. Tabulated data points required for:
Three ideal cases were studied : square and triangular pulses, and the response of a first-order system to an impulse. All cases were studied using a first-order polynomial approximating function. The first-order system response was also studied using a second-order polynomial approximating function. Successive data points were chosen so that $ o k ( t ) represented f ( t ) as accurately as the data justified in the interval between the data points. Only two data points are required to describe the square pulse exactly in the interval 0 to T,. Hence, one would assume zero computation error. This was the observed result using the computational formula given in the Appendix (Equation ,46). Three data points (two of which are the zero end points) will describe a triangular pulse. Again, no computational error resulted from the three-point (two-interval) circulation via Equation A6. Lees and Dougherty (1967) used 100 input intervals to calculate the F I T of a square pulse and 400 intervals to calculate the F I T of a triangular pulse. A more important case is the study of a time response which cannot be represented exactly by a polynomial in every time interval. The case examined is the time response of a first-order system to an impulse disturbance. The system time response is given by:
f ( t ) = (I/T)e-t’r (7 1 The time response of the system was plotted using the analytical expression of Equation 7. The tabulated data, required in the computational algorithms, were taken from this plot using the time intervals indicated in Table 11.
Study of a First-Order System
j ( t ) = (l/r)e-f’i let r = 1 second FIT{f(t)}= l(Cjw 1) w = 100 radians/sec. T, = 7 . 0 sec.
+
No. of Data Points 446 2230 223 22 9
Criterion 1, 2, 3 4 5 This DaDer (Fopla This ,per (SOPj
IV. Comparison T F Exact
Frequency, Rad./Sec.
w
0 .oo
0.01 0.05 0.10 0.50 0.80 1.oo 2.00 3 .OO 5 .OO 10 .oo 50 .OO 100.00 500.00 900.00
--Magnitude 1.0000 1.OoOo 0,9988 0.9950 0.8944 0.7809 0.7071 0.4472 0.3162 0.1961 0.0995 0.0200 0.0100 0.0020 0.0011
Angle 0 .OO -0.57 -2.86 -5.71 -26.56 -38.66 -45.00 -63.43 -71.56 -78.69 -84.29 -88.85 -89.43 -89.89 -89.94
T F Calculated (SOP)
TF Calculated (FOP) ~
3dagnitude
1.0000 1.0000 0,9988 0.9953 0.8961 0,7809 0.7071 0.4475 0.3161 0.1967 0.0995 0.0200 0.0100 0.0020 0.0011
Angle 0 .OO -0.57 -2.84 -5.67 -26.58 -38.68 -44.92 -63.40 -71.56 -78.69 -84.15 -89.74 -89.84 -89.90 -89.94
Magnitude 1 0.9997 0.9988 0.9953 0.8964 0.7841 0.7074 0.4634 0.3444 0.1950 0.1004 0.0200 0.0100 0.0020 0.0011
.oooo
FOP. First-order polynomial approximating function. SOP. Second-order polynomial approximating function. T, arbitrarily selected at 7.0 seconds. even though f(t) is not quite zero a t this point. 746
l & E C
FUNDAMENTALS
Angle 0.00 -0.52 -2.80 -5.67 -26.04 -37.80 -44.89 -61.48 -72.04 -80.59 -85.35 -88.97 -89.49 -89.91 -89.96
Calculations were performed using both first-order and secondorder approximating polynomials. The results of the calculations are shown in Table 11. Only nine data points gave very accurate results with a second-order approximating function, and 22 data points with a first-order approximating function. The sampled-data criteria (Table I j would predict anywhere from 223 to 2330 data points. The sampled data theorem (SDT) is being incorrectly used by several authors to estimate the upper frequency limit. In particular, Lees and Dougherty state that one should use ten times the tabulated data suggested by the SDT. The SDT does not yield correct results, because even though one tabulates discrete time points for a calculation, one chooses these points given the information provided by continuous measurement of the function between the chosen points. This is distinctly different from sampled-data systems, rThere one has no information betm-een successively measured points. The upper frequency limit depends upon both the approximating function, pi: ( t ) , and the pulse function, f ( t ) . Increasing the order of the polynomial approximating function can decrease the number of tabulated data points required for a given error.
Finally the transfer function is written as:
The SDT transfer function has an upper frequency limit as shown in texts on sampled-data systems. This limit is a result of the lack of knowledge of how the function acts between the sampled points. Although Equation 14 is mathematically equivalent t o Equation 10, an important distinction exists-namely, that one uses the information of how functions f ( t ) and r ( t ) act between the tabulated points to choose the tabulated points when these functions are measured continuously. We do not have this information when the functions are measured intermittently. Hence, even though both equations yield identical numerical values for the transfer function, in the case of the sampled system, this calculated transfer function is representative of the system only up to the upper frequency limit. Beyond that limit there is no assurance that the calculated transfer function is reliable. Multiple Sampling Intervals. Equation 10 applies only when the same sampling interval At is used for both the forcing function and the response function. Referring to Figure 1
Additional Considerations
The numerical approximation of the F I T of a function depends upon the shape of the function and the approximating function. The equality of these two fuiictions is the basis for choosing At or the number of tabulated data points required. Possible Point of Confusion. One of the reasons the SDT is being misapplied may be the similarity between the transfer fuiiction commonly calculated numerically and a SDT transfer function. Using a first-order polynomial approximating function, the F I T of a function f ( t ) would be written as Equation 8, assuming f ( t ) has continuous derivatives and constant sampling intervals of length At and that f ( 0 ) = f ( T , ) = 0.
TIME
PROCESS RESPONSE
.I'-l
E f(kAt)e-JwkAtAt
FIT( f ( t ) } = CFf
-
k=O
sin &At
( 3wAt )
CF,=
~
The trZtnsfer function of a system disturbed by r ( t ) and having a response f( t j is given in Equation 9.
TF =
F I T ( f ( t ) ) - CF',xf(kAt,) -
FIT( r ( t ) ) If one chooses Aff can be written as
CF, =
exp (-jwkAtf)At,
cr (kat,) exp (--jwkAf,)At,
(9 1
0
TIME
Figure 1.
-
Typical process functions
At,, then C F , = C F , and Equation 9
The sampled function f " ( t j is given by: TJAt
6 (t
f" ( t ) =
- kAt)f(kAt j
k=O
The Laplace transform of this function is
fw
=
and the F I T is
Cf(kAt)e-ksAf TIME
FIT( f " ( t ) }
xf(kAtje-jukAt
Figure 2. VOL.
8
-
Multiple sampling intervals
NO. 4
NOVEMBER
1969
747
wider spaced data points are sometimes used to describe f ( t ) than r ( t ) , because f ( t ) changes less rapidly than r ( t ) , in a given time span. If a first-order polynomial approximating function is used, f ( t ) could be represented by straight lines over greater periods of time than r ( t ) . Hence, if the same At is used for both functions, many more tabulated points are used than needed to describef(t). When different At’s are used, the transfer function is still given by Equation 9 but C F f # CF,. The use of multiple regions of different At for a typical process response curve (Figure 2 ) should introduce no additional computational error. However, Clements and Schnelle (1963) indicate this to be the case. Apparently the usual algorithm published in the literature (first-order polynomial approximating function) is Equation 8, which is valid only when the function is zero at the end points of the region in which it is being applied. Equation 8 is a modified version of the actual computational algorithm. (Equation A5 or A6 in the Appendix is the complete algorithm.) Most pulse functions are defined as zero at time zero and at some later time when the function returns to its original value. This is not the case at the points which separate regions of different At’s. If Equation 8 is applied without regard to the additional terms when multiple At’s are used, additional computational error may be expected.
analytically to yield :
With constant At values, Equation A5 becomes:
Summary
Choice of the proper number of data points for continuous time records is based on accurately representing the function between any two successive points by the approximating function used in the FIT computational algorithm. The use of multiple regions of constant At is sometimes advantageous, but the complete form of the FIT comgutational algorithm must be used. This can be programmed as easily as the abbreviated form. We have applied these algorithms and the sampling criterion suggested for the analysis of data obtained from the dynamic testing of heat exchangers in our laboratory.
Equation A6 reduces to the commonly used, special case, trapezoidal approximation to the FIT, since i = 0 and
K ( A t ) = Tp. A similar derivation for a second-order representation of f ( t ) in the sampling intervals leads to:
Appendix, Fourier Transform of f(t), First-Order Approximation
The FIT of f ( t ) is given by TP
FIT
(A1 1
f(t)e-jotdt
=
wherej(t) = 0 for t 5 0 and t 2 T p . If the interval (0, T p ) is divided into unequal subintervals of length Atk Equation A1 can be written: N-1
FIT = k-0
\
6x+i
AV-l
j(t)e-jW
= k-0
8k
In
(A2 )
. . ., 2n
(A7)
This equation results from integration of ‘ p k ( t ) over two intervals. The as,are calculated from the end points of adjacent intervals, Equation A7 is valid for unequal sampling intervals. literature Cited
where
IL= l ~ t l j ( t ) e - ~ % t
lrL’ +
A first-order approximation of f ( t ) in each interval yields Ik E
(CYok
where rk =
allrk)e-j.jYtdt
(-43 )
t - ek
The parameters L Y O ~and fflk are computed from the values of = &+I. Equation A2 can be integrated
f ( t ) at t = Ok and t 748
k = 0, 2, 4, 6,
l&EC
FUNDAMENTALS
Clements, W.C., Schnelle, K. B., Ind. Enq. Chem. Process Design Develop. 2, 94 (1963). Draper, C. W., McKay, W., Lees, S., “Instrument Engineering,” ,IlcGraw-Hill,t (New York, 1953. Dreifke, G. E., Effect of Input Pulse Shape and Width on the Accuracy of Dynamic System Analysis from Experimental Pulse Data,” Sc.D. thesis, Washington University, St. Louis, 1961. Filon, L. S. G., Proc. Roy. Soc. Edinburgh 49, 38 (1928). Hougen, J. O., Chem. B i g . Progr. Mono. Ser., e , 60 (1964). Hougen, J . O., Walsh, R. A., Chem. Enq. Proqr. 67, 69 (March 1961). Lees, S., Dougherty, R. C., Basic Eng. 1967,445. RECEIVED for review June 28, 1968 ACCEPTED April 15, 1969
