将Kiss FFT应用于音频样本，为什么会得到NaN输出？

Question

标题已经解释了我的问题。

我要做的事情非常简单：

• 加载MP3音轨（通过libmpg123）
• 读取样本
• 在样本上应用Kiss FFT

到目前为止，我尝试过以下方法：

inline float scale(kiss_fft_scalar val)
{
    int g = 0;
    return val < 0 ? val*(1/32768.0f ) : val*(1/32767.0f);
}

void main()
{
    mpg123_handle *m = NULL;
    int  channels = 0, encoding = 0;
    long rate = 0;
    int err = MPG123_OK;

    err = mpg123_init();        
    m = mpg123_new(NULL, &err);
    mpg123_open(m, "L:\\audio-io\\audio-analysis\\samples\\zero.mp3");
    mpg123_getformat(m, &rate, &channels, &encoding);

    err = mpg123_format_none(m);
    err = mpg123_format(m, rate, channels, encoding);

    // Get 2048 samples
    const int TIME = 2048;

    // 16-bit integer encoded in bytes, hence x2 size
    unsigned char* buffer = new unsigned char[TIME*2];
    size_t done = 0;
    err = mpg123_read(m, buffer, TIME*2, &done);

    short* samples = new short[done/2];
    int index = 0;

    // Iterate 2 bytes at a time
    for (int i = 0; i < done; i += 2)
    {
        unsigned char first = buffer[i];
        unsigned char second = buffer[i + 1];
        samples[index++] = (first | (second << 8));
    }

    // Array to store the calculated data
    int speclen = TIME / 2 + 1;
    float* output = new float[speclen];

    kiss_fftr_cfg config;
    kiss_fft_cpx* spectrum;

    config = kiss_fftr_alloc(TIME, 0, NULL, NULL);
    spectrum = (kiss_fft_cpx*) malloc(sizeof(kiss_fft_cpx) * TIME);

    // Right here...
    kiss_fftr(config, (kiss_fft_scalar*) samples, spectrum);

    for (int i = 0; i < speclen; i++)
    {
        float re = scale(spectrum[i].r) * TIME;
        float im = scale(spectrum[i].i) * TIME;

        output[i] = sqrtf(re*re + im*im);
    }

    return;
}

问题发生在这一行：kiss_fftr(config, (kiss_fft_scalar*) samples, spectrum); 其中samples包含音频样本（16位），spectrum应该保存输出数据。

函数完成后，在调试器窗口中发生了以下情况。 有人能给我一个简单的例子，说明如何在音频（16位编码）样本上应用Kiss FFT函数吗？

- ains

---

KissFFT里面没有文档或者示例代码吗？ - Alexey Frunze

这不是我的问题，但它或许反映了你方面表现出的缺乏努力？ - Alexey Frunze

附上的示例代码展示了我尝试过的。我在谷歌上未能找到任何类似的情况。我相信代码基本上是正确的。我正在努力弄清楚为什么会得到 NaN 值。 - ains

你可以先尝试使用简单的信号进行测试：全零、全一、正弦波等等，以验证FFT本身是否正常工作。我无法相信没有足够的示例代码或文档可以解释KissFFT例程在这些最简单情况下应该如何使用。一旦您完成了这个步骤，就可以开始处理mp3数据。加油！ - Alexey Frunze

是的，我也尝试过了。（在示例代码中，zero.mp3是一个填充了零频率的文件）信不信由你，仍然出现了相同的结果（NaN）。 - ains

你好，我正在与同样的问题苦苦挣扎。我已经使用了全零和全一进行了测试，结果完美无误。但是当我尝试使用实时音频时，它给出了错误的输出... - Sunny Shah

Answer 1

你需要查找代码中的错误。我的测试代码似乎工作得很好。

使用浮点数进行复数正向FFT：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "kiss_fft.h"

#ifndef M_PI
#define M_PI 3.14159265358979324
#endif

#define N 16

