使用 FFTW 计算波长、幅度和相位

Question

我正在尝试使用 C++ 和 FFTW 从一系列数字中获取幅度、波长和相位，但运气不佳，我想知道是否有人可以帮助我？

我已经在互联网上进行了搜索，但由于我是 FFT 新手，所以我无法找到任何我能理解的足够具体的内容

我有 370 个双精度值，正在从文件中读取到 inputRecord 数组中，例如1.11567, 1.11599, 1.11679

这就是我迄今为止在 FFTW 上所做的事情：

double *input;
fftw_complex *output;
fftw_plan p;
int binCount = inputCount/ 2; // inputCount = the number of records in the input file

input = (double*)fftw_malloc(sizeof(double) * inputCount);
output = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * inputCount);

p = fftw_plan_dft_r2c_1d(barCount, input, output, FFTW_ESTIMATE);

// Copy the double numbers into the FFTW input array 
for (int i = 0; i < barCount; i++)
{
    input[i] = inputRecord[i];
}

fftw_execute(p); //performs the fourier transform, fills `output'

然后我循环遍历 FFTW 输出数组，尝试通过此计算获取幅度：

Amplitude = 2 * abs(FFTOutput[i]) / 370

我还想获得波长（不知道如何）和正弦波的相位（也不知道如何做到这一点）

非常感谢任何帮助。

Answer 1

有同样的问题。我希望能够创建一组样本，然后让 FFTW 返回我用于创建该组样本的参数。

这段代码创建了示例。您可以给它频率、偏置、幅度和相位。

    int const N = 32;
    fftw_complex in[N + 1];

    double const freq = 3; // Hz
    double const bias = 1.5;
    double const magnitude = 1.0;
    // double const phase = 0;
    double const phase = M_PI / 2.0;

    double const omega = (2 * M_PI * freq);
    double const increment = omega / N;
    double sample = 0.0;
    for (int time = 0; time < (N + 1); ++time) {
        in[time][0] = bias + (magnitude * ::cos(sample + phase));
        in[time][1] = 0.0;
        printf("%-10s: %3d %+9.5f\n", "input", time, in[time][0]);
        sample += increment;
    }

该块使用 FFTW 来分析并报告这些相同的值：

#include <fftw3.h>


    /* forward Fourier transform, save the result in 'out' */
    fftw_complex out[N];
    fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    printf("\n");
    printf("Expected phase: %9.5f\n", phase);
    printf("%-10s: %3s %9s %9s %9s %9s\n", "FFT", "Hz", "out[0]", "out[1]", "mag", "phase");
    for (int i = 0; i < (N + 1); i++) {
        double const rpt_real = out[i][0];
        double const rpt_imag = out[i][1];
        if (i == 0) {
            printf("%-10s: %3d %+9.5f %9s %9s %9s\n", "bias", i, rpt_real / N, "n/a", "n/a", "n/a");
        } else {
            // phase - 2 * PI => reported phase
            double const calc_phase = std::atan2(rpt_imag, rpt_real);
            double const calc_magnitude = std::sqrt(std::pow(rpt_real, 2) + std::pow(rpt_imag, 2)) * (2.0 / N);
            printf("%-10s: %3d %+9.5f %+9.5f %9.5f %9.5f\n", "freq.", i, rpt_real, rpt_imag, calc_magnitude, calc_phase);
        }
    }
    fftw_destroy_plan(p);

这是示例输出：

FFT       :  Hz    out[0]    out[1]       mag     phase
bias      :   0  +1.50000       n/a       n/a       n/a
freq.     :   1  -0.00000  -0.00000   0.00000  -2.54329
freq.     :   2  -0.00000  -0.00000   0.00000  -2.16145
freq.     :   3  +0.00000 +16.00000   1.00000   1.57080
freq.     :   4  -0.00000  +0.00000   0.00000   2.27161

在 out[0]: 中找到偏差，并显示 +1.5 与

double const bias = 1.5;

匹配

在 out[3] 中找到与

double const freq = 3;

匹配的频率

震级可在“mag”列下找到。 1.00000 匹配

double const magnitude = 1.0;

该相可在“相”列 1.57080 中找到，它与

double const phase = M_PI / 2.0;

相匹配

我尝试过各种频率、偏差、幅度和相位。据我所知，它每次都给出正确的值。

（使用一组 GTest UT 来自动确认所有这些。）