Algorithms for system identification applying throughout Fast Fourier Transform (FFT) to the major calculating operations are introduced. It is shown that by using data of about as twice length as system settling time and by truncating the incorrect correlation functions resulting from them, errors owing to finiteness of data can be avoided. It is shown that so as to suppress the effects owing to statistical fluctuation of input data or output noise, superposition of data in frequency domain is effective, and also the damping terms of poles or zeros can be efficiently evaluated by utilizing the phase change of the spectra of the impulse response sequence. The proposed method can be efficiently applied to relatively higher order systems or relatively rapidly time-variant systems because of high accuracy and high speed processing of FFT. Moreover, it needs not to assume the order of the system a priori, and yields a reasonable lower order approximating system in itself.