void TestFft(const char* title, const kiss_fft_cpx in[N], kiss_fft_cpx out[N])
{
  kiss_fft_cfg cfg;

  printf("%s\n", title);

  if ((cfg = kiss_fft_alloc(N, 0/*is_inverse_fft*/, NULL, NULL)) != NULL)
  {
    size_t i;

    kiss_fft(cfg, in, out);
    free(cfg);

    for (i = 0; i < N; i++)
      printf(" in[%2zu] = %+f , %+f    "
             "out[%2zu] = %+f , %+f\n",
             i, in[i].r, in[i].i,
             i, out[i].r, out[i].i);
  }
  else
  {
    printf("not enough memory?\n");
    exit(-1);
  }
}

int main(void)
{
  kiss_fft_cpx in[N], out[N];
  size_t i;

  for (i = 0; i < N; i++)
    in[i].r = in[i].i = 0;
  TestFft("Zeroes (complex)", in, out);

  for (i = 0; i < N; i++)
    in[i].r = 1, in[i].i = 0;
  TestFft("Ones (complex)", in, out);

  for (i = 0; i < N; i++)
    in[i].r = sin(2 * M_PI * 4 * i / N), in[i].i = 0;
  TestFft("SineWave (complex)", in, out);

  return 0;
}

输出:

Zeroes (complex)
 in[ 0] = +0.000000 , +0.000000    out[ 0] = +0.000000 , +0.000000
 in[ 1] = +0.000000 , +0.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +0.000000 , +0.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = +0.000000 , +0.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +0.000000 , +0.000000    out[ 4] = +0.000000 , +0.000000
 in[ 5] = +0.000000 , +0.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +0.000000 , +0.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = +0.000000 , +0.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +0.000000 , +0.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +0.000000 , +0.000000    out[ 9] = +0.000000 , +0.000000
 in[10] = +0.000000 , +0.000000    out[10] = +0.000000 , +0.000000
 in[11] = +0.000000 , +0.000000    out[11] = +0.000000 , +0.000000
 in[12] = +0.000000 , +0.000000    out[12] = +0.000000 , +0.000000
 in[13] = +0.000000 , +0.000000    out[13] = +0.000000 , +0.000000
 in[14] = +0.000000 , +0.000000    out[14] = +0.000000 , +0.000000
 in[15] = +0.000000 , +0.000000    out[15] = +0.000000 , +0.000000
Ones (complex)
 in[ 0] = +1.000000 , +0.000000    out[ 0] = +16.000000 , +0.000000
 in[ 1] = +1.000000 , +0.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +1.000000 , +0.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = +1.000000 , +0.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +1.000000 , +0.000000    out[ 4] = +0.000000 , +0.000000
 in[ 5] = +1.000000 , +0.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +1.000000 , +0.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = +1.000000 , +0.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +1.000000 , +0.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +1.000000 , +0.000000    out[ 9] = +0.000000 , +0.000000
 in[10] = +1.000000 , +0.000000    out[10] = +0.000000 , +0.000000
 in[11] = +1.000000 , +0.000000    out[11] = +0.000000 , +0.000000
 in[12] = +1.000000 , +0.000000    out[12] = +0.000000 , +0.000000
 in[13] = +1.000000 , +0.000000    out[13] = +0.000000 , +0.000000
 in[14] = +1.000000 , +0.000000    out[14] = +0.000000 , +0.000000
 in[15] = +1.000000 , +0.000000    out[15] = +0.000000 , +0.000000
SineWave (complex)
 in[ 0] = +0.000000 , +0.000000    out[ 0] = +0.000000 , +0.000000
 in[ 1] = +1.000000 , +0.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +0.000000 , +0.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = -1.000000 , +0.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +0.000000 , +0.000000    out[ 4] = +0.000000 , -8.000000
 in[ 5] = +1.000000 , +0.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +0.000000 , +0.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = -1.000000 , +0.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +0.000000 , +0.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +1.000000 , +0.000000    out[ 9] = +0.000000 , +0.000000
 in[10] = +0.000000 , +0.000000    out[10] = +0.000000 , +0.000000
 in[11] = -1.000000 , +0.000000    out[11] = +0.000000 , +0.000000
 in[12] = +0.000000 , +0.000000    out[12] = +0.000000 , +8.000000
 in[13] = +1.000000 , +0.000000    out[13] = +0.000000 , +0.000000
 in[14] = +0.000000 , +0.000000    out[14] = +0.000000 , +0.000000
 in[15] = -1.000000 , +0.000000    out[15] = +0.000000 , +0.000000

