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 <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
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 < 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 < 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 < 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 < 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 < 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<<"\n\nfrequency and amplitude estimation results\n\n n magnitude frequency amplitude\n\n";
for (n = 370; n < 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 < N-2; h++) { // ***** bit reverse starts here ***
if (h < 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 <= j) {
j = j - k; k = k/2;
}
j = j + k;
} //****** bit reverse done ******
for (k = 1; k < N; k = k*2) {
l1 = l2; l2 = l2*2;
u1 = 1.0; u2 = 0.0;
for (j = 0; j < l1; j++) {
for (h = j; h < 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.
