Obtaining mode shapes through the Karhunen-Lo `eve expansion for distributed-parameter linear systems

. The Karhunen-Lo`eve expansion is a powerful spectral technique for the analysis and synthesis of dynamical systems. It consists in decomposing a spatial correlation matrix, which can be obtained through numerical or physical experiments. The decomposition produces orthogonal eigenfunctions or proper orthogonal modes, and eigenvalues that provide a measurement of how much energy is contained in each mode. The relation between KL modes and mode shapes of linear vibrating systems has already been derived and demonstrated for two and three dofs mass-spring-damper systems. The purpose of this paper is to extend this investigation to more complex distributed-parameter linear systems. A plane truss and a simply supported plate subjected to impulsive forces, commonly used in modal analysis are studied. The resulting KL modes are compared to the analytical mode shapes. Damping and random noise effects in the procedure performance are evaluated. Two methods for indirectly obtaining natural frequencies are also presented.


Introduction
The Karhunen-Loève (KL) expansion or decomposition initially appeared in the signal processing literature, where it is customary called principal components analysis (PCA). It is a powerful statistical tool to perform data analysis and compression. It soon found many useful applications such as voice and image recognition. In Mechanical Engineering, where it is also known as the proper orthogonal decomposition (POD), it was firstly employed to uncover coherent structures in turbulent flow fields [7]. Coherent structures can be defined as recurrent spatial forms that are energy dominant [4]. Since then, it has been consistently developed and applied to many fluid dynamics problems [11][12][13].
Interestingly enough, only recently has the KL expansion attracted the attention of structural dynamicists seeking an alternative approach to obtaining reducedorder models of linear and nonlinear dynamical sys-tems [6,8,14]. The objective of reducing a mathematical model is to obtain a simpler one that not only still has a good degree of predictive capability but also is suitable for its intended application. This can be, for instance, the design of a real-time feedback controller that is not too computationally intensive, or simply the performance of a parameter analysis.
The KL method is primarily a statistical procedure. One initially supposes that the observed system dynamics can be modelled as a second-order ergodic stochastic process. The method consists then in constructing a spatial autocorrelation tensor from data obtained through numerical or physical experiments and performing its spectral decomposition. Since it deals only with data, there is no distinction between linear or nonlinear systems and it can even be implemented, in the physically obtained data case, without any previous knowledge about the mechanical characteristics of the system. The autocorrelation tensor is by definition Hermitian and positive semi-definite. Therefore, its decomposition provides a set of orthogonal eigenfunctions (called proper orthogonal modes, POMs, or empirical eigenmodes) and nonnegative real eigenvalues (or proper orthogonal values, POVs, or empirical eigenvalues). These POMs can then be used as a basis for the dynamics projection and in the construction of a reduced-order model through the retention of a finite number of them. An important property of the expansion is that the magnitude of a POV is a measure of the energy contained in the respective POM. Furthermore, the expansion is in a sense optimal, meaning that no other linear decomposition can better reproduce the particular dynamics which generated the POMs with the same number of modes. There are two ways of constructing the expansion: the direct and the snapshot methods. Each one has its domain of application, as will be later discussed, although this paper deals only with the direct method.
The relation between KL modes and mode shapes for general vibrating structural systems described by second-order ordinary differential equations has been derived in [2]. The purpose of this paper is to extend that discussion to more complex linear distributedparameter system. For this purpose, a plane truss and a simply supported plate have been chosen as the objects of study. They have both been excited by means of impulsive forces as usually is done in modal analysis. The POMs obtained were then compared to the intrinsic physical mode shapes and several important questions addressed.

The Karhunen-Loève decomposition
Let a dynamical system be governed by equations whose solutions give a flow, i.e., a function of time and space that describes the evolution of a particular state. Let the flow be denoted by u(x, t) and defined on a spatial domain x ∈ D and a time domain t ∈ [0, ∞).

