Spline-Interpolation-Based FFT Approach to Fast Simulation of Multivariate Stochastic Processes

The spline-interpolation-based fast Fourier transform FFT algorithm, designated as the SFFT algorithm, is proposed in the present paper to further enhance the computational speed of simulating the multivariate stochastic processes. The proposed SFFT algorithm first introduces the spline interpolation technique to reduce the number of the Cholesky decomposition of a spectral density matrix and subsequently uses the FFT algorithm to further enhance the computational speed. In order to highlight the superiority of the SFFT algorithm, the simulations of the multivariate stationary longitudinal wind velocity fluctuations have been carried out, respectively, with resorting to the SFFT-based and FFT-based spectral representation SR methods, taking into consideration that the elements of cross-power spectral density matrix are the complex values. The numerical simulation results show that though introducing the spline interpolation approximation in decomposing the cross-power spectral density matrix, the SFFT algorithm can achieve the results without a loss of precision with reference to the FFT algorithm. In comparison with the FFT algorithm, the SFFT algorithm provides much higher computational efficiency. Likewise, the superiority of the SFFT algorithm is becoming more remarkable with the dividing number of frequency, the number of samples, and the time length of samples going up.


Introduction
Monte Carlo technique has widely been employed for simulating the stochastic processes which are either one-dimensional or multidimensional, univariate or multivariate, homogeneous or nonhomogeneous, stationary or nonstationary, and Gaussian or non-Gaussian. The Monte Carlo simulation methods are able to generate the sample functions that accurately provide the probabilistic characteristics of the corresponding stochastic processes. For the simulation of stochastic processes, the following approaches are now available: 1 Mathematical Problems in Engineering autoregressive AR method, such as Mignolet and Spanos 1, 2 , Iannuzzi and Spinelli 3 , Deodatis and Shinozuka 4 , and Novak et al. 5 ; 2 autoregressive moving average ARMA method, for example, Gersch and Yonemoto 6 , Kozin 7 ,Kareem and Li 8 ,and Rossi et al. 9 ; 3 spectral representation SR method, for instance, Shinozuka and Jan 10 , Grigoriu 11 , Deodatis 12 , Grigoriu 13 , and Chen and Letchford 14 .
It is known that among foregoing simulation methods, the SR method has very high demands on both the computer memory and speed. Notwithstanding this, the SR does not have the problem of model selection, as does the AR and ARMA methods. Likewise, it is easy to implement and has high accuracy. Hence, the SR method has yet been receiving increasing attention in simulating the multivariate stochastic processes with the target cross power spectral density CPSD matrices. The primary concept of the SR method dates back to Rice 15 ,Goto and Toki 16 ,and Borgman 17 . For the multidimensional, multivariate, and nonstationary cases, Shinozuka 18,19 first established the SR approach. Subsequently, Deodatis and Shinozuka 20 extended the application of the SR method to the simulation of stochastic waves. Li and Kareem 21 developed a hybrid discrete Fourier transform and digital filtering approach to simulate the multivariate random process. Likewise, Grigoriu 11 compared two different SR models. But, it is worth pointing out that the generated sample functions with resorting to the early algorithm of the SR method by Shinozuka and Jan 10 are not ergodic. In order to generate the ergodic sample functions, several attempts have been made via modifying the existing algorithms. Shinozuka et al. 22 introduced the idea of double-indexing frequencies. But, the obtained sample functions based on the proposed formula are still not ergodic. Deodatis 12 further extended the SR method to simulate the multivariate ergodic stochastic processes. Likewise, the capabilities and efficiency of the proposed algorithm were demonstrated in detail using a one-dimensional trivariate process as an example. On the other hand, the computational efficiency of the early algorithm of the SR method is yet low. In order to cope with this issue, the fast Fourier transform FFT algorithm was introduced into the SR method. Yang 23,24 showed that the FFT algorithm can remarkably enhance the computational efficiency of the SR method and proposed a formula to simulate the random envelop processes. Shinozuka 25 extended the application of the FFT algorithm to the multidimensional cases. Wittig and Sinha 26 showed that the SR method combined with the FFT algorithm seems to be more than an order of magnitude faster than other simulation methods. Further, Li and Kareem 27 used the FFT-based approach to simulate the multivariate nonstationary Gaussian random processes with the prescribed evolutionary spectral description.
Likewise, the multicorrelated stationary random processes, such as the wind velocity or pressure fluctuations, on structures can be transformed into a set of subprocesses through diagonalizing their covariance or CPSD matrices with resorting to either the Cholesky lower or upper triangular or eigenvector decomposition. The eigenvector decomposition offers physically meaningful insight into the process as each eigenvector eigenmode may be characterized on the basis of its spatial distributions. Theoretically, the eigenvector decomposition is based on the Karhunen-Loeve expansion, consequently referred also to as the proper orthogonal decomposition POD 28 . The POD approach has been widely utilized to reduce the dimensions or variables in the large scale systems to enhance the computational speed, such as Holmes et al. 29 , Di Paola and Gullo 30 , Rathinam and Petzold 31 , Chen and Kareem 32 , and Solari and Carassale 33 . In 34 , Carassale and Solari proposed some strategies such as interpolation of the factorizations of the CPSD matrix, cutoff of the spectral representation based on POD and parallelization of the code, which provide significant computational advantages enabling the simulation of large wind fields. Likewise, Ding et al. Engineering   3 35 also introduced interpolation approximation of the factorizations of the CPSD matrix and cutoff of the spectral representation based on POD to the SR algorithm for the simulation of wind velocity fields on large scale structures. The numerical advantage of the POD technique, akin to the modal analysis of structural dynamics, relies mainly on the reduced-order representation via the truncation of the higher eigenmodes associated with smaller eigenvalues. Of course, this reduced-order representation must warrant that the important characteristics of the stochastic processes and related quantities remain unchanged, or the modification resulting from the approximate representation is acceptable. However, the truncation of higher modes may not necessarily work effectively in the case of local responses, implying that there exists a possibility of underestimating the local wind loads and corresponding effects 36 . In the case of using the CPSD matrix-based POD technique, the similar observations made by Chen and Kareem 37 once again underscore the foregoing phenomenon.

