Page I A STANDARD DEVIATION COMPUTER by Fred V. Brock Daniel J. Provine Technical Report No. 1 ORA Project 03632 NATIONAL SCIENCE FOUNDATION GRANT G-11404 and ORA Project 02515 POWER REACTOR DEVELOPIENT COMPANY Detroit, Michigan College of Engineering -Department of Engineering Mechanics Meteorological Laboratories Institute of Science & Technology Great Lakes Research Division Special Report No. 10 THE UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN January 1961

Page II I i I I

Page III TABLE OF CONTENTS Page LIST OF FIGURES v ABSTRACT vii 1. INTRODUCTION 1 2. MATHEMATICAL DESIGN CRITERIA 5 3. COMPUTING CIRCUITS 11 4. EVALUATION OF THE COMPUTER 17 A. Laboratory Tests 17 B. Field Test 17 5. SUMMARY AND CONCLUSIONS 21 ACKNOWLEDGMENTS 23 REFERENCES 25 iii

Page IV I i I I A 9 4

Page V LIST OF FIGURES Figure Page 1 Plot of x(t - T/2) - i(t - T/2) and the second-order approximation as a function of the input frequency, sampling time product. 6 2 Plot of the integration of sin DT over sampling period T with the second-order approximation. 8 5 Computing circuit for generation of 1.25 E which is an estimate of s. All resistor anid capacitor values are in units of megohms and microfarads, respectively. 12 4 Computing circuit for generation of s. All resistors and capacitors values are in units of megohms and microfarads, respectively. 1 5 Minimum amplifier circuit for computation of 1.25 E. 15 6 Plot of E as a function of CDT with experimental values indicated. 18 v

Page VII ABS TRAC T The standard deviation computer is a small electronic analog device which is designed to calculate the standard deviation of a voltage signal with less than 5% error. The output is a continuous voltage signal proportional to the standard deviation which is generated with a computation time equal to the sampling period. Mathematically, this instrument computes the mean absolute deviation which is related to the standard deviation by a constant factor. It accepts signals in the frequency range of from d-c to one kilocycle over a wide range of amplitudes and computes with sampling times up to 3600 seconds or more. In field tests, the computer operated successfully on the output of a wind vane. vii

Page 1 1. INTRODUCTION Progress in the study of atmospheric turbulence has been hindered by the lack of satisfactory means of measuring directly the significant turbulent characteristics of the wind at any given time and place. Earlier methods, such as those used by Scrase and Best,2 required analysis of the traces made by a sensitive anemometer, wind vane, or bivane. In one method the "mean fluctuation ratio,," which is the mean of the absolute values of the eddy velocity expressed as a percentage of the mean wind speed for a given period, is determined by a time-consuming analysis of wind-speed chart traces. The "gustiness" or "intensity of turbulence," which is the ratio of the root-mean-square value of the oscillations of the wind speed to the mean wind speed, is also cumbersome to evaluate. The gustiness or intensity may also be determined approximately from fluctuating wind direction traces, as shown by Best,2 but this is also time-consuming and difficult to establish on a routine basis so that the evaluation of records may be done by nonscientific personnel. The first instrument designed specifically to give a direct measurement of atmospheric turbulence was the bridled-cup turbulence integrator or bridledcup gust accelerometer developed by Hewson and Gill.-)4 This relatively simple instrument gives the total acceleration, without regard to sign, occurring in the fluctuating wind speed during any specified time interval. For example, the instrument may show a total wind acceleration of 800 m sec-2 during a period of 50 min, and this value may be obtained directly from the chart by unskilled workers without any exercise of scientific judgment on their part. The instrument thus gives a measure of a fundamental property of turbulent flow, its acceleration, and of the turbulent diffusion associated with it. Its general usefulness has been limited, however, by the fact that the theories of atmospheric turbulence and diffusion have not yet developed substantially in terms of the accelerations of turbulent flow, although a beginning in this direction has been made by Cramer.5 There has been a recent trend toward the specification of the turbulent characteristics of the atmosphere in terms of the statistical properties, such as the standard deviation, of the wind speed and direction. This approach was developed by Hewson, Gill, Cramer, and Record,6 and used by Friedman,7 Falk et al.,8 and Record and Cramer9 in field studies of atmospheric diffusion. Further development and comparison with field results have been undertaken by Hay and Pasquill;10 and Panofsky and McCormick1 have suggested a semi-empirical expression for the variation with height of the standard deviation of the vertical component of the wind velocity. A direct method of computing an estimate of the standard deviations of the wind direction from some arbitrarily fixed bearing has been presented by Jones and Pasquill.-12 The experimental system reported by them was designed to com 1

