Most accurate way to find RMS of signals fundamental
This is the start of an answer containing my initial thoughts that may be helpful (I wrote quickly due to limited time so may have errors, I will review in more detail later but can at least offer some direction):

One approach is to use this resonator that was developed in this post to isolate your frequency of interest: FFT analysis for Vibration Signal. The issue you will have is allowing for enough time for the filter to settle to it's final amplitude on the output (which is based on the bandwidth of the filter, translated to time constants in time, and depending on your accuracy needed is how many taus you will need to wait until you have signal that you can take the rms value of. I have not confirmed this but it also seems that you could be clever and use the known information of settling time so that you can utilize the initial settling signal's information as part of your answer (weighting the rms calculation by the result in an "optimal combining" fashion. 

You could also do FFT's with windowing in which case refer to this paper by fred harris regarding the effect the window will have on your rms calculation: fred harris On the Use of Windows

Note with an FFT approach that each bin in the FFT is a filter that is responsive to the signal anywhere within its filtering bandwidth. Without any further windowing, it is a "rectangular" window which has the highest frequency selectivity but lowest dynamic range, as evidenced in the plot below of a simple 4 point DFT. Any signal power that "leaks" into the bandwidth of an FFT bin will become part of your RMS calculation. The filter responses as shown below approach Sinc functions as the number of points increase, with the null to null spacing of the main lobe for each filter as 1/N where N is the number of points in the FFT. When you use a window, the sidelobes significantly decrease, but at the expense of the main lobe. 



Given you have spurs that you are concerned with that are only 6 Hz away from a signal that is 16KHz, if you assume you used a window that increased your main lobe 3x (for example) to get your desired sidelobe attuenuation, meaning was 3/N null to null, Then you would need a minmumum FFT length on the order of $3\times 16000/6=8000$ which is not so bad.  (you can take the FT of various windows to see for yourself what each would do- my favorite is the Kaiser window because you can adjust selectivity and bandwidth with the additional $\beta$ factor).

Note these particular items that will effect you if using a window: coherent gain, noise gain, scalloping loss.

Coherent Gain: The standard deviation of a single tone will be modified by the coherent gain of the window. The coherent gain is found by summing all the samples in the window itself:

$$G_c=\sum_0^{N-1}w_n $$ 

For instance, the coherent gain of a rectangular window is N, so the net gain of a window compared to the rectangular window is $G_c/N$, which if you try yourself with common windows, you will matches what is tabulated in table 1 in fred harris' paper.

The noise gain however is how the standard deviation of the noise components in your signal will scale, where noise is anything that is uncorrelated from sample to sample. The noise gain due to a window increases as the sum of the squares:

$$G_n=\sqrt{\sum_0^{N-1}w_n^2}$$

So note the processing gain that results, for the rectangular window the signal will increase by $20log_{10}(N)$ while the noise will increase by $10log_{10}(N)$ for a net gain of $10log_{10}(N)$ in SNR. You can use this and the measurement of the variance of your entire signal for a more accurate estimate of the signal component since the variance measurement is S+N, and your post FFT measurement after scaling will be S+N/K where K is the resulting processing gain for your window used.

The final item to mention is "scalloping loss" which occurs if your signal of interest is not in the center of the bin (observe the immediate roll-off of a sinc function). You can either account for this mathematically in your coherent gain if you have specific knowledge of your signal, or you can minimize the effect of scalloping loss by using overlap-add techniques (detailed in fred harris' paper).

Having written this out, between the FFT appraoch and my first suggestion of using the tuned 2nd order resonator, I favor the FFT approach as being simplest and most straightforward. You also have significantly more samples than the minimum needed to reject your interference signals, which are sufficiently low so you should be able to extract a good estimate of the standard deviation of your signal. Note the overall SNR is also of interest as you can use that to establish the confidence interval for your result. For this see Estimating confidence intervals for the mean value