Mathematical Problems in
In order to illuminate solely the effectiveness of interpolation in fast simulation of multivariate stochastic processes without POD technique, the spline function employing the nature of both the smoothness and continuance is utilized to reduce the number of the Cholesky decomposition of power spectral matrix by means of interpolation. The splineinterpolation-based FFT algorithm, designated as the SFFT algorithm, then is proposed in the present paper. Therefore, different from 34, 35 , the main purpose of the present study is to evaluate the computational efficiency and accuracy of the SFFT-based SR method with reference to the FFT-based SR method through simulating the multivariate stationary longitudinal wind velocity fluctuations with the phase angles.

Simulation of Multivariate Stochastic Processes
According to Deodatis 12 and Shinozuka 19 , the multivariate stationary stochastic process f 0 j t j 1, 2, . . . , n with the mean value equal to zero can be simulated by the following series: H jm ω ml cos ω ml t − ϑ jm ω ml Φ ml j 1, 2, . . . , n , 2 where H ω is the Cholesky decomposition of CPSD matrix S 0 ω see the appendix and a lower triangular matrix; N is the sufficiently large dividing number of frequency; Δω ω up /N is the circular frequency increment; ω up refers to the upper cutoff circular frequency, with the condition that, when ω > ω up , the value of S 0 ω is trivial and negligible; Φ ml is the sequences of independent random phase angles, distributed uniformly over the interval 0, 2π ; the double-indexing frequency ω ml lΔω − N −m /N Δω l −1 Δω m/N Δω. It is noted that the simulated stochastic process using the SR method is asymptotically Gaussian as N → ∞ Deodatis 12 . Shinozuka and Deodatis 38 provided rigorous derivations and elaborations about asymptotic Gaussian of the simulated stochastic process according to the central limit theorem. It can be considered approximately Gaussian for most practical applications if N is greater than approximately 100 Deodatis and Micaletti 39 .
In order to avoid aliasing, the time interval Δt has to obey the condition as follows:

Mathematical Problems in Engineering
Likewise, the entire period of the sample functions in 2.1 can be calculated as follows: It has been proved by Deodatis 12 that the obtained results based on 2.1 possesses the ergodicity. Apparently, once the CPSD matrix is determined and the associated parameters, such as N, ω up , and Δω, are properly chosen, the stationary one-dimensional multivariate Gaussian stochastic process can then be simulated quite well with resorting to 2.1 .

Spline-Interpolation-Based FFT (SFFT) Algorithm
For the simulation of the multivariate stationary stochastic processes, the spline-interpolation-based FFT algorithm, referred to as the SFFT algorithm, is proposed in the present paper.

