In this paper, we propose a novel algorithm for the identification of Hammerstein systems. Adopting a Bayesian approach, we model the impulse response of the unknown linear dynamic system as a realization of a zero-mean Gaussian process. The covariance matrix (or kernel) of this process is given by the recently introduced stable-spline kernel, which encodes information on the stability and regularity of the impulse response. The static non-linearity of the model is identified using an Empirical Bayes approach, i.e. by maximizing the output marginal likelihood, which is obtained by integrating out the unknown impulse response. The related optimization problem is solved adopting a novel iterative scheme based on the Expectation-Maximization (EM) method, where each iteration consists in a simple sequence of update rules. Numerical experiments show that the proposed method compares favorably with a standard algorithm for Hammerstein system identification.