Main hypothesis
In order to define the autocorrelation tensor for the flow, one must model it as a second-order stochastic process. However, it is desirable to avoid the mathematical description of the sample space, σ-algebra, and probability measure associated with the flow. A great advantage of this methodology is that this description is unnecessary, though two additional assumptions are needed: the flow is supposed to be strict-sense timestationary and ergodic [9].
Let v(x, t) define the deviation from the mean flow, Hence, v(x, t) is a stochastic process with zero mean and consequently its autocorrelation tensor equals its autocovariance tensor [10]. If v(x, t) is real, then the two-point spatial autocorrelation function is defined by the dyadic product A final assumption regarding the flows u(x, t) and v(x, t) is that they are continuous in quadratic mean, implying the continuity of the autocorrelation tensor R(x, x ) in its spatial domain [9].

Model reduction
In order to find a reduced-order flow model that still reveals the main features contained in the dynamics, one can search for an expansion of the form with i.e., the modes are uncorrelated, and meaning that the set {ψ n } is orthonormal and N = dim v. Inserting Eq. (3) into Eq. (2) and using the relation given by Eq. (4), one obtains Since by definition and according to our assumptions, R(x, x ) is positive, semi-definite, Hermitian and continuous on D, Mercer's Theorem [9] guarantees the existence and uniqueness of the spectral representation of R(x, x ) given by Eq. (6), where the elements belonging to the orthonormal set {ψ n } are the eigenfunctions of the integral operator with kernel R(x, x ) and the set {λ n } is formed by the corresponding real and nonnegative eigenvalues so that Then, the Karhunen-Loève Theorem [9] states that a continuous second-order stochastic process with auto-covariance tensor K(x, x ) can be expanded in a series analogous to (3) where {ψ n } are the eigenfunctions of the integral operator with kernel K(x, x ) and {λ n } are the corresponding eigenvalues. Since for v(x, t) the autocovariance equals its autocorrelation function, we have thus proved that the expansion stated in Eq. (3) is realizable. The set {ψ n } is formed by the POMs, also called coherent structures.
The original flow u(x, t) can, therefore, be reconstructed with reduced dimension through the truncation of the series Eq. (3) and addition of the mean flow: where the temporal coefficients A n are easily found by projecting the flow onto each POM ψ n , i.e., Finally, the eigenvalue may be written, using the ergodic hypothesis, as indicating that it is a measure of the mean energy contained in each mode. Besides, it can be shown that the total mean energy equals the sum of all eigenvalues, that is, E = n λ n [11].

Pratical construction of POMs
As previously mentioned, there are two practical methods available for the construction of the KL expansion, namely the original direct and the more recent snapshot methods. Both will be briefly presented and the respective advantages discussed, so that later on it will be clear why the direct method constitutes itself in the best choice for this work.

Direct method
In this method, the displacements of a dynamical system are measured or calculated at N locations and labelled u 1 (t), u 2 (t), . . . , u N (t). Sampling these displacements M times, we can form the following M ×N ensemble matrix: . . . . . . . . . . . .
Thus, using the time-stationarity and ergodicity hypothesis, the variation from the mean can be calculated as and the spatial correlation matrix of dimension N × N formed as The POMs are then given by the eigenvectors of R which are orthogonal due to its symmetry. Eigenvalues will provide the POVs. Clearly, the matrix dimension is determined by the number of sampling points N . For a three-dimensional flow, i.e. v(x, t) ∈ R 3 , the number of operations necessary for the diagonalization of R is

Snapshot method
This method was firstly introduced in [11]. It uses the fact that due to the assumed ergodicity, the spatial autocorrelation tensor can be expressed as = lim where v (m) (x) = v(x, mτ) is referred to as a snapshot and τ is the sampling time which should be greater than the correlation time. However, in practice, one would have to deal with a finite number of snapshots. Hence, the kernel R(x, x ) in Eq. (15) would become degenerated and, therefore, would have eigenfunctions which are a linear combination of the snapshots [3,11]: where the coefficients A km are still to be determined. Introducing Eqs (15) and (15) in Eq. (7) these coef-ficients would present themselves as solutions to the eigenequation defined by In the above expression, ·, · defines an inner product on an L p 2 (D) space, where p = dim v. Thus, the determination of the POMs requires the spectral decomposition of a matrix whose dimension is determined by the number of snapshots M . Actually, this procedure requires O(M 3 ) operations [1]. The number of sampling points N would only indirectly enter the calculation through the evaluation of the inner products.
It is thus clear that the direct method should be applied to experimentally gathered data where there is, in general, a rich time history obtained in a relatively small number of locations or to numerically generated data with moderate spatial resolution. Otherwise, in the case of multidimensional simulated flows with high spatial resolutions, the snapshot method is to be preferred.

