Since many years, the pioneering work of Hu and Tonder is used to generate rough surfaces with prescribed statistical moments (skewness and kurtosis) along with spatial properties (correlation lengths). The present work enlightens the drawbacks of this method and it proposes an original approach based on a hybrid analytical/numerical method. Simulations are conducted on very different surface specimens and the method is validated over a wide range of statistical moments. The results are obtained with high accuracy (beyond what is usually needed) and very short computing times (the order of a second).