A numerical procedure for the simulation of non-Gaussian surfaces has been developed. Fast Fourier Transform is used to generate the surfaces. The amplitude part is constructed from given spectral density or auto-correlation function. The phase part is obtained by using Johnson translator system. By iteration, this method can simulate surfaces with given skewness and kurtosis and with spectral density or auto-correlation function. (C) 2003 Published by Elsevier Ltd.