Modelling signals as being periodic is common in many applications. Such periodic signals can be represented by a weighted sum of sinusoids with frequencies being an integer multiple of the fundamental frequency. Due to its widespread use, numerous methods have been proposed to estimate the fundamental frequency, and the maximum likelihood (ML) estimator is the most accurate estimator in statistical terms. When the noise is assumed to be white and Gaussian, the ML estimator is identical to the non-linear least squares (NLS) estimator. Despite being optimal in a statistical sense, the NLS estimator has a high computational complexity. In this paper, we propose an algorithm for lowering this complexity significantly by showing that the NLS estimator can be computed efficiently by solving two Toeplitz-plus-Hankel systems of equations and by exploiting the recursive-in-order matrix structures of these systems. Specifically, the proposed algorithm reduces the time complexity to the same order as that of the popular harmonic summation method which is an approximate NLS estimator. The performance of the proposed algorithm is assessed via Monte Carlo and timing studies. These show that the proposed algorithm speeds up the evaluation of the NLS estimator by a factor of 50–100 for typical scenarios.