best window for close frequency components
I'm presuming that you have some or all 9 of the input frequencies available at one time, and you just want to find out which of them is present (you didn't specify in your problem statement).

Two methods come to mind: 1) compute DFTs at the frequencies of interest, or 2) use the Gaussian ratio technique, as described in reference_1 and reference _2.

Regarding 1) above: When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer.  In other words, you can compute a DFT bin centered at ANY frequency, not just those that correspond to the usual k = integer of your FFT.  For example, given your sample rate and N, a frequency of F = 16351 would correspond to a DFT bin computed at k = (N*F)/sample_rate = 1024*16351/44100 = 379.669 (ie: between bins 379 and 380 of your FFT). 

So for the DFT, X(F) = (sum over n){x[n]*e[-j*twopi*n*k/N]}, where X[F] is the complex result, N = 1024, F = 16351, k = 379.669, and n = 0 to N-1.

Some of the other values for your problem are:

Frequency   and value of   k

16,566  ==> 384.662

16,781  ==> 389.659

16,996  ==> 394.646

I'll leave it to you to figure out the rest.  And, as you can see, you've got about 5 bins of separation between the input frequencies.  That may, or may not be enough based on your SNR and method used.  And you might want to use N = 1025 instead – it makes your input frequencies much more 'bin-centered.'

Method 2) above uses an FFT, which, given the right circumstances, can be extraordinarily accurate when computing frequencies.  Note that in the first reference above, by using a Gaussian window, and then taking the log ratio of amplitudes from adjacent bins, you are actually solving for frequency, subject to the requirements of the technique.  Some code (DevC++ compiler) for your specific problem, followed by the programs' output is shown below (explanation to follow):

#include &lt;cstdio&gt;
#include &lt;cstdlib&gt;
#include &lt;iostream&gt;
#include &lt;cmath&gt;
using namespace std;
void fft_recur(long, double* r, double* i);
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O
const long N = 1024;   double sample_rate = 44100., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ;
long n;      double  t, twopi = 2*3.141592653589793238;
double amp_sq[N], F, fbin, amp, power, c = 2.3;

for (n = 0; n &lt; N; n++)  {  // generate test data
    t = twopi*n/sample_rate;
    r[n] = 1.*cos(16351.*t)+1.*cos(16996.*t)+1.*cos(17211.*t)+1.*cos(17426.*t) ;
} // end for

//Gaussian window as computed in [MCEACH94] program
double z;
for (n = 0; n &lt; N/2; n++) {
    z = (.5+(N/2)-1.-n)*(twopi/2.)/N;
    w[n] = exp(-c*z*z);
    w[N-1-n] = w[n];
}

//multiply input data by window points
for (n = 0; n &lt; N; n++)  {
    r[n] = w[n]*r[n];    // we're overwriting r[n] with the windowed r[n]
} //end for

fft_recur(N, r, i);       // note: fft is not scaled

// compute amplitude squared
for (n = 0; n &lt; N/2+1; n++)  {
    amp_sq[n] = r[n]*r[n]+i[n]*i[n];
} //end for

// compute natural log of FFT outputs
for (n = 0; n &lt; N/2+1; n++)  {
    r[n] = .5*log( amp_sq[n] ) ; // r[n] is now the natural log of the amplitude of FFT bin 'n'
} // end for                     // .5*log(A squared) = log(A)

cout&lt;&lt;"\n\nfrequency and amplitude estimation results\n\n n        magnitude      frequency       amplitude\n\n";
for (n = 370; n &lt; 415; n++)  {
    fbin = n+.5+.5*c*(r[n+1]-r[n]);   // r[n+1] and r[n] are ln(b) and ln(a)
    F = (sample_rate/N)*fbin;
    amp = sqrt(c)*exp(r[n]+((fbin-n)*(fbin-n))/c)*(2.0*sqrt(twopi/2.)/N) ;
    printf("%2d\t%9.5f\t%9.5f\t%9.5f\n",n,amp_sq[n],F,amp);
} //end for

