真正的 FFT 输出
Real FFT output
我已经使用 kiss fft 库将 fft 实现到 at32ucb 系列 ucontroller 中,目前正在努力处理 fft 的输出。
我的目的是分析来自压电扬声器的声音。
目前,发声器的频率是 420Hz,这是我从 fft 输出成功获得的(用示波器交叉检查)。但是,如果我将函数发生器波形放入系统中,输出频率只是预期的一半。
我怀疑是我弄错的频率仓计算公式;当前使用,fft_peak_magnitude_index*采样频率/fft_size。
我的输入是真实的并且是真实的 fft。 (输出样本 = 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 %04d\n", (int)peak, (int)freq_bin, (int)frequency);
sprintf(filtResult, "%04d %04d %04d\n", (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
with:
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
我已经使用 kiss fft 库将 fft 实现到 at32ucb 系列 ucontroller 中,目前正在努力处理 fft 的输出。 我的目的是分析来自压电扬声器的声音。 目前,发声器的频率是 420Hz,这是我从 fft 输出成功获得的(用示波器交叉检查)。但是,如果我将函数发生器波形放入系统中,输出频率只是预期的一半。 我怀疑是我弄错的频率仓计算公式;当前使用,fft_peak_magnitude_index*采样频率/fft_size。 我的输入是真实的并且是真实的 fft。 (输出样本 = 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 %04d\n", (int)peak, (int)freq_bin, (int)frequency);
sprintf(filtResult, "%04d %04d %04d\n", (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
with:
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