## 3 Enhanced Interpolation Approach
### Time-Frequency Decoupling with POD Interpolation  
Some approaches, including DFT [26], POD [21; 28] and non-negative matrix factorization [29], were utilized to factorize the time-frequency coupled functions \(H_{jk}(\omega,t)\) or \(\sqrt{S_{jj}(\omega,t)}\) for invoking the FFT into the nonstationary simulation. However, these schemes would be rather time-consuming for the simulations with extended durations and high sample frequencies. The POD interpolation was utilized for time-frequency decoupling in this section. Due to the auto-EPSDs \(S_{jj}(\omega,t)\) being continuous functions of \(\omega\) and \(t\) for the ground motions, the uniformly distributed time-frequency interpolation points \((\widetilde{\omega},\widetilde{t})\) were selected and determined as  
\[\widetilde{\omega}_{r}=(\omega_{u}-\omega_{1})\bigg{(}\frac{r-1}{N_{\widetilde {\omega}}-1}\bigg{)}+\omega_{1}\ r=1,2,\cdots,N_{\widetilde{\omega}} \tag{26}\]  
\[\widetilde{t}_{i}=\frac{i-1}{N_{\widetilde{t}}-1}T_{0}\ i=1,2,\cdots,N_{ \widetilde{t}} \tag{27}\]  
where \(N_{\widetilde{\omega}}\) and \(N_{\widetilde{t}}\) are, respectively, the number of the time-frequency interpolation points along frequency and time directions, \(N_{\widetilde{\omega}}\ll N\) and \(N_{\widetilde{t}}\ll M\) (\(M\) is the amount of time instants of samples). In a matrix, the square root of EPSD at the time-frequency interpolation points \(\sqrt{S_{jj}(\widetilde{\omega}_{1},\widetilde{t}_{1})}\) is given as  
\[\mathbf{L}= \begin{bmatrix}\sqrt{S_{jj}(\widetilde{\omega}_{1},\widetilde{t}_ {1})}&\sqrt{S_{jj}(\widetilde{\omega}_{1},\widetilde{t}_{2})}&\cdots&\sqrt{S _{jj}(\widetilde{\omega}_{1},\widetilde{t}_{N_{\widetilde{t}}})}\\ \sqrt{S_{jj}(\widetilde{\omega}_{2},\widetilde{t}_{1})}&\sqrt{S_{jj}( \widetilde{\omega}_{2},\widetilde{t}_{2})}&\cdots&\sqrt{S_{jj}(\widetilde{ \omega}_{2},\widetilde{t}_{N_{\widetilde{t}}})}\\ \vdots&\vdots&\ddots&\vdots\\ \sqrt{S_{jj}(\widetilde{\omega}_{N_{\widetilde{\omega}}},\widetilde{t}_{1})} &\sqrt{S_{jj}(\widetilde{\omega}_{N_{\widetilde{\omega}}},\widetilde{t}_{2}) }&\cdots&\sqrt{S_{jj}(\widetilde{\omega}_{N_{\widetilde{\omega}}},\widetilde{ t}_{N_{\widetilde{t}}})}\\ \end{bmatrix}  \tag{28}\]This matrix's columns can be thought of as column vectors. A set of optimal orthogonal bases \(\Phi=\left[\mathbf{\Phi}_{1},\mathbf{\Phi}_{2},\cdots,\mathbf{\Phi}_{N_{\bar{ \omega}}}\right]_{N_{\bar{\omega}}\times N_{\bar{\omega}}}\) for these vectors can be obtained by the following eigenvector decomposition  
\[\mathbf{R}\mathbf{\Phi}_{q}=\lambda_{q}\mathbf{\Phi}_{q}\,,\,q=1,2,\cdots,N_{ \bar{\omega}} \tag{29}\]  
where \(\mathbf{\Phi}_{q}\) and \(\lambda_{q}\) are, respectively, the \(q\)th eigenvector and eigenvalue; and \(\mathbf{R}\) is the correlation matrix of these column vectors, which could be calculated using  
\[\mathbf{R}=\frac{1}{N_{\bar{t}}}\mathbf{L}\mathbf{L}^{T} \tag{30}\]  
Additionally, the principal coordinates could be determined by  
\[\mathbf{a}_{q}^{T}=\mathbf{\Phi}_{q}^{T}\mathbf{L},\,q=1,2,\cdots,N_{\bar{t}} \tag{31}\]  
where \(\mathbf{a}_{q}\) is the \(q\)th projection vector of the column vectors. By reorganizing the eigenvalues in descending order, the eigenvalues at lower orders containing more energy were reserved. \(\mathbf{L}\) could be roughly recreated by  
\[\mathbf{L}\approx\sum_{q=1}^{N_{\Phi}}\mathbf{\Phi}_{q}\mathbf{a}_{q}^{T}\,,\,N _{\Phi}\ll N_{\bar{t}} \tag{32}\]  
where the first \(N_{\Phi}\) items were chosen to meet the accuracy requirement. Therefore, those coupled functions were factorized as  
\[\sqrt{S_{jj}(\widetilde{\omega},\widetilde{t})}\approx\sum_{q=1}^{N_{\Phi}}a_ {q}^{jj}(\widetilde{t})\mathbf{\Phi}_{q}^{jj}(\widetilde{\omega}),\,j=1,2, \cdots,n \tag{33}\]  
where \(a_{q}^{jj}(\widetilde{t})\) and \(\Phi_{q}^{jj}(\widetilde{\omega})\) are, respectively, the discrete functions corresponding to \(\mathbf{a}_{q}\) and \(\mathbf{\Phi}_{q}\).