## 1 Introduction  
The seismic response of large-scale and long-span structures is seriously affected by the spatially varying characteristics of ground motions [1; 2]. Therefore, the seismic response analysis of such structures must adopt multicorrelation seismic inputs. Despite the fact that numerous real earthquake recordings have been gathered, they are not sufficient for a seismic response analysis of structures under multicorrelation seismic excitation due to the limitations caused by the number of strong-motion arrays, seismic environments and soil conditions [2]. As a result, it is of great significance to develop an acceptable multicorrelation ground motion simulation approach for the seismic response analysis of large and long-span structures. Multicorrelation ground motions can be simulated with a variety of methods [2]. Among them, the spectral representation method (SRM) is simple and accurate, so it is widely used. [PERSON] et al. [3] proposed a spectral decomposition approach for simulating the spatially varying ground motions based on the SRM. [PERSON] [4] proposed a nonstationary SRM for the simulation of multivariate nonstationary ground motions described by the evolutionary power spectral density (EPSD). [PERSON] and [PERSON] [5] proposed an SRM-based method to simulate multipoint correlated ground motions with complete nonstationarity and spectral compatibility. Based on the well-known uniform ground motion model, [PERSON] et al. [6] developed a new model for simulating multicorrelation ground motions for engineering purposes. [PERSON] and [PERSON] [7] utilized the SRM for the simulation of multicorrelation ground motions considering the site effect. However, the simulation of spatially varying ground motions may involve a large number of simulation points. In this case, the SRM is computationally expensive in memory and time. Therefore, it is necessary to accelerate the simulation with the SRM to adapt to the time history response analysis of large-scale and long-span structures under multicorrelation ground motions. The low efficiency of the SRM lies in the decomposition of the spectral matrix and the summation of cosine terms. In recent years, there has been a large number of studies devoted to improving the simulation efficiency from these two aspects. In order to accelerate the decomposition of the spectral matrix, the closed-form formation of Cholesky's decomposition is proposed to prevent the repeating decomposition at each frequency in the wind field simulation [8; 9]. [PERSON] and [PERSON] [10] further provided an enhanced closed-form solution to relax the limiting conditions and enlarge its range of application. However, it is still limited in the exponential coherence function, which may not be satisfied for seismic ground motions. To reduce the number of calculations of direct decomposition, a phase separation method is developed for avoiding complex number operations [11; 12]. Some interpolation approaches are also presented to reduce the decomposition amount in the wind field simulation [13; 14; 15] and the ground motions simulation [11; 12]. On the other hand, [PERSON] [16] introduced the fast Fourier transform (FFT) technique to expedite the summation of cosine items in a stationary simulation. [PERSON] and [PERSON] [17] further reduced the amount of the FFT in the simulation. In a nonstationary simulation, the uniformly modulated EPSD is widely used for simplification purposes and the FFT technique can be directly utilized. To better reflect the amplitude and frequency of the nonstationary feature, the time-frequency nonseparable EPSD is recommended [18]. Due to the time-frequency coupling function being involved in the summation of cosine items, the FFT is directly unavailable for the non-separable nonstationary simulation. In order to introduce the FFT, [PERSON] and [PERSON] [19] decoupled time-frequency coupling functions with the discrete Fourier transform (DFT) and digital filtering technique.