One of the most advanced classes of techniques for ground response analysis is based on the use of Transfer Functions. They represent the ratio of Fourier spectrum of amplitude motion at the free surface to the corresponding spectrum of the bedrock motion and they are applied in frequency domain usually by FFT method. However, Fourier spectrum only shows the dominant frequency in each time step and is unable to represent all frequency contents in every time step and this drawback leads to inaccurate results. In this research, this process is optimized by decomposing the input motion into different frequency sub-bands using Wavelet Multi-level Decomposition. Each component is then processed with transfer Function relating to the corresponding component frequency. Taking inverse FFT from all components, the ground motion can be recovered by summing up the results. The nonlinear behavior is approximated using an iterative procedure with nonlinear soil properties. The results of this procedure show better accuracy with respect to field observations than does the Conventional method. The proposed method can also be applied to other engineering disciplines with similar procedure.