The substructure technique based on the flexibility formulation in the time domain, which takes into account linear unbounded soils including the local nonlinearity of near field, has been proposed by Wolf, Hayashi and Motosaka. These formulations involve the input motion and dynamic flexibility coefficients, which are not formulated in the time domain but are calculated as the inverse Fourier transform of the corresponding values in the frequency domain for the sake of simplicity in numerical treatments. Thus these methods depend on the numerical treatments of FFT algorithm which may affect the stability of response analyses and the accuracy of the results. In this paper, we discuss the numerical treatments of FFT algorithm based on distribution theory considering causality condition in calculating the inverse Fourier transform of flexibility coefficients and input motion. Then we propose the response analysis method which is advantageous in computational stability, accuracy and efficiency compared with the methods of previous studies.