我假设您一次有部分或全部 9 个输入频率可用，而您只想找出其中存在哪些（您没有在问题陈述中指定）。
想到了两种方法：1) 在感兴趣的频率处计算 DFT，或 2) 使用高斯比技术，如参考 1 和参考 _2 中所述。
关于上述1）：当摩西下山时，他没有携带石板，上面命令在公式中：F = k *（df）= k * sample_rate / N，k必须是整数。换句话说，您可以计算以任何频率为中心的 DFT bin，而不仅仅是那些对应于通常 k = FFT 整数的 DFT bin。例如，给定您的采样率和 N，频率 F = 16351 将对应于在 k = (N*F)/sample_rate = 1024*16351/44100 = 379.669 计算的 DFT bin（即：在 bin 379 和 380 之间你的 FFT）。
所以对于 DFT，X(F) = (sum over n){x[n]*e[-j*twopi*n*k/N]}，其中 X[F] 是复数结果，N = 1024， F = 16351，k = 379.669，n = 0 到 N-1。
您的问题的其他一些值是：
k的频率和值
16,566 ==> 384.662
16,781 ==> 389.659
16,996 ==> 394.646
我会把它留给你来弄清楚其余的。而且，如您所见，输入频率之间有大约 5 个间隔。根据您的 SNR 和使用的方法，这可能不够，也可能不够。您可能想改用 N = 1025 - 它使您的输入频率更加“以 bin 为中心”。
上面的方法 2) 使用 FFT，在适当的情况下，它在计算频率时可以非常准确。请注意，在上面的第一个参考文献中，通过使用高斯窗口，然后从相邻的 bin 中获取幅度的对数比，您实际上是在求解频率，这取决于技术的要求。针对您的特定问题的一些代码（DevC++ 编译器），以及程序的输出如下所示（解释如下）：
#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
首先，生成 4 个音调的输入（我在计算时间时使用 twopi，因为它比将其写入余弦 4 次更容易）。然后，生成高斯窗口，然后将其应用于输入（顺便说一下，此处显示的窗口来自 reference_1，有点错误——我不想使用它的“校正”版本，因为我不想写一个页面来描述我是如何得出它的）。然后我们进行 1024 点 FFT，然后进行幅度平方的（非归一化）计算。然后是计算每个 FFT bin 的自然对数。最后，我通过上面参考文献 1和参考文献 2中的方程计算频率和幅度。
从代码后面的结果可以看出，我们在第一个输入频率的 bin 379 和 380 之间有一个孤立的峰值，估计为 16350.90700，幅度为 0.99982。我们还可以看到，其他 3 个输入从 bin 392 到 407 表现出一点“峰值合并”，但峰值和相关的估计频率和幅度非常准确。
如果我有：1）调整参数，2）使用“正确”窗口，2）使用非递归 FFT，或 3）使用该技术的一种变体（我不是将在这里展示）。
您可能还注意到，随着参数 'c' 的值变小（'c' 与高斯窗口的 'beamwidth' 相关），您实际上可以改进一些估计，但它可能会导致一些结果由于数字敏感性而爆炸。高斯比率技术的有效“Q”或“增益”约为 5,000 – 10,000，如果您知道如何正确调整参数，它可能会更高。
当然，如果你的 SNR 太低，你可能不得不增加 N 以获得更多的积分时间，但受制于信号稳定性的限制。