Spline Interpolation Technique
Since the Cholesky decomposition has to be conducted separately for each frequency ω ml , the number of performing the Cholesky decomposition with respect to 2.1 then is n × N. The computational effort of the SR method is tremendous in the large scale system, though the FFT algorithm can significantly enhance the computational speed.
For the multivariate stationary stochastic processes with the phase angles, the CPSD matrix S 0 ω is the complex-valued matrix. Then, the element H jm ω |H jm ω |e iϑ jm of the obtained lower triangular matrix H ω based on the Cholesky decomposition is the complex value. It is noted that |H jm ω | varies continuously with regard to the circular frequency. Therefore, as long as |H jm ω | at some appropriate circular frequency points are calculated, |H jm ω | at other circular frequency points can then be obtained by the cubic spline interpolation. Since the spline function has the nature of both the smoothness and continuance, it has widely been utilized to interpolate and fit data in engineering.
In the circular frequency interval ω i−1 , ω i , H jm ω is expressed as cubic polynomial, and then H jm ω is a linear function: Then, the continuous double integral of H jm ω can be expressed as ♦ Setting up (3.9) according to the continuation of divided into r subintervals by r − 1 frequencies .
in ω 1 , ω 2 , . . ., ω r−1 .  in which A i and B i are integral constants. According to the boundary conditions

6
Mathematical Problems in Engineering 3.5 Then, the derivative of H jm ω can be written as

3.6
In the interval ω i−1 , ω i , and, in next interval ω i , ω i 1 , According to 3.7 and 3.8 , we can obtain Mathematical Problems in Engineering 7 In order to determine P i i 0, 1, 2, . . . , r , the boundary conditions are assumed as H jm ω 0 0 and H jm ω r 0. Then, we can obtain the following functions: According to 3.9 , 3.10 and 3.11 , P i i 0, 1, 2, . . . , r can be determined through the following equation system

8 Mathematical Problems in Engineering
Introducing the spline interpolation into |H jm ω |, 2.1 for the simulation of the multivariate stationary stochastic processes can then be rewritten as follows:

Fast Fourier Transform (FFT) Algorithm
Taking into account that introducing the FFT algorithm renders very high computation efficiency in simulating the stationary processes Li and Kareem 27 , with resorting to the FFT algorithm, 3.14 can then be derived as follows:

3.16
It is noteworthy that g jm pΔt can be obtained in terms of the inverse FFT of B jml .
Mathematical Problems in Engineering

Numerical Results
In the following, two cases of performing the simulation of multivariate stationary fluctuating wind velocity are taken into consideration in order to demonstrate the capabilities and computational efficiency of the proposed SFFT algorithm.

Description of Case 1
A hundred velocity points to be simulated are evenly distributed from 10 m to 208 m heights along a vertical line see Case 1 in Figure 2 . The process corresponding to these hundred components is denoted by f 0 j t f 0 1 t , f 0 2 t , . . . , f 0 100 t . It is assumed that the mean value of the process is equal to zero. Then, the elements of its CPSD matrix may be given as in which S j ω represents the power spectral density function of f 0 j t ; Γ jk ω denotes the coherence function between f 0 j t and f 0 k t ; θ jk ω refers to the phase angle between S j ω and S k ω .
It is emphasized here that the stochastic process described by 4.1 and 4.2 is nonhomogeneous in space, since the longitudinal velocity fluctuations f 0 1 t , f 0 2 t , . . . , f 0 100 t have different frequency contents, namely, S 1 ω / S 2 ω / · · · / S 100 ω . where z is the elevation above the ground, in meters; ω is the circular frequency, in rad/s; u * is the shear velocity of the flow, in m/s; U z is the mean wind speed at height z, in m/s. Taking into account the following coherence function between the velocity fluctuations at two different heights z j and z k suggested by Davenport 41 , in which U z j and U z k are the mean wind speeds at heights z j and z k , respectively; Δz |z j −z k |; C z is a constant that may be set equal to 10 for structural design purposes Kristensen and Jensen 42 and Simiu and Scanlan 43 . Furthermore, the expression of the phase angle θ jk ω is taken from Di Paola 44 and Simiu and Scanlan 43 ,

