我已经使用kiss fft库将fft实现到32ub系列ucontroller中,目前正在努力处理fft的输出。我的目的是分析来自压电扬声器的声音。目前,音响的频率是420Hz,我成功地从fft输出(与示波器交叉检查)。然而,如果我将函数生成器波形放入系统中,输出频率只是预期的一半。我怀疑是频率仓的计算公式搞错了;当前使用fft_peak_magnitude_index*采样频率/fft_size。我的输入是真实的,做的是真实的。(输出样本= N/2)在fft之前也要做iir滤波和加窗。任何建议都会有很大的帮助!
// IIR filter calculation, n = 256 fft points
for (ctr=0; ctr<n; ctr++)
{
// filter calculation
y[ctr] = num_coef[0]*x[ctr];
y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
y1[ctr] = y[ctr] - 510; //eliminate dc offset
// hamming window
hamming[ctr] = (0.54-((0.46) * cos(2*M_PI*ctr/n)));
window[ctr] = hamming[ctr]*y1[ctr];
fft_input[ctr].r = window[ctr];
fft_input[ctr].i = 0;
fft_output[ctr].r = 0;
fft_output[ctr].i = 0;
}
kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);
peak = 0;
freq_bin = 0;
for (ctr=0; ctr<n1; ctr++)
{
fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);
if(fft_mag[ctr] > peak)
{
peak = fft_mag[ctr];
freq_bin = ctr;
}
frequency = (freq_bin*(10989/n)); // 10989 is the sampling freq
//************************************
//Usart write
char filtResult[10];
//sprintf(filtResult, "%04d %04d %04dn", (int)peak, (int)freq_bin, (int)frequency);
sprintf(filtResult, "%04d %04d %04dn", (int)x[ctr], (int)fft_mag[ctr], (int)frequency);
char c;
char *ptr = &filtResult[0];
do
{
c = *ptr;
ptr++;
usart_bw_write_char(&AVR32_USART2, (int)c);
// sendByte(c);
} while (c != 'n');
}
主要问题可能是您如何声明fft_input
。根据您之前的问题,您将fft_input
分配为kiss_fft_cpx
数组。另一方面,函数kiss_fftr
期望一个标量数组。通过使用:
kiss_fft_scalar
kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);
KissFFT本质上看到的是一个实值数据数组,其中包含每秒钟采样的零(您填充的虚部)。这实际上是原始信号的上采样版本(尽管没有插值),即具有有效两倍采样率的信号(这在freq_bin
到frequency
转换中没有考虑)。为了解决这个问题,我建议您将数据打包到kiss_fft_scalar
数组中:
kiss_fft_scalar fft_input[n];
...
for (ctr=0; ctr<n; ctr++)
{
...
fft_input[ctr] = window[ctr];
...
}
kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, fft_input, fft_output);
还请注意,在寻找峰值大小时,您可能只对最终最大的峰值感兴趣,而不是运行最大值。因此,您可以将循环限制为只计算峰值(如果需要,在以下sprintf
语句中使用freq_bin
而不是ctr
作为数组索引):
for (ctr=0; ctr<n1; ctr++)
{
fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);
if(fft_mag[ctr] > peak)
{
peak = fft_mag[ctr];
freq_bin = ctr;
}
} // close the loop here before computing "frequency"
最后,当计算与最大幅度的bin相关的频率时,您需要确保使用浮点运算完成计算。如果我怀疑n
是一个整数,那么您的公式将使用整数算术执行10989/n
因子,从而导致截断。这可以简单地使用:
frequency = (freq_bin*(10989.0/n)); // 10989 is the sampling freq