Relationship between POMs and mode shapes
The first demonstration that the Karhunen-Loève expansion could be used in order to obtain the eigenfunctions of linear operators appeared in [1]. This would imply obtaining mode shapes of a distributed-parameter vibrating system since these are the eigenfunctions of the spatial linear operator subjected to boundary conditions. A more straightforward demonstration was presented in [2] for discrete systems and is briefly reviewed here.
Consider a discretized dynamical system described by Mẅ + Kw = 0. The motion for this system can be expressed as where the w i and ω i are the mass normalized mode shapes and respective natural frequencies [5] and a i , φ i are constants to be determined from initial conditions. Since for an oscillatory system, the mean flow is zero, using the direct method, one can promptly form the ensemble matrix: . . . where Hence, the spatial correlation matrix is Creating an adjusted autocorrelation matrix given bŷ R = RM, it is possible to show that its eigenvectors converge to the modal vectors: Using the orthogonality condition w T i Mw j = δ ij , the above expression simplifies tô Therefore, each term (1/M )w i c T i c j will tend to zero as M → ∞, except for the term (1/M )w i c T i c i which is proportional to w i . This means that the mode shapes are eigenvectors of the adjusted autocorrelation matrix. Furthermore, in the case of a global mass matrix which is equal to or proportional to the identity matrix, the above result would imply that the mode shapes are eigenvectors of the autocorrelation matrix and hence POMs or, alternatively, that the POMs would be converging to the mode shapes when M → ∞. When damping is present, all functions c i (t) → 0 as t → ∞. However, it would still be possible to apply this procedure to structures with very small damping where a great number of oscillations would be observed [2]. The fact that the convergence of the adjusted matrix eigenvectors to the mode shapes occurs as M → ∞ was determinant for the choice of the direct method in the rest of this work. Another interesting point is that in many experimental cases, the mass matrix may be unknown and it would be impossible to calculate the adjusted matrix. This issue will be addressed in the next section.

Modal analysis of a plane truss through the KL expansion
In order to explore the ideas expressed earlier, we have chosen to perform the KL expansion on a plane truss depicted in Fig. 1.

