My approach to this is as follows (just due to both its simplicity, high accuracy and applicability to the measurement of any tones regardless of sampling rate and frequency relations):
Yes use windowing. I like the Kaiser window due to its excellent time-bandwidth properties, and that you can easily trade resolution bandwidth (how closely we can resolve individual tones) with dynamic range (how much difference in dB you can see a weaker tone next to a stronger tone that is beyond the resolution bandwidth).
To measure the tone, both magnitude and phase (relative to your sampling clock which serves as the phase reference) with very high accuracy, do the following:
Ensure the capture duration in time is long enough to meet the resolution bandwidth requirement. If there was no further windowing applied (native rectangular window), the resolution bandwidth in Hz (RBW) otherwise referred to as equivalent noise bandwidth (ENBW),is simply the inverse of the total capture time in seconds - windowing stretches this out (see my answer to DSP.SE# 94702 for the exact formula for RBW for a given window) but a 1.5 to 2x increase is a reasonable assumption covering most practical cases. Therefore, if you want to resolve tones that are spaced 1 KHz apart, the capture duration needs to be sufficiently larger than 2 ms. (The longer the better as long as the signal is sufficiently stationary over this duration, this is just setting the minimum requirement). Sampling at 192 KHz to get 10 ms of data (for example) is only 1920 samples, so this could easily be much longer.
Window the captured signal over that total time duration.
Zero pad the resulting windowed signal out to 10x its length (or longer) and then take the FFT (basically in MATLAB, Octave or Python scipy.signal, just specify a longer length for the FFT and it will do the zero padding. This results in minimizing the scalloping loss and any narrow band tone that is individually within the resolution bandwidth of the measurement will be accurately captured both in frequency and magnitude.
Divide the FFT by the sum of the window, and the resulting magnitude of the two peaks in the FFT for every sinusoid will be accurately represented according to the two exponential tones in Euler's formula as given below (representing the magnitude and starting phase of the individual exponential tones):
$$A\cos(\omega t + \phi) = \frac{A}{2}e^{(j\omega t + \phi)} + \frac{A}{2}e^{(-j\omega t - \phi)}$$
I don't think it gets easier than this, and the accuracy is only limited by how much you zero pad (which does the interpolation) and what the SNR is for that tone within the resolution bandwidth of the measurement.
Note: I want to stress how important the windowing is with this approach (and of little consequence to the SNR of the result). Try what I suggest without windowing when measuring low frequency tones or tones near Nyquist. Without windowing the positive and negative frequency tones will interact and degrade the result and windowing with the sidelobe reduction significantly reduces that error, quickly to negligible results. The resulting method is then robust and independent of frequency for generalized measurement cases.
For further details on how to accurately measure tones versus how to accurately measure noise densities see DSP.SE #87723 . Also I will be presenting this very topic in much more detail at the 2024 DSP Online Conference that will take place at the end of next month!