The problem of identifying vibrating systems from frequency response data is considered. The unknown model of a vibrating system is parameterized as a sum of second order transfer functions. The model can be guaranteed to be stable if the denominator coefficients in the second order transfer functions are constrained to be positive. The nonlinear least squares cost function is adopted and the numerator and denominator coefficients are estimated in a separable manner. The denominator coefficients are estimated using simulated annealing and the numerator coefficients are obtained by solving the linear least squares problem. A numerical example is given to demonstrate the validity of the proposed algorithm.