Mode shapes determination
The truss was modelled by the finite element method using simple bar elements and an inconsistent mass matrix. The vertical and horizontal bars have lengths l 1 = 15 cm and l 2 = 20 cm, respectively, as shown in Fig. 1. Cross-sectional area for the elements is 1 cm 2 , mass density is 7860 kg/m 3 and elastic modulus is 200 GPa. Since the truss is clamped in its left extremity, implying that nodes 1 and 6 are fixed, we arrive at a 16 dof system, for other nodes are capable of moving in both horizontal and vertical directions. Fig. 2 presents the truss together with its mass normalized analytical mode shapes and respective natural frequencies.
The KL expansion was then performed for unitary impulsive forces acting on node 10 in both horizontal and vertical directions. The response was measured at all 8 nodes in both directions. The time sampling rate was 0.1 s which is slightly greater than the first natural period for a reason that will be later explained, and the time span was from 0 to 100 s,so that 1000 time samples were available. Since the FEM discretization generated a system with a small number of degrees of freedom, it is ideal for the application of the direct method. The eigenvalue problem was solved through the singular value decomposition, for this is a robust technique that avoids numerical instabilities in the eigenvalue problem solution that may appear due to the optimality property of the expansion. Figure 3 presents the resulting POMs obtained from the correlation matrix R as in Eq. (20), as well as the respective POVs expressed as a percentage of the total energy.
The inconsistent global mass matrix obtained with the FEM discretization though diagonal is not proportional to the identity matrix. In this case and in many others where the mass matrix is block diagonal, it is not possible to assure that the POMs are converging to the analytical mode shapes, since the hypothesis formulated in section 4 will not be true. However, it is clear that at least qualitatively, the POMs represent these mode shapes very well. Naturally, we could have multiplied the correlation matrix by the global mass matrix in order to obtain the adjusted correlation matrixR and then calculate its eigenvectors. This was not performed because one of the objectives here is to verify whether the POMs obtained directly from R can be taken as a good approximation for the discretized model mode shapes. Moreover, the mode shapes for a physical structure are the eigenfunctions of an Hermitian operator (the stiffness operator), and therefore are mutually orthogonal. Hence, the POMs from R are probably a better representation for the real mode shapes than those calculated from the discretized model. Figure 4 presents the resulting POMs from a vertical impulse applied at the same node 10. In this case, there was an increase in the energy contained in the first POM and in the remaining pure flexural modes such as the second, third, and tenth POMs, when compared to the previous case. Furthermore, there were some difficulties in the representation of some mode shapes. The ninth POM, for example, is a poor representation of the fourteenth mode shape and the fifth POM is completely distinct from any of the mode shapes. This result is a consequence of the small excitation provided to some modes due to the restriction of the excitation force to one single direction.
The optimality property of the KL expansion is clearly demonstrated in Fig. 5 where the reconstruction of the horizontal displacement of node 10 is performed using the first one and the first two POMs presented in Fig. 4. As stated previously, a great advantage of the expansion is the fact that the eigenvalues are a direct measure of the energy contained in the related POMs.
Here, the first one and first two eigenvalues account for 97.5% and 99.5%, respectively, of the total energy. One can observe that these first two modes suffice for a good dynamics representation.
On the other hand, this same optimality property means that this is a noncausal procedure where the empirical modes cannot be determined before the system response is known. Moreover, changing the impulse excitation direction will result in a completely different energy distribution among the POMs. Figure 6 shows the POMs obtained for a horizontal impulse applied at node 10.
In this case, we have a more homogeneous energy distribution where the first POM accounts for only 37.8% of it. This results in a necessity of using a greater number of POMs to accomplish an accurate dynamics representation. The first graph in Fig. 7 depicts the same horizontal displacement of node 10 used in the former case. One can note that even using 4 POMs, a good representation is not attained. This happens because for the dynamics resulting from the horizontal impulse application, the first four POMs respond for only 87% of the energy, far from the 99% threshold usually recommended [11]. The dynamics reconstruction is even worse if one attempts to use the first four POMs obtained from the vertical impulse excitation (Fig. 4), as shown in the second graph. Although in that former case these four POMs accounted for more than 99% of  the energy, this is not true for the dynamics resulting from the horizontal impulse. Comparing Figs 4 and 6, one notes that for the present case, those four POMs represent less than 50% of the energy and thus are un-able to accurately represent the horizontal displacement of node 10. Therefore, it is possible to conclude that the KL expansion may not be robust enough to generate a reduced-order model when excitation conditions change dramatically. The third graph also presents a modal reconstruction using the first four analytical mode shapes. The result is not bad because the first three mode shapes correspond qualitatively to the first 3 POMs obtained for the horizontal impulse. Another reconstruction was then performed in the same graph using POMs one, two, three and six which are related to the first four mode shapes. One can observe that this representation is practically indistinguishable from the modal one. However, the evaluation of the error norm between both reconstructions and the original response yields 4, 90 × 10 −3 for the modal case and 4, 61 × 10 −3 for these POMs. This is an interesting result and demonstrates once again how the optimality property works in order to provide empirical modes that better represent the sampled dynamics.
Finally, the effect of structural damping was investigated. Damping was introduced in the system of second-order ODEs through a damping matrix which was proportional to the stiffness matrix, that is C = βK, where β = 1 × 10 −4 . This proportional damping corresponds to modal damping factors ranging from ζ 1 = 0.0057 to ζ 16 = 0.1930 according to the relation ζ i = βω i /2. The truss was once again excited by means of impulsive forces acting on node 10 in both directions and the same sampling conditions were ob-served. The POMs obtained with the KL expansion for this case are presented in Fig. 8.
It is clear that the presence of proportional damping had the effect of concentrating the energy in the first POM because the first mode shape has the smallest damping factor. Moreover, damping made it more difficult for some POMs to approximate the respective mode shapes.

