Inspired by the recent promising developments of Bayesian learning techniques in the context of system identification, this paper proposes a Transfer Function estimator, based on Gaussian process regression. Contrary to existing kernel-based impulse response estimators, a frequency domain approach is adopted. This leads to a formulation and implementation which is seamlessly valid for both continuous- and discrete-time systems, and which conveniently enables the selection of the frequency band of interest. A pragmatic approach is proposed in an output error framework, from sampled input and output data. The transient is dealt with by estimating it simultaneously with the transfer function. Modelling the transfer function and the transient as Gaussian processes allows for the incorporation of relevant prior knowledge on the system, in the form of suitably designed kernels. The SS (Stable Spline) and DC (Diagonal Correlated) kernels from the literature are translated to the frequency domain, and are proven to impose the stability of the estimated transfer function. Specifically, the estimates are shown to be stable rational functions in the frequency variable. The hyperparameters of the kernel are tuned via marginal likelihood maximisation. The comparison between the proposed method and three existing methods from the literature - the regularised finite impulse response (RFIR) estimator, the Local Polynomial Method (LPM), and the Local Rational Method for Frequency Response Function estimation - is illustrated on simulations on multiple case studies. (C) 2016 Elsevier Ltd. All rights reserved.