4.5
In 4.2 , e −iθ jk ω is a measure of the wave passage delay due to the apparent velocity of waves υ j,k app . Likewise, the apparent velocity of waves can be assumed to be the following form suggested by Simiu and Scanlan 43 : in which C θ is an appropriate coefficient that has to be determined from experimental data. In the present paper, the expression by Peil and Telljohann 45 is taken into consideration, which has the form υ j,k app U z j z k /2 /5.5. The presence of the phase angle turns into a time shift of the peak of the cross-correlation functions from the data of the measured fluctuations at different heights.
Let it be supposed that the mean wind velocity at the first point 1 z 10 m is U 10 26 m/s and that the surface rough roughness length is  to T 8192 s. The mean wind speeds at other points are computed in terms of the logarithmic law i.e., the vertical profile , which has the form as follows: U z j ln z j /z 0 ln z s /z 0 U z s . 4.7

Description of Case 2
Fifty velocity points to be simulated are evenly distributed from 10 m to 255 m heights along a vertical line see Case 2 in Figure 2 . The process corresponding to these fifty components is

Numerical Analyses
The comparisons among several different schemes are made and listed in Table 1. The Cholesky decomposition with the SFFT algorithm is implemented on 128 circular frequency points. It can be seen that in Case 1, the elapsed time for the simulation with the FFT algorithm is 215.9 min, while that for the simulation with the SFFT algorithm is only 10.2 min; in the Case 2, the elapsed time for the simulation with the FFT algorithm is 5.9 min, while that for the simulation with the SFFT algorithm is only 0.5 min. Apparently, in the calculation speed the SFFT algorithm has considerable advantage over the FFT algorithm. Likewise, the superiority becomes more remarkable with the increasing of N, n, and T . This is because the number of the Cholesky decomposition with the FFT algorithm equals n × N, while that with the SFFT algorithm keeps constant with the value equal to 128. For a determined SFFT or FFT algorithm, however, the computational efficiency becomes lower with the increasing of the parameters N, n, and T .
Since in Case 1 the number of the stimulated time histories is very large, without any loss of generality, arbitrary six time histories of all are chosen to demonstrate the effort of simulating the longitudinal wind velocity fluctuations by resorting to the proposed SFFT algorithm. The generated samples of longitudinal wind velocity fluctuations at six points 1, 11, 26, 46, 71, and 100 see Case 1 in Figure 2 , denoted by f 1 t , f 11 t , f 26 t , f 46 t , f 71 t , and f 100 t , respectively, are displayed in Figures 3 a , 3 b , 3 c , 3  extracted out, respectively, as shown Figure 4, in order to better visualize the differences and the similarities among these six time histories.
As far as the correlation degree among these six time histories is concerned, it is clearly seen from Figure 4 that with the decreasing of the distance between arbitrary two points, the loss of coherence between two samples corresponding, respectively, to these two points will  Figure 2 . This behavior is controlled by the coherence functions see 4.4 . Furthermore, the phase angle can also be detected clearly at 100 s time instant between arbitrary two samples of these six time histories. Likewise, the phase angle of these six time histories will become larger with the increasing of the distance between arbitrary two time histories. Similarly, this phenomenon can be controlled by the phase angle functions see 4.5 .
The forgoing elucidations indicate that the proposed SFFT algorithm is able to simulate the longitudinal wind velocity fluctuations that are spatially correlated according to a prescribed coherence function and that possess the phase angles following a prescribed phase angle function.
Plotted in Figure 5 is the temporal autocorrelation function R jj τ j 1, 11, 26, 46, 71, 100 of the generated sample function shown in Figure 3, along with the target autocorrelation function R 0 jj τ j 1, 11, 26, 46, 71, 100 . It is worth pointing out that the target autocorrelation function R 0 jj τ is computed with resorting to the Wiener-Khintchine transformation. As can be seen in Figure 5, the simulated autocorrelation function R jj τ practically coincides with the target autocorrelation function R 0 jj τ . Displayed in Figure 6 is the temporal cross-correlation function R jk τ j 1, k 11, 26, 46, 71, 100 of the generated sample function presented in Figure 3, along with the target cross-correlation function R 0 jk τ j 1, k 11, 26, 46, 71, 100 . The same, the target cross-correlation function R 0 jk τ is calculated according to the Wiener-Khintchine transformation. As expected, the simulated cross-correlation function R jk τ j / k practically coincides    with the target cross-correlation function R 0 jk τ j / k . Apparently, the temporal auto correlation and cross-correlation functions of any sample function f j t j 1, 2, . . . , 100 are, respectively, identical to the target autocorrelation and cross-correlation functions based on the SFFT algorithm. Moreover, Figure 6 demonstrates that the presence of the phase angle turns into a time shift of the peak of the cross-correlation functions. Likewise, the time shift will become larger with the increasing of the distance between arbitrary two points. It is proved once again that the proposed SFFT algorithm can efficiently carry out the simulation of longitudinal wind velocity fluctuations with the phase angles.
From Figure 7, it is clear that the estimated spectra derived from these sample time histories and the target spectra are in remarkably good agreement with each other. Likewise, it is seen from Figures 8 and 9 that the mean and variance values of these 100 time histories generated with the SFFT algorithm approach those with the FFT algorithm, which implies that the SFFT and FFT algorithms practically achieve the same computation accuracy. More specifically, Figure 8 shows that the variance values of the generated time histories with resorting to both the SFFT and FFT are identical to the target values, respectively. It is observed in Figure 9 that the mean values of the generated time histories with resorting to both the SFFT and FFT are not more than 0.2 m/s, effectively meaning that these mean values are considerably small with respect to the zero target value.
In order to demonstrate Gaussian of the simulated wind velocity fluctuations, Figure 10 shows the probability density functions of generated samples displayed in Figure 3 versus corresponding target Gaussian . It is clear that the probability density functions estimated from these samples match closely Gaussian probability density function.

