This paper describes a new algorithm for initializing and estimating Wiener-Hammerstein models which consist of two linear parts with a static nonlinearity in between. The algorithm makes use of the best linear model of the system, which is a consistent estimate of the systems dynamics for Gaussian excitations. The linear model is split in all possible ways into two sub-models. For all possible splits, a Wiener-Hammerstein model is initialized which means that a nonlinearity is introduced in between the two sub-models. The linear parameters of this nonlinearity can be estimated using least-squares. All initialized models can then be ranked depending on the fit. Typically, one is only interested in the best one, for which all parameters are fitted using prediction error minimization. The paper explains the algorithm in detail and consistency of the initialization is proven. Computational aspects are investigated, showing that in most realistic cases, the number of splits of the initial linear model remains low enough to make the algorithm useful. The algorithm is illustrated on an example where it is shown that the initialization is a tool to avoid many local minima. (C) 2011 Elsevier Ltd. All rights reserved.