Kiss_fftr 的 KissFFT 输出答案

【问题标题】：KissFFT output of kiss_fftrKiss_fftr 的 KissFFT 输出
【发布时间】：2012-06-05 13:33:01
【问题描述】：

我正在通过套接字连接接收包含 320 个样本的数据包中的 PCM 数据。声音的采样率为每秒 8000 个样本。我正在做这样的事情：

int size = 160 * 2;//160;
int isinverse = 1;
kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_in[size];
kiss_fft_cpx fft_out[size];
kiss_fft_cpx fft_reconstructed[size];

kiss_fftr_cfg fft = kiss_fftr_alloc(size*2 ,0 ,0,0);
kiss_fftr_cfg ifft = kiss_fftr_alloc(size*2,isinverse,0,0);

for (int i = 0; i < size; i++) {
    fft_in[i].r = zero;
    fft_in[i].i = zero;
    fft_out[i].r = zero;
    fft_out[i].i = zero;
    fft_reconstructed[i].r = zero;
    fft_reconstructed[i].i = zero;
}

// got my data through socket connection

for (int i = 0; i < size; i++) {
     // samples are type of short
     fft_in[i].r = samples[i];
     fft_in[i].i = zero;
     fft_out[i].r = zero;
     fft_out[i].i = zero;
 }

 kiss_fftr(fft, (kiss_fft_scalar*) fft_in, fft_out);
 kiss_fftri(ifft, fft_out, (kiss_fft_scalar*)fft_reconstructed);

 // lets normalize samples
 for (int i = 0; i < size; i++) {
     short* samples = (short*) bufTmp1;
     samples[i] = rint(fft_reconstructed[i].r/(size*2));
 }

之后，我填充 OpenAL 缓冲区并播放它们。一切正常，但我想对kiss_fftr 和kiss_fftri 之间的音频进行一些过滤。我认为的起点是将声音从时域转换为频域，但我真的不明白我从kiss_fftr 函数接收到什么样的数据。每个复数中存储了哪些信息，它的实部和虚部可以告诉我有关频率的信息。而且我不知道fft_out 中覆盖了哪些频率（什么频率跨度） - 哪些索引对应于哪些频率。

我是信号处理和傅里叶变换主题的新手。

有什么帮助吗？

【问题讨论】：

这是题外话；可能更适合例如dsp.stackexchange.com。但实际上，您可能应该为自己准备一本不错的 DSP 入门书籍（例如 this），因为这个主题太大，无法在 StackOverflow 风格的答案中充分涵盖。

标签： c audio filter fft kissfft

【解决方案1】：

在双脚投入 C 实现之前，请熟悉数字滤波器，尤其是 FIR filters。

您可以使用GNU Octave's 信号工具箱之类的工具设计 FIR 滤波器。查看命令 fir1（最简单的）、firls 或 remez。或者，您可以通过网页设计 FIR 滤波器。快速网络搜索“online fir filter design”找到this（我没用过，但它似乎使用了remez或firpm命令中使用的等波纹设计）

尝试首先使用直接卷积（不带 FFT）来实现您的过滤器，看看速度是否可以接受 - 这是一条更简单的路径。如果您需要基于 FFT 的方法，kissfft/tools/kiss_fastfir.c 文件中有一个重叠保存的示例实现。

【讨论】：

【解决方案2】：

我将尝试直接回答您的问题。

// a) the real and imaginary components of the output need to be combined to calculate the amplitude at each frequency. 

float ar,ai,scaling; 

scaling=1.0/(float)size;

// then for each output [i] from the FFT...

ar = fft_out[i].r; 
ai = fft_out[i].i; 
amplitude[i] = 2.0 * sqrtf( ar*ar + ai*ai ) * scaling ;

// b) which index refers to which frequency? This can be calculated as follows. Only the first half of the FFT results are needed (assuming your 8KHz sampling rate) 

for(i=1;i<(size/2);i++) freq = (float)i / (1/8000) / (float)size ; 

// c)  phase (range +/- PI) for each frequency is calculated like this: 

phase[i] = phase = atan2(fft_out[i].i / fft_out[i].r);

【讨论】：

很好的解释，我有 2 个问题。 1）你能解释更多关于阶段是什么吗？在 b) 如果 i=1 给出 freq = 25; i=2 给出 50。 Bin1 是 [0,25]hz 和 bin2 是 [26,50]hz 吗？灵态是在那儿进来的吗？ 2）频率的另一半会发生什么？如果我们以 8Khz 采样，在 i=max(160) 时我们只有 4000hz。 4400hz+呢？我们应该在垃圾箱中生成 8000hz 的频率吗？
"奈奎斯特-香农采样定理说，当采样率是被采样信号最大频率的两倍以上时，可以重建信号。这条规则告诉我们，要再现高达人类听到 20,000 赫兹，我们必须每秒至少采集 40,001 个样本”daugroupassignment.wikia.com/wiki/…

【解决方案3】：

您可能想要研究的是使用重叠添加或重叠保存算法的 FFT 快速卷积。您需要将每个 FFT 的长度扩展为所需滤波器的脉冲长度。这是因为 (1) FFT/IFFT 卷积是圆形的，并且 (2) FFT 阵列结果中的每个索引对应于几乎所有频率（Sinc 形响应），而不仅仅是一个（即使大部分接近一个），所以任何一个bin 修改将在整个频率响应中泄漏（某些精确的周期性频率除外）。

【讨论】：