FFTW 信号频谱答案

【问题标题】：FFTW signal spectrumFFTW 信号频谱
【发布时间】：2017-04-04 13:45:30
【问题描述】：

我目前正在实现信号频谱图。我的输入数据：1.051,1.365,1.837,2.334,2.486,2.68,2.878,2.579,2.11,1.605,0.936,0.73,1.036,1.336,1.796,2.043,2.462,2.935,2.892,2.609,2.151,1.641,0.961， 0.738,0.767,1.303,1.503,2.004,2.435,2.931,2.908,2.641,1.935,1.423,0.988,0.74,0.755,1.018,1.465,1.966,2.662,2.919,2.918,2.669,1.976,1.462,1.017,0.757， 0.746,0.989,1.426,1.927,2.637,2.908,2.927,2.441,2.015,1.501,1.302,1.025,0.739,0.962,1.644,2.144,205,2.895,2.935,205,2.895,2.935,2.467,2.304,1.797,1.331,1.039,0.732， 0.936,1.608,2.15,2.575,2.88,2.688,2.493,2.343,1.835,1.365,1.054,0.725,0.913,1.568,2.067,2.288,2.867,2.69,2.516,2.379,1.877,1.397,1.072,0.721,1.144， 1.277, 1.77

您可以在我的屏幕截图中看到这些数据。这是底部图表。它只包含 100 个点。并且不可能增加它的数量。

输入信号频率为 1000 Hz。并且可以更改。

我使用 FFTW 获取频谱。

输入信号频率为 1000 Hz。在顶部图表中，它仅显示约 200 Hz。这是主要问题。我想可能是我的代码错误或分数不够。

我的数据分析代码：

    QVectorDouble points(100);
    points = this->reader->ReadCOM(100);
    double timePassed = this->reader->timePassed;
    unsigned int n = points.count();
    double timeShift = timePassed / n;

    QVectorDouble signalX(n), signalY(n);

    for (unsigned int i = 0; i < n; i++)
    {
        signalX[i] = i*timeShift; // x goes from 0 to  timePassed to take points amount
        signalY[i] = points[i];
    }

    fftw::maxthreads = get_max_threads();
    unsigned int np = n / 2 + 1;
    size_t align = sizeof(Complex);

    array1<Complex> F(np, align);
    array1<double> f(n, align);

    rcfft1d Forward(n, f, F);

    for(unsigned int i = 0; i < n; i++) {
        f[i] = points[i];
    }

    Forward.fft(f, F);

    QVectorDouble spectrX(np), spectrY(np);
    for (int i = 0; i < np; i++)
    {
      spectrX[i] = i * 20; //multiply by 20 because np is (100/2 + 1) and chart maximum xOrigin is 1000
      spectrY[i] = abs(F[i]) / np;
    }

【问题讨论】：

“我做得对吗”范围很广 - 您是否有证据表明存在/不存在特定问题？
输入信号频率为 1000 Hz。在顶部图表上，它仅显示约 200 Hz。这是主要问题）我想可能是我的代码错误或分数不够。
最好更改 fft 图 x 轴上的单位。现在看起来好像您正在绘制 fft bin 编号，实际上您希望绘制该特定 bin 的频率。
@TarasNikulin 看来 OP 的运行长度为 1024 fft 并使用 bin 编号作为 x 轴
@AlexZywicki 1000 这是我自己设置的范围

标签： c++ signals fft fftw spectrum

【解决方案1】：

以下是一些不太正确的地方：

在时域图中，1000Hz 采样率下的 100 个点应该给你 0.1 秒的总持续时间。假设 timePassed 为 0.1 秒，用于 signalX 的公式是正确的。你应该看看为什么这个变量实际上是 0.4 秒左右。
由于您在 0.1 秒内有 8 个完整的周期，因此您应该获得 8*1/0.1 = 80Hz 的峰值频率。因此，显然频谱图的 x 轴缩放有问题。更具体地说，频率增量应设置为1000/np 或np == 10 为10Hz。这将生成一个高达 500Hz 的图，这是采样率的一半（符合 Nyquist 定理），其中峰值将出现在正确的 80Hz 频率处。

【讨论】：