Spectral analysis with high selectivity but low number of points I'm trying to create spectral plots similar to one below (it's from the datasheet of an AD1986 sample rate converter, page 8):

In the comments below the plot, it says &quot;Wideband FFT plot, 16K points&quot;.
I'm particularly interested in having high spectral selectivity (if that's the right term), with the 1kHz peak being very narrow, with 'only' having 16K FFT points.
I'm doing my exercise on a 1000.621 Hz sine wave, quantized to 16 bits with a 44.1kHz sample rate.
This code to create all the plots below can be found here.
Step 1: straightforward FFT of 1s of samples, with and without window:

There are 2 issues with this:

there's a lot of variability in the noise
there's still a lot of spreading around the peak, despite using a Blackman window

Step 2: Average 100 FFTs over 100s of samples

This averages the noise, and it overlays perfectly with the original FFT in terms of signal power.
Step 3: Do a single FFT over 100s of samples

Increasing the number of samples by 100, has reduced the width of the FFT bins.
Of course, the noise isn't averaged, so I'm back to square one on that one. What's more, the noise floor due to quantization will be lower than for step 2, due to the increased number of bins.
Step 4: Merge 4410000 FFT bins of step 3 back to 44100 FFT bins
I could average a bunch of these to reduce the noise, just like in step 2, but the AD1896 image says it's a 16K point FFT. Which means that they used a different method.
So I think that I should also be able to merge 100 FFT bins back to 1 bin, so that I get the combination of both a narrow peak and a noise floor that that is similar in value to the one in step 2.
Here's the result after adding the absolute values of each bucket of 100 FFTs bin into 1 FFT bins:

This is starting to look good, but the peak value is now well above 0dB. And the noise floor around ~-128dB is much higher than the one that was around ~-145dB of step 2.
What I really want is the result of step 2 with the same level of noise floor, but without the wide peak of step 3.
I have the feeling that adding the absolute values of each bin is not the right approach, but I don't know which would be the right approach, and why?
While I'd still like to know the answer to my question about how to merge multiple FFT bins into one, the real solution to increasing the selectivity was much easier to solve: I just had to use a different windowing function.
Here's the result of step 2 for 3 cases: no windowing function (blue), using the Blackmann window that I used earlier (green), and using a Kaiser window with a beta factor of 20 (orange):

In my numpy code, the change is as simple as this:
    w = np.blackman(len(y)/num_ffts); corr = len(w)/sum(w)

by this:
    w = np.kaiser(len(y)/num_ffts, 20); corr = len(w)/sum(w)