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

Question

我目前正在尝试弄清楚如何使用库的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，而我正在寻找的是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;

    for(int u = 0; u < N; ++u){
        freq(u) = Fs * u / double(N);
    }

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

这导致了以下绘图结果： 这看起来有点不对劲..缩放显然没有太多意义，而且存在两个值尖峰的事实让我有点怀疑..

- Sito

---

应该是类似于 k * SampleRate/N 的式子，其中 N 是点数， -N/2 <=k < N/2。需要检查 N 是奇数还是偶数。此外，根据实现方式，频率可以存储在 -fmax 到 fmax 或 0 到 fmax 然后是 -fmax 到 0，这意味着绘图前需要重新组织数据。我知道这是一个 C++ 的问题，但是你可以查看 python 的文档 numpy.fft 和 numpy.fft.fftfreq 的实现。 - mikuszefski

这个有帮助吗？https://stackoverflow.com/questions/49676830/eigen-fft-library - Jerry Jeremiah

总的来说，你的FFT只知道信号的点数。因此，答案以总时间间隔为单位。如果你想将其映射到实际值，你需要知道两个点之间的时间。这通常以“每秒点数”的形式表示，即采样率。其倒数是两个点之间的时间，这就是numpy实现中使用的内容。 - mikuszefski

@Sito：请注意，与Ts相比，您的正弦波频率非常低。也许尝试另一个Ts以使尖峰更靠近中心... Ts = 10 /（2 * pi） 对我来说似乎是合理的。按照您现在的设置，尖峰应该在第一个箱中：对我来说看起来还不错。 - Fusho

@Sito：有两个峰是因为Eigen返回一个双频谱（如果信号是实数，则对称）。 - Fusho

Answer 1

从你计算的time(u)来看，我会说你的采样周期Ts是2*pi/1000 [s]，这导致Fs = 1/Ts = 1000/(2*pi) [Hz]。你计算的正弦波的模拟频率f0将会是：

1*t = 2*pi*f0*t [radians]
f0 = 1/(2*pi) [Hz]

请注意，Fs >> f0。

在数字领域中，频率总是跨越2*pi [弧度]（它可以是[-pi,pi)或[0,2*pi) ，但Eigen返回后者）。因此，您需要一致地将范围[0,2*pi)分成N个bin。例如，如果索引为k，则相关的归一化频率为f=2*pi*k/N [弧度]。

要知道每个归一化频率bin对应的模拟频率f，请计算f = (fs*k/N) [赫兹]，其中fs是采样频率。

关于Eigen FFT doc的缩放和全频谱特性：

1) 缩放：其他库（FFTW、IMKL、KISSFFT）不执行缩放，因此在正反变换后会产生恒定的增益，所以 IFFT(FFT(x)) = Kx；这样做是为了避免进行向量乘法。缺点是在 Matlab/octave 中正确工作的算法在 C++ 中实现后不会表现出相同的行为。Eigen/FFT 的区别在于：执行可逆缩放，因此 IFFT(FFT(x)) = x。

2) 实数 FFT 半频谱：其他库仅使用频率谱的一半（加上一个额外的样本用于奈奎斯特 bin），另一半是第一半的共轭对称。这节省了一些内存和复制。缺点是调用者需要针对复杂 vs 实数的 bin 数量具有特殊逻辑。Eigen/FFT 的区别在于：从正向变换返回完整的频谱。这通过消除实数 vs 复杂分开专门化来促进通用模板编程。在逆变换中，如果输出类型为实数，则实际上只使用了一半的频谱。

因此，您应该期望获得收益，只需进行测试ifft(fft(x)) == x（测试为“误差功率”<<“信号功率”）。 您可以除以N以获得归一化版本。

另一方面，您看到的两个峰是由点2引起的。您上面发布的图只是变换的一侧，如果信号是实数，则另一侧是对称的。您可以删除输出的高半部分。

---

这段代码：

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

unsigned const N = 1000;  //
double const Fs  = 32;    // [Hz]
double const Ts  = 1./Fs; // [s] 
const double f0  = 5;     // [Hz]
std::complex<double> f(std::complex<double> const & t){
    return std::sin(2*M_PI*f0*t);
}

int main(){
    std::ofstream xrec("xrec.txt");
    Eigen::VectorXcd time(N);
    Eigen::VectorXcd f_values(N);
    Eigen::VectorXd freq(N);
    for(int u = 0; u < N; ++u){
        time(u) = u * Ts;
        f_values(u) = f(time(u));
        freq(u) = Fs * u / double(N);
    }

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

    for(int u = 0; u < N; ++u){
        xrec << freq(u) << " " << std::abs(f_freq(u)) << "\n"; 
    }
}

生成 xrec.txt 文件。然后，您可以使用此 gnuplot 脚本生成图形：

set key off
set grid
set output "figure.png"
set xlabel "Frequency [Hz]"
plot [-1:34] [-10:500] "xrec.txt" with impulses, "xrec.txt" with points pt 4

在图中，您可以看到5和27 Hz处有两个尖峰，正如此代码所预期的那样。我更改了这些值以更好地理解发生了什么，只需尝试其他任何值即可。

在您展示的绘图风格中，X 轴范围为 [0,16)，而不是此绘图中的 [0,32)，但由于您的信号是实数，谱是对称的，因此您可以舍弃其中一半。

- Fusho

---

我根据您建议的计算方法修改了我的问题，并绘制了结果。但是这似乎很奇怪，有两个峰值，而且与我试图复制的图片相比，缩放似乎有点奇怪... - Sito

Answer 2

通常情况下，图书馆使用以下公式计算DFT：

X[k] = sum_n(x[n] * exp(-2*pi * i * k * n/N)

其中，

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

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

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

在k>N/2的索引处，由于混叠，频率非常高，实际上对应于较低的频率。

这是N=8的示例：

我没有特别查询过Eigen，但我认为它并没有什么不同。

- CygnusX1

---

快速傅里叶变换难道不是和离散傅里叶变换不一样吗？你的回答在两种情境下都有意义吗？因为我很确定我的例子中使用的算法是FFT而不是DFT。 - Sito

离散傅里叶变换是一种数学运算。它解释了输出和输入之间的关系，但并没有说明如何高效地计算它。快速傅里叶变换是执行该操作的算法。由于输入和输出都是数字化的，因此我们谈论的是离散而不是连续的傅里叶变换。 - CygnusX1

首先，感谢您的回答！不幸的是，这并不是我想要的。我需要/想要的是一个使用 Eigen 的 FFT 算法工作示例，而不是理论描述。无论如何，感谢您花费时间撰写此内容。 - Sito

我认为答案就在这里？freq = k，其中k是您的f_freq矩阵中的列索引。 - CygnusX1