Natural frequencies determination
So far nothing has been mentioned about how to obtain the system natural frequencies from the KL expansion. At a first glance, the expansion provides only the mode shapes, in the case of linear systems. However, this drawback can be overcome by using the POMs as initial conditions and performing a Fast Fourier Transform of the resulting dynamical response. Hence, the first three POMs depicted in Fig. 4 and that correspond to the first, second, and fourth mode shapes were given as initial conditions to the truss in three distinct occasions. The vertical displacement of node 10 was then sampled in 0.001s steps until 1s and the FFT of this signal performed. The result is presented in Fig. 9 where the natural frequencies characteristic peaks for these mode shapes are clearly visible.
Nonetheless, in the real world, it is usually not realizable to give a certain spatial form to a structure.  A more straightforward strategy for finding the natural frequencies is to perform a FFT of the temporal coefficients (also called principal coordinates) of the expan-sion Eq. (8) since, in the linear case, they correspond to the modal amplitudes. Figure 10 presents the results of the FFT performed on the temporal coefficients   calculated from the dynamics projection onto the first three POMs, according to Eq. (9). These POMs appear in Fig. 3 and correspond to the first, third and second mode shapes. Once again, the respective natural fre-quencies peaks are clearly visible. Table 1 compares the natural frequencies obtained using the FFT through both methods with the analytical natural frequencies.
An excellent agreement is observed.   Fig. 11. First twelve POMs for impulsive force applied at the middle of plate.

Dynamics of a uniform rectangular simply supported plate
The equation of free motion for such a plate of dimensions a × b is [5]: where x = (x, y) and with the following mode shapes and natural frequencies: satisfying the orthogonality condition Thus, the plate dynamics for zero initial conditions and an unitary impulsive force applied at sin ω mn t.

