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

Authors: Jinhua Li; Chunxiang Li; Shuisheng Chen
Journal: Mathematical Problems in Engineering
 (2011)

---

## Abstract

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.


---

## Body

## 1. 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) 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. [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.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 spline-interpolation-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.


## 2. Simulation of Multivariate Stochastic Processes

According to Deodatis [12] and Shinozuka [19], the multivariate stationary stochastic process fj0(t)(j=1,2,…,n) with the mean value equal to zero can be simulated by the following series:(2.1)fj(t)=2Δω∑m=1n∑l=1N|Hjm(ωml)|cos⁡[ωmlt-ϑjm(ωml)+Φml](j=1,2,…,n),
where H(ω) is the Cholesky decomposition of CPSD matrix S0(ω) (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 S0(ω) 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 asN→∞ (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:(2.2)Δt≤2π2ωup.
Likewise, the entire period of the sample functions in (2.1) can be calculated as follows:(2.3)T0=2πnΔω=2πnNωup.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).


## 3. 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.


### 3.1. 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 matrixS0(ω) is the complex-valued matrix. Then, the element Hjm(ω)=|Hjm(ω)|eiϑjm of the obtained lower triangular matrix H(ω) based on the Cholesky decomposition is the complex value. It is noted that |Hjm(ω)| varies continuously with regard to the circular frequency. Therefore, as long as |Hjm(ω)| at some appropriate circular frequency points are calculated, |Hjm(ω)| 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.The circular frequency interval[0,ωu] is evenly divided into r subintervals by r-1 frequencies  ω1,ω2,…,ωr-1(0=ω0<ω1<ω2<⋯<ωr-1<ωr=ωu). The corresponding |Hjm(ωi)|(i=0,1,2,…,r) can then be calculated. It is assumed that H̃jm(ωi)=|Hjm(ωi)|(i=0,1,2,…,r) and second continuous derivative H̃jm′′(ωi)=Pi(i=0,1,2,…,r), then H̃jm(ω) at other frequency points can be obtained by the cubic spline interpolation. The flow scheme of building spline interpolation H̃jm(ω) is shown in Figure 1.Figure 1
Flow scheme of building spline interpolation functionH̃jm(ω).In the circular frequency interval[ωi-1,ωi],H̃jm(ω) is expressed as cubic polynomial, and then H̃jm′′(ω) is a linear function:(3.1)H̃jm′′(ω)=Pi-1ωi-ωωi-ωi-1+Piω-ωi-1ωi-ωi-1=rωu(Pi-1(ωi-ω)+Pi(ω-ωi-1)).
Then, the continuous double integral of H̃jm′′(ω) can be expressed as  (3.2)H̃jm(ω)=r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+Ai(ωi-ω)+Bi(ω-ωi-1),
in which Ai and Bi are integral constants. According to the boundary conditions(3.3)H̃jm(ωi-1)=ωu26r2Pi-1+Aiωur,H̃jm(ωi)=ωu26r2Pi+Biωur,
the integral constants Ai and Bi can be expressed as(3.4)Ai=(H̃jm(ωi-1)-ωu26r2Pi-1)rωu,Bi=(H̃jm(ωi)-ωu26r2Pi)rωu.
Substituting (3.4) into (3.2), (3.5)H̃jm(ω)=r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+(H̃jm(ωi-1)-ωu26r2Pi-1)rωu(ωi-ω)+(H̃jm(ωi)-ωu26r2Pi)rωu(ω-ωi-1)(ωi-1≤ω≤ωi;i=1,2,3,…,r).
Then, the derivative of H̃jm(ω) can be written  as(3.6)H̃jm′(ω)=r2ωu(ω-ωi-1)2Pi-r2ωu(ωi-ω)2Pi-1+ωu6r(Pi-1-Pi)+rωu(H̃jm(ωi)-H̃jm(ωi-1))(ωi-1≤ω≤ωi;i=1,2,3,…,r).
In the interval [ωi-1,ωi],(3.7)H̃jm′(ωi)=ωu3rPi+ωu6rPi-1+rωu(H̃jm(ωi)-H̃jm(ωi-1)),
and, in next interval [ωi,ωi+1],(3.8)H̃jm′(ωi)=-ωu3rPi-ωu6rPi+1+rωu(H̃jm(ωi+1)-H̃jm(ωi)).
According to (3.7) and (3.8), we can obtain  (3.9)Pi+1+4Pi+Pi-1=6r2ωu2(H̃jm(ωi+1)-2H̃jm(ωi)+H̃jm(ωi-1))(i=1,2,3,…,r-1).
In order to determine Pi(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:(3.10)2P0+P1=6r2ωu2(H̃jm(ω1)-H̃jm(ω0)),(3.11)2Pr+Pr-1=6r2ωu2(H̃jm(ωr-1)-H̃jm(ωr)).
According to (3.9), (3.10) and (3.11),  Pi(i=0,1,2,…,r) can be determined through the following equation system(3.12)2P0+P1=6r2ωu2(H̃jm(ω1)-H̃jm(ω0)),P2+4P1+P0=6r2ωu2(H̃jm(ω2)-2H̃jm(ω1)+H̃jm(ω0)),P3+4P2+P1=6r2ωu2(H̃jm(ω3)-2H̃jm(ω2)+H̃jm(ω1)),⋮Pr+4Pr-1+Pr-2=6r2ωu2(H̃jm(ωr)-2H̃jm(ωr-1)+H̃jm(ωr-2)),2Pr+Pr-1=6r2ωu2(H̃jm(ωr-1)-H̃jm(ωr)).
Then, substituting Pi(i=0,1,2,…,r)  into (3.5), we can obtain (3.13)H̃jm(ω)={r6ωu[(ω1-ω)3P0+(ω-ω0)3P1]+(H̃jm(ω0)-ωu26r2P0)rωu(ω1-ω)+(H̃jm(ω1)-ωu26r2P1)rωu(ω-ω0)(ω0≤ω≤ω1),r6ωu[(ω2-ω)3P1+(ω-ω1)3P2]+(H̃jm(ω1)-ωu26r2P1)rωu(ω2-ω)+(H̃jm(ω2)-ωu26r2P2)rωu(ω-ω1)(ω1≤ω≤ω2),⋮r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+(H̃jm(ωi-1)-ωu26r2Pi-1)rωu(ωi-ω)+(H̃jm(ωi)-ωu26r2Pi)rωu(ω-ωi-1)(ωr-1≤ω≤ωr).Introducing the spline interpolation into|Hjm(ω)|, (2.1) for the simulation of the multivariate stationary stochastic processes can then be rewritten as follows:(3.14)fj(t)=2Δω∑m=1n∑l=1NH̃jm(ωml)cos⁡[ωmlt-ϑjm(ωml)+Φml](j=1,2,…,n).


### 3.2. 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.15)fj(pΔt)=Re{∑m=1jh̃jm(pΔt)exp⁡[i(mΔωn)(pΔt)]},
where(3.16)p=0,1,…,n×M-1;M=T0nΔt;T0=2πnΔω=2πnNωup;h̃jm(pΔt)={g̃jm(pΔt),p=0,1,…,M-1,g̃jm[(p-M)Δt],p=M,M+1,…,2M-1,⋮⋮g̃jm[(p-(n-1)M)Δt],p=(n-1)M,(n-1)M+1,…,nM-1,g̃jm(pΔt)=∑l=1M-1B̃jmlexp⁡[ilp2πM](p=0,1,…,M-1),B̃jml={2ΔωH̃jm(lΔω+mΔωn)exp⁡[-iϑjm(lΔω+mΔωn)]exp⁡(iΦml),0≤l≤N,0,N<l≤M-1.It is noteworthy thatg̃jm(pΔt) can be obtained in terms of the inverse FFT of B̃jml.


## 3.1. 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 matrixS0(ω) is the complex-valued matrix. Then, the element Hjm(ω)=|Hjm(ω)|eiϑjm of the obtained lower triangular matrix H(ω) based on the Cholesky decomposition is the complex value. It is noted that |Hjm(ω)| varies continuously with regard to the circular frequency. Therefore, as long as |Hjm(ω)| at some appropriate circular frequency points are calculated, |Hjm(ω)| 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.The circular frequency interval[0,ωu] is evenly divided into r subintervals by r-1 frequencies  ω1,ω2,…,ωr-1(0=ω0<ω1<ω2<⋯<ωr-1<ωr=ωu). The corresponding |Hjm(ωi)|(i=0,1,2,…,r) can then be calculated. It is assumed that H̃jm(ωi)=|Hjm(ωi)|(i=0,1,2,…,r) and second continuous derivative H̃jm′′(ωi)=Pi(i=0,1,2,…,r), then H̃jm(ω) at other frequency points can be obtained by the cubic spline interpolation. The flow scheme of building spline interpolation H̃jm(ω) is shown in Figure 1.Figure 1
Flow scheme of building spline interpolation functionH̃jm(ω).In the circular frequency interval[ωi-1,ωi],H̃jm(ω) is expressed as cubic polynomial, and then H̃jm′′(ω) is a linear function:(3.1)H̃jm′′(ω)=Pi-1ωi-ωωi-ωi-1+Piω-ωi-1ωi-ωi-1=rωu(Pi-1(ωi-ω)+Pi(ω-ωi-1)).
Then, the continuous double integral of H̃jm′′(ω) can be expressed as  (3.2)H̃jm(ω)=r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+Ai(ωi-ω)+Bi(ω-ωi-1),
in which Ai and Bi are integral constants. According to the boundary conditions(3.3)H̃jm(ωi-1)=ωu26r2Pi-1+Aiωur,H̃jm(ωi)=ωu26r2Pi+Biωur,
the integral constants Ai and Bi can be expressed as(3.4)Ai=(H̃jm(ωi-1)-ωu26r2Pi-1)rωu,Bi=(H̃jm(ωi)-ωu26r2Pi)rωu.
Substituting (3.4) into (3.2), (3.5)H̃jm(ω)=r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+(H̃jm(ωi-1)-ωu26r2Pi-1)rωu(ωi-ω)+(H̃jm(ωi)-ωu26r2Pi)rωu(ω-ωi-1)(ωi-1≤ω≤ωi;i=1,2,3,…,r).
Then, the derivative of H̃jm(ω) can be written  as(3.6)H̃jm′(ω)=r2ωu(ω-ωi-1)2Pi-r2ωu(ωi-ω)2Pi-1+ωu6r(Pi-1-Pi)+rωu(H̃jm(ωi)-H̃jm(ωi-1))(ωi-1≤ω≤ωi;i=1,2,3,…,r).
In the interval [ωi-1,ωi],(3.7)H̃jm′(ωi)=ωu3rPi+ωu6rPi-1+rωu(H̃jm(ωi)-H̃jm(ωi-1)),
and, in next interval [ωi,ωi+1],(3.8)H̃jm′(ωi)=-ωu3rPi-ωu6rPi+1+rωu(H̃jm(ωi+1)-H̃jm(ωi)).
According to (3.7) and (3.8), we can obtain  (3.9)Pi+1+4Pi+Pi-1=6r2ωu2(H̃jm(ωi+1)-2H̃jm(ωi)+H̃jm(ωi-1))(i=1,2,3,…,r-1).
In order to determine Pi(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:(3.10)2P0+P1=6r2ωu2(H̃jm(ω1)-H̃jm(ω0)),(3.11)2Pr+Pr-1=6r2ωu2(H̃jm(ωr-1)-H̃jm(ωr)).
According to (3.9), (3.10) and (3.11),  Pi(i=0,1,2,…,r) can be determined through the following equation system(3.12)2P0+P1=6r2ωu2(H̃jm(ω1)-H̃jm(ω0)),P2+4P1+P0=6r2ωu2(H̃jm(ω2)-2H̃jm(ω1)+H̃jm(ω0)),P3+4P2+P1=6r2ωu2(H̃jm(ω3)-2H̃jm(ω2)+H̃jm(ω1)),⋮Pr+4Pr-1+Pr-2=6r2ωu2(H̃jm(ωr)-2H̃jm(ωr-1)+H̃jm(ωr-2)),2Pr+Pr-1=6r2ωu2(H̃jm(ωr-1)-H̃jm(ωr)).
Then, substituting Pi(i=0,1,2,…,r)  into (3.5), we can obtain (3.13)H̃jm(ω)={r6ωu[(ω1-ω)3P0+(ω-ω0)3P1]+(H̃jm(ω0)-ωu26r2P0)rωu(ω1-ω)+(H̃jm(ω1)-ωu26r2P1)rωu(ω-ω0)(ω0≤ω≤ω1),r6ωu[(ω2-ω)3P1+(ω-ω1)3P2]+(H̃jm(ω1)-ωu26r2P1)rωu(ω2-ω)+(H̃jm(ω2)-ωu26r2P2)rωu(ω-ω1)(ω1≤ω≤ω2),⋮r6ωu[(ωi-ω)3Pi-1+(ω-ωi-1)3Pi]+(H̃jm(ωi-1)-ωu26r2Pi-1)rωu(ωi-ω)+(H̃jm(ωi)-ωu26r2Pi)rωu(ω-ωi-1)(ωr-1≤ω≤ωr).Introducing the spline interpolation into|Hjm(ω)|, (2.1) for the simulation of the multivariate stationary stochastic processes can then be rewritten as follows:(3.14)fj(t)=2Δω∑m=1n∑l=1NH̃jm(ωml)cos⁡[ωmlt-ϑjm(ωml)+Φml](j=1,2,…,n).


## 3.2. 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.15)fj(pΔt)=Re{∑m=1jh̃jm(pΔt)exp⁡[i(mΔωn)(pΔt)]},
where(3.16)p=0,1,…,n×M-1;M=T0nΔt;T0=2πnΔω=2πnNωup;h̃jm(pΔt)={g̃jm(pΔt),p=0,1,…,M-1,g̃jm[(p-M)Δt],p=M,M+1,…,2M-1,⋮⋮g̃jm[(p-(n-1)M)Δt],p=(n-1)M,(n-1)M+1,…,nM-1,g̃jm(pΔt)=∑l=1M-1B̃jmlexp⁡[ilp2πM](p=0,1,…,M-1),B̃jml={2ΔωH̃jm(lΔω+mΔωn)exp⁡[-iϑjm(lΔω+mΔωn)]exp⁡(iΦml),0≤l≤N,0,N<l≤M-1.It is noteworthy thatg̃jm(pΔt) can be obtained in terms of the inverse FFT of B̃jml.


## 4. 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.


### 4.1. 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 Figure2). The process corresponding to these hundred components is denoted by  fj0(t)[f10(t),f20(t),…,f1000(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 (4.1)Sjj0(ω)=Sj(ω)(j=1,2,…,100),(4.2)Sjk0(ω)=Γjk(ω)Sj(ω)Sk(ω)e-iθjk(ω)(j,k=1,2,…,100,j≠k),
in which Sj(ω) represents the power spectral density function of fj0(t);Γjk(ω) denotes the coherence function between fj0(t) and fk0(t); θjk(ω) refers to the phase angle between Sj(ω) and Sk(ω).Figure 2
Configuration of points along a vertical line in Case 1 and Case 2.It is emphasized here that the stochastic process described by (4.1) and (4.2) is nonhomogeneous in space, since the longitudinal velocity fluctuations f10(t),f20(t),…,f1000(t) have different frequency contents, namely, S1(ω)≠S2(ω)≠⋯≠S100(ω).The following expression proposed by Kaimal et al. [40] is selected to model the two-sided power spectral density (PSD) function of the longitudinal wind velocity fluctuations at different height:(4.3)S(z,ω)=122002πu*2zU(z)1[1+50(ωz/2πU(z))]5/3,
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 heightszj and zk suggested by Davenport [41],(4.4)Γ(Δz,ω)=exp⁡[-ω2πCzΔz(1/2)[U(zj)+U(zk)]]
in which U(zj) and U(zk) are the mean wind speeds at heights zj and zk, respectively; Δz=|zj-zk|; Cz 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)θjk(ω)=ω(zj-zk)υapp(j,k).
In (4.2),  e-iθjk(ω) is a measure of the wave passage delay due to the apparent velocity of waves υapp(j,k).Likewise, the apparent velocity of waves can be assumed to be the following form suggested by Simiu and Scanlan [43]:(4.6)υapp(j,k)=π[U(zj)+U(zk)]Cθ,
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 υapp(j,k)=U[(zj+zk)/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=10m) is U(10)=26m/s and that the surface rough roughness length is z0=0.001266 m (corresponding shear velocity of the flow u*=1.76 m/s). It is worth mentioning that the values z0=0.001266mand u*=1.76m/s are taken from Simiu and Scanlan [43]. The upper cutoff frequency is ωup=2π rad/s. The dividing number of frequency is N=8192. The simulation is carried out at 1×M=1×8192×2=16384 time instants with a time step Δt=0.5s, over a time length equal to T=8192s. 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:(4.7)U(zj)=ln⁡(zj/z0)ln⁡(zs/z0)U(zs).


### 4.2. 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 Figure2). The process corresponding to these fifty components is represented by fj0(t)[f10(t),f20(t),…,f500(t)]. The dividing number of frequency is N=2048. The numerical simulation is implemented at 1×M=1×2048×2=4096 time instants with a time step Δt=0.5s, over a time length equal to T=2048s. The other conditions are the same as Case  1.


### 4.3. Numerical Analyses

The comparisons among several different schemes are made and listed in Table1. 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.Table 1
Comparisons among several different schemes for wind velocity simulation.
CasesAlgorithmsDividing number of frequency(N)Number of points(n)T (sec)Elapsed time (min)Case 1SFFT8192100819210.2FFT81921008192215.9Case 2SFFT20485020480.5FFT20485020485.9Since 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 Figure2), denoted by f1(t),f11(t),f26(t),f46(t),f71(t), and f100(t), respectively, are displayed in Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f) for a time length T=8192s. Likewise, the parts in the first 600 s of these six time histories are further extracted out, respectively, as shown Figure 4, in order to better visualize the differences and the similarities among these six time histories.Generated sample functions for longitudinal wind velocity fluctuations at six different heights in Case 1.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mGenerated sample functions for longitudinal wind velocity fluctuations at six different heights, over first 600 s of generated sample functions displayed in Figure3.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mAs far as the correlation degree among these six time histories is concerned, it is clearly seen from Figure4 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 become smaller. For example, since the distance between the points 1 and 11 is the smallest in these six points, consequently, there is the smallest loss of coherence between f1(t) and f11(t). On the other hand, there is the biggest loss of coherence between f1(t) and f100(t), as the point 100 is located at 198 m height from the point 1 (see Case 1 in 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 Figure5 is the temporal autocorrelation function [Rjj(τ)(j=1,11,26,46,71,100)] of the generated sample function shown in Figure 3, along with the target autocorrelation function [Rjj0(τ)(j=1,11,26,46,71,100)]. It is worth pointing out that the target autocorrelation function [Rjj0(τ)] is computed with resorting to the Wiener-Khintchine transformation. As can be seen in Figure 5, the simulated autocorrelation function [Rjj(τ)] practically coincides with the target autocorrelation function [Rjj0(τ)]. Displayed in Figure 6 is the temporal cross-correlation function [Rjk(τ)(j=1,k=11,26,46,71,100)] of the generated sample function presented in Figure 3, along with the target cross-correlation function [Rjk0(τ)(j=1,k=11,26,46,71,100)]. The same, the target cross-correlation function [Rjk0(τ)] is calculated according to the Wiener-Khintchine transformation. As expected, the simulated cross-correlation function [Rjk(τ)(j≠k)] practically coincides with the target cross-correlation function [Rjk0(τ)(j≠k)]. Apparently, the temporal auto correlation and cross-correlation functions of any sample function [fj(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.Temporal autocorrelation functions of generated sample function displayed in Figure3 versus corresponding targets.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mTemporal cross-correlation functions of generated sample functions displayed in Figure3 versus corresponding targets.
(a)
z=10m and z=30m(b)
z=10m and z=60m(c)
z=10m and z=100m(d)
z=10m and z=150m(e)
z=10m and z=208mFrom Figure7, 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.Power spectrum of generated sample functions displayed in Figure3 versus corresponding targets.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mFigure 8
Variance values of sample generated by resorting to the SFFT or the FFT with respect to corresponding targets in Case 1.Figure 9
Mean values of sample generated with resorting to the SFFT or the FFT with reference to corresponding targets in Case 1.In order to demonstrate Gaussian of the simulated wind velocity fluctuations, Figure10 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.Probability density functions of generated sample functions displayed in Figure3 versus corresponding target (Gaussian).
(a)
z=10m(b)
z=30m(c)
z=60m(d)
z=100m(e)
z=150m(f)
z=208m


### 4.4. 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:(4.8)RMSE=1N∑i=1N(yi-Yi)2,EF=1-(∑i=1N|yi-Yi|∑i=1N|Yi|),
in which, N is the total number of time points in the generated sample; Yi and yi 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.Table 2
Comparison of time expense among different simulation schemes (noted that the data come from [46]).
Different schemesTime expense (sec)10 location points20 location pointsSFFT8.3529.84CL-FFT12.7647.80SL-FFT10.8040.36NN-FFT20.2474.98FFT43.03301.55Table 3
Comparisons of MRSE and EF among different simulation schemes (noted that the data come from [46]).
Locatio pointMRSEEFSFFTCL-FFTSL-FFTNN-FFTSFFTCL-FFTSL-FFTNN-FFT10.16930.18960.23420.16550.98120.98040.97470.981820.16880.17910.23300.16480.98100.98020.97450.981730.16810.17820.23130.16430.98140.98060.97510.982140.16750.17740.22950.16360.98130.98050.97510.981950.16670.17640.22780.16330.98080.98010.97460.981560.16620.17580.22630.16250.98090.98020.97470.981670.16540.17480.22430.16210.98110.98040.97510.981780.16450.17370.22220.16170.98100.98030.97520.981690.16380.17270.22040.16070.98130.98060.97560.9819100.16270.17130.21770.15970.98200.98130.97650.9825


## 4.1. 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 Figure2). The process corresponding to these hundred components is denoted by  fj0(t)[f10(t),f20(t),…,f1000(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 (4.1)Sjj0(ω)=Sj(ω)(j=1,2,…,100),(4.2)Sjk0(ω)=Γjk(ω)Sj(ω)Sk(ω)e-iθjk(ω)(j,k=1,2,…,100,j≠k),
in which Sj(ω) represents the power spectral density function of fj0(t);Γjk(ω) denotes the coherence function between fj0(t) and fk0(t); θjk(ω) refers to the phase angle between Sj(ω) and Sk(ω).Figure 2
Configuration of points along a vertical line in Case 1 and Case 2.It is emphasized here that the stochastic process described by (4.1) and (4.2) is nonhomogeneous in space, since the longitudinal velocity fluctuations f10(t),f20(t),…,f1000(t) have different frequency contents, namely, S1(ω)≠S2(ω)≠⋯≠S100(ω).The following expression proposed by Kaimal et al. [40] is selected to model the two-sided power spectral density (PSD) function of the longitudinal wind velocity fluctuations at different height:(4.3)S(z,ω)=122002πu*2zU(z)1[1+50(ωz/2πU(z))]5/3,
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 heightszj and zk suggested by Davenport [41],(4.4)Γ(Δz,ω)=exp⁡[-ω2πCzΔz(1/2)[U(zj)+U(zk)]]
in which U(zj) and U(zk) are the mean wind speeds at heights zj and zk, respectively; Δz=|zj-zk|; Cz 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)θjk(ω)=ω(zj-zk)υapp(j,k).
In (4.2),  e-iθjk(ω) is a measure of the wave passage delay due to the apparent velocity of waves υapp(j,k).Likewise, the apparent velocity of waves can be assumed to be the following form suggested by Simiu and Scanlan [43]:(4.6)υapp(j,k)=π[U(zj)+U(zk)]Cθ,
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 υapp(j,k)=U[(zj+zk)/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=10m) is U(10)=26m/s and that the surface rough roughness length is z0=0.001266 m (corresponding shear velocity of the flow u*=1.76 m/s). It is worth mentioning that the values z0=0.001266mand u*=1.76m/s are taken from Simiu and Scanlan [43]. The upper cutoff frequency is ωup=2π rad/s. The dividing number of frequency is N=8192. The simulation is carried out at 1×M=1×8192×2=16384 time instants with a time step Δt=0.5s, over a time length equal to T=8192s. 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:(4.7)U(zj)=ln⁡(zj/z0)ln⁡(zs/z0)U(zs).


## 4.2. 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 Figure2). The process corresponding to these fifty components is represented by fj0(t)[f10(t),f20(t),…,f500(t)]. The dividing number of frequency is N=2048. The numerical simulation is implemented at 1×M=1×2048×2=4096 time instants with a time step Δt=0.5s, over a time length equal to T=2048s. The other conditions are the same as Case  1.


## 4.3. Numerical Analyses

The comparisons among several different schemes are made and listed in Table1. 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.Table 1
Comparisons among several different schemes for wind velocity simulation.
CasesAlgorithmsDividing number of frequency(N)Number of points(n)T (sec)Elapsed time (min)Case 1SFFT8192100819210.2FFT81921008192215.9Case 2SFFT20485020480.5FFT20485020485.9Since 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 Figure2), denoted by f1(t),f11(t),f26(t),f46(t),f71(t), and f100(t), respectively, are displayed in Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f) for a time length T=8192s. Likewise, the parts in the first 600 s of these six time histories are further extracted out, respectively, as shown Figure 4, in order to better visualize the differences and the similarities among these six time histories.Generated sample functions for longitudinal wind velocity fluctuations at six different heights in Case 1.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mGenerated sample functions for longitudinal wind velocity fluctuations at six different heights, over first 600 s of generated sample functions displayed in Figure3.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mAs far as the correlation degree among these six time histories is concerned, it is clearly seen from Figure4 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 become smaller. For example, since the distance between the points 1 and 11 is the smallest in these six points, consequently, there is the smallest loss of coherence between f1(t) and f11(t). On the other hand, there is the biggest loss of coherence between f1(t) and f100(t), as the point 100 is located at 198 m height from the point 1 (see Case 1 in 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 Figure5 is the temporal autocorrelation function [Rjj(τ)(j=1,11,26,46,71,100)] of the generated sample function shown in Figure 3, along with the target autocorrelation function [Rjj0(τ)(j=1,11,26,46,71,100)]. It is worth pointing out that the target autocorrelation function [Rjj0(τ)] is computed with resorting to the Wiener-Khintchine transformation. As can be seen in Figure 5, the simulated autocorrelation function [Rjj(τ)] practically coincides with the target autocorrelation function [Rjj0(τ)]. Displayed in Figure 6 is the temporal cross-correlation function [Rjk(τ)(j=1,k=11,26,46,71,100)] of the generated sample function presented in Figure 3, along with the target cross-correlation function [Rjk0(τ)(j=1,k=11,26,46,71,100)]. The same, the target cross-correlation function [Rjk0(τ)] is calculated according to the Wiener-Khintchine transformation. As expected, the simulated cross-correlation function [Rjk(τ)(j≠k)] practically coincides with the target cross-correlation function [Rjk0(τ)(j≠k)]. Apparently, the temporal auto correlation and cross-correlation functions of any sample function [fj(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.Temporal autocorrelation functions of generated sample function displayed in Figure3 versus corresponding targets.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mTemporal cross-correlation functions of generated sample functions displayed in Figure3 versus corresponding targets.
(a)
z=10m and z=30m(b)
z=10m and z=60m(c)
z=10m and z=100m(d)
z=10m and z=150m(e)
z=10m and z=208mFrom Figure7, 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.Power spectrum of generated sample functions displayed in Figure3 versus corresponding targets.
(a)
Heightz=10m(b)
Heightz=30m(c)
Heightz=60m(d)
Heightz=100m(e)
Heightz=150m(f)
Heightz=208mFigure 8
Variance values of sample generated by resorting to the SFFT or the FFT with respect to corresponding targets in Case 1.Figure 9
Mean values of sample generated with resorting to the SFFT or the FFT with reference to corresponding targets in Case 1.In order to demonstrate Gaussian of the simulated wind velocity fluctuations, Figure10 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.Probability density functions of generated sample functions displayed in Figure3 versus corresponding target (Gaussian).
(a)
z=10m(b)
z=30m(c)
z=60m(d)
z=100m(e)
z=150m(f)
z=208m


## 4.4. 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:(4.8)RMSE=1N∑i=1N(yi-Yi)2,EF=1-(∑i=1N|yi-Yi|∑i=1N|Yi|),
in which, N is the total number of time points in the generated sample; Yi and yi 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.Table 2
Comparison of time expense among different simulation schemes (noted that the data come from [46]).
Different schemesTime expense (sec)10 location points20 location pointsSFFT8.3529.84CL-FFT12.7647.80SL-FFT10.8040.36NN-FFT20.2474.98FFT43.03301.55Table 3
Comparisons of MRSE and EF among different simulation schemes (noted that the data come from [46]).
Locatio pointMRSEEFSFFTCL-FFTSL-FFTNN-FFTSFFTCL-FFTSL-FFTNN-FFT10.16930.18960.23420.16550.98120.98040.97470.981820.16880.17910.23300.16480.98100.98020.97450.981730.16810.17820.23130.16430.98140.98060.97510.982140.16750.17740.22950.16360.98130.98050.97510.981950.16670.17640.22780.16330.98080.98010.97460.981560.16620.17580.22630.16250.98090.98020.97470.981670.16540.17480.22430.16210.98110.98040.97510.981780.16450.17370.22220.16170.98100.98030.97520.981690.16380.17270.22040.16070.98130.98060.97560.9819100.16270.17130.21770.15970.98200.98130.97650.9825


## 5. 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 parametersN,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.



---

*Source: 842183-2011-11-01.xml*