This paper considers black- and grey-box continuous-time transfer function estimation from frequency response measurements. The first contribution is a bilinear mapping of the original problem from the imaginary axis onto the unitdisk. This improves the numerics of the underlying Sanathanan-Koerner iterations and the more recent instrumental-variable iterations. Orthonormal rational basis functions on the unit disk are utilized. Each iteration step necessitates a minimal state-space realization with these basis functions. One such derivation is the second contribution. System identification with these basis functions yield zero-pole-gain models. The third contribution is an efficient method to express transfer function coefficient constraints in terms of the orthonormal rational basis functions. This allows for estimating transfer function models with arbitrary relative degrees (including improper models), along with other fixed and bounded parameter values. The algorithm is implemented in the tfest function in System Identification Toolbox (Release 2016b, for use with MATLAB) for frequency domain data. Two examples are presented to demonstrate the algorithm performance.