system ("PAUSE");
return 0;
} // end main
//******************** fft_recur ***********************
void fft_recur(long N, double *r, double *i)  {
long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1;
double c = -1.0, s = 0.0, t1, t2, u1, u2;

for (h = 0; h &lt; N-2; h++) {    // ***** bit reverse starts here ***
    if (h &lt; j) {
       t1 = r[h]; r[h] = r[j]; r[j] = t1; 
       t2 = i[h]; i[h] = i[j]; i[j] = t2;
    }
    k = i2;
    while (k &lt;= j) {
          j = j - k;      k = k/2;
    }
    j = j + k;
}   //****** bit reverse done ******

for (k = 1; k &lt; N; k = k*2) {
    l1 = l2;     l2 = l2*2;
    u1 = 1.0;    u2 = 0.0;
    for (j = 0; j &lt; l1; j++) {
        for (h = j; h &lt; N; h = h + l2) {
            i1 = h + l1;
            t2 = (r[i1] - i[i1])*u2 ;
            t1 = t2 + r[i1]*(u1 - u2) ;
            t2 = t2 + i[i1]*(u1 + u2) ;
            r[i1] = r[h] - t1;
            i[i1] = i[h] - t2;
            r[h]  = r[h] + t1;
            i[h]  = i[h] + t2;              
        } // end for over h
        t1 = u1 * c - u2 * s;
        u2 = u1 * s + u2 * c;
        u1 = t1;  //x = u1 - u2;    y = u1 + u2;
    } // end for over j
    s = - sqrt((1.0 - c) / 2.0);
    c =   sqrt((1.0 + c) / 2.0);
} // end for over k

} // end fft


frequency and amplitude estimation results

 n        magnitude      frequency       amplitude

370       0.01156       15959.53142       0.00065
371       0.01328       16002.86947       0.00070
372       0.01542       16046.20441       0.00076
373       0.01810       16089.42371       0.00082
374       0.02138       16136.92108       0.00095
375       0.03020       16209.62311       0.00211
376       0.14116       16370.55303       3.20360
377      77.00512       16349.81449       0.95808
378     3193.91932      16351.21160       1.00396
379     24622.71009     16350.90700       0.99982
380     32939.04014     16351.10559       0.99850
381     7803.22610      16350.67285       1.01020
382     319.11953       16353.81879       0.87280
383       2.60322       16324.16039       7.57917
384       0.00113       16605.10228       0.00051
385       0.00724       16606.09917       0.00052
386       0.00851       16650.89175       0.00058
387       0.01072       16694.51462       0.00065
388       0.01381       16737.66957       0.00074
389       0.01787       16785.61197       0.00090
390       0.02814       16862.78585       0.00251
391       0.17568       17013.74935       2.83597
392      85.58946       16994.84074       0.95935
393     3413.62331      16996.22071       1.00408
394     25288.20406     16995.88948       0.99974
395     32472.56144     16995.79357       1.00043
396     7296.89833      16957.17339       4.09340
397      60.55471       17218.13259       1.24775
398     3330.74781      17212.33274       1.01810
399     25463.82887     17210.95691       0.99978
400     32350.54927     17210.69761       1.00169
401     7144.91966      17170.91626       4.35710
402      55.60789       17436.25628       1.40301
403     3423.65695      17427.15709       1.01474
404     25642.98140     17426.03425       1.00002
405     32244.80235     17425.93874       1.00073
406     7095.46971      17426.35995       0.98918
407     278.99401       17421.20777       1.27242
408       1.56515       17494.25360       0.02622
409       0.02946       17619.31213       0.00091
410       0.01521       17677.62446       0.00071
411       0.01452       17719.83432       0.00069
412       0.01340       17762.64014       0.00066
413       0.01224       17805.61294       0.00063
414       0.01113       17848.67454       0.00060


First, 4 tones of the input are generated (I use twopi when computing time because it's easier than writing it 4 times into the cosines).  Then, the Gaussian window is generated, after which it is applied to the input (by the way, the window shown here, which is from reference_1, is slightly wrong – I don't want to use the 'corrected' version of it because I don't want to write a page describing how I derived it).  Then we take a 1024 point FFT, followed by the (non-normalized) calculation of amplitude squared.  Then comes the calculation of the natural log of each FFT bin.  Finally, I compute the frequencies and amplitudes via the equations in reference 1 and reference 2above.

As you can see from the results following the code, we have an isolated peak between bin 379 and 380 from the first input frequency, which is estimated as 16350.90700 and an amplitude of .99982.  We can also see that the other 3 inputs exhibit a little 'peak merging' from bins 392 to 407, but the peaks and associated estimated frequencies and amplitudes are quite accurate.

The results could be even more accurate if I had: 1) adjusted parameters, 2) used the 'correct' window, 2) used a non-recursive FFT, or 3) used one of the variations of the technique (which I am not going to show here).

You might also notice that, as your value of the parameter 'c' gets smaller ('c' is related to the 'beamwidth' of the Gaussian window), you can actually improve some of your estimates, but it can cause some results to blow up due to numerical sensitivities.  The effective 'Q' or 'gain' of the Gaussian ratio technique is on the order of 5,000 – 10,000, and it can be much higher if you know how to adjust the parameters properly.

Of course, if your SNR is too low, you may have to increase N to get more integration time, subject to the limits of signal stability.