Accuracy and Time Expense of Various Interpolations
Four kinds of interpolation techniques, such as SFFT, the cubic Lagrangian interpolation based FFT CL-FFT , square Lagrangian interpolation based FFT SL-FFT , and neural network interpolation based FFT NN-FFT , are taken into consideration to simulate the wind velocity field for the evaluation of accuracy and time expense of various interpolations 46 . The root mean square error RMSE and error factor EF 47 are introduced to evaluate the accuracy of these four interpolation-based FFT techniques: in which, N is the total number of time points in the generated sample; Y i and y i signify generated wind velocities at a certain instant employing the FFT and the interpolation-based FFT, respectively. It is noted that, the interpolation based FFT is more precise if the value of RMSE is closer to zero; the interpolation-based FFT is more precise if the value of EF is closer to one. The data coming from 46 are listed in Tables 2 and 3. Table 2 shows the elapsed time in the simulation using the SL-FFT, CL-FFT, SFFT, NN-FFT, and FFT, respectively, and Table 3 records both the RMSE and EF values. It is from Table 2 seen that the efficiency of the SFFT  is the highest due to the minimal time expense. It is from Table 3 observed that the NN-FFT presents the highest accuracy. Worth noting, the accuracy by the SFFT is a little lower than that by the NN-FFT. Forasmuch, the SFFT, generally, is the best.

Conclusions
This study presented the SFFT algorithm to further enhance the computational efficiency of simulating the multivariate stochastic processes. More specifically, the main purpose of the proposed SFFT algorithm, simultaneously with resorting to the spline interpolation and FFT algorithms, is to reduce the number of the Cholesky decomposition. The main conclusions of the present research effort on the SFFT algorithm can be drawn as follows.
1 The proposed SFFT algorithm is able to efficiently generate the longitudinal wind velocity fluctuations that are spatially correlated according to a prescribed coherence function and that possess the phase angles following a prescribed phase angle function.   2 The proposed SFFT algorithm has considerable advantage over the FFT algorithm. Likewise, this superiority will become more and more remarkable with the increasing of the parameters N, n, and T .
3 Although introducing the spline interpolation approximation in decomposing the cross-power spectral density matrix, the SFFT and FFT algorithms practically achieve the same computational accuracy.
4 Employing the SFFT algorithm, the temporal autocorrelation and cross-correlation functions of any sample function are identical to the target autocorrelation and cross-correlation functions, respectively. Likewise, the estimated spectra derived from these sample time histories and the target spectra are in remarkably good agreement with each other. 5 The SFFT, generally, is the best in many interpolation methods.
Eventually, it is worth emphasizing that since the spline interpolation is able to significantly reduce the number of the Cholesky decomposition of time-varying power spectral matrix, the SFFT algorithm can improve the performance of the SR method for the simulation of the multivariate nonstationary stochastic processes. It is known that the power spectrum of nonstationary stochastic process is time-varying; thus, meaning that the FFT algorithm is very hard to be used in this case. Only when the time-varying power spectrum satisfies some special conditions, the FFT can be utilized. However, it is anticipated that the SFFT algorithm can be used to enhance the computational efficiency of simulating the multivariate nonstationary stochastic processes, since the spline interpolation technique is not affected by the time-varying power spectrum.