Page 2 pute the statistical mean of the standard deviation of wind-direction fluctuations when averaged over 5 seconds and sampled over 3 minutes. Study of the fundamental properties of atmospheric turbulence and diffusion would be facilitated by a similar device capable of calculating the standard deviation with less than 5% error and reliable enough for continuous operation. This accuracy requirement is imposed to match the accuracy of the sensor and so that the data need not be recomputed. The computer should provide several pre-set values of the sampling period up to 3600 seconds. Since the computer is to perform mathematical operations on the signal, it should be possible to convert it readily to other uses, such as computing the correlation function of two signals. On-site operation demands that the computer be compact and rugged with simple but effective controls. These requirements suggest the use of an electronic analog circuit which would consist of high-quality passive components, e.g., resistors and capacitors, and of low-cost, commercially available operational amplifiers. Using these components, precise mathematical operations can be performed and this type of computer will accept any direct analog voltage signal and will drive any conventional recorder or display device. 2

Page 3 2. MATHEMATICAL DESIGN CRITERIA The mean and standard deviation of a random signal, x(t), over a fixed period of time are given by = -fC x(t) dt (1) T 2 s = | f J+~ [x(t) - x]2 dt (2) where x(t) = input signal t = time T = sampling period x = mean over period T s = standard deviation over period T. Another measure of the dispersion which will be useful is the mean absolute deviation given by E = |x(t) - X| dt (3) L T where E is the mean absolute deviation over the period T. If x(t) is normally distributed, the standard deviation and mean absolute are related by s = 1.25 E (4) This relation may also hold for some non-normal distributions but if (4) does not hold, the distribution is not normal.15 The above formulae may be modified for continuous computation by changing the limits of integration to t - (T/2) and t + (T/2) instead of -T/2 and +T/2, respectively. At this point, a difficulty arises because no physical system can 5

Page 4 perform integration in future time as required by the limit t + (T/2). This difficulty can be overcome by delaying the computation by T/2 seconds so that the mean is given by t iare x(t - T/2) = 2 x(t) dt (5) t-T The expressions for the standard deviation and the mean absolute deviation are t E(t - T) 1 x(t - T/2) - Y(- - T/2)1 dt (7) T -T Equations (6) and (7) require comparison of the input with the mean of the in-1 put computed over a period symmetrical about the input so that if the mean is given by (5), the instantaneous value to be compared with it is x(t - T/2) rather than x(t). Delay of the input by T/2 seconds for large T requires expensive equipment not compatible with the demand for low initial cost. What is needed is an approximation to x(t - T/2) - x(t - r/2) which is sufficiently accurate but which can be implemented with a simple electrical analog. If the approximation can be stated as a ratio of two polynomials in the Laplace transform domain, in the form \2 n Yo(S) a + a1S + a2S +.. + aSn n (8) Y.(S) bo + blS + b2S +... + bnS where S is the transform variable, then there exists a direct electrical analog which can be built up with relatively inexpensive components. The Laplace transform of x(t - T/2) - 7(t - T/2) is Lx(t - T/2)- (t - T/2)] = Lx(t - T/2) - L j x(t) d = X(S) e - - e ) t-T -rS -*,-TS X(S - T/2) - X(S - T/2) 2 1 e X() T x(s) TS 4