使用浮点数进行实值正向FFT：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "kiss_fftr.h"

#ifndef M_PI
#define M_PI 3.14159265358979324
#endif

#define N 16

void TestFftReal(const char* title, const kiss_fft_scalar in[N], kiss_fft_cpx out[N / 2 + 1])
{
  kiss_fftr_cfg cfg;

  printf("%s\n", title);

  if ((cfg = kiss_fftr_alloc(N, 0/*is_inverse_fft*/, NULL, NULL)) != NULL)
  {
    size_t i;

    kiss_fftr(cfg, in, out);
    free(cfg);

    for (i = 0; i < N; i++)
    {
      printf(" in[%2zu] = %+f    ",
             i, in[i]);
      if (i < N / 2 + 1)
        printf("out[%2zu] = %+f , %+f",
               i, out[i].r, out[i].i);
      printf("\n");
    }
  }
  else
  {
    printf("not enough memory?\n");
    exit(-1);
  }
}

int main(void)
{
  kiss_fft_scalar in[N];
  kiss_fft_cpx out[N / 2 + 1];
  size_t i;

  for (i = 0; i < N; i++)
    in[i] = 0;
  TestFftReal("Zeroes (real)", in, out);

  for (i = 0; i < N; i++)
    in[i] = 1;
  TestFftReal("Ones (real)", in, out);

  for (i = 0; i < N; i++)
    in[i] = sin(2 * M_PI * 4 * i / N);
  TestFftReal("SineWave (real)", in, out);

  return 0;
}

输出：

Zeroes (real)
 in[ 0] = +0.000000    out[ 0] = +0.000000 , +0.000000
 in[ 1] = +0.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +0.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = +0.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +0.000000    out[ 4] = +0.000000 , +0.000000
 in[ 5] = +0.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +0.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = +0.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +0.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +0.000000    
 in[10] = +0.000000    
 in[11] = +0.000000    
 in[12] = +0.000000    
 in[13] = +0.000000    
 in[14] = +0.000000    
 in[15] = +0.000000    
Ones (real)
 in[ 0] = +1.000000    out[ 0] = +16.000000 , +0.000000
 in[ 1] = +1.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +1.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = +1.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +1.000000    out[ 4] = +0.000000 , +0.000000
 in[ 5] = +1.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +1.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = +1.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +1.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +1.000000    
 in[10] = +1.000000    
 in[11] = +1.000000    
 in[12] = +1.000000    
 in[13] = +1.000000    
 in[14] = +1.000000    
 in[15] = +1.000000    
SineWave (real)
 in[ 0] = +0.000000    out[ 0] = +0.000000 , +0.000000
 in[ 1] = +1.000000    out[ 1] = +0.000000 , +0.000000
 in[ 2] = +0.000000    out[ 2] = +0.000000 , +0.000000
 in[ 3] = -1.000000    out[ 3] = +0.000000 , +0.000000
 in[ 4] = +0.000000    out[ 4] = +0.000000 , -8.000000
 in[ 5] = +1.000000    out[ 5] = +0.000000 , +0.000000
 in[ 6] = +0.000000    out[ 6] = +0.000000 , +0.000000
 in[ 7] = -1.000000    out[ 7] = +0.000000 , +0.000000
 in[ 8] = +0.000000    out[ 8] = +0.000000 , +0.000000
 in[ 9] = +1.000000    
 in[10] = +0.000000    
 in[11] = -1.000000    
 in[12] = +0.000000    
 in[13] = +1.000000    
 in[14] = +0.000000    
 in[15] = -1.000000    

- Alexey Frunze

---

你好，我正在与同样的问题苦苦挣扎。我已经使用了全零和全一进行了测试，结果完美无误。但是当我尝试使用实时音频时，它给出了错误的输出... - Sunny Shah

