In phase-shifting fringe projection profilometry, fringe harmonics have been recognized as one of the main error-inducing factors. Generally, the response of a phase-shifting algorithm to fringe harmonics strongly depends on the used phase shifts and is usually unpredictable, especially when using nonuniform phase shifts. For this reason, it is difficult to eliminate the phase-measuring errors caused by fringe harmonics in a general case, for example, when the phase shifts are not uniform but arbitrarily valued. To overcome this problem, this paper analyzes the phase error function related to each fringe harmonic under the condition of using arbitrary phase shifts, reveals the partial orthogonality of these functions, and then derives an algorithm for decomposing the harmonics-induced phase errors. In the implementation, this algorithm calculates a coarse phase map first in the least squares sense from captured fringe patterns, and then estimates the coefficients of fringe harmonics from this calculated phase map by use of the partial orthogonality of the error functions. By using the estimated harmonic coefficients, the phase map is updated, thus having improved accuracies so that the harmonics-induced phase errors are compensated for through an iterative procedure. The simulation and experimental results demonstrate this method to be effective and efficient in measuring fringe phases.