使用Eigen :: FFT执行FFT时的频率

Question

bounty到期在6小时内。该问题的答案有资格获得+ 50声望赏金。Sito正在寻找规范答案：

由于Eigen库在科学计算中得到了广泛使用，因此FFT的简单示例及其使用方法可能对广大受众有益。如果提供的示例在某种程度上不足，希望有人可以提供一些说明性的示例。

我目前正在试图弄清楚如何正确使用Eigen的FFT算法。让我们假设我有一个功能

std::complex<double> f(std::complex<double> const & t){
    return std::sin(t);
}

然后我使用此函数进行计算

Eigen::VectorXcd time(1000);
Eigen::VectorXcd f_values(1000);
for(int u = 0; u < 1000; ++u){
    time(u) = u* 2. * M_PI / 1000;
    f_values(u) = f(time(u));
}

我现在想计算f_values的傅立叶变换，所以我这样做

Eigen::FFT<double> fft;
Eigen::VectorXcd f_freq(1000);
fft.fwd(f_freq, f_values);

现在，我想对此作图，但是要做到这一点，我需要评估f_freq的频率，但是我真的不知道如何获得这些频率。所以我的问题归结为找到Eigen::VectorXcd，其中包含绘制此类内容的频率（对不起，我使用图片作为描述，但是我认为这样更清晰了，如果我尝试用文字来描述它的话...情节中的amplitude应该对应于我的f_freq和我所看到的for是图片中freq的值...）。

以下是上面的代码段，放入单个文件中：

#include <eigen3/Eigen/Dense>
#include <eigen3/unsupported/Eigen/FFT>
#include <complex>
#include <cmath>

std::complex<double> f(std::complex<double> const & t){
     return std::sin(t);
}

int main(){
    Eigen::VectorXcd time(1000);
    Eigen::VectorXcd f_values(1000);
    for(int u = 0; u < 1000; ++u){
        time(u) = u* 2. * M_PI / 1000;
        f_values(u) = f(time(u));
    }

    Eigen::FFT<double> fft;
    Eigen::VectorXcd f_freq(1000);
    fft.fwd(f_freq, f_values);
    //freq = ....
}

---

我实现了以下建议的答案之一：

#include <eigen3/Eigen/Dense>
#include <eigen3/unsupported/Eigen/FFT>
#include <complex>
#include <cmath>
#include <iostream>
#include <fstream>

std::complex<double> f(std::complex<double> const & t){
     return std::sin(1.*t);
}

int main(){
    std::ofstream freq_out("frequencies.txt");
    std::ofstream f_freq_out("f_freq.txt");

    unsigned const N = 1000.;
    Eigen::VectorXcd time(N);
    Eigen::VectorXcd f_values(N);
    for(int u = 0; u < N; ++u){
        time(u) = u* 2. * M_PI / double(N);
        f_values(u) = f(time(u));
    }

    Eigen::FFT<double> fft;
    Eigen::VectorXcd f_freq(N);
    Eigen::VectorXd freq(N);
    fft.fwd(f_freq, f_values);

    double const Ts = 2. * M_PI/double(N);
    double const Fs = 1./Ts;

    freq = time.real() * Fs/2.;

    freq_out << freq; 
    f_freq_out << f_freq.cwiseAbs();
}

这将导致下图似乎有点不对。缩放比例当然没有多大意义，但事实上，存在两个峰值会导致我有点怀疑。.]

我目前正在试图弄清楚如何精确使用本征的FFT算法。让我们假设我有一个函数std :: complex f（std :: complex const＆t）{return std :: ...

Answer 1

通常，库使用以下公式计算DFT：

X [k] = sum_n（x [n] * exp（-2 * pi * i * k * n / N）] >>

where

• X-是傅立叶域数组。值表示相应的正弦/余弦函数的振幅。
• k-是傅立叶域中的索引，它也唯一地定义了频率
• i-只是数学i-复数（0 + 1i）
• N-是数组的大小

因此，在索引k，您的频率的长度是整个信号输入的1/k。特别是：

• X[0]是您的平均值
• X[1]对应于在整个域中恰好适合一次的正弦/余弦函数
• X[2]对应于适合您的域两次的正弦/余弦函数...依此类推...

在索引k> N / 2处，频率是如此之高，以至于由于混叠，它实际上对应于较低的频率。

这里是N = 8的示例：

我没有特别询问Eigen，但我认为没有什么不同。

Answer 2

在数字域上，频率始终跨越2*pi [radians]（根据需要，它可以从-pi到pi或从0到2*pi）。因此，您需要将范围[0,2*pi)划分为N bin。例如，如果index为k，则关联的归一化频率为f=2*pi*k/N [radians]。