Bottom Line: The OP is experiencing the effects of spectral leakage and the primary motivation for windowing with the Discrete Fourier Transform. The reason the windowing approach achieves much better estimation accuracy is explained as follows:
The OP's equation is a sum of real sinusoids, and his estimator to determine those sinusoids is a form of the Discrete Fourier Transform (DFT). Without additional windowing (which is a rectangular window) it is well understood that the DFT suffers from a considerable amount of spectral leakage, in that the Kernel of the rectangular window (the Discrete-Time Fourier Transform or DTFT of the rectangular window which is of the form of an aliased Sinc function) convolves in frequency with the underlying two impulses in frequency of each real sinusoid. Since the aliased Sinc function decreases relatively gradually with frequency offset from the location of the impulses representing the tones, and has such high side-lobes, many other frequency bins are affected causing errors in their own measurements (and a significant decrease in dynamic range or the ability to detect a lower level tone in the presence of a strong one). Windows are selected that have the property of much lower side-lobe levels, and roll-off much faster in frequency, which therefore causes any given tone to have much less impact on the estimation of other tones present.
This is made much clearer if we observe the DTFT instead of the DFT, as the DTFT is a continuous function of frequency and the DFT is samples on this continuous function. By observing the DTFT we can more easily get the full intuitive picture of what is being described above. Further we also see how approximating the DTFT (through zero-padding prior to doing the DFT) can also be used to approximate the frequencies). I show this with plots of the OP's original frequencies where we can see the effect of the windowing directly:
First consider the OP's example of 1024 samples with just F1 alone with and without additional windowing. First without additional windowing (rectangular window) The log magnitude plot of the DTFT with the corresponding 1024 DFT samples of just F1 is given in the plot below:
This is specifically the DTFT of
$$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$$
for n = 0 to 1023. And f1 = 0.123
(Note: the above is just a single tone at $f1$ with additional phase shift as can be derived from the relationship $a\cos(x)+b\cos(x) = \sqrt{a^2+b^2}\cos(x-\tan^{-1}(b/a))$
Magnitude Spectrum of F1 Signal Only - Rectangular Window
If we zoom in on one of the spectral peaks, we can clearly see the Sinc shaped Kernel (in this case given the number of samples the aliasing would be very small so this is a very close approximation of an actual Sinc) that has convolved with the underlying impulse that would be at the F1 frequency location, and we see the spectral leakage and how it affects the other frequency locations. F1 as 0.123 is very close to bin number 126 as given by $1024 \times .123 = 125.952$.
Zoom In on F1 Peak Location
Below I show the two frequency components as the OP has used, but still as separate signals to demonstrate how the spectral leakage of each will effect the measurement of the other.
This is specifically
$$sig1 = 1.0\cos(2\pi f_1 n) + 0.5\sin(2\pi f_1 n)$$
$$sig2 = 0.05\cos(2\pi f_2 n) + 0.01\sin(2\pi f_2 n)$$
For n = 0 to 1023 and f1 = 0.123 and f2 = 0.312874
Magnitude Spectrum of F1 and F2 Signals (Separate) - Rectangular Window
Notice if we zoom in on the peak for F2 how significant the leakage from F1 is will effect the estimate of the magnitude for F2, given that the leakage from F1 is only 29 dB lower! If the magnitude of F2 was made lower, the resulting error would be even more pronounced. (Also notice how with the DTFT we can locate the F2 frequency precisely at its fractional bin location of $1024 \times 0.312874 \approx 320.3830$ which is where the peak of the orange curve in the plot below would be located, but once combined with the F2 signal given by the blue curve, this would cause a frequency error as well since the peak of the blue curve happens slightly to the right of the orange F1 curve.
Zoom In on F2 Peak Location - With Rectangular Window
If we then compute the DTFT of the composite waveform with both F1 and F2, we see the resulting amplitude and frequency error in vicinity of the F2 Peak location as shown in the plot below (where Green is the DTFT of the composite waveform summing both F1 and F2 components).
This error is predicted by the coherent sum of the two waveforms and thus can be predicted from rectangular window Kernel and knowledge of the relative amplitude and frequency offsets of the two signals. Given the waveforms are complex the phase of the components will impact the summation result (thus the inteference could feasibly reduce or increase the amplitude of the affected component).
In this example the following were the signal contributions at the F2 peak location (not in dB scale as above to show the summation directly).
F2 alone at bin 320.38: magnitude = 52.10 (given as 52.00 -j3.26)
F1 alone at bin 320.38: magnitude = 1.70 (given as 1.42-j0.94)
F1 + F2 at bin 320.38: magnitude = 53.58 (given as 53.42-j4.19)
Where F1 + F2 as measured is predicted by the sum of F1 and F2 alone when the complex components are added.
Finally if we observe the above plot with the use of windowing (we'll use the same Hann window the OP used) we can see the significant reduction in spectral leakage and then exactly how this would contributed to a better estimation of both the frequency and magnitude of the underlying signals.
Magnitude Spectrum of F1 and F2 Signals (Separate) - Hann Window
Zooming in on the F2 spectral peak now shows that the leakage from F1 is now over 100 dB lower, but we also see the drawback of windowing- the main lobe of the Kernel is wider leading to reduced frequency selectivity. (The rectangular window has the best selectivity with the narrowest possible main lobe, but the worst dynamic range due to the high sidelobe level and slow roll-off).
Zoom In on F2 Peak Location - With Hann Window
For more info on DTFT vs DFT and approximating the DTFT with zero padding see Getting the DTFT from the DFT samples and For 2D signals can it be said that the frequency response is the same as the Fourier transform?
For additional details on windowing see:
Why would one use a Hann or Bartlett window?
For additional details on spectral leakage see: Intuition for sidelobes in FFT
Suggestions for Efficient Processing
Assuming the objective is to only find the magnitudes for each of the tones, if 0.01 dB accuracy is sufficient, a suggestion is to use a flattop window, which provides nearly 100 dB of sidelobe rejection while maintaining a very flat main lobe in frequency extending out to $\pm 0.5$ bins. This means no further interpolation would be necessary to estimate the magnitude of a tone located at a fractional bin frequency to get an estimate within 0.01 dB of the actual magnitude.
Below shows an example of the Kernel (DTFT) of a flattop window with just 128 samples, showing nearly 100 dB of sidelobe rejection reached at a 5 bin offset.
Zooming in to see the flatness, we see that with 128 points we maintain less than 0.01 dB flatness over $\pm 0.5$ bins.
If the $0.01$ dB magnitude accuracy was acceptable, the solution to minimize processing would be to choose the maximum DFT size required to ensure no two tones are spaced closer than 5 bins based on the minimum tone proximimity (a 58.7 dB leakage which occurs at approximately 4.6 bins would cause a maximum amplitude error of 0.01 dB). This would set the minimum number of samples to process, after which can be multiplied by the flattop window and fft'd where the maximum bins would be sufficient to provide the amplitude estimates.
Here is an interesting paper with further information on designing the flattop window to any passband error specification (at the expense of sidelobe performance), as well as other optimized flattop solutions: Design and Improvement of Flattop Windows with Semi Infinite Optimization.