Application of the KL decomposition to the dynamics of a simply supported plate
A distributed-parameter system such as this can be viewed as a natural extension of a discrete system with the number of modes tending to infinity. Indeed, since the independent mode shapes for the plate are known, the orthogonality relation Eq. (25) is equivalent to a discrete system with a mass matrix proportional to the  identity matrix. Hence, it is expected that the spatial correlation matrix eigenvectors found through the direct method will converge to the analytical mode shapes.
As an example, an impulsive force was applied to the middle (f 1 = f 1 = 1/2) of a retangular (a = 1.9m, b = 1m) steel plate 5mm thick. In order to compute the analytical response, 256 mode shapes were employed. These particular dimensions assure us that there are no repeated natural frequencies. The spatial domain was discretized in a grid with points spaced by 0.05m. The response was sampled using a 1s time interval (this will be later discussed) which is slightly greater than the first mode oscillation period. The resulting POMs were computed and are shown in Fig. 11. One thousand time samples were needed to attain these results. Notice also that the POVs are normalized, so that their sum is one.
Although the results are good, it is clear that no mode shape with m or n even was encountered. That is, naturally, due to the position of the impulsive force. Mode shapes with a node in the middle of the plate would not be excited and, thus, do not appear in the statistically determined POMs. Shifting the force application point to f 1 = f 2 = 1/4 yields the new POMs depicted in Fig. 12. It is possible to observe some new even mode shapes and the fact that the energy is less concentrated in the first POM.
Shifting the force far away from the center affects negatively the results since the most excited mode shapes are high frequency ones that are seldom of interest. Furthermore, the number of time samples to achieve convergence increases a lot. Fortunately, in a physical experiment, the singularity present in the application of a numerical impulse would be avoided and there should be no difficulty in obtaining low frequencies mode shapes through the use of an impact hammer around the center of the plate.
The robustness of the procedure to noise was also investigated. The data were contaminated by addition of gaussian white noise. The procedure consisted in summing the exact response to numeric generated gaussian white noise previously to calculating the POMs. However, these proved to be quite insensitive to noise. Only when the signal to noise ratio was of O(10) was there noticeable qualitative change in the POMs as shown in     Fig. 13. Even then, the mode shapes are clearly recognizable in the POMs. However, this level of noise is quite high and not to be expected in experimental setups. In this example, the impulsive force was also applied at f 1 = f 2 = 1/4.
The time span of 1000s used so far is rather unrealistic since structural damping will always be present and because of that, free plate vibrations will not be sustained for so long. As mentioned previously, if a structure is lightly damped, then it should be possible to find its mode shapes through the KL decomposition. The presence of damping and consequent decrease in available time span and the necessity of a great number of snapshots or time samples to attain convergence would mean decreasing dramatically the sampling interval. On the other hand, the ergodicity assumption used so far implies the use of uncorrelated snapshots. But we are actually dealing with a deterministic system, meaning that this would, indeed, be impossible. Choosing the sampling time of 1 s (above the first period) was a strategy to minimize this problem. Results showed that increasing the sampling time above this threshold made no difference but decreasing it did. Introducing a modal damping factor ζ mn = 0.01 and sampling the plate dynamics in 0.001 steps until 1s (giving the same 1000 samples) produces the POMs depicted in Fig. 14. Table 2 also presents the Frobenius norm of the difference between the first POM and the first plate mode shape for different sampling intervals and considering the undamped case. It is clear that the error monotonically decreases when the interval is increased until past the first natural period. From then on, there is no established pattern and the error seems to oscillate around the value associated with the first period. Therefore, it is logical to conclude that decreasing sampling time to deal with damping is not a good strategy. Perhaps a better option, specially in an experimental case, would be to apply another impulsive force after some time, thus increasing the available time span for measurements.
Finally, the behavior of a square plate was investigated. This system has infinite pairs of equal frequencies. It was excited by means of an unitary impulse applied to its center. The time span was 1000s sampled at every second. The first four resulting POMs are presented in Fig. 15. It is interesting to note that the KL decomposition was unable to differentiate between mode shapes having the same frequency. Because they are vibrating in phase, statistically they cannot be distinguished and ended up added by the procedure. One result, for instance, is the second POM which is, in reality, the sum of mode shapes w 13 (x, y)+w 31 (x, y). And the fourth POM is the result of w 15 (x, y) + w 51 (x, y) (note that the even modes do not appear due to the applied impulse position, as discussed earlier). This is a great disadvantage for performing modal analysis in unknown systems via the KL decomposition. However, if the objective is to construct a reducedorder model, the method continues to be optimal. Using only the first two POMs, it is possible to practically reconstruct the original dynamics for the plate central point as can be seen in Fig. 16.

Conclusions
This work discussed the use of the Karhunen-Loève expansion as a tool to perform modal analysis of distributed-parameter linear system. It was shown that the proper orthogonal modes obtained with the procedure are a good approximation for the analytical mode shapes of a discretized model with diagonal global mass matrix. A great advantage of using the KL expansion for modal analysis is that the proper orthogonal values provide the contribution of the associated POMs in the global dynamics. Moreover, it was also demonstrated that the convergence capacity of the expansion is little affected by the existence of gaussian noise.
Nevertheless, there are also some disadvantages in its application. First of all, the sampling time should be at least close to the first oscillation period, what possibly will not be known a priori in an experimental case. Moreover, only lightly damped structures can be the object of analysis by this method. Besides that, the expansion cannot differentiate mode shapes with the same frequency which end up added. Finally, there is also the necessity of measuring the structure dynamical response in a fairly large number of points. Therefore, although it is possible that the Karhunen-Loève expansion may replace traditional modal analysis in some applications, it is for now more likely that it will be used as an additional tool, specially when choosing a modal basis for a reduced order model, due to its optimality property.