This paper gives an algorithm for identifying spectral densities using orthonormal basis functions. Mathematically, this amounts to identifying a time-invariant linear SISO system with the additional constraint that the transfer function should be positive-real. Thus, we solve the long-standing problem of how to incorporate this positivity constraint while using orthonormal basis functions. The procedure is a variant of the THREE algorithm introduced by Byrnes, Georgiou and Lindquist. The relation between and numerical properties of the proposed and the THREE algorithms are discussed. The orthonormal basis functions are better scaled for a concentrated pole selection in the basis, which increases the accuracy of the estimates. A numerical example which highlights this phenomenon and illustrates the algorithm is given.