Page 5 The infinite series expression for (9) is X(S - T/2) - X(S - T/2) X( S)., 1)n+l )n+2 Ln+2 n () (TS) + n=O n+2 2n+2 (n + 3)' - (n + 3) -.. - (10) which can be approximated by the ratio of two polynomials in S of the form of (8). For a second-order polynomial, the approximation is X(S - T/2) - X(S - T/2) X(S) 0.067T2 s2 0.067T2 S2 + 0.50TS + 1 (11) This was obtained by graphical fitting. The result may be seen in Fig. 1. The exact solution, needed for comparison, is more readily obtained by direct integration than by inversion of the transform given in (9). For a sinusoidal input of unity amplitude, t. x(t) = lt sin ct dt t-T x(t - T/2) - x(t - T/2) and =.[ sin 2 CT sin c (t - /2) =-.. sin n (t - T/2) - T (12) (13) The term in brackets in (13) is the amplitude of the variations about the mean and the phase shift is just T/2 seconds. Equation (11) may also be put converting to a Fourier transform. the amplitude term is given by the in the amplitude form similar to (13) by This is done by replacing S by jo. Then real part: o.o672 T2 (1 -0.0672 T2)2 + (0.50 T 12 i + 00.50WT) (~io~i (lla) 5

Page 6 I Ix91 I. -N x-x for x=sin ut eq.(13) / -- 2 nd order approximation eq. (1l) / 0 Experimental points 0.01o/ 1.. II I I....I... I L111L 0.2 1 10 100 OT Fig. 1. Plot of x(t - T/2) - X(t - T/2) and the second-ordexapproximation as a function of the input frequency, sampling time product. * * If

