UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Implementation of a modal filtering procedure Fraser, David Raye 1988

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1988_A7 F72.pdf [ 4.1MB ]
Metadata
JSON: 831-1.0064893.json
JSON-LD: 831-1.0064893-ld.json
RDF/XML (Pretty): 831-1.0064893-rdf.xml
RDF/JSON: 831-1.0064893-rdf.json
Turtle: 831-1.0064893-turtle.txt
N-Triples: 831-1.0064893-rdf-ntriples.txt
Original Record: 831-1.0064893-source.json
Full Text
831-1.0064893-fulltext.txt
Citation
831-1.0064893.ris

Full Text

I M P L E M E N T A T I O N OF A M O D A L F I L T E R I N G P R O C E D U R E David Raye Fraser B. Sc. H. (App. Sc.) Queen's University at Kingston, 1986 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F M A S T E R O F A P P L I E D S C I E N C E in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F E L E C T R I C A L E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A Apr i l 1988 © David Raye Fraser In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Electrical Engineering The University of Brit ish Columbia 1956 Ma in Mal l Vancouver, Canada Date: ;/ 20 , im Abstract A F O R T R A N program has been developed in order to investigate the process of modal parameter estimation and non-parametric system identification. The theory underlying the process of modal parameter estimation is reviewed and the decoupling of a M I M O system into several SISO systems is demonstrated. Modal filtering is shown to be useful in the field of non-parametric system identification and it is shown that it may also be of some use in the field of signal processing. The program is documented. It simulates the output of a n-th order system from which a smaller order subsystem can be decoupled. The modal parameters of a subsys-tem output signal and its first two derivatives and the modal parameters of a second subsystem output and its first derivative are calculated. The unit step response of the theoretical system and the subsystem are then calculated. The signals are then modal filtered to produce the periodic unit step response and the periodic unit square wave response. Finally, the discrete Fourier coefficients of the periodic unit step response are calculated. ii Table of Contents Abstract ii Acknowledgement viii 1 Introduction 1 2 The Modal Function 4 2.1 Identification of Modal Parameters 9 2.1.1 Relationship Between the Modal Parameters of System States . . 9 2.1.2 Identification of the M x Parameters 12 2.1.3 Numerical Modal Parameter Estimation 13 2.1.4 Justification of Analytical Integration Technique 16 2.2 Solving for the Modal Parameters 19 2.3 Solving for the Modal Constant 20 3 Non-Parametric System Identification 22 3.1 The Unit Step Response 22 3.2 The Periodic Unit Step Response 24 3.3 The Periodic Unit Square Wave Response 25 4 The Modal Program 28 4.1 Data Input 28 4.2 Noise Simulation . 31 4.3 Modal Parameter Calculation 32 iii 4.3.1 Setting up the Moda l System 32 4.3.2 Obtaining the Modal Parameters 35 4.3.3 Evaluating the Modal Constant 35 4.3.4 Finding the Modal Parameters for the Second and Th i rd States . 36 4.4 Calculating the System Responses 36 4.4.1 Unit Step Response 36 4.4.2 Periodic Unit Step Response 37 4.4.3 Periodic Unit Square Wave Response 37 4.4.4 Calculation of the Finite Fourier Transform Coefficients 38 5 Modal Parameter Estimation with Noise 39 5.1 Sensitivity Analysis of Modal Parameters 39 5.2 Statistical Model for the Simulation of Noise 42 5.3 Noise Reduction Through Numerical Integration 44 5.4 Further Concerns Associated with Modal Functions 45 6 Discussion and Conclusions 48 Bibliography 49 A A Trial Run of Program Modal 53 B A Method of Reducing the Bias of Noisy Estimates 69 C Listing of Program Modal 72 C . l Program Moda l Listing 72 C.2 Subroutine Vandm Listing 101 D Listings of Other Subroutines 103 iv D . l Listing of the Subroutines D S L I M P and F I N V 103 D . l . l Listing of Subroutine D S L I M P 103 D.l.2 Listing of Subroutine F I N V 109 D.2 Listing of Subroutine F F T 112 v L i s t o f T a b l e s A . l Summary of the modal parameter estimates for example 1 64 A.2 Summary of non-parametric system responses for subsystem SI 65 A.3 Summary of non-parametric system responses for Subsystem SI 65 A . 4 Summary of F F T coefficients 68 B. 5 Modal parameter estimates in the prescence of additive measurement noise; VAR=0.0002 69 B.6 Modal parameter estimates in the prescence of additive measurement noise; VAR=0.0004 71 B.7 Moda l parameter estimates in the prescence of additive measurement noise; VAR=0.0008 71 vi L i s t o f F i g u r e s 2.1 Interval, w(t), used to identify modal parameters 14 2.2 Interval, w(t), used to calculate the M 0 constant 20 4.3 Intervals, «;,-, used to evaluate the modal parameters for the two subsys-tems 33 5.4 Interval, w(t) 44 A.5 Control input, c 0 for example A 54 A.6 Intervals for example A 55 A.7 System responses for SI 66 A . 8 System responses for S2 67 B. 9 Graph indicating the linear relationship between modal parameter esti-mate biases and the noise standard deviation 70 vi i Acknowledgement First and foremost, I would like to acknowledge the invaluable help of my supervisor, Dr . E . V . Bohn , for his academic and financial aid without which this thesis would not have been possible. Also, I would like to thank the many people at U . B . C . who have made my stay here a pleasurable one. viii Chapter 1 Introduction Early efforts to develop system models relied primarily on determining the response of the system to sinusoidal and transient signals. These methods were continuous in nature and usually involved differential equations or transfer function representations of the system. As most realistic, physical systems are of a continuous nature, this was the obvious approach. As digital computers came into more widespread use and became more powerful, the trend developed to use discrete methods of system models with a corresponding decrease in the use of continuous model identification (CMI) . One of the difficulties involved with C M I was the need to know the time-derivatives of the signal, usually in the presence of noise. Whether the techniques involved the measurement or the actual calculation of the time-derivatives, they tended to only work reliably in cases where the input/output noise was small. In practice these algorithms were useful only in deterministic, low noise situations. The difficulty in expressing time-derivatives in terms of input-output data lead, in the early nineteen sixties, to the development of new methods that avoided evaluation of derivatives. One of the offsprings of this work was the 'method function' or 'modulating function' method. In the continuous single-input single-output (SISO) case the output and input are related and expressed as a weighted sum of output and input signal sets. Usually the signal sets took the form of time-derivatives or delayed time-derivatives of the output or input. Since the calculation of these signal set components involve time-derivatives 1 Chapter 1. Introduction 2 that are usually unavailable to the observer, it was desirable to perform some sort of linear operation on the signal bases so as to obtain terms that could be measured. This is the basis for most method function algorithms. Many different methods of method functions have evolved. Shinbrot [21] multi-plied the signal set by known functions and then integrated over the data sampling period. Others (Perdreauville and Goodson [15]) involved characterizing the system by partial differential equations and identifying the system parameters by the use of two dimensional method functions. One problem with method functions is the requirement for evaluation of integral terms. To overcome this obstacle, method functions of time that are the output of linear, time-invariant niters are chosen so that the definite integration terms may be measured directly from the filter output. B y constructing a chain of successive cas-caded filters all the required definite integrals can be evaluated. The Poisson moment functional ( P M F ) method (Saha and Rao [19]) is one example of such a technique. Recently, the use of Walsh functions and block pulse functions has been investigated (Rao and Palanisamy [17], Bohn [1]). This method involves transforming the signal into a Walsh series, or sum of weighted Walsh functions, similar to a Fourier series. In this case a 'persistently exciting' Walsh function input such as a square wave is used. Ideally the input is periodic, which allows the use of simple averaging to eliminate random, measurement noise from the estimated data. A n important advantage of the Walsh function approach is that the evaluation of integral terms that arise can be reduced to simple algebra (Bohn [2], [3]). In the case of a multiple-input multiple-output (MIMO) system, the problem of estimating a large number of unknown parameters arises for which most of the above algorithms start to loose their usefulness. Also there is the problem of initial conditions of the system, which not only add extra parameters to the estimation problem, but Chapter 1. Introduction 3 also leads to uncertainty in the calculation of multiple integrals. The method of 'modal functions' (Bohn and El-Shafey [7]) has been developed to alleviate these probems. The obstacle of initial conditions is avoided because the modal function representation of a signal is independent of initial conditions. It is shown in this thesis that for the case of a M I M O system, several SISO subsystems can be decoupled from the system and the parameters solved for relatively easily since the number of unknowns in the subsystems is significantly less than in the original system. It is also shown that modal function techniques can be used in the field of non-parametric system identification and signal processing. The modal representation of a signal can be used to identify the unit step response of the unknown system from sampled input-output data. Also the non-periodic unit step response can be filtered by a modal filter to produce a periodic unit step response that is an exact representation of the unit step response. Furthermore, the output to a unit square wave input to the system can be calculated from the non-periodic data. From the periodic unit step and square wave responses the fast Fourier transform ( F F T ) of the signals can be calculated exactly from non-periodic data. This eliminates the problem of windows and the associated spectral smearing that arises in standard F F T methods. In this thesis the concepts of a modal function and the derivation of the modal unit step responses are reviewed. A computer program is described that will simulate the analytic modal filter relations on an IBM PC. Emphasis in the thesis is placed on the analysis of modal function properties that are new and that are useful for parametric and non-parametric system identification. The PC program that is developed is meant to investigate the properties in a simple and convenient manner for the user. C h a p t e r 2 T h e M o d a l F u n c t i o n Through the use of modal functions it is possible to decouple a MIMO system into several SISO subsystems of smaller order. Consider a MIMO linear time-invariant system with a strictly proper transfer function and a state-space representation x = Ax + Bu, y = Cx, (2.1) where x is the state, u is the control, and y is the output. For illustrative purposes we will consider a fourth order system consisting of two subsystems. The system to be discussed is a one input-two output system with phase variable canonical form. 0 1 0 0 Xi bi -au —a 12 --«13 - a 1 4 ^2 b2 *3 0 0 0 1 X3 bs 0 0 -a-23 - a 2 4 X4 y i 1 0 0 0 0 0 1 0 The most convenient representation for system identification is to rewrite the system in input-output form. To do this we must obtain an input-output differential equation without reference to the internal system states. Start by rewriting (2.2) as = 41) = x2 + 6jU (2.3) 4 Chapter 2. The Modal Function 5 = - o n ^ i — a12x2 — al3x3 — aux4 + b2u (2-4) y[ l) = 41) = 3 4 + * 3 U (2.5) x{^ = —a23x3 — a 2 4 x 4 + b4u (2.6) where denotes the k-th derivative of y with respect to t. B y rearranging (2.3) to get x2 explicitely we obtain X2 = y[1] - ftiu. (2.7) B y performing a similar operation on (2.4) , y ( 2 ) - + anyi + a12y[^ - a u 6 x u + al3y2 + a 1 4 x 4 = b2u (2.8) results. Performing the same procedure on the third and fourth states of the system gives Xi = y 2 1 } - b3u (2.9) Therefore, we obtain l4 2 ) - 6 3 u ( 1 ) + a23y2 + a24y^ - a2ib3u = b4u. (2.10) Substituting (2.9) into (2.8) to eliminate the x4 term gives y j 2 ) - & ! U ( 1 ) + anyi + a12y[^ - anbiu + al3y2 + auy^ - a 1 4 6 3 u = b2u. (2.11) Note that we have now written expressions that relate the system outputs to the system inputs without using any of the internal system states. From (2.11) and (2.10) we want to write the system (2.2) in what is known as the 'input-output' form. To do this we must group the inputs on the right hand side of an equation and the outputs on the left hand side. This gives {y[2) + any? + anS/i) + [a^y^ + a13y2) = 6 x u ( 1 ) + {anbi + b2 + a 1 4 6 3 )u y{2] + a24y {2 ] + a23y2 = b3uw + (a 2 4 6 3 + h)u (2.12) Chapter 2. The Modal Function 6 If we express (2.12) as a matr ix equation we obtain the input-output form for the system as follows, d n ( y i ) + du{y2) _ Fi{u) — l^-lcJJ 0 + d22{y2) J [ F2{u) where dn(yi) di 2 (y 2 ) d22(y2) F2{u) (2) , (1) , (!) , (2) , (1) , y2 + °24y2 + °23y2 6 X « ( 1 ) + (au&x + 62 + aiib3)u 6 3 « ( 1 ) + (a24&3 + bi)u. (2.14) In the case where y\(t) and y2(t) are completely decoupled, dl2(y2) = 0, and if yi{t) and y2(t) should have a l l modes in common then rf2i(yi) would be non-zero. In practice, the input-output representation of a system is the most convenient way for dealing w i t h modal functions. Also note that by taking Laplace transforms of (2.13) a transfer function representation of the system is obtained. To obtain the characteristic equation of the system we observe that for the zero input case, dn(yi) + <*i2(y2) o + d 2 2 (y 2 ) s + a12s + au aus + ai3 0 s2 + a2is + a23 Yi{s) Let us define the following notation; £{dij(yk)} = di](s)Yk(s). = 0. (2.15) (2.16) Chapter 2. The Modal Function 7 Now we can rewrite (2.15) more compactly as dn{s) dn{s) 0 d22(s) Yi(s) Y2(s) (2.17) det dn{s)d22(s) = 0. (2.18) The characteristic equation of the system is then given by dn{s) d12{s) 0 d22 (s) In the i l lustrative example being used here we obtain the following characteristic equation. 0 = dn{s)d22(s) = (s2 + a12s + an)(s2 + a2is + a2S) s4 + dss3 + d2s2 + dis + d0 (2.19) More generally the characteristic equation of the system is as follows: sn + a ^ i s " - 1 + • • • + do = 0 (2.20) and yields the system poles Sk, k = 1 , . . . , n. For simplicity, it is assumed that the poles are distinct. The modal output function for the system is defined to be M(y(t), n) = y(t) + miy(£ - td) -j (- mny(t - ntd), where are the modal parameters, td is the modal time delay, and (2.21) zk = e'kt'1k = 1 , . . . , are the zeros of the modal characteristic equation, (2.22) c(z, n) = zn + mxzn~Y H V mn = 0 (2.23) Chapter 2. The Modal Function 8 For distinct poles, the output response for zero input has the form yQ{t) = CleSlt + ••• + cnts"1. (2.24) Since c(zk,n) —0, it is seen that M{y0{t),n) =0. (2.25) To obtain persistent excitation, general piecewise constant inputs are chosen and (2.25) is used to eliminate unknown initial conditions. The unit-step response vector gT(t) = [ffi(£)> 02(01 is denned as the solution of (2.13) with input u = 1, t > 0, satisfying the initial conditions g{0) =0,^(0) =0 (2.26) Derivatives of the unit-step on the right-hand side of (2.12) result in the unit-impulse function 6(t), where 6(0) — 0, 8(0+) = 0. Initial conditions for t = 0 + are given by g{0+) = 0 (2.27) ff^V) = h (2.28) ffl'V) = bi-auh . (2.29) When t > 0+, the unit step response function has the following representation gi{t) = gi{oo) + cneait + ••• + cineSnt. (2.30) where and where (OO) = a i l b l + b* + Q l 4 6 3 ~ a13g2(°°) ^ 9i an , x cbubz + 64 , o 2(oo) = . (2.32) <*23 Chapter 2. The Modal Function 9 W h e n t > riitd, it follows from (2.23) that for each subsystem M{gi{t),ni) = gi{t) + ma</,(i - td) + • • • + mingi{t - n,i d) = M ; »0 (2.33) where (see 2.23) Mi0 = c ( l ,n , )y , (oo) . (2.34) 2.1 Identification of Modal Parameters 2.1.1 Relationship Between the Modal Parameters of System States A n analysis in this section is given of modal parameter estimation from input-output data and the use of modal functions for estimating output derivatives. Let z(t) be defined as the modal vector and let y{t) = CleSlt + ••• + c n e f l " f = ez{t) where e = [ 1 , 1 , . . . , 1], and / ( t ) = [ d ^ , . . . , ^ " ' ] . The modal output equation for the zero-input response (2.35) is defined by y(t) + MiYit) = 0 where M i = [mi , m 2 , . . . , ran], and YT{t) = [ y ( < - r i ) , . . . , y ( i - r „ ) ] . where rx < r2 < • • • < rn are modal time-shift parameters. If we define (2.35) (2.36) (2.37) (2.38) (2.39) (2.40) HT2 «-S2T2 (2.41) Chapter 2. The Modal Function then Vmz(t) = Y(t), and from (2.38) it follows that e + M1Vm = 0. We can s imilar ly define the modal equation for output derivatives: y^(t)+Mi+1Y{t)=0. F r o m (2.35) it is seen that yW(t) = [s\,...,s in}z(t) = eD t(sk)z(t), where D(sk) s i 0 0 52 0 0 0 0 • • • sn Using the definition of Vm we can rewrite (2.43) as eD i{sk)+Mi+1Vm = 0. Define the following Vandermonde matr ix , Va, as 1 1 • • 1 _ -Si 52 ••• 5n V a — '1 "2 and define e, = [ 0 , . . . , 0 , 1 , 0 , . . . ,0], w i t h the 1 in the i-th posit ion. Using the above definition it is seen that, eD l(sk) - ei+1Va, Chapter 2. The Modal Function 11 and (2.46) can be written as ei+1Va + Mi+1Vm = 0. (2.50) Note that eD {{sk) = ei+1Va = [s[,.. .,s\]. Summarizing the above for y,y^ l\y^2\ we can see that MxVm = -e M2Vm = -eD(sk) = MxVmD{sk) M3Vm = -eD2(sk) = MxVmD\sk). Now we can state that M2Vm = MxVmD(sk) M2 = MtVnDMV-1 M2 = M~A. Furthermore, MsVm = MxVmD\sk) Ms = MyVmD2(sk)V-1 M 3 = M2VmD{Sk)V~ l M3 — M2A. where m ' A = VmD{sk)Vn It is seen from (2.56) and (2.60) that modal output derivative parameters Mk, can be determined from M i and the system matrix A. The analytic expressions and (2.61) are useful for the numerical evaluation of Mk when the system poles are known. Chapter 2. The Modal Function 12 2.1.2 Identification of the Mi Parameters For piecewise constant inputs Uj over a data interval Tj{tj < t < tj+i} it is seen from (2.33) that y(t) + MiY(t) = M0Uj,tj + ntd<t< ti+1 (2.62) To identify the M i parameters we note that equation (2.62) yields (wy(«),y(0> + M1(wj{t),Y{t)) = M o(u; ;(0,u ;(t)). (2.63) where (u>j(t),y(t)) represents the inner product of the two functions and w(t) is a suitable function denned on Tj. Now, choose Tj,Wj such that {wj{t),Uj{t))=0. (2.64) To evaluate the left hand side terms of (2.63) we must evaluate the following {wj{t),e !t} = f'1 eetdt - f'2 estdt + • • • + f'" eHdt 1 N = - E « i * e " t j f c (2-65) 5 k=i Using (2.37) and (2.65) yields N = D-\8k)Zj. (2.66) The data vector (2.66) is useful for the analysis of the identification procedure. Using (2.66) it is seen that the processed output data has the representation K-(*),y(0> = eiT^wMzit)) = [tD-\ak)Z,]y (2.67) {wj{t),Y{t)) = \vmD-\ak)Z,] . (2.68) Chapter 2. The Modal Function 13 Substituting (2.67) and (2.68) into (2.63) and choosing n data intervals gives [eD-^sjz] + M1[vmD-1{sk)z]=0 (2.69) where Z = [Zu Z2,..., Zn) (2.70) Matr ix inversion gives [eD- l{sk)z] \vmD- l(sk)z]~X + Mi = 0 (2.71) or, eV-i + Mi = 0 Mi = -eV~\ (2.72) The result given by (2.72) is the defining equation for Mi (see 2.51). 2.1.3 N u m e r i c a l M o d a l P a r a m e t e r E s t i m a t i o n For an n — th order system we must identify n modal parameters as well as the M 0 constant in (2.34). Given the definitions of a modal output function for output yi{i) yi(t) + m n y 1 ( i - * * ) + • • • + m l n y i ( i - ntd) = M10Ui, (2.73) we can take the inner product of the above equation with respect to some interval, / , , and obtain (w,y(t)) + mn(w,y{t - td)) H h min(w,y(t - ntd)) = (w,M10Ui) (2.74) where /+0O w{t)f{t)dt. (2.75) -oo Chapter 2. The Modal Function 14 1 1 -1 a 1 1 ( t i n e Figure 2.1: Interval, w(t), used to identify modal parameters. In order to identify the n parameters we choose the interval , I, given by the function w(t), as (see figure 2.1) w(t) = 1 , atd < t < btd - 1 ,btd <t < ctd 0 , otherwise w here b = a + c and where td is a suitable sampling interval. Evaluating (w(t),y(t)} we obtain: (w{t),y{t)) = {w{t),c0 + CleSit + . . . + C n e»» ' ) = (w{t),c0) + {w{t),c1e~'' t) + ----r(w{t),cne f" t) = co<u/(0,l) + ^(^(t)^"' ) + •• • + c„(«;(0,e'"f) (2.76) (2.77) (2.78) Chapter 2. The Modal Function 15 The inner product terms of (2.78) can be evaluated analytical ly by {w{t),eSit} = / eSitdt- / eSitdt J atd J bt,i = I(^|M-_I( e'.-«| r f' (2.79) Si V / _ A f2e S< bt<l — e s i a t d _ e s i c t d \ Si^ )' Clearly, ( w ( t ) , c ) = 0 (2.80) where c is any constant. Therefore (2.74) reduces to {w{t),yi(t)) + mn(w(t), yi{t - td)) + ••• + mln{w{t),yi(t - ntd)) = 0 (2.81) In order to identify the n modal parameters we must evaluate (2.81) over n distinct intervals, Wi(t), that satisfy (2.76), and choose a different value for u , , in (2.73), to simulate the persistant excitation of the system modes. To identify and y^2\ the following equations are used: y[1](t) + m2lyl{t) + ••• + m2(n_1)yi(t - (n - l)td) = M20Ui, (2.82) y l 2 ) ( i ) + msiy i ( t ) + • • • + m S ( „_ i )y i ( t - (n - l)td) = M 3 0 u , - . (2.83) B y taking the inner products of (2.82) and (2.83) over an interval Wi as in (2.76), the following is obtained {w{t),yP{t)) + m21(w{t),yi{t)) + ••• + m 2 ( n _ 1 ) ( « ; ( i ) , y i ( i - (n - l)td)) = 0, (2.84) {w{t),y{2){t)) + m31{w{t),yi{t)) + ••• + m 3 ( n_i)(u;(0,yi(* - (n - l)td)) = 0. (2.85) Since y[ 1 } (t) = C l 5 i e 3 1 ' + • • • + c n 5 n e s " J (2.86) Chapter 2. The Modal Function 16 and y i 2 ) ( i ) =Clsle°Lt + ----rCnSy»t (2.87) it follows that {w{t),y[ l){t)) = C l 6 ! ( u ; ( t ) , e S l < ) + • • • + cnsn(w{t), e3"1), (2.88) (w(t),y[2\t)) =clSl(w(t),e^) + -.- + cns2n(w(t),e^). (2.89) It is also necessary to calculate the {w(t),yi(t — jtd)) terms. Th i s is done by the following equation (w(t),yi(t - jtd)) = Cl(w{t), e3^-^) + ••• + cn(w(t),es'^-^) = cie-3litd{w{t),eSit) + ••• + cne-a" jtd{w{t),e3"t) (2.90) As in the case of the first state of the system, the modal parameters for yW and yW states can be identified by evaluating (2.84) and (2.85) over n distinct intervals. 2.1.4 Justification of Analytical Integration Technique In practice we cannot analytical ly evaluate the {w(t), y(t)) terms of (2.74) but must use some method of numerical integration. However for the computer s imulation of (2.74) a process of analytical integration has been used. It is neccessary to justify the validity of this choice. The theorem below shows that given a signal that may be expressed as a linear combinat ion of exponential terms, there exists an exact analytic relationship between the result obtained through analytic integration and through a trapezoidal approximation to the integration. Theorem 2.1.1 If a signal, y(t), can be expressed as y(t) = CleSlt + c 2 e S 2 t + • • • + c n e 5 » f , Chapter 2. The Modal Function 17 and defining "+00 / w{t)y{t)dt. (2.91) -oo then, if (Tj,y(t)) and (Ij,y{t)) refer to evaluating (2.91) by trapezoidal approximation and by analytic integration respectively, the following relationship holds: (Tj,?*) = F ( S ) ( J ; - , e * < ) , (2.92) sA esA + 1 2~ ' e a A - 1' where F(s) = — . ( 2 . 9 3 ) Proof: Given that y{t) = C l e a i t + c 2 e S 2 t + • • • + c „ e s " t , we can write (w{t),y{t)) = Cl{w{t),e3lt) + c2(w{t), e3*) + • • • + cn(w{t), e a"<). (2.94) Therefore to evaluate (2.94)  w e must only evaluate {w{t),e3it),i = l , . . . , n . (2.95) Using analytic integration, and assuming for now that (Ij,w(t)) = l,tj <t< tj+1, <Ii>e") = / ' + > <ft = — s < y + 1 _ e 3 t ' ( e s A r > A - 1) (2.96) where tj+\ — tj = NjA. Evaluating the same expression in (2.95) but using a trapezoidal approximation, of N partitions, we get {Tj,e3t) = ^ [e3'1 + 2e 3^ +A) + • • • + 2e 3^ N~ 1^  + c"''+'] = ^ [ l + 2(o + --- + oN-1)+oN], (2.97) where o = esA and t ; + 1 = tj + NA Chapter 2. The Modal Function 18 Now if we let x = o + a 2 + • • • + oN~1 then ox a2 + --- + oN~1+oN and subtracting the two equations, we get x ( l - a ) = o - o N o - a N x = 1 - a Substituting (2.98) into (2.97) we get Therefore 1 + g » + » ( - < r » ) 1 - a est' A 2(1 So) V + °N ~ ° ~ °N+1 + 2(7 ~ 2<jN] c"'' A{1 - v){l - o») 2(1-a) 2 o - 1 L J A e s A + 1 2 e s A  2 C . A _ ! ^ ' C / (2.98) (2.99) (2.100) <r,-,e"> = F ( « ) ( I , . , e r f ) sA esA + 1 where F(s) = 2 e a A - 1 (2.101) (2.102) and we can easily relate the approximate integration of (2.94) to the exact analytic integration. I Chapter 2. The Modal Function 19 Theorem 2.1.1 shows that it is valid to use the exact analytic integration formulas for the PC Moda l program since we can easily relate this to the result that we would get using the trapezoidal approximation that would be used in practice if we were estimating the modal parameters of an actual physical system using sampled output data. 2.2 Solving for the Modal Parameters Once the inner product modal equations (2.81), (2.84), and (2.85), have been evalu-ated over n intervals the following system of equations can be set up and the modal parameters solved for: (™i,y{t ( f i , y ( < -2**)) • • («>i,y(* - ntd)) m n -(wi , y ( 0 ) - * - ) > - 2td)> • • («>a,y(t - ntd)) mu = - ( « > 2 , y ( 0 ) (™n,y{t - * - ) > (wn,y{t -2*«)> • • {wn,y{t - ntd)) _ - < W n , y ( 0 ) (2.103) and (u/i ,y( t )> (wi,y(t - *<)> • • {wuy{t- (n - mn - ( w i , y ( * >(0> ' <^2,y(0) (w2,y(t - *-)> • • {w2,y(t- (n ~ ! ) « - ) > mi2 = - ( ^ 2 , y( ' ]{t)) ( w „ , y ( « ) > (wn,y{t - *-)> • • {u>n,y{t- (n " ! ) « - ) > . _ ~ ( w „ , y ( ]{t)) . (2.104) where i — 2, 3. To solve (2.103) and (2.104) Gaussian elimination is used. Chapter 2. The Modal Function 20 i 1 w(t) a b tine Figure 2.2: Interval, w( t ) , used to calculate the M 0 constant. 2.3 Solving for the Modal Constant To solve for M j o , a new interval, w is defined of the following form (see figure 2.2) 1 ,atd <t < btd w{t) = 0 , otherwise (2.105) Revaluating (2.74) we see that (w(t),Mw) = / Mwdt = (6 - a)tdMw Jatd (2.106) and that {w{t),e") = / estdt Jat,i (2.107) Therefore, from (2.74) and (2.79) (vj{t),yi{t))+mn(w{t),y1{t-td)) + --- + mln(w{t),yl{t-ntd)} = (b-a)tdMl0 (2.108) or M ™ = 7 7 i — T (M0>yi(0) + m n (M0>y i ( ' - '<*)) + • • • + min(w(<) ,yi ( ' - ntd)}) td(o — a) (2.109) Chapter 2. The Modal Function 21 where ( iu(t) ,yi(t - itd)),i = 0 , 1 , . . . ,n is given by (2.90) and (2.107). Similarly, M 2 0 = q ) ((w{t),y[1}{t)) + m21(w{t),yi{t)) + ••• + m2n{w{t),yi{t - (n - l)t d)>) (2.110) and M » = 777~ 7 ( M O . ^ O ) + m31(w{t),yi{t)) + • • • + m 3 n (u ; (<) ,y i ( ' - (n - 1)^))) td\o — a)  v 7 (2.111) where (w(t),y'(t)) is given by (2.88) and (2.90) and {w(t),y"(t)) is given by (2.89) and (2.90). The above process w i l l identify the modal parameters of a system for y, and yW. It has been shown that through the use of modal parameters a signal may be repre-sented by a continuous-time A R M A type of equation independent of in i t ia l conditions. For a M I M O system, each output can be represented separately, which decreases the number of parameters that must be estimated and decreases the time required to obtain accurate estimates. The parameters can be identified from input-output data or from system poles, and the knowledge of the system poles enables us to obtain modal function representations of the time-derivatives of the signal w i t h very little extra computational effort. The estimates of the parameters can be easily carried out using simple matr ix algebra which is easily implemented on a digital computer. Chapter 3 Non-Parametric System Identification Non-parametric system identification can be made in the time-domain or in the fre-quency domain. In the time-domain this usually involves identification of the impulse response or the unit-step response. In the frequency domain this involves determining the frequency response of the system (Wellstead [26], Young [27]). Through the use of modal parameters and the modal output function it is possible to obtain the non-parametric unit step response directly from input-output data. It is also possible to convert non-periodic data to periodic data and determine the frequency response of the system. This may have applications in the field of spectral analysis. 3.1 The Unit Step Response To simplify the notation let g(t) represent one of the elements of the unit-step response matrix . If the system input, u, is piecewise continuous over a sequence of time intervals of length Tk > ntd then the output response y{t), associated with g{t), can be given , for t > T _ i as where u(k) is the input over the interval 0 + < t — (2\ + • • • + Tk) < T^+i, Vo = y(T_\), Au(k) = u(k) — u(k — 1), and v(t) = y0 + g(t + r_i)tt(0). Note that yQ represents the initial conditions of the system. y(0 y0 + g{t + r_i)u(0) + (/(i)Au(l) + • • • , = v[t) +g{t)Au{l) + (3.112) 22 Chapter 3. Non-Parametric System Identification 23 The modal function for y(t) is given by M{y{t),n) = M{y0(t),n)+M{g(t + T-1),n)u{0)+M{g{t),n)Au{l) + ---Mou(O) + M(g(t),n)Au{l) ,t>0+ (3.113) M(v{t),n) ,t<0 where M{v(t),n) = M 0 u(0), t > 0 +. This can be generalized for the k-th interval and for t > 0 + as M ( y ( * + r i + --- + r f c _ 1 ) ,n ) = MQu{k - 1) + M{g{t),n)Au{k). (3.114) Rewriting M(y( t ) ,n ) as Y(t,n) we obtain F((* + Ti + • •• + Tk-i),n) = M0u[k - 1) + G ( « , n ) A u ( f c ) . (3.115) B y using the simplified case for the first interval, k = 1, the following is obtained, Y(t,n) = M0u(0) + G(t,n)Au{l). (3.116) Solving (3.116) for G(t,n) we obtain, G ( t ' n ) - A ^ l j ' ( 3 ' 1 1 7 ) where for the unit step response A « ( l ) = 1. From the definition of the modal function: G{t,n) - g(t) + rmgit -td) + --- + mng{t - ntd), (3.118) solving for g(t) yields g{t) = G{t,n) - ml9{t -td) mng(t - ntd). (3.119) It is seen from (3.117) and (3.119) that, since g(t) = 0, t < 0, g(t) can be recursively evaluated for t > 0 given only input-output data. Chapter 3. Non-Parametric System Identification 2 4 The modal function G(t,n) can be identified by use of (3.117). The non-parametric identification of g(t) is then obtained by recursive use of (3.119). It is seen from (2.33) that G(t, n) = Mo when t > ntd. It is of interest to investigate the properties of G(t, n) in the transition interval 0 < t < ntd. Consider the generalization of (2.14) to an n-th order differential equation of the form d(y(t)) = F(u(t)) having a unit-step input solution g(t). It follows from the definition of the unit-step response that d{G{t,n)) =F{a{t)), (3.120) where a(t) is piecewise constant in the transition interval. 1 ,0+<t<td l + m i ,td<t<2td a(t) = < ( 3 . 1 2 1 ) 1 + m i -I H mn ,ntd <t 3.2 T h e P e r i o d i c U n i t S t e p R e s p o n s e The piecewise constant function a(t) and the solution of (3.120) can be expressed in terms of Walsh functions. In using Walsh functions it is convenient to use a normalized time x — such that the data interval corresponds to — | < x < | . For notational convenience, o = | is used. Since G(ntd,n) — M0, G(0, n) = 0 it is possible to define Gp(x,n), the periodic modal unit step response on the interval —a < x < o by Gp(x,n) = -oM0 + G(x,n), = - 0 . 5 M 0 -t- G(x,n). (3.122) The function Gp(x, n) satisfies the periodicity condition Gp{x,n) = - G p ( ( x - o),n) = -Gp{{x + o),n) (3.123) Chapter 3. Non-Parametric System Identification 25 for an odd symmetric function of period 2o — 1. It is seen from (3.117) and (3.122) that the periodic function Gp(x,n) can be iden-tified from non-periodic output data. It is shown in the next section that the system frequency response can be obtained from Gp(x,n). This is an important result since it allows the use of non-periodic data in determining the frequency response. Classi-cal spectral analysis methods require window functions and the use of approximations when non-periodic data is sampled to identify spectral components. For analysis and computer simulation the equation Gp(x, n) = -oM0 + g{t) + mig(t - td) + • • • + mng{t - ntd) (3.124) obtained by substituting (3.118) into (3.122) is useful. 3.3 T h e P e r i o d i c U n i t S q u a r e W a v e R e s p o n s e Traditional spectral analyses methods use periodic inputs and outputs to identify the system transfer function. In this case, by averaging over each cycle of the input-output , we can reduce the effect of measurement noise. However, for efficient parameter estimation we usually desire a persistently excited input and in most cases, this is a non-periodic signal. In the previous section it was shown that a periodic modal function, Gp(x,n) can be identified from non-periodic data. This implies that there should be a relationship between Gp(x, n) and H(ju), the transfer function. To obtain this relationship we consider a periodic unit square wave input and let R(x, n) represent the periodic output response. Since the system input is the Walsh function sal{l,x), we use the notation R(x,n) = RSl(x,n). The following theorem is proven in [4]. Chapter 3. Non-Parametric System Identification 26 T h e o r e m 3.3.1 Let (n + l)Td — o = 0.5. The modal function Gp(x,n) and response function RSl (x, n) satisfy the modal equation Gp(x,n) = a (R3l(x,n) + miRSl(x - r,n) H h mnRSl(x - nr,n)) (3.125) and the matrix equation L[Rtl{x,n)] = L[Gp{x,n)\C-1 where L[Gp{x,n)} = [Gp(x,n) Gp(x - r,n) • • • Gp(x - nr,n)} (3.126) C where c 0 = o, Cjt = oirik, k = 1, Proof: It follows from (8.120), (3.121), and (8.119) that D(Gp(x,n)) = N{a(x)), where a(x) satisfies the modal equation a(x) = o [sa/(l, x) + misal(l, x — r) + • • • + mnsal(l, x — nr)}. (3.128) It is seen from (3.128) that the input is a sequence of sal(l,x — jr), j = 0,... ,n, terms. Equation (8.125) follows from the definition of Rgi(x,n) and by noting that Gp(x,n) and sal(l,x) satisfy condition (3.122). The existance proof for C~l follows by contradiction. If det(C) = 0, then constants ^k, k — 1,... n, exist such that F(x) = RSi (x, n) + -)!RSi (x - r, n) + • • • + ^nRSl[x - nr, n) = 0. (3.129) Equation (8.129), however, contradicts the fact that since D[RSl(x,n)} = sal(l,x) ^ 0, then D{F{x)) ^ 0. i r \ ' / r \ ' / f \ ' / j \ / Co ~c„Ci c 0 C l -c 2 (3.127) ^ c n c n - l ' ' ' Co j ...,n. Chapter 3. Non-Parametric System Identification 27 It follows from a Fourier series input-output representation that the transfer function matrix is given by H{juk) = <77T2fc r R3l{x,n)e j2*kxdx (3.130) Jo where wk = 2nk. Since R3l(x — a) = —R3l(x) and since (n + l)rd = o, it follows from (3.125) that 1 rax m 2 -m„ 1 mi mn m„_i Rsi{x) = 2 Gp(x,n) Gp(x - Td,n) -mi — m 2 • • • —mn 1 RH(x — nrj) Gp(x — nrd,n) The matrix in (3.131) is a Toeplitz matrix. Solving for RSl(x) gives RSl{x) = ioGp(x,n) + iiGv(x — rd,n) H h ^ nGp(x - nrd,n), where . (3.131) (3.132) which can be used in (3.130) to determine the frequency response. Once again the periodicity conditions that R3l(x,n) = -R3l(x - o,n) = -R3l(x + o,n) (3.133) holds as in the case of Gp(x,n) I From (3.130) and (3.131) we can see that H(juk) can now be identified from non-periodic data without having to resort to the use of window functions and their resulting spectral smearing errors. For the modal function approach, no assumptions need be placed on the normalized output data, y(x). Chapter 4 The Modal Program The program Moda l is written in F O R T R A N and is stored on floppy diskette for use with an IBM P. C. or any compatible micro computer. The program simulates a system that can be decoupled into a subsystem of order Ni and a subsystem of order N2, N2 < N\. Any or all of the poles of the two sub-systems may be common to both. This is to simulate mode coupling between the two subsystems. The modal parameters of the first 3 states ,yi[t), y[(t), and y"{t), of the first sub-system as well as the modal parameters of the first two states ,2/2(') and y'2(t) , of the second subsystem are calculated. The modal parameters for yi(t) and y2(<) are then used to calculate the unit step response, g(t), of the system and the subsystem. From that, the periodic unit step response and the periodic square wave response are calculated. Finally the finite Fourier coefficients are calculated for the periodic unit step response. 4.1 Data Input The input that the user must supply simulates the noise-free output of a system and its subsystem in the form of (2.24) y{t) = c0 + cxeSit + ••• + c „ e 3 " < 28 Chapter 4. The Modal Program 29 This is justified in that the output samples of any system may be processed using Prony's method to yield the c coefficients and the s exponents of (2.24). It should be noted that complex data entries must be entered in the form The user is first prompted for the order of the first subsystem , then for the order of the second subsystem, and finally for the number of modes of subsystem S i that are coupled with subsystem S2. After that the uncoupled coefficients of yi(t) and their corresponding complex exponential values are entered from the keyboard. The coeffi-cients and exponential terms of subsystem S2 are then entered. For S2, the coefficients and exponents for the coupled modes must be entered before those of the uncoupled modes. For the coupled modes of the subsystems we can look at the coupled modes of S i as driving inputs to the subsystem. The coefficients for the coupled modes of S i , cc, are therefore related to the corresponding coefficients of S2 by a constant term. If SI has Ni modes, of which Ns are states, and Nc are modes coupled with S2, then For S2, all JV2 modes are states. (For S2, Ns = N2.) We can now rewrite the output yi(t) to separate the states from the coupled modes as (4.134) N^Ns+Nc- (4.135) Vi{t) = c 1 0 + c u e a i l t + • • • + c 1 J V s e" ws« + yc{t), (4.136) where ye{t) ^ cCleXlt + • • • + ccNoeXNc\ (4.137) Chapter 4. The Modal Program 30 In (4.137), Aj and cc< are the exponents and coefficients of the coupled modes. As before, we can express 2/2 (') as y2(0 = c 2 0 + c 2 i e S 2 1 < + • • • + c 2 J V 9e' 9 W»*. (4.138) In this case, the following modes are coupled: 52(1) = Sl{Ns + 1) 52(2) =Sl(Ns + 2) (4.139) S2{NC) = S1{NS + Nc) = 5 l (JV a ) The relationship between the coefficients of the coupled modes can now be deter-mined, ( 4 l 4 0 ) where d\2 and d 2 2 a r e from (2.13) and are given by M A ; ) = ( A y - S i i ) ( A , ~su)---{\3:-s1Ns) (4.141) d«(Ay) = rjv 0 + 1 ( A , - - r i ) • • • (Ay - r ^ ) . (4.142) The ry values must also be entered from the keyboard when prompted for. In order to evaluate the modal parameters of an n-th order system the system must be persistently excited over n intervals. To simulate this, the user must specify Ni values of c 0 , the steady-state output in (2.24). This has the effect of simulating a piecewise constant input to the system. For the estimation of the modal parameters the values of CQ have little effect. However, the calculation of the system responses use the c 0 value for the first interval. Therefore, to obtain realistic 'unit' input responses, c 0 , should be of order 1. As well, the duration of that input must be specified for each subsystem. This is done by choosing the duration of a half cycle of the interval w in (2.74). Finally, the modal delay td , must be specified. Chapter 4. The Modal Program 31 4 . 2 Noise Simulation In order to simulate the presence of measurement noise in the system a pseudo-random Gaussian noise is produced and added to the noise-free output (4.136, 4.138). The following algorithm [8] was used to simulate the noise (~N(0,1)) and is implemented in the subroutine Random. • Choose Ro (between 0 and 1) • i?i = 9821/2 + 0.211327 • Ri = fractional part of i? i • y = yJ-2logR1sm(2nRl) • w=z + iy • RQ = Ri The value W represents a complex random noise. The user is prompted for the noise variance and the initializing value of Ro, thus producing Gaussian noise of any variance by the following formula. If X ~ N (0,1) and Y~N{0,o2) (4.143) then y - xo. Chapter 4. The Modal Program 32 4.3 Modal Parameter Calculation 4.3.1 Setting up the Modal System The modal parameters are calculated as described in section 2.1. A loop is set up to calculate the rows of the systems of (2.103) and (2.104). The number of rows equals the order of the system or the subsystem, that is Ni or N2. The (wi,es>'),i — 1,...,N terms are calculated by the function Ex . Though some of these values may be complex, once they are scaled and summed to produce (wi(t),yW(t)),i = 1 , . . . ,N,j = 0,1,2 the result is real. A n additive Gaussian noise is added to all the inner product terms to simulate the effects of noise in the system. The elements of the matrix of delayed inner product terms (2.103) are evaluated using (2.90) . This almost completes the calculations for the first interval u>j. What remains to be done is to make a time shift to account for the modal time constraint given in (2.62). The beginning of the second interval is taken to be the order of the subsystem plus the length of the sampling interval, iu,-, measured in units of time delays, td. The delay corresponding to the number of time delays times the order of the subsystem is chosen to meet the modal time constraint. The total interval delay, i l t- = (Ni + 2t2i)td, where f u is the beginning of the next interval, and t2i is half the duration of «;,-, is evaluated for both subsystems, and the maximum is used as tw This ensures greater sampled signal levels for y(t) (see (2.24)), as y(t) decays exponentially with time. Before we start the coefficient updating we first save a copy of the original coef-ficients to be used in evaluating the modal constant. The updating itself is done by Chapter 4. The Modal Program 33 t ime Figure 4.3: Intervals, «;,-, used to evaluate the modal parameters for the two subsystems. Chapter 4. The Modal Program 34 making the following observations. If over the first interval 0 < t < t 1 ? y{t) = c0 + cneSit + ••• + c l ne 3"' then (4.144) y{t) - c0 = c n e a i t + • • • + c l n e S n t and f = clSle^ + • • • + cnsnes"1 , H - l . S n l (4.145) by evaluating y[t), ^ , . . . , d"tn-l at t = ti we obtain the updated coefficients for the second interval £\ -\ h c N s = cio + C i e S l t l H h C j v e ^ 1 1 - c2o - cci sici + • • • + sN„c^l = s i C i e ' 1 ' 1 H (- SNCrteSNtl — AiccJ — • C C J V 0 »1 * i + • • •+*#rlg^ = s,1"s-vi" + -+4s"v - A f s _ 1 c c l A ^ - ' c c ^ . where the cc",- for SI are calculated before the actual updating from the updated coeffi-cients of S2. The S2 coefficients are updated before the S i coefficients. Note, that for S2, Ns = N2 and Nc = 0. Therefore we can set up the following Vandermonde system to solve for the £ coefficients for the next interval 1 1 Sl Si _ n - l Q r i - 1 6 1 6 1 1 „n-l = dt y{ti) - c0 - ccx C C J V ( 7 AiccJ - • • • - \NCCCN0 t=t. t=tL (4.146) Chapter 4. The Modal Program 35 where n = Ns. The coefficient updating is done by the subroutine Update and the Vandermonde system is solved by the subroutine Vandm which is based on algorithm 5.6-2 in [9]. 4.3.2 Obtaining the Modal Parameters The modal parameters for the first states are obtained by solving the matrix equation (2.103) by the use of the subroutine Dslimp [14]. The subroutine performs Gaussian elimination in double precision arithmetic by the use of an iterative routine. It continues the iterations until one of three events occur. • the error vector calculated by the subroutine is less than a preset tolerance, the variable E P S . • the iterations start to diverge. • until a preset maximum number of iterations is performed. The subroutine Check prints out the modal parameters, the value of E P S , which indicates the accuracy of the solution returned by Dslimp, and evaluates the equations of (2.81) to ensure that the equality holds. 4.3.3 Evaluating the Modal Constant The evaluation of the Mi0 term is done by the subroutine Const which evaluated (2.109) to obtain M.o-The (w(t),eat) terms are calculated by the function Wo where w(t) is given by (2.105). The function Woy in turn, uses Wo to evaluate the (w(t),y(t — itd)),i = 1 , . . . , N terms of (2.109) and to calculate Mi0. Chapter 4. The Modal Program 36 4.3.4 Finding the Modal Parameters for the Second and Third States With the work that had to be done to calculate the modal parameters of the first state of the system and the subsystem the calculation of the modal parameters of the second and third state is relatively trivial. The matrix of (2.104) is the matrix of (2.103) but shifted one column to the right. The ntd delayed terms of (2.103) drop out and the first column of (2.104) is simply the non-delayed inner product terms of (2.103). Once this shifting is complete the modal parameters are obtained by the same procedure of the subroutines Dslimp, Check, and Const described above. 4.4 Calculating the System Responses Program Moda l calculated three different responses of the system and the subsystem: • the unit step response. and through the process of modal filtering: • the periodic unit step response. • the periodic unit square wave, sal(l,x), response. 4.4.1 Unit Step Response The unit step responses for the system and the subsystem are calculated using the theory outlined in chapter 3. The user is prompted for the number of values of the unit step response desired, up to fifty points. The response is sampled at intervals td. The function Y is used to evaluate y(t) and Y(t,n) is calculated using (2.33) which gives G(t,n) using (3.117). For purposes of convenience and without loss of generality Chapter 4. The Modal Program 37 it has been decided to let v(t) = y(i), t < 0. Then since y{t) = v(t)+g{t) Au{l) (4.147) due to the fact that g(t) = 0, t < 0. Then g(t) is evaluated using the recursive relationship of (3.119). 4 . 4 . 2 Periodic Unit Step Response The periodic unit step response, an odd symmetric function of period 2o = 2(N + is calculated recursively from (3.118). 4 . 4 . 3 Periodic Unit Square Wave Response The response to an input of the Walsh function, sal(l ,x), is also odd, symmetric, and of period 2o. The response is given by (3.132). In order to generate the RSl(x,n) values we must know the 7J parameters. The subroutine Incoef derives the 7 parameters from the modal parameters. Incoef sets up the matrix C in (3.127) and then calls upon the subroutine Finv [14] that inverts C. From C - 1 the 7, parameters can be read off the first column of C - 1 . For both the periodic unit step response and the periodic unit square wave response only the N + l values from t = 0 to t = o are calculated. The other N + l values are generated form these. Chapter 4. The Modal Program 38 4 . 4 . 4 Calculation of the Finite Fourier Transform Coefficients To calculate the finite Fourier transform coefficients of the periodic unit step response we use the subroutines Prefft and Fft [12] to implement a fast Fourier transform al-gorithm on one period of Gp(t,n). Since the algorithm only works for a number of data points equal to 2k, we must pass to the subroutine more than one period of data values. In the subroutine Setfft, which calls Prefft and Fft , the values of one complete period of Gp(t, n) is stored in an array. Then subsequent values from the next period are added to the array until the total number of values is 2k. These values are then used by the fast Fourier transform subroutines. Chapter 5 Modal Parameter Estimation with Noise In any reasonable, physical system there is always some uncertainty in readings or measurements made on the outputs of the system. In the program Modal , it is possible to add a zero-mean Gaussian noise to the system under consideration and observe the effects of the noise on the modal parameter estimation and non-parametric system identification. However, it is beyond the scope of the present program to implement any process of statistical averaging or other technique to attempt to reduce the effect of the noise. What the PC program does is to give a measure of the 'sensitivity' of the various modal function techniques to an additive noise. This can be of assistance in estimating the number of data intervals required to achieve a given accuracy in the presence of noise. The use of more sophisticated techniques to reduce or eliminate the effects of noise will have to be dealt with at a future time. 5.1 Sensitivity Analysis of Modal Parameters The addition of noise to the systems modelled by the program Moda l has the effect of giving the user an indication of the sensitivity of the modal parameter estimates to an additive Gaussian noise. During the process of modal parameter estimation we are required to sample the outputs of the system. The measurement of output will contain an uncertainty due in part to the inherent uncertainty of the measurement device. This uncertainty can be best modelled as an additive, random noise with uniform distribution, centered at zero and with range [ — w h e r e W is the accuracy of the 39 Chapter 5. Modal Parameter Estimation with Noise 40 measurement device. fx(x) = I 0 , otherwise k ,-W < x < W (5.148) In practice, the inner product terms used to estimate the modal parameters must be evaluated by some technique of numerical integration such as the trapezoidal rule (see section 5.3). In evaluating the trapezoidal approximation to the integral, a number of data samples must be summed. Since all of the data samples are corrupted by an additive uniformly distributed noise, the central limit theorem [18], says that the final approximate integration will be corrupted by an additive Gaussian noise. To illustrate the effects of the addition of additive Gaussian noise to a signal it will be demonstrated what happens when noise is added to a first order signal. (The result may be easily generalized to the case of a n-th order signal by the use of matrix notation.) Consider the case of Taking the inner products of (5.149) with the interval w(t) as in (2.74) we obtain y{t) = c 0 + cxest. (5.149) (5.150) where E = - [2e3btd - esati - esct'1] . s (5.151) From (5.150) we obtain the modal function representation of (5.149), y{t) + Myy{t - td) = M 0 . (5.152) For the noise free case it is observed that ( W ( 0 , y ( i ) ) + M 1 ( « ; ( f ) , y ( ' - i ( i ) ) = 0 Chapter 5. Modal Parameter Estimation with Noise 41 CiE + c^EM^-*1'1 = 0 l + Mle-3ti = 0 (5.153) where (5.153) is the definition of M i . If noise is added to the system, the following is seen, y(t) = y{t)+n(t) (5.154) (w{t),y(t)) = {w{t),y{t))+n, (5.155) and the following data equation is obtained, ctE + n + Mi [cie~stdE + n] = 0 . (5.156) If (5.156) is solved for the change in modal parameter estimates caused by noise, M i — M i , it is seen that ( M i - Mx)Cle-at' lE = [l + M i ] n M i - M i = n [l + M i ] E^e^cT1 (5.157) where — =the noise to signal ratio. To use a numerical example choose a = l , b=2, c=3, £<j=0.1, and s=-3. Then = I [2c-°-6 - c-°-3 - c0-9] = -0.01659 = JETV 0 ' 3 = -44.659 M i = - e " 0 3 = -0.74082 Therefore the relative change in M i due to noise is M i - M i (n\ r 1 "I . , fn\ , ^MT S U) I1 + Mil <-44'659' = I5'62 k) <5'158> where M i has been used instead of M i in the [l + term. Chapter 5. Modal Parameter Estimation with Noise 42 It can be seen that the term that most influences the sensitivity of Mx is the E~1eati< term that is quite large in this case. This demonstrates the importance of choosing the factors that make up the E~ lestd term; the sample interval, td, and the actual form of the interval w(t). For these preliminary investigations, w(t) was chosen simply to eliminate the constant Mo terms and td was completely arbitrary. Clearly the choice of these parameters is an area that deserves closer study. Equation (5.157) is suitable for sensitivity studies and for non-statistical averaging of the minimum number of intervals. In (5.158) it is seen that if sensitivity of less than 1% is desired, we require and if cx « 1, then 15.62 = 0.01 0.01 4 n = = 6.40 x 10~4 15.62 5.2 Statistical Model for the Simulation of Noise The general vector form of (5.157) is given by M-M = n(e + M){VmEC)~ l = 7^(e + M ) r (5.159) where T = ( V m £ ' C ) _ 1 and the factor 7 is chosen such that the elements of matrix T are of order unity. It is clear from (5.158) and (5.159) that some form of averaging must be used to reduce the effect of yn, on the estimate. Statistical ensemble averaging yields, since E[n] = 0, E\M\ — M - 7 E [ n M ] r . (5.160) Mult iplying (5.159) by n yields - E\nM) = io\e + M)T, (5.161) Chapter 5. Modal Parameter Estimation with Noise 43 where o 2 — E\n 2) and where the approximation E[n2M] — o 2M is used. Substituting (5.161) into (5.160) yields E[M] =M+ (-ya)2(e + M)r 2. (5.162) Equation (5.162) clearly illustrates the parameter bias term that is common to all least-square estimation procedures. The standard methods of eliminating bias is to use two-stage least-squares or maximum-likelihood estimation, or instrumental variables. In the modal function approach, due to the elimination of transients, the control input is a natural instrumental variable. Furthermore, since continuous-time data is sampled, integration can be used as a form of data compression and noise reduction. This is of importance since 7 may be large (see (5.158)). Equation (5.162) is a standard least-squares form showing parameter bias and can be implemented in a computer program. It follows from the discussion given in section 5.3 that the reduced effective noise variance is determined by o — (5.163) where /? is an attenuation factor associated with a high-frequency analog noise filter, Ni =number of samples used for integration, and iV 2 =number of intervals used. Example: Let , o = 0.1, and -7 = 100.0. Then (io) = io~ 2 If it is desired that PyjNiNt = 3.1 x 102, then one possible choice of Nx, iV 2 , and f3 is that Ni = N2 = 100, and /? = 3.1. Chapter 5. Modal Parameter Estimation with Noise 44 i 1 (t) r time a i Figure 5.4: Interval, w(t). 5.3 N o i s e R e d u c t i o n T h r o u g h N u m e r i c a l In tegra t ion To reduce the sensitivity of modal parameter estimates to the factor 7 a it is desirable to use preprocessing reduce the variance of n by as much as possible. It can be shown that by using the process of trapezoidal approximation that the noise variance may be reduced to an arbitrari ly small level. Assume that the system output , y(t), is corrupted by an additive uniform noise so that y{t) = y{t)+ n{t). (5.164) Therefore, for an interval, w(t), as shown in figure 5.4, (w{t),y{t)) = {w{t),y{t)) + (w(t),n{t)) = (w{t),y{t)) + 77. (5.165) The first term of (5.165) can be evaluated either analytical ly or by trapezoidal ap-proximat ion as shown previously but the V term must be evaluated by an approximation method since we cannot model n(t) as a sum of weighted exponentials. Chapter 5. Modal Parameter Estimation with Noise 45 Evaluating u we get A v + 2{u{tx + A ) + h + [N - 1) A) ) + - \u(t) + 2(i/(t + A ) + • • • + u(t + (N — 1) A)) + i/(t2)] (5.166) — + 2{u{h + A ) + • • • + u(tx + {N - 1) A ) + 0 +u{t + A ( + • • • + u(t + {N - 1)A)) + u(t2)}. (5.167) Since each u(t) term is a uniformly distributed random variable and we are adding up 2N of them, we can, for sufficiently large N, say that V ~ iV(0,<72) due to the central limit theorem [18]. The next step is to evaluate E [ F ( « 1 ) P t ( < 1 ) ] . Let E [i/(ti)^ r(ti)] = o2I. Then E [jTu1 2 A c r ^ 2 c r 2 + 2 ( c r 2 + • • • + CT2 + 0 + c r 2 H h c r 2 ) + C T ' 2(w-i) terms = ( ^ ) [2 + 4 ( i V - l ) ] J = A V 2 (5.168) From (5.168) we can see that if we choose enough partitions for the trapezoidal approzimation, A will become as small as we like and E \vvT] will be as close to zero as we wish. So we can, with a sufficient number of samples, reduce the noise variance in estimating modal parameters. 5.4 Further Concerns Associated with Modal Functions There has been much work in recent years dealing with the fitting of exponentials and specifically dealing with Prony's method. Since the method of modal functions is derived from Prony's method, these works have a direct bearing on this thesis. Chapter 5. Modal Parameter Estimation with Noise 46 A basic consideration is whether techniques based on Prony's method perform any better than the older, more traditional F F T approach. Bucker [5] and Spitznogle and Quazi [22] compared the two techniques and found that in some cases, Prony's method was superior to the F F T approach. In particular, Prony's method gave better resolution in the estimation of poles and behaved better in low-noise situations. Van Blar icum and Mi t t ra [24] showed how the poles and residues of a system can be extracted directly from its transient response when the number of poles are known. However, in the case where the number of poles or the order of the system is unknown it is necessary to estimate the order of the system before applying either Prony's method or modal parameter estimation. Holt and Anti l l [10] used singular value decomposition (SVD) to determine the number of terms to use in Prony's method by observing the singular values of an equally spaced data matrix. If a system is of order n, then the singular values of a data matrix of order greater than n should become quite small after the n-th singular value. Van Blar icum and Mi t t ra [25] used standard deviations to determine the number of poles of a system when noise was present. Another area of interest concerning modal parameter estimation is what an optimum sampling period would be. Crittenden et al. [6] investigated the effect of sampling period on Prony's method when identifying state parameters of a continuous time system. The problem of what constitutes 'persistently excited inputs' has been looked at by Schaubert [20] in the case of step-like excitation. When using modal functions, the Co coefficients and the intervals, w(t), represent a form of pseudo-random step-like excitation but just what would be an optimum input must be investigated further. To deal with noise two procedures seem to have been put forward. The first involves the modelling of a system by a higher order model. B y using forward-backward linear prediction ( B F L P ) the actual poles of a system may be distinguished from those poles Chapter 5. Modal Parameter Estimation with Noise 47 due to noise by superimposing the results on a unit circle. While the legitimate poles are mirrored across the unit circle for the two cases, the noise poles remain constant (Kumaresan and Tufts [11] and Tufts and Kumaresan [23]). Poggio et al. [16] also investigated the use of rank-overspecification and filtering on noise. The second method consists of using S V D to decompose a noisy data matrix into a data matrix and a noise matrix. This is done by observing the singular values of the noisy data matrix and constructing a new data matrix by excluding the smaller singular values that are assumed to be caused by noise (Holt and Anti l l [10], Kumaresan and Tufts [11], and Tufts and Kumaresan [23]). Both of these methods could be useful for the estimation of modal parameters and this should be investigated further. C h a p t e r 6 D i s c u s s i o n a n d C o n c l u s i o n s In this thesis, it has been shown that the method of modal functions offers a new and very useful tool for both parametric and non-parametric system identification. Modal functions are independent of system initial conditions and allow the representation of not only a signal, but also its time-derivatives by a simple continuous-time A R M A model. For the case of a M I M O system, the number of parameters that must be estimated is greatly reduced by being able to decouple the larger system into smaller subsystems. This increases the speed and simplifies the parameter estimation problem. The actual estimation of the parameters is carried out through the use of simple matrix algebra and, in practical cases, by using trapezoidal approximation to calculate integral terms, it is possible to reduce the noise variance. Through the use of modal functions we can directly identify the system unit step re-sponse and the square wave response directly from input-output data. As well, modal filtering can produce a periodic unit step response which may be used to find F F T coefficients without the use of window functions. The identification of the periodic square-wave response from non-periodic data may be of use in the field of signal pro-cessing. The PC program Modal will enable the user to investigate the above processes while varying the different external parameters in an easy and efficient manner and will further the understanding of modal functions. 48 B i b l i o g r a p h y [1] Bohn , Erik V . "Measurement of continuous-time linear system parameters via Walsh functions." IEEE Transactions on Industrial Electronics 29 (1982): 38-46. [2] Bohn , Er ik V . "Recurssive evaluation of Walsh coefficients for multiple integrals of Walsh series." Automatica 20 (1984): 243. [3] B o h n , Er ik V . "Walsh function decoupled parameter estimation equations for dy-namic continuous-time models with time-delay." 7th IFAC/IFORS Symposium. Identification and System Parameter Estimation York, U .K. : 799. [4] B o h n , Erik V . " Non-Parametric and Parametric Identification of M I M O Systems via Moda l Functions". (unpublished). [5] Bucker, H.P. "Comparaison of F F T and Prony algorithms for bearing estimation of narrow-band signals in a realistic ocean environment." Journal of the Acoustical Society of America 61 (1977): 756-762. [6] Crittenden, J . L . , R . J . Mulhol land, J . Hill IV, and E . F . Martinez. "Sampling and data analysis of Prony's method applied to model identification." International Journal of Systems Science 14 (1983): 571-584. [7] El-Shafey, M . H . and Erik V . Bohn. " The use of modal functions for continuous time system identification and state observer design." International Journal of Control 45 (1987): 1723-1736. 49 Bibliography 50 [8] Froberg, Car l -Er ik . Numerical- Mathematics: Theory and Computer Applications. Don Mil ls, Ontario: Benjamin/Cummings Publishing Company, 1985. [9] Golub, Gene H . and Charles F . Van Loan. Matrix Computations. Baltimore, Mary-land: The Johns Hopkins University Press, 1985. [10] Holt , J . N . and R . J . Anti l l . "Determining the number of terms in a Prony algorithm exponential fit." Mathematical Biosciences 36 (1977): 319-332. [11] Kumaresan, R. and D .W. Tufts. "Estimating the parameters of exponentially damped sinusoids and pole-zero modeling in noise." IEEE Transactions on Acous-tics, Speech, and Signal Processing 30 (1982): 833-840. [12] Marple, S .L . , Jr. Digital Spectral Analyses with Applications. Englewood Cliffs, New Jersey: Prentice-Hall, 1987. [13] Merchant, Michael J . Applied FORTRAN Programming with Standard FORTRAN, WATFOR, WATFIV and Structured WATFIV. Belmont, California: Wadsworth Publishing Company, Inc., 1976. [14] Nicol, Tom. UBC MATRIX: A Guide to Solving Matrix Problems. Vancouver, Brit ish Columbia: U B C Computing Centre, 1982. [15] Perdreauville, F . J . and R . E . Goodson. " Identification of systems described by partial differential equations." Jou. of Basic, Engg. Trans. ASME 88, Series-D, No.2 (1966): 463-468. [16] Poggio, A . J . , M . L . Van Blaricum, E . K . Miller, and R. Mit tra. "Evaluation of a processing technique for transient data." IEEE Transaction on Antennas and Propagation 26 (1978): 165. Bibliography 51 [17] Rao, G.P . and K .R . Palanisamy. "Improved algorithms for parameter identification in continuous systems via Walsh functions." IEE Proceedings Part D 130 (1983): 9-16. [18] Ross, Sheldon. A First Course in Probability. New York, N.Y.: Macmil lan Pub-lishing Company, 1984. [19] Saha, D . C . and G.P . Rao. Lecture Notes in Control and Informational Sciences: Vol. 56, Identification of Continuous Dynamical Systems. The Poisson Moment Functional (PMF) Approach.. E d . A . V . Balakrishnan and M . Thoma. New York, N.Y.: Springer-Verlag, 1983. [20] Schaubert, D .H . "Application of Prony's method to time-domain reflectometer data and equivalent circuit synthesis." IEEE Transactions on Antennas and Prop-agation 27 (1979): 180. [21] Shinbrot, M . " O n the analysis of linear and nonlinear systems." Transactions of the ASME 79 (1957): 547-552. [22] Spitznogle, F .R . and A . H . Quazi . "Representation and analysis of time-limited signals using a complex exponential algorithm." Journal of the Acoustical Society of America 47 (1969): 1150-1155. [23] Tufts, D . W . and R. Kumaresan. "Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood." Proceedings of the IEEE 70 (1982): 975-989. [24] Van Blar icum, M . L . and R. Mit tra. " A Technique for extracting the poles and residues of a system directly from its transient response." IEEE Transactions on Antennas and Propagation 23 (1975): 777-781. Bibliography 52 [25] Van Blar icum, M . L . and R. Mit tra. "Problems and solutions associated with Prony's method for processing transient data." IEEE Transactions on Antennas and Propagation 26 (1978): 174-182. [26] Wellstead, P .E . " Non-parametric methods for system identification." Automatica 17 (1981): 55. [27] Young, P. "Parameter estimation for continuous- time models - a survey." Auto-matica 17 (1981): 23. Appendix A A Trial Run of Program Modal To demonstrate the use of the program Modal the following example is given. The following system of one input and two outputs was simulated. yi(0 = c 1 0 + 100e~ f - 100e~ 2 t + 100e _ 3 t y2(0 = Cio - 100e- 3 < + lOOe" 2 ' , For this system, it can be seen that subsystem SI consists of the modes (—1, —2, —3) and that subsystem S2 consists of the modes (—3, —2) (ie. the third mode of SI is coupled with the first mode of S2). The factored roots of d^iui) were given as ( — 1, — 1). For values of c i 0 , the following were chosen for each interval: 100, —100,0. This cor-responds to a control input of a unit-impulse for interval 1, and unit-steps for intervals 2 and 3 (see figure A.5) . The integration intervals, Wi(t), t = 1,2,3, were defined by the duration of half of the interval, and were chosen to be (see figure A.6) • for SI: 8, 8, 4. • for S2: 8, 8. The sampling period was chosen to be 0.1 seconds and a value of 0.0001 was chosen for V A R , the reduced effective noise variance (see 5.163). What follows next are the prompts and messages displayed on the P C monitor ( in t y p e w r i t e r p r i n t ) and the proper responses from the user. Modal parameter e s t i m a t i o n program. 53 Appendix A. A Trial Run of Program Modal 54 1 100 \ i C O • • . _ ! — 1 •100 1 1 1 1 B 40 time i i i i t • i Figure A.5: Control input, c 0 for example A . F i r s t , sample data must be generated. Each of the 2 outputs i s assumed t o be a sum of complex e x p o n e n t i a l s . Subsystems; SI has 3 s t a t e s and S2 has 2 s t a t e s . E n t e r # of modes of subsystem S I . 3 Ent e r # of modes of subsystem S2. 2 Ent e r the number of uncoupled modes of S I . 2 A l l complex e n t r i e s must be of the form ( x , y ) . Subsystem S I , uncoupled modes o n l y ! E n t e r r e a l c o e f f i c i e n t 1 100 Appendix A. A Trial Run of Program Modal 55 rf(t) SI w(t) 32 311 41 ^5 19 time 19 21 29 37 88 49 time Figure A.6: Intervals for example A . Appendix A. A Trial Run of Program Modal E n t e r r e a l c o e f f i c i e n t 2 -100 E n t e r complex e x p o n e n t i a l 1 (-i.,o.) E n t e r complex e x p o n e n t i a l 2 (-2..0.) Subsystem S2: E n t e r c o u p l e d modes b e f o r e uncoupled! E n t e r r e a l c o e f f i c i e n t 1 -100 E n t e r r e a l c o e f f i c i e n t 2 100 E n t e r complex e x p o n e n t i a l 1 (-3..0.) En t e r complex e x p o n e n t i a l 2 (-2..0.) En t e r the r o o t s of the char, e q u a t i o n , r ( i ) . E n t e r complex r o o t , R ( l ) . (-1..0.) En t e r complex r o o t , R ( 2 ) . (-i.,o.) To i d e n t i f y the modal parameters of S I , S2 you must s p e c i f y 3 i n t e r v a l s . S p e c i f y the l e n g t h of each h a l f - c y c l e i n t e r v a l as a m u l t i p l e of T. En t e r CO v a l u e f o r i n t e r v a l 1 100 E n t e r CO v a l u e f o r i n t e r v a l 2 Appendix A. A Trial Run of Program Modal 57 -100 E n t e r CO v a l u e f o r i n t e r v a l 3 0 Subsystem 1 En t e r the h a l f - c y c l e d u r a t i o n of i n t e r v a l 1 8 En t e r the h a l f - c y c l e d u r a t i o n of i n t e r v a l 2 8 En t e r the h a l f - c y c l e d u r a t i o n of i n t e r v a l 3 4 Subsystem 2 En t e r the h a l f - c y c l e d u r a t i o n of i n t e r v a l 1 8 En t e r the h a l f - c y c l e d u r a t i o n of i n t e r v a l 2 8 En t e r the sampling p e r i o d , T, i n seconds. .1 E n t e r the n o i s e v a r i a n c e . 0.0001 Random number gen e r a t o r i n i t i a l i z a t i o n . E n t e r R ( 0 ) ; (between 0 and 1) .78478143 This completes the information that the user must enter from the PC keyboard. The output is displayed on the screen and stored in a file named M O D A L . L I S . For this example the file M O D A L . L I S is shown below. Appendix A. A Trial Run of Program Modal 58 Modal Parameter E s t i m a t i o n Program. Data g e n e r a t i o n parameters: Subsystem 1 has 3 s t a t e s and 3 modes. Subsystem 2 has 2 s t a t e s and 2 modes. C l ( 1)= .100000E+03 C l ( 3)= .100000E+03 C2( 1)= -.100000E+03 S l ( 1)= -.100000E+01 S l ( 3)= -.300000E+01 S2( 1)= -.300000E+01 R( 1)= -.100000E+01 CO f o r i n t e r v a l 1= CO f o r i n t e r v a l 2= CO f o r i n t e r v a l 3= Subsystem 1 I n t e r v a l # 1 has a I n t e r v a l # 2 has a I n t e r v a l # 3 has a Subsystem 2 I n t e r v a l # 1 has a I n t e r v a l # 2 has a .100000E+03 .000000E+00 -.200000E+01 .OOOOOOE+OO -.200000E+01 .000000E+00 100000E+01 .000000E+00 .000000E+00 C l ( 2)= -.100000E+03 .000000E+00 .000000E+00 C l ( .000000E+00 C2( 2)= .OOOOOOE+OO S l ( 2)= .000000E+00 S l ( .000000E+00 S2( 2)= .000000E+00 R( 2)= .100000E+03 -.100000E+03 .000000E+00 h a l f c y c l e of 8*T s. h a l f c y c l e of 8*T s. h a l f c y c l e of 4*T s. h a l f c y c l e of 8*T s. h a l f c y c l e of 8*T s. The sampling p e r i o d i s .100E+00 seconds. Appendix A. A Trial Run of Program Modal 59 The a d d i t i v e n o i s e ~N(0, .0001) New C parameters f o r i n t e r v a l 2 C l ( 1)= .814957E+03 .000000E+00 C l ( 2)= C l ( 3)= .40033BE+03 .000000E+00 C l ( C2( 1)= -.400335E+03 .000000E+00 C2( 2)= New C parameters f o r i n t e r v a l 3 C l ( 1)= -.278108E+03 .000000E+00 C l ( 2)= C l ( 3)= -.198660E+03 .000000E+00 C l ( C2( 1)= .198660E+03 .000000E+00 C2( 2)= The modal parameters, M 1 a r e : ( 1)= -.252332E+01 ( 2)= .212513E+01 ( 3)= -.596610E+00 ( Checking! EPS adds up t o .185583E-14 The modal i d e n t i t y adds up t o , -.313777E-05 The modal i d e n t i t y adds up t o , -.370271E-05 The modal i d e n t i t y adds up t o , -.124320E-05 The M 1 c o n s t a n t i s .548101E+00. .100224E+04 .000000E+00 .602237E+03 .000000E+00 .477579E+03 .000000E+00 .286527E+03 .000000E+00 The modal parameters, M 4 a r e : ( 1)= -.156070E+01 ( 2)= .607591E+00 Checking! Appendix A. A Trial Run of Program Modal EPS adds up t o .000000E+00 The modal i d e n t i t y adds up t o , .613459E-07 The modal i d e n t i t y adds up t o , .931981E-06 The M 4 co n s t a n t i s .468830E+01. 0***ITERATIVE IMPROVEMENT IS DIVERGING The modal parameters, M 2 a r e : ( 1)= -.128752E+02 ( 2)= .161173E+02 ( 3)= -.322496E+01 ( Checking! EPS adds up t o -.346528E-13 The modal i d e n t i t y adds up t o , -.315092E-05 The modal i d e n t i t y adds up t o , -.301753E-04 The modal i d e n t i t y adds up t o , -.474889E-05 The M 2 co n s t a n t i s .169810E+01. Appendix A. A Trial Run of Program Modal 0***ITERATTVE IMPROVEMENT IS DIVERGING The modal parameters, M 3 a r e : ( 1)= -.422736E+02 ( 2)= .942608E+02 ( 3)= -.515318E+02 ( Checking! EPS adds up t o -.568457E-12 The modal i d e n t i t y adds up t o , -.705720E-05 The modal i d e n t i t y adds up t o , .573695E-06 The modal i d e n t i t y adds up t o , -.385358E-05 The M 3 c o n s t a n t i s .449344E+02. 0***ITERATIVE IMPROVEMENT IS DIVERGING The modal parameters, M 5 a r e : ( 1)= -.753411E+01 ( 2)= .780764E+01 Checking! EPS adds up t o -.743869E-14 The modal i d e n t i t y adds up t o , -.689855E-06 The modal i d e n t i t y adds up t o , -.111486E-04 Appendix A. A Trial Run of Program Modal The M 5 co n s t a n t i s .273397E+02. The f i r s t 10 v a l u e s of the s t e p response are: Subsystem 1 G( 1) = .199878E+03 G( 2) = .182303E+03 G( 3) = .168947E+03 G( 4) = .158627E+03 G( 5) = .150512E+03 G( 6) = .144020E+03 G( 7) = .138739E+03 G( 8) = .134377E+03 G( 9) = .130724E+03 G(10)= .127627E+03 Subsystem 2 G( 1) = .999805E+02 G( 2) = .107748E+03 G( 3) = .112087E+03 G( 4) = .114145E+03 G( 5) = .114727E+03 G( 6) = .114386E+03 G( 7) = .113503E+03 G( 8) = .112334E+03 G( 9) = .111046E+03 G(10)= .109747E+03 The 4 v a l u e s of the p e r i o d i c s t e p response are Subsystem 1 GP( 1)= .199604E+03 GP( 2)= -.322326E+03 GP( 3)= .133432E+03 GP( 4)= .213424E+00 Subsystem 2 GP( 1)= .976363E+02 GP( 2)= -.506359E+02 GP( 3)= .232760E+01 GP( The P e r i o d i c Square Wave Response i s g i v e n by: Subsystem 1 Appendix A. A Trial Run of Program Modal 63 GAMMA( 1)= -.253115E+01 GAMMA( 2)= -.138587E+01 GAMMA( 3)= .409452E+00 GAMMA( 4)= .246822E+01 Subsystem 2 GAMMA( 1)= .474376E-01 GAMMA( 2)= .109263E+01 GAMMA( 3)= .167645E+01 GAMMA( The 4 v a l u e s of the per. square wave r e s p . are Subsystem 1 RS( 1)= .236009E+03 RS( 2)= .209802E+03 RS( 3)= .190169E+03 RS( 4)= .175230E+03 Subsystem 2 RS( 1)= .869770E+02 RS( 2)= .100376E+03 RS( 3)= .108466E+03 RS( FFT c o e f f i c i e n t s of G p ( t ) . Subsystem 1 F o u r i e r C o e f f i c i e n t s : Modal f u n c t i o n . ( 1)= .000000E+00 .000000E+00 ( 2) = -.569322E+01 .188673E+02 ( 3) = .OOOOOOE+OO .OOOOOOE+OO ( 4) = .855347E+02 .722399E+02 ( 5) = .OOOOOOE+OO .OOOOOOE+OO ( 6) = .855347E+02 -.722399E+02 ( 7) = .000000E+00 .000000E+00 ( 8) = -.569321E+01 -.188673E+02 Subsystem 2 F o u r i e r C o e f f i c i e n t s : Modal f u n c t i o n . ( 1)= .470004E+01 .OOOOOOE+OO ( 2) = .460756E+01 . 162702E+02 ( 3) = .483083E+01 -.953087E+01 ( 4) = .479253E+01 -.279153E+01 ( 5) = .449472E+02 .000000E+00 ( 6) = .479253E+01 .279153E+01 ( 7) = .483083E+01 .953087E+01 ( 8) = .460756E+01 -.162702E+02 Appendix A. A Trial Run of Program Modal 64 Modal Parameter Estimates Subsystem S i Subsystem S2 m a -2.52332 m 4 1 -1.56070 m 1 2 2.12513 m 4 2 0.607591 mis -0.596610 MlQ 0.548101 M 4 0 4.68830 m2i -12.8752 m 5 1 -7.53411 m22 16.1173 m52 7.80764 m23 -3.22496 M20 1.69810 27.3397 m3i -42.2736 ™32 94.2608 m33 -51.5318 M30 44.9344 Table A . l : Summary of the modal parameter estimates for example 1. The modal parameter estimates are summarized in table A . l . The non-parametric system responses for subsystem SI are summarized in table A.2 and for subsystem S2 in table A .3 . The F F T coefficients of the periodic unit step response are summarized in table A .3 . The system responses for SI and S2 are graphed in figures A.7 and A.8 respectively. endix A. A Trial Run of Program Modal System Responses, Subsystem SI Time G{t,n) Gp[t,n) R3l(t,n) 0.0 199.878 199.604 236.009 0.1 182.303 -322.326 209.802 0.2 168.947 133.432 190.169 0.3 158.627 0.213424 175.230 0.4 150.512 -199.604 -236.009 0.5 144.020 322.326 -209.802 0.6 138.739 -133.432 -190.169 0.7 134.377 -0.213424 -175.230 0.8 130.724 0.9 127.627 Table A.2: Summary of non-parametric system responses for subsystem SI. System Responses, Subsystem S2 Time G{t,n) Gp{t,n) RSl(t,n) 0.0 99.9805 97.6363 86.9770 0.1 107.748 -50.6359 100.376 0.2 112.087 2.32760 108.466 0.3 114.145 -97.6363 -86.9770 0.4 114.727 50.6359 -100.376 0.5 114.386 -2.32760 -108.466 0.6 113.503 0.7 112.334 0.8 111.046 0.9 109.747 Table A.3: Summary of non-parametric system responses for Subsystem S i . Appendix A. A Trial Run of Program Modal 66 Unit Step Response of SI aoo 1 B O -1 6 0 -1 4 0 - 1 P e r i o d i c Unit Step Response of SI « 2 0 0 - . 0-* 0 0 -aoo*-* ' . l \ 0 .2 /0 .3 0.\^0.JS 0.6 X(.J/0.B 0.9 lne Square Wave Response of SI 8 0 0 -1 0 0 -- 1 0 0 -- 2 0 0 -0 0.1 0.2 0.3 0 .4 \0 .S 0.6 0.7 0L8 0.9 Rime Figure A.7: System responses for S i . Appendix A. A Trial Run of Program Modal Uni t Step Response of S2 1 1 9 - 1 i « H 0.1 0.2 0.3 0.4 O.S 0.6 0.7 0.8 0.9 time P e r i o d i c Uni t Step Response of S2 too-, - 1 0 C M 0 . 4 / 0.5 0.6 0.7 0.8 0.9 time lOO-i Square Wave Response of S2 .« 0 0.1 0.2 0.3 \ 0 . 4 0.5 0 - 1 0 0 - 1 .6 0.7 0.8 0.9 time F i g u r e A . 8 : S y s t e m responses for S2 . Appendix A. A Trial Run of Program Modal 68 F F T Coefficients Subsystem S i Subsystem S2 0.000000 +J 0.000000 4.70004 +3 0.000000 -5.69322 +3 18.8673 4.60756 +3 16.2702 0.000000 +3 0.000000 4.83083 -3 9.53087 85.5347 +3 72.2399 4.79253 -3 2.79153 0.000000 +3 0.000000 44.9472 +3 0.000000 85.5347 -3 72.2399 4.79253 +3 2.79153 0.000000 +3 0.000000 4.83083 +3 9.53087 -5.69321 -3 18.8673 4.60756 -3 16.2702 Table A . 4 : Summary of F F T coefficients. Appendix B A Method of Reducing the Bias of Noisy Estimates It was shown in equation (5.162) that a bias term arises when estimating modal pa-rameters in the presence of noise. B y varying the variables V A R and R(0) in the PC program it is possible to obtain enough random number samples to simulate this bias effect. It was observed that M — M appeared to be a linear function of \ / V A R scaled. B y plott ing the bias against VVAR the slope of the resulting line may be obtained and then used to obtain an unbiased estimate of the modal parameters. The example in appendix A was run w i t h variances of 0.0, 0.0002, 0.0004, and 0.0008, and w i t h R(0) values of 0.25, 0.5, and 0.75 for each value of V A R . The modal parameter estimates for SI are shown in table B .5 , B.6, and B.7. In figure B.9 the mean of M — M is plotted against v / V A R and it can be seen that the bias effect appears to be linear in y/VAR and offer a promissing technique to obtain unbiased parameter estimates. M o d a l parameter estimates, VAR=0 .0002 M i - M i M 2 - M 2 M 3 - M 3 R(0) -0.062660 -0.066270 -0.068750 0.113200 0.124710 0.126440 -0.050261 -0.056922 -0.054151 0.25 0.50 0.75 Table B.5 : M o d a l parameter estimates in the prescence of additive measurement noise; V A R = 0 . 0 0 0 2 . 69 Appendix B. A Method of Reducing the Bias of Noisy Estimates 70 e r r o r in Ml o.oo. - 0 . 0 2 -0 0.014 0 . 0 2 - 0 . 0 4 H -0 .06 - 0 . 0 8 -- 0 . 1 0 -- 0 . 1 2 -0.028 1 — J l / A R e r r o r In M2 0.25-1 0.20 0.15 0.10 0 .05- I 0.00 -U 0 1 1 0.014 0 . 0 2 l l / A R 1—>• 0.028 e r r o r i n M3 o O . Q O -0 -0 .02 • -0 .04 ' -0 .06 --0 .08 --0 .10 -0.014 • 0.02 0.028 1 — > i V A R Figure B.9: Graph indicating the linear relationship between modal parameter estimate biases and the noise standard deviation. Appendix B. A Method of Reducing the Bias of Noisy Estimates 71 Modal parameter estimates, VAR=0.0004 Mi-Mi M 2 - M 2 M 3 - M 3 R(0) -0.088690 -0.091520 -0.095560 0.160400 0.172300 0.175720 -0.071286 -0.078662 -0.075221 0.25 0.50 0.75 Table B.6: Modal parameter estimates in the prescence of additive measurement noise; VAR=0.0004. Modal parameter estimates, VAR=0.0008 Af i - Mi . M 2 - M 2 M 3 - M 3 R(0) -0.12559 -0.12528 -0.13193 0.227460 0.235980 0.242560 -0.101229 -0.107764 -0.103764 0.25 0.50 0.75 Table B.7: Modal parameter estimates in the prescence of additive measurement noise; VAR=0.0008. Appendix C Listing of Program Modal Below is a listing of the program Modal and the subroutine Vandm, both written by the author. A l l the program code used was written in F O R T R A N and was compiled using Microsoft F O R T R A N 7 7 V3.31, August 1985. During the compilation of the program, the following object ( * .OBJ) files and libraries (*.LIB) were linked: The object files M O D A L , D S L I M P , F I N V , V A N D M , and F F T , and the libraries F O R T R A N , 8087, and M A T H . C . l Program Modal Listing. PROGRAM MODAL C C T h i s program w i l l c a l c u l a t e the modal parameters of a C 2 subsystem system. The f i r s t subsystem c o n s i s t s of NI C modes and 3 s t a t e s w h i l e the second subsystem has N2 C modes and 2 s t a t e s . C C The output of t h i s program i s s t o r e d i n the f i l e C MODAL.LIS. C C V a r i a b l e Legend: 72 Appendix C. Listing of Program Modal C C S I , S2- the complex exponents. C C I , C2- the complex c o e f f i c i e n t s ; these are updated C f o r every i n t e r v a l . C C01, C02- the o r i g i n a l complex c o e f f i c i e n t s . C INNER1, INNER2- the v a l u e s of the i n n e r p r o d u c t s of C the e x p o n e n t i a l terms and the i n t e r v a l s . C W - the complex random number ~N(0,1). C EPS- the t o l e r a n c e of the s u b r o u t i n e DSLIMP, used t o C f o r the modal parameters. C INNY1, INNY2- the va l u e s of the i n n e r products of C y ( t ) and the i n t e r v a l , w. C INNYP1, INNYP2- the v a l u e s of the i n n e r p r o d u c t s of C y ' ( t ) and w. C INYPP1- the v a l u e s of the i n n e r p r o d u c t s of y ' ' ( t ) and C w. C A l , A2- the m a t r i x of the v a l u e s of the del a y e d output C and w. C M l , M2, M3, M4, M5- the modal parameters. C T- the sampling p e r i o d . C CO- the v a l u e s of Co f o r each i n t e r v a l ; s i m u l a t e s the C change i n i n p u t . C A- dummy v a r i a b l e f o r r e a d i n g i n a r e a l v a l u e and making C i t complex. C MO- the modal c o n s t a n t s . C GI, G2- the v a l u e s of the u n i t s t e p response. Appendix C. Listing of Program Modal C GP1, GP2- the v a l u e s of the p e r i o d i c u n i t s t e p response. C GG1, GG2- temporary s t o r a g e f o r the v a l u e of the modal C f u n c t i o n of the u n i t s t e p response, G(x,n). C ADD1, ADD2- dummy v a r i a b l e s . C GAM1, GAM4- the gamma parameters. C SIGMA- the h a l f p e r i o d of the p e r i o d i c s i g n a l s . C RSI, RS2- the v a l u e s of the p e r i o d i c u n i t square wave C response. C R, RR- used t o i n i t i a l i z e the random number g e n e r a t o r . C VAR- the v a r i a n c e of the n o i s e . C I , J , K- l o o p c o u n t e r s . C T2- the d u r a t i o n s of the h a l f i n t e r v a l s of w. C N l , N2- the o r d e r of the system and the subsystem. C T l - d u r a t i o n of each i n t e r v a l , w. C N3- the number of v a l u e s of the u n i t s t e p response C d e s i r e d . C N4- the number of d i s t i n c t v a l u e s i n the p e r i o d i c responses. C COMPLEX Sl(10),S2(10).INNERl(lO),INNER2(10),TEMP(5),C1(10),C2(10) COMPLEX COl(lO),C02(10),W,CC(10),D1(10),D2(10),D3(10),AC(10) DOUBLE PRECISION EPS REAL INNYl(lO),INNY2(10).Al(lO.lO),A2(10.10),M1(10),M2(10) REAL M3(10),M4(10),M5(10),INNYP1(10).INNYP2(10),INYPP1(10) REAL T,C0(10),A,M0(10),G1(50),G2(50),GP1(10),GP2(10),ADD1 REAL ADD2.GG1,GG2.GAM1(11),GAM4(11).SIGMA.RS1(10),RS2(10) REAL R,RR,VAR.TNT Appendix C. Listing of Program Modal INTEGER I , J.K.T2K10) ,T22(10) ,N1 .N2,N3.N4,NC.NS.N,T1 OPEN(3,FILE='MODAL.LIS'.STATUS='OLD') EPS=1.0D-14 WRITE(*,100) 100 FORMAT(IX,'Modal parameter e s t i m a t i o n program.') WRITE(*,101) 101 FORMAT(IX,"First, sample d a t a must be generated.') WRITE(*,102) 102 FORMAT(IX,'Each of the 2 outputs i s assumed t o be a sum ') WRITE(*,202) 202 FORMAT(IX,'of complex e x p o n e n t i a l s . ' ) WRITE(*,103) 103 FORMAT(IX,'Subsystems; SI has 3 s t a t e s and S2 has 2 s t a t e s . WRITE(*.203) 203 FORMAT(IX,'Enter # of modes of subsystem SI.*) READ(*,*)N1 WRITE(*,204) 204 FORMAT(IX,'Enter # of modes of subsystem S2.') READ(*,*)N2 WRITE(*,290) 290 FORMAT(IX,'Enter the number of uncoupled modes of SI.') READ(*,*)NS WRITE(*,104) 104 FORMAT(IX,'All complex e n t r i e s must be of the form ( x . y ) . ' ) WRITE(*,105) 105 FORMAT(IX,'Subsystem S I , uncoupled modes o n l y ! ' ) Appendix C. Listing of Program Modal C C E n t e r the d a t a . C DO 1 1=1,NS WRITE(*,106)I 106 FORMAT(IX,'Enter r e a l c o e f f i c i e n t ' , 1 3 ) READ(*,*)A 1 C01(I)=CMPLX(A,0.) DO 2 1=1,NS WRITE(*.107)I 107 FORMAT(IX,'Enter complex e x p o n e n t i a l . ' , 1 3 ) 2 READ(*.*)S1(I) WRITE(*,108) 108 FORMAT(IX,'Subsystem S2: En t e r coupled modes b e f o r e uncoupled!') DO 3 1=1.N2 WRITE(*.106)I READ(*,*)A 3 C02(I)=CMPLX(A,0.) DO 903 1=1,N2 WRITE(*.107)I 903 READ(*,*)S2(I) WRITE(*.330) 330 FORMAT(IX,'Enter the r o o t s of the char, e q u a t i o n , r ( i ) . ' ) NC=N1-NS DO 331 I=1,NC+1 WRITE(*.332)1 Appendix C. Listing of Program Modal 77 332 FORMAT(IX,'Enter complex r o o t , R ( ' , I 2 , ' ) . ' ) 331 READ(*,*)AC(I) DO 800 J=1,NC 800 S1(NS+J)=S2(J) DO 805 J=1,NC D2(J)=(1..0.) DO 801 K=1,NS 801 D2(J)=(S2(J)-S1(K))*D2(J) 805 CONTINUE DO 806 J=1.NC D1(J)=AC(NC+1) DO 802 K=1,NC 802 D1(J)=(S2(J)-AC(K))*D1(J) 806 CONTINUE DO 803 J=1,NC 803 D3(J)=D1(J)/D2(J) DO 804 1=1,NC 804 C01(NS+I)=-D3(I)*C02(I) WRITE(*,109) 109 FORMAT(IX,'To i d e n t i f y the modal parameters of S I , S2 you must') WRITE(*,110)N1 110 FORMAT ( I X . ' s p e c i f y ',12,' i n t e r v a l s . S p e c i f y the l e n g t h o f ) WRITE(*,210) 210 FORMAT(IX.'each h a l f - c y c l e i n t e r v a l as a m u l t i p l e of T.') DO 5 1=1,NI WRITE(*,111)1 Appendix C. Listing of Program Modal 111 FORMAT(IX,'Enter CO v a l u e f o r i n t e r v a l ',12) 5 READ(*,*)CO(I) WRITE(*,58) 58 FORMAT(IX,'Subsystem 1') DO 55 1=1,NI 211 FORMAT(13) WRITE(*.112)1 112 FORMAT(IX,'Enter the h a l f - c y c l e d u r a t i o n of i n t e r v a l ',12) 55 READ(*,211)T21(I) WRITE(*,57) 57 FORMAT(IX,'Subsystem 2') DO 56 1=1.N2 WRITE(*.112)1 56 READ(*,211)T22(I) WRITE(*,113) 113 FORMAT(IX.'Enter the sampling p e r i o d . T, i n seconds.*) READ(*,*)T C C W r i t e the system parameters t o the output f i l e , MODAL.LIS. C WRITE(3,300) 300 FORMAT(IX,'Modal Parameter E s t i m a t i o n Program.',/) WRITE(3.301) 301 FORMAT(IX.'Data g e n e r a t i o n parameters:'./) WRITE(3,302)N1 302 FORMAT(IX,'Subsystem 1 has 3 s t a t e s and '.12,' modes.') Appendix C. Listing of Program Modal 79 WRITE(3.303)N2 303 FORMAT(IX,'Subsystem 2 has 2 s t a t e s and ',12,' modes.',/) WRITE(3.304)(I,C01(I),1=1,N1) 304 F0RMAT(2(2X,'Cl(',12.')='.2E14.6)) WRITE(3,305)(I,C02(I),I=1,N2) 305 F0RMAT(2(2X,'C2(*.12,*)='.2E14.6)) WRITE(3.306)(I.S1(I),I=1.N1) 306 F0RMAT(2(2X.'Sl('.12,')='.2E14.6)) WRITE(3.307)(I,S2(I),I=1.N2) 307 FORMAT(2(2X,'S2(',12,')='.2E14.6)) WRITE(3,700)(I.AC(I),I=1.NC+1) 700 FORMAT(2(2X.*R('.12.')='.2E14.6)) WRITE(3,308)(I.C0(I),I=1.N1) 308 FORMAT(5X.'CO f o r i n t e r v a l '.12.'='.E14.6) WRITE(3.58) WRITE(3.309)(I.T21(I),I=1.N1) 309 FORMAT(5X,'Interval #',12.' has a h a l f c y c l e of '.I3,'*T s.') WRITE(3.57) WRITE(3.309)(I,T22(I),1=1,N2) WRITE(3.310)T 310 FORMAT(/.IX,'The sampling p e r i o d i s '.E10.3,' seconds.') WRITE(*.312) C C Random number g e n e r a t o r i n i t i a l i z a t i o n . C 312 FORMAT(IX,'Enter the n o i s e v a r i a n c e . ' ) Appendix C. Listing of Program Modal 80 READ(*,*)VAR WRITE(3,316) VAR 316 FORMAT(IX,'The a d d i t i v e n o i s e ~N(0.',F7.4.')') WRITE(*,313) 313 FORMAT(IX,'Random number g e n e r a t o r i n i t i a l i z a t i o n . ' ) WRITE(*,315) 315 FORMAT(IX.'Enter R(O); (between 0 and 1)') READ(*,314)R 314 F0RMAT(F12.9) RR=9821.*R+0.211327 R=RR-AINT(RR) C C Save the o r i g i n a l c o e f f i c i e n t s s i n c e the v a l u e s w i l l get C updated every i n t e r v a l . C DO 29 1=1,Nl C1(I)=C01(I) IF(I.LE.N2) C2(I)=C02(I) 29 CONTINUE C C C a l c u l a t e <w,exp(st)>. C DO 6 1=1,Nl DO 7 J=1.N1 INNER1(J)=EX(T.S1,T21,J.I.N1) IF(J.GT.N2.0R.I.GT.N2)G0 TO 7 Appendix C. Listing of Program Modal 81 INNER2(J)=EX(T,S2,T22,J,I,N2) 7 CONTINUE C C C a l c u l a t e <w,y(t)>, <w.y'(t)>, and <w,y"(t)>. C DO 18 J=l,5 18 TEMP(J)=(0.,0.) DO 8 K=1,N1 TEMP(1)=TEMP(1)-CI(K)*INNER1(K) TEMP(2)=TEMP(2)-CI(K)*INNER1(K)*S1(K) TEMPO) =TEMP(3) -CI (K)*INNER1 (K) *S1(K)**2 IF(K.GT.N2.0R.I.GT.N2) GO TO 8 TEMP(4)=TEMP(4)-C2(K)*INNER2(K) TEMP(5)=TEMP(5)-C2(K)*INNER2(K)*S2(K) 8 CONTINUE CALL RANDOM(R.W) INNY1(I)=REAL(TEMP(1)+W*SQRT(VAR)) CALL RANDOM(R.W) INNYP1(I)=REAL(TEMP(2)+W*SQRT(VAR)) CALL RANDOM(R.W) INYPP1(I)=REAL(TEMP(3)+W*SQRT(VAR)) CALL RANDOM(R.W) IF(I.LE.N2) INNY2(I)=REAL(TEMP(4)+W*SQRT(VAR)) CALL RANDOM(R.W) IF(I.LE.N2) INNYP2(I)=REAL(TEMP(5)+W*SQRT(VAR)) Appendix C. Listing of Program Modal C C a l c u l a t e the <w,y(t-iT)> terms and s t o r e them i n A C DO 9 K=1,N1 TEMP(1)=(0.,0.) TEMP(4)=(0.,0.) DO 10 J=1,N1 TEMP(1)=TEMP(1)+C1(J)*CEXP(-FL0AT(K)*S1(J)*T)*INNER1(J) IF(I.GT.N2.0R.K.GT.N2.0R.J.GT.N2)G0 TO 10 TEMP(4)=TEMP(4)+C2(J)*CEXP(-FL0AT(K)*S2(J)*T)*INNER2(J) 10 CONTINUE CALL RANDOM(R.W) Al(I,K)=REAL(TEMP(1)+W*SQRT(VAR)) CALL RANDOM(R.W) IF(K.LE.N2.AND.I.LE.N2) A2(I,K)=REAL(TEMP(4)+W*SQRT(VAR)) 9 CONTINUE C C Update the c o e f f i c i e n t s f o r the next i n t e r v a l . C NT=N1+2*T21(I) N=N2+2*T22(I) IF(N.GT.NT) NT=N I F ( I . E Q . l ) T1=NT TNT=NT*T CALL UPDATE(S2,C2,CC,CO.TNT.N2,0,I) DO 420 J=1,NC 420 CC(J)=-D3(J)*C2(J) Appendix C. Listing of Program Modal 83 CALL UPDATECSl,CI,CC,CO,TNT,NS,NC,I) DO 421 J=1,NC 421 C1(NS+J)=CC(J) K=I + 1 IF(K.LE.Nl) THEN WRITE(3,311) K 311 FORMAT(IX,'New C parameters f o r i n t e r v a l ',13) WRITE(3,304)(J,C1(J),J=1,N1) WRITE(3,305)(J.C2(J),J=l,N2) END IF 6 CONTINUE C C Sol v e the m a t r i x e q u a t i o n s f o r the modal parameters and C check t o see t h a t the modal i d e n t i t y adds up t o z e r o . C CALL DSLIMP(A1,INNY1,M1.N1.EPS) CALL CHECK(N1,A1,Ml,1.INNY1.EPS) CALL CONST(M0(1).Ml,N1,T,T1.1.CO.COl,S1.0) CALL DSLIMP(A2,INNY2.M4.N2.EPS) CALL CHECK(N2.A2.M4.4,INNY2.EPS) CALL CONST(M0(4),M4,N2.T,T1,4,CO.C02,S2.0) C C Update A m a t r i x f o r f u r t h e r c a l c u l a t i o n s . C DO 11 1=1.NI DO 12 J=N1,2,-1 Appendix C. Listing of Program Modal A 1 ( I , J ) = A 1 ( I , J - l ) 12 IFCI.LE.N2.AND.J.LE.N2) A 2 ( I , J ) = A 2 ( I . J - l ) IFQ.LE.N2) A2(I.1)=-INNY2(I) 11 A1(I,1)=-INNY1(I) C C Solv e f o r the modal parameters of the second and t h i r d C s t a t e s and check the modal i d e n t i t i e s . C CALL DSLIMPCA1.INNYP1,M2.N1.EPS) CALL CHECK(N1.Al,M2.2,INNYP1,EPS) CALL C0NST(MO(2).M2,Nl,T.Tl.2,CO,C01.SI,1) CALL DSLIMP(A1,INYPP1.M3.N1.EPS) CALL CHECK(Nl.Al.M3.3,INYPP1,EPS) CALL C0NST(MO(3).M3.N1,T,T1,3,C0,C01,S1,2) CALL DSLIMP(A2,INNYP2,M5,N2,EPS) CALL CHECK(N2,A2,M5.5.INNYP2.EPS) CALL CONST(M0(5),M5,N2.T,T1,5.CO,C02,S2.1) C C U n i t s t e p response. C WRITE(*,400) 400 FORMAT(IX,"Unit s t e p response g e n e r a t i o n . ' ) WRITE(*,401) 401 FORMAT(IX,'Enter the number of data p o i n t s d e s i r e d . ' ) READ(*,*)N3 DO 19 1=1,N3 Appendix C. Listing of Program Modal 85 GG1=-M0(1) GG2=-M0(4) DEL=(I-1)*T ADD1=Y(C0(1),C01,S1,N1.DEL) ADD2=Y(C0(1).C02.S2.N2.DEL) DO 21 K=1,N1 DEL=(I-1-K)*T ADD1=ADD1+(M1(K)*Y(CO(l),C01,S1,N1.DEL)) IF(K.LE.N2)ADD2=ADD2+(M4(K)*Y(CO(1) ,C02,S2,N2,DEL)) 21 CONTINUE GG1=GG1+ADD1 GG2=GG2+ADD2 G1(I)=GG1 G2(I)=GG2 J=MINO(Nl+l.I) DO 22 K=1.J-1 G1(I)=G1(I)-(M1(K)*G1(I-K)) IF(K.LE.N2)G2(I)=G2(I)-(M4(K)*G2(I-K)) 22 CONTINUE 19 CONTINUE C C P r i n t out C WRITE(*,402)N3 WRITE(3.402)N3 402 FORMAT(IX,'The f i r s t ',13.' v a l u e s of the ste p response are:') Appendix C. Listing of Program Modal 86 WRITE(*,404) WRITE(3,404) 404 FORMAT(IX,'Subsystem 1*) WRITE(*,403)(I,G1(I),I=1,N3) WRITE(3,403) ( I . G K I ) ,I=1,N3) 403 FORMAT(2(2X,'G(',12,')='.E14.6)) WRITE(*,405) WRITE(3.405) 405 FORMAT(IX,'Subsystem 2') WRITE(*.403)(I,G2(I),I=1.N3) WRITE(3,403)(I,G2(I),I=1,N3) PAUSE C C P e r i o d i c u n i t s t e p response. C DO 23 I=1,N1+1 GP1(I)=-0.5*M0(1)+G1(I) GP2(I)=-0.5*M0(4)+G2(I) DO 24 J=1,I-1 GP1(I)=GP1(I)+(M1(J)*G1(I-J)) IF(J.LE.N2)GP2(I)=GP2(I)+(M4(J)*G2(I-J)) 24 CONTINUE 23 CONTINUE C C P r i n t out p e r i o d i c s t e p response. C Appendix C. Listing of Program Modal 87 N4=Nl+l WRITE(*,406)N4 WRITE(3,406)N4 406 FORMAT(IX.'The '.13,' v a l u e s of the p e r i o d i c s t e p response are') WRITE(*i404) WRITE(3,404) WRITE(*.407)(I.GPl(I),I=1,N1+1) WRITE(3.407)(I,GP1(I),I=1,N1+1) 407 F0RMAT(2(2X,'GP(',12.')='.E14.6)) WRITE(*,405) WRITE(3,405) WRITE(*,407)(I.GP2CI),I=1.N2+1) WRITE(3,407)(I.GP2(I).I=1.N2+1) PAUSE C C Square wave response. C WRITE(*,500) WRITE(3.500) 500 FORMAT(IX,'The P e r i o d i c Square Wave Response i s g i v e n by:') WRITE(*,404) WRITE(3.404) C C F i g u r e out the gamma parameters. C SIGMA=0.5 Appendix C. Listing of Program Modal 88 CALL INC0EF(M1.N1.GAM1,SIGMA) WRITE(*,405) WRITE(3.405) CALL INC0EF(M4.N2,GAM4,SIGMA) C C C a l c u l a t e the p e r i o d i c square wave response. C DO 25 I=1,N1+1 RS1(I)=GAM1(1)*GP1(I) DO 26 J=1,N1 IFU.GT.J) RS1(I)=RS1(I)+GAM1(J+1)*GP1(I-J) 26 I F ( I . L E . J ) RS1(I)=RS1(I)-GAM1(J+1)*GP1(N1+1+I-J) IFCI.GT.N2+1) GOTO 25 RS2(I)=GAM4(1)*GP2(I) DO 27 J=1,N2 IFQ.GT. J) RS2(I)=RS2(I)+GAM4(J+1)*GP2(I-J) 27 I F ( I . L E . J ) RS2(I)=RS2(I)-GAM4(J+1)*GP2(N2+1+I-J) 25 CONTINUE C C P r i n t out the p e r i o d i c square wave response. C WRITE(*,408)N4 WRITE(3,408)N4 408 FORMAT(IX,'The ',13,' v a l u e s of the per. square wave re s p . are*) WRITE(*,404) WRITE(3.404) Appendix C. Listing of Program Modal 89 WRITE(*,409)(I.RSl(I),I=1.N1+1) WRITE(3,409)(I,RS1(I),1=1,N1+1) 409 F0RMAT(2(2X, 'RSC ,12. ') = ' .E14.6)) WRITE(*,405) WRITE(3,405) WRITE(*,409)(I,RS2(I),I=1,N2+1) WRITE(3,409)(I,RS2(I),1=1,N2+1) PAUSE C C C a l c u l a t e the F o u r i e r c o e f f i c i e n t s C WRITE(3.410) WRITE(*.410) 410 FORMAT(IX.*FFT c o e f f i c i e n t s of G p ( t ) . ' ) WRITE(3.404) WRITE(*.404) CALL SETFFT(GP1.N1+1,T) WRITE(3,405) WRITE(*,405) CALL SETFFTCGP2.N2+1.T) STOP END C C SUBROUTINE INCOEF(M.N,GAM,SIGMA) C Appendix C. Listing of Program Modal C T h i s s u b r o u t i n e s e t s up a m a t r i x of s c a l e d modal parameters C and i n v e r t s i t t o f i n d the gamma parameters t h a t r e l a t e the C u n i t p e r i o d i c square wave response t o the p e r i o d i c u n i t s t e p C response. C C V a r i a b l e s . C C M- the modal parameters. C GAM- the gamma parameters. C MAT- the s c a l e d modal m a t r i x . C MAT2- the i n v e r t e d m a t r i x . C SIGMA- sigma=(N+l)T. C COND- checks the success of the m a t r i x i n v e r s i o n . C I . J - c o u n t e r s . C N- the o r d e r of the system. C REAL M(10),GAM(11).MAT(ll.ll),MAT2(11,11) REAL SIGMA,COND INTEGER N . I . J C C Set up the s c a l e d modal m a t r i x . C DO 1 1=1,N+l DO 2 J=1,N+1 IF(I.EQ.J) MAT(I,J)=SIGMA I F U . L T . J) MAT ( I ,J)=-M(N+1 + I-J) * SIGMA Appendix C. Listing of Program Modal IF(I.GT.J) MAT(I,J)=M(I-J)*SIGMA 2 CONTINUE 1 CONTINUE C C I n v e r t the m a t r i x . C CALL FINV(N+1.11,MAT,11,MAT2,COND) IF(COND.EQ.O.O) WRITE(3,100) 100 FORMAT(IX,'Inversion has f a i l e d ! ! ! ! ! ' ) C C E x t r a c t the gamma parameters from the i n v e r t e d m a t r i x . C DO 3 1=1.N+l 3 GAM(I)=MAT2(I,1) WRITE(*,101)(I.GAM(I),1=1.N+l) WRITE(3,101)(I,GAM(I),1=1.N+l) 101 F0RMAT(2(2X.'GAMMA(',12,')=',E14.6)) RETURN END C C SUBROUTINE CHECK(N,A,MOD,M,INNY.EPS) C C T h i s s u b r o u t i n e e v a l u a t e s the modal i d e n t i t y and sees t h a t C i t adds up t o z e r o . In t h i s way you can be sure t h a t the C parameters are c o r r e c t . Appendix C. Listing of Program Modal C C V a r i a b l e s . C C A- the m a t r i x of de l a y e d output i n n e r product terms. C MOD- the modal parameters. C INNY- the <w,y(t)>, <w,y'(t)>, or <w,y"(t)> terms. C CHECK- the c a l c u l a t e d v a l u e of the modal i d e n t i t y . C EPS- the degree of accuracy of the s o l u t i o n of the modal C parameters. C I - a counter C N- the degree of the system. C M- an i n d i c a t o r of what s t a t e of what subsystem. C REAL A(10.10),MOD(10).INNY(IO).CHECK DOUBLE PRECISION EPS INTEGER I.M.N WRITE(3.100)M WRITE(3.101)(I.MOD(I).I=1.N) WRITEC3.102) WRITE(*,100)M 100 FORMAT(IX,'The modal parameters, M'.I2,' are:') WRITE(*,101)(I.MOD(I).1=1,N) 101 FORMAT(2(2X.'(',12,')='.E14.6)) WRITE(*,102) 102 FORMAT(IX,'Checking!') WRITE(*,600)EPS Appendix C. Listing of Program Modal WRITE(3.600)EPS 600 FORMAT(IX,'EPS adds up t o \E14.6) C C E v a l u a t e the modal i d e n t i t y . C DO 2 J=1,N CHECK=INNY(J) DO 1 1=1,N 1 CHECK=CHECK-(MOD(I)*A(J,I)) WRITE(3,103)CHECK WRITE(*,103)CHECK 103 FORMAT(IX,'The modal i d e n t i t y adds up t o , '.E14.6,/) 2 CONTINUE RETURN END C C FUNCTION EX(T,S,T2,J.I.N) C C T h i s f u n c t i o n e v a l u a t e s the i n n e r p r o d u c t , <w,exp(st)> C when w(t) i s a s a l ( l . x ) f u n c t i o n of d u r a t i o n 2 * T 2 ( I ) . C C V a r i a b l e s . C C S- the a r r a y of complex e x p o n e n t i a l s C EX- dummy v a r i a b l e Appendix C. Listing of Program Modal 94 C T- sampling p e r i o d C MID- the midpoint of w(t) C BEGIN- the b e g i n n i n g of w(t) C END- the endpoint of w(t) C T2- a r r a y of the h a l f - c y c l e d u r a t i o n s of w(t) C I.J.N- c o u n t e r s C COMPLEX S(5).EX REAL T.MID,BEGIN.END INTEGER T2(5) . I . J . N BEGIN=FLOAT(N) MID=FL0AT(N+T2(I)) END=FL0AT(N+2*T2(I)) EX=2.*CEXP(MID*S(J)*T)-CEXP(BEGIN*S(J)*T) EX=(EX-CEXP(END*S(J)*T))/S(J) RETURN END C C SUBROUTINE UPDATE(S,C,CC.CO.TNT.N.NC.M) C C T h i s s u b r o u t i n e w i l l update the C c o e f f i c i e n t s a f t e r C an i n p u t change. C C V a r i a b l e s . C Appendix C. Listing of Program Modal c S- complex e x p o n t i a l s c B- ac c u m u l a t i n g v a r i a b l e a r r a y c X- dummy v a r i a b l e f o r s t o r i n g v a l u e s of powers of S c C- c o e f f i c i e n t s c CO- CO c o e f f i c i e n t s c T- sampling p e r i o d c I , J - c o u n t ers c N- o r d e r of system c T2- h a l f - c y c l e d u r a t i o n of i n t e r v a l s c M- i n t e r v a l number c T l - time e l a p s e d s i n c e the b e g i n n i n g of the f i r s t i n t e r v a l c COMPLEX S(10),B(10),X.C(10),CC(10) REAL CO(10).TNT INTEGER I.J.N.M.NC INTEGER CI.N3 C C Set up the Vandermonde system t o be s o l v e d f o r the C new c o e f f i c i e n t s . C DO 1 1=1,N B(I)=(0.,0.) I F ( I . E Q . l ) B(I)=C0(M) DO 2 J=1.N+NC X=S(J)**(I-1) IF(AIMAG(S(J)).EQ.O.) X=CMPLX(REAL(X),0.) Appendix C. Listing of Program Modal 96 2 B(I)=B(I)+X*C(J)*CEXP(S(J)*TNT) IF(NC.EQ.O) GO TO 1 DO 4 J=1.NC X=S(N+J)**(I-1) IF(AIMAG(S(N+J)).EQ.O.) X=CMPLX(REAL(X).0.) 4 B ( I ) = B ( I ) - C C ( J ) * X 1 I F ( I . E Q . l ) B(I)=B(I)-C0(M+1) CALL VANDM(S,B,N) C C E l i m i n a t e any unwanted imaginary p a r t s of the new C c o e f f i c i e n t s . C DO 3 1=1.N IF(AIMAG(S(I)).EQ.O.) B(I)=CMPLX(REAL(B(I)),0.) 3 C ( I ) = B ( I ) RETURN END C C SUBROUTINE CONST(MO.MOD,N.T.Tl.I.CO,C.S.DERIV) C C T h i s s u b r o u t i n e c a l c u l a t e s the modal o f f s e t c o n s t a n t , C MO. C C V a r i a b l e s . C Appendix C. Listing of Program Modal 97 C C- c o e f f i c i e n t s C S- complex e x p o n e n t i a l s C MO- the modal co n s t a n t term, MO C MOD- a r r a y of modal parameters C T- the sampling p e r i o d C CO- CO c o n s t a n t s C N- the o r d e r of the system C T l - i n t e r v a l d u r a t i o n C I . J - c o u n t e r s C DERIV- the s t a t e of the system C DEL- the d e l a y of each term C COMPLEX C(10),S(10) REAL MO.MOD(IO),T,C0(10) INTEGER N.Tl,I.J,DERIV,DEL M0=W0Y(C0,N,T1,T,C,S,0,DERIV) DO 1 J=1.N DEL=J IFU.NE.1.AND.I.NE.4) DEL=J-1 1 MO=MO+M0D(J)*WOY(CO.N.Tl,T.C.S.DEL,0) MO=MO/(FLOAT(T1-N)*T) WRITE(*,100)I.M0 WRITE(3.100)I.M0 100 FORMAT(IX,'The M',I2,' co n s t a n t i s '.E14.6,'.'.//) PAUSE RETURN Appendix C. Listing of Program Modal 98 END C C FUNCTION WOY(CO,N,Tl,T,C,S,DEL,DERIV) C C T h i s f u n c t i o n c a l c u l a t e s the <w,y(t-itd)> terms f o r C the MO c o n s t a n t c a l c u l a t i o n . C COMPLEX C(10),S(10).WHY.X REAL CO(10).T.WOY INTEGER N.Tl.DEL.I.DERIV WHY=CO(1)*FLOAT(T1-N)*T IF(DERIV.NE.O) WHY=(0..0.) DO 1 1=1.N X=C(I)*(S(I)**DERIV)*WO(S(I),N,T1,T)*CEXP(-FL0AT(DEL)*S(I)*T) 1 WHY=WHY+X WOY=REAL(WHY) RETURN END C C FUNCTION WO(S.N.Tl.T) C C T h i s f u n c t i o n c a l c u l a t e s the <w,exp(st)> terms f o r C the MO term c a l c u l a t i o . C Appendix C. Listing of Program Modal COMPLEX S.WO REAL T INTEGER N.Tl WO=(CEXP(FLOAT(T1)*S*T)-CEXP(FLOAT(N)*S*T)) /S RETURN END C C FUNCTION Y(CO,C,S,N,TAU) C C T h i s f u n c t i o n e v a l u a t e s y ( t ) . C COMPLEX C(10),S(10),YC REAL TAU.CO INTEGER N,I YC=CMPLX(CO,0.) DO 1 1=1.N 1 YC=YC+(C(I)*CEXP(S(I)*TAU)) Y=REAL(YC) IF(TAU.LT.O.) RETURN Y=2.*Y RETURN END C C SUBROUTINE RANDOM(R.W) Appendix C. Listing of Program Modal 100 C C Th i s s u b r o u t i n e generates a complex random number C where i t s r e a l and imaginary p a r t s b o t h ~N(0,1). C REAL R.Rl.Z.Y COMPLEX W Rl=9821.*R+0.211327 R1=R1-AINT(R1) Z=SQRT(-2.*AL0G(R))*C0S(2.*3.14*R1) Y=SQRT(-2.*AL0G(R))*SIN(2.*3.14*R1) W=CMPLX(Z,Y) R=R1 RETURN END C C SUBROUTINE SETFFT(DATA,N,T) C C T h i s s u b r o u t i n e prepares d a t a v a l u e s f o r use by C a FFT r o u t i n e . T h i s s u b r o u t i n e was adapted from C a book by Marple. C COMPLEX X(64),W(64) REAL DATA(10),T INTEGER NF.NEXP.I.J DO 4 1=1.5 Appendix C. Listing of Program Modal 101 4 I F ( ( 2 * N ) . L E . ( 2 * * 1 ) ) GO TO 5 5 NEXP=I DO 1 1=0,2,2 DO 2 J=1,N 2 X((I*N)+J)=CMPLX(DATA(J),0.) DO 3 J=1,N 3 X(((I+1)*N)+J)=CMPLX(-DATA(J),0.) 1 CONTINUE NF=2**NEXP WRITE(*,102)(I.X(I),I=1.NF) 102 F0RMAT(2(2X,'X('.12.')=',2E14.6)) CALL PREFFT(NF.O.NEXP.W) CALL FFT(NF,0,T.NEXP.W.X) WRITE(3,100) WRITE(*,100) 100 FORMAT(IX,'Fourier C o e f f i c i e n t s : Modal f u n c t i o n . ' ) WRITE(3,101)(I,X(I),I=1.NF) WRITE(*,101)(I,X(I),I=1,NF) 101 FORMAT(2(2X,'('.12.')=',2E14.6)) RETURN END C.2 S u b r o u t i n e V a n d m L i s t i n g . SUBROUTINE VANDM(X.B.N) C Appendix C. Listing of Program Modal 102 C S o l v e s the system of l i n e a r e q uations Vz=b C where V i s a Van der Monde m a t r i x . C COMPLEX X(10),B(10) INTEGER I.K.N DO 1 K=1.N-1 DO 2 I=N,K+1.-1 2 B ( I ) = B ( I ) - X ( K ) * B ( I - 1 ) 1 CONTINUE DO 3 K=N-1,1.-1 DO 4 I=K+1.N 4 B ( I ) = B ( I ) / ( X ( I ) - X ( I - K ) ) DO 5 I=K.N-1 5 B(I)=B(I)-B(I+1) 3 CONTINUE RETURN END Appendix D Listings of Other Subroutines D . l Listing of the Subroutines D S L I M P and F I N V The following two subroutines were taken from the University of British Columbia, M T S - G mainframe system and are described [14]. Some changes were made to both subroutines. D . l . l Listing of Subroutine D S L I M P SUBROUTINE DSLIMP(AA,BB,XX,N.EPS) IMPLICIT REAL*8 (A-H.O-Z) DIMENSION A(10.10).T(10.10),B(10),X(10),RZ(10),IPS(10) REAL AA(IO.IO),BB(10),XX(10) DOUBLE PRECISION EPS DO 13 1=1,N DO 11 J=1,N 11 A(I,J)=DBLE(AA(I.J)) 13 B(I)=DBLE(BB(I)) L=10 LT=1 ITMAX=10 IF(LT.NE.l) GO TO 10 103 Appendix D. Listings of Other Subroutines 104 CALL LRD(N,N,L,A,IPS.L,T) CALL DETM(N,IPS,L,T,DET,JEXP) 10 CALL DBS(N,1,L.B,X,IPS,L,T) XN0RM=O.DO DO 1 1=1.N 1 XN0RM=DMAX1(XN0RM.DABS(X(I))) IF(XNORM.LE.O.DO) RETURN EPS=EPS*XNORM ZX=l.D+60 LD=0 DO 2 LL=1.ITMAX DO 3 1=1.N DSUMM=0.0 DO 4 K=1,N DA=A(I,K) DX=X(K) 4 DSUMM=DSUMM+DA*DX DSUMM=B(I)-DSUMM 3 RZ(I)=DSUMM CALL DBS(N.l.L.RZ.RZ.IPS.L.T) ZN0RM=O.DO DO 5 1=1.N ZN0RM=DMAX1(ZNORM.DABS(RZ(I))) 5 X(I)=X(I)+RZ(I) DO 12 1=1,N 12 XX(I)=SNGL(X(I)) Appendix D. Listings of Other Subroutines 105 IF(ZNORM.GT.ZX) GO TO 50 IF((ZNORM-EPS).LT.O.DO) GO TO 60 ZX=ZNORM GO TO 2 50 IF(ZNORM.GT.10.D0*ZX) GO TO 70 LD=LD+1 IF(LD.GE.3) GO TO 70 2 CONTINUE LL=LL-1 WRITE(*,250) WRITE(3,250) GO TO 71 70 WRITE(*,247) WRITE(3,247) 71 EPS=-ZNORM NITER=LL RETURN 60 EPS=ZNORM NITER=LL RETURN 250 FORMAT(/'0***ITERATIVE IMPROVEMENT DID NOT CONVERGE'/) 247 F0RMAT(/'O***ITERATIVE IMPROVEMENT IS DIVERGING'/) END SUBROUTINE LRD(N,N1,NDIMA,A.IP,NDIMT,T) IMPLICIT REAL*8(A-H.0-Z) DIMENSION A(IO.IO),T(10.10),IP(10) Appendix D. Listings of Other Subroutines 106 DO 7 J=1.N DO 7 1=1.N 7 T ( I . J ) = A ( I . J ) IP(N)=1 DO 6 K=1.N IF(K.EQ.N) GO TO 5 KP1=K+1 M=K DO 1 I=KP1,N IF(DABS(T(I,K) ) .GT.DABS(T(M,K))) M=I 1 CONTINUE IP(K)=M IF(M.NE.K) IP(N)=-IP(N) TEMP=T(M.K) T(M,K)=T(K.K) T(K.K)=TEMP IF(TEMP.EQ.O.DO) GO TO 5 DO 2 I=KP1.N 2 T(I,K)=-T(I.K)/TEMP DO 4 J=KP1.N TEMP=T(M.J) T(M.J)=T(K.J) T(K,J)=TEMP IF(TEMP.EQ.O.DO) GO TO 4 DO 3 I=KP1.N 3 T(I.J)=T(I.J)+T(I,K)*TEMP Appendix D. Listings of Other Subroutines 4 CONTINUE 5 IF(T(K.K):EQ.O.DO) IP(N)=0 6 CONTINUE RETURN END SUBROUTINE DETM(N,IP,NDIMT,T,DET,JEXP) IMPLICIT REAL*8 (A-H.O-Z) DIMENSION T(IO.IO),IP(10) DET=IP(N) JEXP=0 DO 1 1=1.N TEMP=T(I.I) IF(DABS(DET).LE.1.D15) GO TO 2 DET=DET*1.D-15 JEXP=JEXP+15 GO TO 3 2 IF(DABS(DET).GE.1D-15) GO TO 3 DET=DET*1.D+15 JEXP=JEXP-15 3 DET=DET*TEMP 1 CONTINUE RETURN END SUBROUTINE DBS(N.NSOL.NDIMBX.B.X.IP,NDIMT,T) IMPLICIT REAL*8(A-H.0-Z) DIMENSION T(IO.IO),X(10),B(10) Appendix D. Listings of Other Subroutines 108 INTEGER IP(IO).K.N.NMl.KPl.KB NM1=N-1 DO 2 1=1.N 2 X ( I ) = B ( I ) IF(N.EQ.1)G0 TO 9 DO 10 K=1.NM1 KP1=K+1 M=IP(K) TEMP=X(M) X(M)=X(K) X(K)=TEMP DO 7 I=KP1.N 7 X(I)=X(I)+T(I.K)*TEMP 10 CONTINUE DO 8 KB=1.NM1 KM1=N-KB K=KM1+1 X(K)=X(K)/T(K.K) TEMP=-X(K) DO 8 1=1,KM1 8 X(I)=X(I)+T(I,K)*TEMP 9 X(1)=X(1)/T(1.1) RETURN END Appendix D. Listings of Other Subroutines 109 D.1.2 Listing of Subroutine FINV SUBROUTINE FINV(N, NDIMT.Tl,NDIMA,A,COND) C C C a l c u l a t e s the i n v e r s e of a m a t r i x . C DIMENSION A ( 1 1 , 1 1 ) . I P ( l l ) , T 1 ( 1 1 , 1 1 ) REAL*8 DSQRT.CSUMA.CSUMB.TA DO 30 J=1.N DO 30 1=1,N 30 A ( I , J ) = T 1 ( I , J ) IEXP=0 IF(N.EQ.l) GO TO 1991 DEt=1.0 KSW=1 GO TO 260 45 CSUMA=CSUMB DO 199 K=1,N AMAX=ABS(A(K,K)) IMAX=K IF(K.EQ.N) GO TO 65 50 KP=K+1 DO 60 I=KP,N AIK=ABS(A(I,K)) IF(AIK.LE.AMAX) GO TO 60 55 AMAX=AIK Appendix D. Listings of Other Subroutines 110 IMAX=I 60 CONTINUE 65 IFCAMAX.EQ.O.) GO TO 300 IP(K)=IMAX IF(K.EQ.IMAX) GO TO 100 DET=-DET 100 DET=DET*A(IMAX,K) QDET=ABS(DET) IF(QDET.LT.1.E15) GO TO 701 DET=DET*1.E-15 IEXP=IEXP+15 701 IF(QDET.GT.1.E-15) GO TO 750 DET=DET*1.E15 IEXP=IEXP-15 750 T=1./A(IMAX,K) A(IMAX.K)=A(K.K) A(K.K)=-1.0 DO 1999 1=1,N A(I,K)=-A(I.K)*T 1999 CONTINUE DO 144 J=1.N IF(J.EQ.K) GO TO 144 TEMP=A(IMAX,J) IF(K.EQ.IMAX) GO TO 140 A(IMAX,J)=A(K.J) 75 A(K,J)=TEMP Appendix D. Listings of Other Subroutines 140 A(K,J)=TEMP*T DO 109 1=1,N IF(I.EQ.K) GO TO 109 A(I,J)=A(I.J)+TEMP*A(I,K) 109 CONTINUE 144 CONTINUE 199 CONTINUE NM1=N-1 DO 250 KK=1.NM1 K=N-KK 210 J=IP(K) IF(J.EQ.K) GO TO 250 220 DO 225 1=1,N T=A(I,J) A(I,J)=A(I,K) 225 A(I,K)=T 250 CONTINUE KSW=2 260 CSUMB=O.ODO DO 271 J=1,N DO 270 1=1,N TA=A(I,J) CSUMB=CSUMB+TA*TA 270 CONTINUE 271 CONTINUE GO TO (45,275), KSW Appendix D. Listings of Other Subroutines 112 275 FN=N COND=DSQRT(CSUMA*CSUMB)/FN RETURN 300 WRITE(*,310) K.AMAX 310 FORMAT(IX,'STEP',13,'PIVOT =',E18.8,'INVERSION STOPPED',/) DET=0.0 IEXP=0 COND=0.0 RETURN 1991 IF(A(1,1).EQ.O.) GO TO 1992 DET=A(1.1) A ( l . l ) = l . / A ( l , l ) C0ND=1. RETURN 1992 K=l AMAX=0. GO TO 300 END D .2 L i s t i n g o f S u b r o u t i n e F F T The subroutine F F T was obtained from [12]. The subroutine as cited here combines the two subroutines P R E F F T and F F T . C These two programs s e t up the complex e x p o n e n t i a l t a b l e (PREFFT) and C compute the d i s c r e t e - t i m e F o u r i e r s e r i e s of an a r r a y of complex data C samples u s i n g a d e c i m a t i o n - i n - f r e q u e n c y f a s t F o u r i e r t r a n s f o r m (FFT) Appendix D. Listings of Other Subroutines 113 C a l g o r i t h m . C SUBROUTINE PREFFT (N.MODE.NEXP,W) C C Input Parameters: C C N Number of data samples t o be pro c e s s e d (integer-must be a C power of two) C MODE - Set t o 0 f o r d i s c r e t e - t i m e F o u r i e r s e r i e s (Eq. 2.C.l) o r C 1 f o r i n v e r s e (Eq. 2.C.2) C C Output Parameters: C C NEXP - I n d i c a t e s power-of-2 exponent such t h a t N=2**NEXP . C W i l l be s e t t o -1 t o i n d i c a t e e r r o r c o n d i t i o n i f N C i s not a power of 2 ( t h i s i n t e g e r used by sub. FFT) C W - Complex e x p o n e n t i a l a r r a y C C Notes: C C E x t e r n a l a r r a y W must be dimensioned .GE. N by c a l l i n g program. C COMPLEX W(l).C1.C2 NEXP=1 NT=2**NEXP IF (NT .GE. N) GO TO 10 Appendix D. Listings of Other Subroutines 114 NEXP=NEXP+1 GO TO 5 10 IF (NT .EQ. N) GO TO 15 NEXP=-1 RETURN 15 S=8.*ATAN(1.)/FLOAT(NT) C1=CMPLX(C0S(S).-SIN(S)) IF (MODE .NE. 0) C1=C0NJG(C1) C2=(l.,0.) DO 20 K=1,NT W(K)=C2 20 C2=C2*C1 RETURN END SUBROUTINE FFT (N.MODE.T.NEXP.W.X) C C Input Parameters: C C N.MODE.NEXP.W - See parameter l i s t f o r s u b r o u t i n e PREFFT C T - Sample i n t e r v a l i n seconds C X - A r r a y of N complex data samples, X ( l ) t o X(N) C C Output Parameters: C Appendix D. Listings of Other Subroutines 115 C X - N complex t r a n s f o r m v a l u e s r e p l a c e o r i g i n a l d a t a samples C indexed from k=l t o k=N, r e p r e s e n t i n g the f r e q u e n c i e s C ( k - l ) / N T h e r t z C C Notes: C C E x t e r n a l a r r a y X must be dimensioned .GE. N by c a l l i n g program. C COMPLEX X(l),W(1),C1.C2 MM=1 LL=N DO 70 K=1,NEXP NN=LL/2 JJ=MM+1 DO 40 I=1.N,LL KK=I+NN C1=X(I)+X(KK) X(KK)=X(I)-X(KK) 40 X(I)=C1 IF (NN .EQ. 1) GO TO 70 DO 60 J=2,NN C2=W(JJ) DO 50 I=J,N,LL KK=I+NN C1=X(I)+X(KK) X(KK)=(X(I)-X(KK))*C2 Appendix D. Listings of Other Subroutines 50 X(I)=C1 60 JJ=JJ+MM LL=NN MM=MM*2 70 CONTINUE NV2=N/2 NM1=N-1 J=l DO 90 1=1,NM1 IF ( I .GE. J) GO TO 80 C1=X(J) X(J)=X(I) X(I)=C1 80 K=NV2 85 IF (K .GE. J) GO TO 90 J=J-K K=K/2 GO TO 85 90 J=J+K IF (MODE .EQ. 0) S=T IF (MODE .NE. 0) S=l./(T*FLOAT(N)) DO 100 1=1,N 100 X ( I ) = X ( I ) * S RETURN END 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064893/manifest

Comment

Related Items