OP did allow the base frequency to vary when he did the fit. So if you found the first residual and then used the same varying-frequency fit, you might not get

OP did allow the base frequency to vary when he did the fit. So if you found the first residual and then used the same varying-frequency fit, you might not get an exact harmonic of the base frequency. That would not be the Fourier transform! But suppose you fixed the base frequency (which wouldn't be a bad idea). You still seem to have a nonlinear fit, because the phase (the "f" in the model equation in OP) is buried inside sin(). Why are we needing a nonlinear function-fitting process when the Fourier transform is linear? Of course, you can bring the phase outside by adding in a cos() term with its own amplitude. Now instead of the phase "f" you have an interplay between the amplitude of the sin and cos terms, and those amplitudes are linear in the data. The resulting fit (or recursive sequence of fits) would indeed be the Fourier transform. The key property is the orthogonality of the various harmonics. That's what allows the sequence of single-frequency fits to not step on each other.