In fringe projection profilometry, the purpose of using two- or multi-frequency fringe patterns is to unwrap the measured phase maps temporally. Using the same patterns, this paper presents a least squares algorithm for, simultaneously with phase-unwrapping, eliminating the influences of fringe harmonics induced by various adverse factors. It is demonstrated that, for most of the points over the measured surface, projecting two sequences of phase-shifting fringe patterns having different frequencies enables providing sufficiently many equations for determining the coefficient of a high order fringe harmonic. As a result, solving these equations in the least squares sense results in a phase map having higher accuracy than that depending only on the fringe patterns of a single frequency. For the other few points which have special phases related to the two frequencies, this system of equations becomes under-determined. For coping with this case, this paper suggests an interpolation-based solution which has a low sensitivity to the variations of reflectivity and slope of the measured surface. Simulation and experimental results verify that the proposed method significantly suppresses the ripple-like artifacts in phase maps induced by fringe harmonics without capturing extra many fringe patterns or correcting the non-sinusoidal profiles of fringes. In addition, this method involves a quasi-pointwise operation, enabling correcting position-dependent phase errors and being helpful for protecting the edges and details of the measurement results from being blurred.