@Sunnyshah 你要么输入了错误的数据，要么期望了错误的结果。只要 InverseFFT(FFT(signal))=signal，问题肯定出在你的代码中而不是 FFT 中。但如果 InverseFFT(FFT(signal))≠signal，则仍然很可能是你没有正确使用(I)FFT例程。恐怕你需要调试你的代码。 - Alexey Frunze

谢谢，我解决了我的问题，是因为我输入的数据有误。但现在我又遇到了另一个问题，请看一下我的代码峰值问题。 - Sunny Shah

Answer 2

当我开始查看这个答案时，我一直在想为什么-8.0出现在虚部而不是实部。在重新阅读FFT的印刷文章时，我意识到我一直在思考幅度问题。

因此，我在复数代码中对printf进行了调整，如下所示。

for (i = 0; i < N; i++)
    printf(" in[%02i]=%+f, %+f  out[%02i]=%+f, %+f M[%02i]=%+f\n",
         i, in[i].r, in[i].i,
         i, out[i].r, out[i].i,
         i, sqrt((out[i].r * out[i].r) + (out[i].i * out[i].i)));

这将显示出大小的答案。

...
SineWave (complex)
 in[00]=+0.000000, +0.000000  out[00]=+0.000000, +0.000000 M[00]=+0.000000
 in[01]=+1.000000, +0.000000  out[01]=+0.000000, +0.000000 M[01]=+0.000000
 in[02]=+0.000000, +0.000000  out[02]=+0.000000, +0.000000 M[02]=+0.000000
 in[03]=-1.000000, +0.000000  out[03]=+0.000000, +0.000000 M[03]=+0.000000
 in[04]=-0.000000, +0.000000  out[04]=-0.000000, -8.000000 M[04]=+8.000000
 in[05]=+1.000000, +0.000000  out[05]=+0.000000, -0.000000 M[05]=+0.000000
 in[06]=+0.000000, +0.000000  out[06]=+0.000000, -0.000000 M[06]=+0.000000
 in[07]=-1.000000, +0.000000  out[07]=+0.000000, -0.000000 M[07]=+0.000000
 in[08]=-0.000000, +0.000000  out[08]=+0.000000, +0.000000 M[08]=+0.000000
 in[09]=+1.000000, +0.000000  out[09]=+0.000000, +0.000000 M[09]=+0.000000
 in[10]=+0.000000, +0.000000  out[10]=+0.000000, +0.000000 M[10]=+0.000000
 in[11]=-1.000000, +0.000000  out[11]=+0.000000, +0.000000 M[11]=+0.000000
 in[12]=-0.000000, +0.000000  out[12]=-0.000000, +8.000000 M[12]=+8.000000
 in[13]=+1.000000, +0.000000  out[13]=+0.000000, -0.000000 M[13]=+0.000000
 in[14]=+0.000000, +0.000000  out[14]=+0.000000, -0.000000 M[14]=+0.000000
 in[15]=-1.000000, +0.000000  out[15]=+0.000000, -0.000000 M[15]=+0.000000

我也尝试过在生成正弦波的for循环中改变频率。

float freq;
...
freq = 6.0;
for (i = 0; i < N; i++)
    in[i].r = sin(2 * M_PI * freq * i / N), in[i].i = 0;

只要我保持在1.0的倍数及其以下，并且在奈奎斯特频率16/2 = 8之内，结果就会很好地从一个bin转移到另一个bin。当然，将频率设置为分数值会使其幅度在bins上分布，并且如果不应用窗函数，则会出现泄漏。如果你仍然像我一样对FFT感到困惑，请尝试使用这样的代码玩耍一段时间，可以在单个屏幕上看到所有结果，事情开始变得更加清晰。

最后，感谢Alexey提供的答案，它帮助我开始使用Kiss FFT。

- TJA

Answer 3

试试这个：

in[i].r = sin(2 * M_PI * freq * (i / N*1.00)), in[i].i = 0;

- Michal Rudnicki

---

如果您能解释一下为什么这个解决方案有效，那将非常有帮助。 - Robert Columbia

变量“i”和“N”都是整数，结果将是相同的。如果您更改代码以添加数组中的乘法数据，则应该是正确的。 - Michal Rudnicki