Page 7 The correct amplitude function for x(t - T/2) - R(t - T/2) given by the term in brackets in Eq. (13) is plotted in Fig. 1 as a function of coT along with the corresponding amplitude function for the approximation given by Eq. (lla). The second integration to be performed, the one indicated in Eqs. (6) and (7), must also be put in the form of (8). If this operation is O 1tatt T o = l t Yi dt then the Laplace transform is Y0 (S) Yi (S) -TS,, i -- e 7S (14) An approximation to (13) can be developed from the series expansion: 1 - e-rS -TS =o ( l)n (TS)n n= (n —+l n=O (n + 1); (15) The second-order approximation to (15) is 12 T2 S2 + 6 TS + 12 (16) The desired amplitude function for a unity sinusoidal input is given by the term in brackets in Eq. (12). This is plotted in Fig. 2 as a function of oT along with the real part of the Fourier transform of (16), which is 2 2) 12 ~-1/2 22 12 j (16a) The mean absolute deviation for the unity amplitude sinusoidal input can be found by combining Eqs(',(7). an4 (153) 7

Page 8 00.6 oJ 0.5 a 0.3 -I St tror ig. 2. Plot of the integration of sin ver sampling period with2 12the second-order approximation. 0.1 0.0 0.1 -- -- -- ------ - — I -- ------- 0.1 I 10 100 Fig. 2. Plot of the integration of sin CT 6ver sampling period T with the second-order approximation. 8

Page 9 E(t - T) = CT - 2 sin 1/2 WT t oT2 t-T Isin w(t - T/2)I dt. (17) The absolute value term can be expressed as a Fourier series |sin w(t - T/2) I 00 _ 2 4 z cos 2 n a(t - T/2) rt "r:n=l 4 n2 - 1 Thus E(t.,- T) t CDT - 2 sin 1/2 wT / CDT2 t-T 2 4 it 00 Z cos 2 n w(t - /2) dt n=l 4 n2 - 1 OT - 2 sin 1/2 WT 22 T2 (t - 2 1T - 2 sin 1/2 T 4 E(t - - - n )T3 00 sin n CT. cos 2 n C(t - T) n=l n(4 n2 - l) (18) This expression will be used to compare with the the computer. frequency-response data of 9

Page 10 I I I rT 1

Page 11 3. COMPUTING CIRCUITS The computing circuits required for generation of the mean absolute deviation and the standard deviation were set up on a small, general-purpose analog computer to demonstrate the feasibility of the system. This is a practical test since this computer has the same quality components and amplifiers as the proposed special purpose computer. The circuit diagrams are shown in Figs. 3 and 4. All the indicated resistors and capacitors are -l quality and the multipliers are the diode quarter-square type. The circuit diagram of Fig. 5 has been divided into two major mathematical operations by the dashed lines. The operations are the generation of |x(t - T/2) - x(t - T/2) in the top half of Fig. 3, and the subsequent integration required for computation of E at the bottom. The coefficients are determined by three circuit components, a voltage divider a, f, 6, or y followed by an RC combination as follows: 3.86 R1 Ci 7.46 R2 C1 T (19) a ( from Eq. (12), and 3.47 R3 C2 6 R4 C2 (20) 5 7 from Eq. (16). Equations (19) and (20) must be satisfied for the concept of sampling time to have meaning. Given that the computer can handle a voltage range from -100 to +100 volts, then, for optimum use of the computer, the voltage input x(t) should have extreme values of +100 volts. If this is not practical, the input resistor Rg shown in Figs. 3 and 4 provides an adjustable gain. If the voltage extremes of x(t) are, say, ~10 volts, then Rg should be 0.1 megohm which makes the effective voltage range ~100 volts since the input gain is 1/Rg. The relation between the voltage at any point in the circuit and the magnitude of the corresponding physical variable is determined by the amplitude scaling factor A which is defined by the relation X(t) = A x(t) (21) where x(t) is the physical quantity and X(t) is the corresponding voltage. This relation holds at every point in the circuit in Figs. 3 and 4. Ordinarily, the 11

Page 12 x LL1 i/ W R / 1.25E - - HIGH-GAIN d-c AMPLIFIER --- VOLTAGE DIVIDER Fig. 3. Computing circuit for generation of 1.25 E which is an estimate of s. All resistor and capacitor values are in units of megohms and microfarads, respectively. 12

Page 13 I A - R 'A S2 100 (O - ---- (X-X). 100 S Diode Multiplier -- HIGH-GAIN d -c AMPLIFIER - VOLTAGE DIVIDER Fig. 4. Computing circuit for generation of s. All resistors and capacitors values are in units of megohms and microfarads, respectively. 13

Page 14 squaring operation changes this relation but the circuit points affected have been relabeled so that (21) holds throughout. The circuit of Fig. 4, for computing the standard deviation, is similar to that of Fig. 3 but has a squaring device instead of an absolute value network and another multiplier set up to perform the square root operation at the end of the circuit. To compute the standard deviation, 8 amplifiers and 2 multipliers are needed but only 7 amplifiers are required for the mean absolute deviation. Simply on the basis of the amount of equipment required, it is least expensive to compute the mean absolute deviation. It is possible to design a circuit to perform the same functions with fewer amplifiers. A 4-amplifier circuit for computing the mean absolute deviation is shown in Fig. 5 where the transfer function of the first amplifier is eo C1 C2 R1 R2 S2 (22) ei R1 R2 C2 C3 S2 + R1(C1 + C2 + Cs)S + 1 and C1 = C3 for unity gain. To determine parameter values, Eq. (22) must be matched with Eq. (11). The transfer function of the last amplifier is eo R5 1 ei R3 R4 R5 C4 C5 S2 + X5s (R4 R5 + R3 R4 + R3 Rs)S + 1 R3 and R5 = 1.25 R3 for a gain of 1.25. Equation (23) must be matched with Eq.(16). This circuit has the disadvantage that it does not work well for values of the sampling period greater than about 10 seconds. An amplifier with an output capacity of 20 milliamp at ~100 volts could not drive the complex feedback circuits for T greater than 12 seconds. Also, the magnitude of the resistors and capacitors required is much greater than for the 7-amplifier circuit. For comparison, the 7-amplifier circuit, using the same amplifiers, can operate with a sampling period of 3600 seconds and with reasonable component values. For example, a T of 3600 seconds can be obtained with 10-megohm resistors, 10 -microfarad capacitors, and voltage divider settings between 0.09 and 0.17. 14

Page 15 0.1 -(X-X) 0.1 Ci 1.25 E C4 a HIGH-GAIN d-c AMPLIFIER Fig. 5. Minimum amplifier circuit for computation of 1.25 E. 15

Page 16 I I I 0 1 i i 0 I 4 I A

Page 17 4. EVALUATION OF THE COMPUTER Evaluation of the computer was directed toward justifying the use of the mean absolute deviation as an approximation to the standard deviation. The computer performance was evaluated with respect to Eq. (7) and an attempt was made to evaluate Eq. (4). The procedure involved both laboratory and field tests.. A. LABORATORY TESTS Since any input can be made up of a combination of sine waves of various frequencies, it is valid to study the frequency response of the computer. A sine wave generator was used as the input, and to facilitate the work, a sampling period of 2 seconds was used. The result was plotted in Fig. 6 in terms of Tr so that the data apply to any sampling period. The output of the computer was compared with the theoretical response given by Eq. (18), which is the integrated form of Eq. (7), where the input is a unity magnitude sine wave. The mean absolute deviation of a sine wave consists of a steady component plus a component which oscillates in time. This is indicated in Fig. 6 by plotting the mean or steady value and the bounds of the oscillating components. For 0.1 orT _ 100, the computer error relative to Eqs. (11) and (16) depended upon the amplitude of the output so that the error was less than 5% at the 5-volt level, less than 35 at the 10-volt level, and less than 2% at the 50-volt level. Error depending on the amplitude level is a phenomenon common to all such electrical computing devices. B. FIELD TEST As a further demonstration of the computer performance, it was given a field test. A potentiometer-type wind vane was mounted on the roof of a building. Voltages applied to the potentiometer were ~10 volts from the computer, and the potentiometer output was fed directly to the computer with a gain of 10 so that ~180 degrees of vane rotation equals ~100 volts. Using the circuits shown in Figs. 3 and 4, the standard deviation and the mean absolute deviation were computed. The sampling times used were 24 and 60 seconds. The amplitude scale factor A, was 1/1.8. This was a fairly severe test as the wind vane was in a region of mechanical turbulence, although complete revolutions of the vane were uncommon. The standard deviation was in the range of from 15 to 20 degrees. No provision was made to keep the vane from rotating the potentiometer arm across the ends 17

Page 18 1. 0.1/ 0 Experimental average E 0a1 A Experimental oscillating E / --------- 0.0 1 l l-I I I ---— I -I- I I I I 10 100 Fig. 6. Plot of E as a function of caT with experimental values indic ated.

Page 19 of the potentiometer winding which causes the amplified computer input to jump from +100 volts to -100 volts or vice versa. To the computer, this looks like part of the signal upon which it is to operate, and this causes a large perturbation in both measures of the deviation for about 3T seconds. Some of the results were recorded with a chart speed of 10 millimeters per second, which is fast enough to permit good resolution of the input signal. From these, E was hand-calculated to check the computer. This calculation, over a 45-second portion of the record where T = 24 seconds, showed an average 4% error in the mean absolute deviation. Instantaeous values of the error were as great as 8%, but this is at least partly attributable to abstracting error. It was stated above, Eq. (4), that E and s should be related by the form factor 1.25 when the input is normally distributed. Having computed both E and s, this can be readily checked. To do this, E and s/E were determined at a number of points in an interval of the record. The form factor and the mean absolute deviation were averaged for the points to eliminate random measurement errors as shown in Table I TABLE I FORM FACTOR FOR SEVERAL VALUES OF THE SAMPLING TIME Sampling Number Length of time, of record, E, Form see points sec deg Factor 24 20 20 14.4 1.21 24 20 20 14.3 1.21 60 20 40 12.7 1.23 60 20 40 14.6 1.24 19

Page 20 j A A I I i I

Page 21 5. SUMMARY AND CONCLUSIONS There were two key steps involved in this method of computing the standard deviation. The first was to find a mathematical model which is an adequate approximation to the defining model and which could be readily implemented as an electrical analog. The second was to find this electrical analog which has all the desired properties, e.g., accuracy, simplicity, and reliability. The defining model was expressed in Eqs. (4) and (7) and the approximate model was partially stated in Eqs. (11) and (16). The complete statement is that the input is operated on by (11); the absolute value is taken; this is operated on by (16); and the result is multiplied by 1.25 to obtain the standard deviation. Computer circuits have been shown (Figso 5 and 4) for obtaining the standard deviation by two methods. The first used the concept of the mean absolute deviation and its relation to the standard deviation, while the second circuit obtained the standard deviation directly with the variance as a byproduct. The second circuit has the advantage of being somewhat more exact, but it is also more expensive to implement. The decision to use the first circuit was made solely on the basis of initial cost. As another measure of economy, the circuit of Fig. 5 can be used, provided that only small values of the sampling time are required. This computer can accept any direct analog voltage signal at any amplitude level. It can attenuate or amplify that signal to its own optimum range. The upper limit in frequency of input is the frequency at which the amplifier band width interferes, so that the upper limit should be at about one kilocycle. The maximum sampling time that can be used is ultimately set by the capacity of the amplifiers to drive the computing networks and by the difficulty in avoiding leakage paths around very large resistors. With amplifiers having a forward gain of 200,000 and an output capacity of 20 milliamp at ~100 volts, the accuracy will be unimpaired with sampling periods up to 3600 seconds but will begin to fall off for larger sampling periods. The primary limit on accuracy is the quality of the mathematical approximations used. The computing circuits were implemented with 1% quality resistors and capacitors, but if more accuracy in the computer were needed, one could easily use 0.05% components. Stability of operation is an inherent feature, even though the d-c amplifiers used are not chopper-stabilized. The computer can operate unattended, with its accuracy unimpaired by amplifier drift, for periods of up to one week. 21

Page 22 The applications for this computer are numerous. In addition to computing the standard deviation for air-pollution studies, it could be used as an aid in the controlled emission of plant effluents, the study of time averaging of the standard deviation of turbulent eddies versus spatial averages, quantifying the turbulent structure of air flow around buildings, and the study of the decay of turbulence with height. 22

Page 23 ACKNOWLEDGMENTS The authors wish to express their grateful appreciation to Professors E. Wendell Hewson and Gerald C. Gill and Mr. Albert W. Stohrer and Mr. Eugene W. Bierly for their consistent interest and constructive advice. The wholehearted cooperation and assistance of Mrs., Anne C. Rivette is also appreciated. 23

Page 24 0aSw

Page 25 REFERENCES 1. Scrase, F. J., 1930, "Some characteristics of eddy motion in the atmosphere," Geophys. Mem. No. 52, 16 p. 2. Best, A. C., 1935, "Transfer of heat and momentum in the lowest layers of the atmosphere," Geophys. Mem. No 65, 66 p. 35 Hewson, E. W., and G. C. Gill, 1944, "Meteorological investigations in the Columbia River valley near Trail, B. C.," in Report Submitted to the Trail Smelter Arbitral Tribunal, U. S. Bureau of Mines Bulletin 453, 23-228. 4. Hewson, E. W., 1956, "Meteorological measurements in air pollution studies," in: C. D. Yaffe, D. H. Byers, and A. D. Hosey, eds., Encyclopedia of Instrumentation for Industrial Hygiene, Ann Arbor, Univ. of Michigan Institute of Industrial Health, pp. 521-540. 5. Cramer, H. E., 1953, "A new approach to the problem of turbulent mixing," J. Meteor., 10, 46-53. 6. 'Hewson, E. W., G. C. Gill, H. E. Cramer, and F. A. Record, 1951, Research on turbulence and diffusion of particulate matter in the lower layers of the atmosphere, Final Rep., Contract No. AF 28(099)-7, M.I.T., 80 p. 7. Friedman, D. G., 1953, "The height variation of lateral gustiness and its effect on lateral diffusion," J. Meteor., 10, 372-379. 8. Falk, L. L., W. R. Chalker, J. A. Greene, C. Bo Cave, and C. W. Thorngate, 1954, J. Air Poll. Control Assoc., 4, 87-96. 9. Record, F. A., and H. E. Cramer, 1958, "Preliminary Analysis of Project Prairie Grass Diffusion Experiments," J. Air Poll. Control Assoc., 8, 240-248. 10. Hay, Jo S., and F. Pasquill, 1959, "Diffusion from a continuous source in relation to the spectrum and scale of turbulence," Advances in Geophysics, 6, 345-365 -11. Panofsky, H. A., and R. A. McCormick, 1960, "The spectrum of vertical velocity near the surface," Quart. J. r. Meteor. Soc., 86, 495-5035 12. Jones, Jo I. P., and F. Pasquill, 1959, "An experimental system for directly recording statistics of the intensity of atmospheric turbulence," Quart. J. r. Meteor. Soc., 85, 225-236. 13. Brooks, C. E. P., and N. Carruthers, 1953, Handbook of Statistical Methods in Meteorology, London, Her Majesty's Stationery Office, 412 pp. Orc