保存在文件中的双样本的 FFTW

Posted 2023-02-22

【中文标题】保存在文件中的双样本的 FFTW

【英文标题】：FFTW of double samples saved in file

【发布时间】：2016-04-06 11:45:31

【问题描述】：

我正在尝试使用 FFTW 库计算 53k 双样本的 FFT，并在此基础上猜测信号的基频是多少。样本由 sndfile 库在 wav 输入文件的基础上生成（程序加载到 wav 文件中，生成双数据样本并保存到文本文件中）。每个样本的数量级范围从 -0.0009 到 +0.0009。使用下面介绍的函数计算 FFT 后，我收到 in[0][0] = -2.142。这是检索数据的不正确方法，或者输入文件不正确。我究竟做错了什么？是传递给 FFTW 函数的数据错误，是文件中存储的数据不正确还是我编写的函数不正确？ Buffer 是一个包含样本的浮点数组。

static fftw_complex* calculateFourier(float* buffer, const unsigned int bufferLen)


    unsigned int i;
    //const unsigned short windowSize = 1024;
    fftw_plan result;
    fftw_complex *out, *in;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * bufferLen);
    out = (fftw_complex*) fftw_malloc( sizeof(fftw_complex) * bufferLen);

    double* doubleBuffer = castToDouble(buffer, bufferLen);

    for(i = 0; i < bufferLen; i++)
    
        in[i][0] = doubleBuffer[i];
        in[i][1] = 0.000000000000;
    

    result = fftw_plan_dft_1d(bufferLen, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    //result = fftw_plan_dft_r2c_1d(fileSize, castToDouble(buffer, bufferLen), out, FFTW_MEASURE);

    fftw_execute(result);

    fftw_destroy_plan(result);

    return out;

static double* castToDouble(float* buffer, const unsigned int bufferLen)

    unsigned int i;
    double* doubleBuffer = (double*)malloc( sizeof(double) * bufferLen); 

    for(i = 0; i < bufferLen; i++)
        doubleBuffer[i] = (double)buffer[i];

    return doubleBuffer;

【问题讨论】：

castToDouble 是做什么的？

由于 sndfile 库将样本保存为浮点数，而 fftw 需要 double 类型的样本，它只是将每个样本强制转换为 double。

fftw 的 fft 未标准化。参见fftw.org/doc/…：如果缓冲区长度为 53000，则 out[0][0]（“零频率”）是缓冲区平均值的 53000 倍，即缓冲区所有项的总和。

【参考方案1】：

我认为问题在于您没有将window function 应用于您的原始数据。因此，傅里叶变换具有与样本缓冲区长度相对应的强频率分量。

尝试通过替换此行来应用汉宁窗：

in[i][0] = doubleBuffer[i];

用这个：

in[i][0] = doubleBuffer[i] * 0.5 * (1.0 - cos(2.0 * M_PI * i / (bufferLen - 1)));

（如果需要重复使用，显然这个窗口函数可以缓存。）

【讨论】：

在这种情况下，输出数组 (out[i][0]) 中的每个实分量都等于 0

汉宁窗在cos参数中缺少i因子，所以它应该是0.5 * (1.0 - cos(2.0 * i * M_PI / (bufferLen - 1)))

@SleuthEye 感谢您发现这一点

@Kokos34 我错过了汉宁窗公式中的一个因素。现在再试一次。

@squeamishossifrage 它仍然没有给我我期望的答案。 out[0][0] 不是信号的基频吗？答案是 -1.87，在我将幅度计算为 sqrt(re^2 + im^2) 之后，它给出的结果是 1.87。我还需要做什么才能得到信号基频的结果？

【参考方案2】：

仔细看看下面的代码有什么问题：

#include <stdio.h>

void prn (float *a, int n);

int main (void)

    float floatarray[] =  0.0, 5.0, 10.0, 15.0, 20.0, 25.0 ;

    prn (floatarray, sizeof floatarray/sizeof *floatarray);

    return 0;


void prn (float *a, int n)

    int i;
    double *d = (double *)a;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);

运行它，会发生什么？为什么？

sizeof 和 float 是什么？ sizeof 和 double 是什么？当演员double *d = (double *)a; 发生时会发生什么？ d 的每个元素需要多少字节？那里有多少？你的

double* doubleBuffer = castToDouble(buffer, bufferLen);

导致完全相同的问题。

需要逐个元素的演员/分配

无论您使用分配有malloc（或calloc）的新double 数组，您都必须执行逐个元素 强制转换/赋值来完成数组的强制转换.下面是相同的代码，但添加了一些使用可变长度数组以避免动态分配内存：

#include <stdio.h>

void castfloatdouble (double *d, float *f, int n);
void prn (float *a, int n);
void prndouble (double *d, int n);

int main (void)

    float floatarray[] =  0.0, 5.0, 10.0, 15.0, 20.0, 25.0 ;
    int n = (int)(sizeof floatarray/sizeof *floatarray);
    double doublearray[n];  /* memset is a good idea to zero a VLA */

    prn (floatarray, n);  /* incorrect results due to typesize mismatch */
    putchar ('\n');

    castfloatdouble (doublearray, floatarray, n);

    prndouble (doublearray, n);  /* correct results after deep cast/assign */
    putchar ('\n');

    return 0;


void castfloatdouble (double *d, float *f, int n)

    int i;

    for (i = 0; i < n; i++)
        d[i] = (double)f[i];


void prn (float *a, int n)

    int i;
    double *d = (double *)a;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);


void prndouble (double *d, int n)

    int i;

    for (i = 0; i < n; i++)
        printf (" a[%d] : %01.0lf\n", i, d[i]);

输出显示不正确和正确的结果

$ ./bin/castfloatdouble2
 a[0] : 2048
 a[1] : 16777220
 a[2] : 805306499
 a[3] : 0
 a[4] : 0
 a[5] : 0

 a[0] : 0
 a[1] : 5
 a[2] : 10
 a[3] : 15
 a[4] : 20
 a[5] : 25

【讨论】：

但我的 castToDouble 分配了适当大小的内存，我猜：double* doubleBuffer = (double*)malloc( sizeof(double) * bufferLen);

在没有看到代码的情况下很难 猜测 它在做什么，但即使你分配了内存，如果你没有进行逐个元素的转换，并且将每个buffer[i] 值分配到新的double 存储中，那么将存在相同的问题。您必须注意如何将 4 字节浮点数组转换为 8 字节双精度数组。记住array 意味着顺序存储在内存中。

我编辑了帖子以包含 castToDouble 函数实现。这是正确的吗？你想让我发布一些额外的代码吗？

是的，应该可以。我将放弃另一种使用静态数组的方式来编辑我的答案。