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

Question

我目前正试图弄清楚如何使用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和我寻找的是freqin的值图片……）。

以下是放在单个文件中的上述代码片段：

#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();
}

这导致以下情节 这似乎有点不对劲..缩放当然没有多大意义，但有两个值飙升的事实让我有点怀疑..

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 [radians]（可以是[-pi,pi)或[0,2*pi)，但Eigen返回后者）。因此，您需要一致地将范围划分[0,2*pi)为Nbin。例如，如果索引为k，则关联的归一化频率为f=2*pi*k/N [radians]。

要了解哪个模拟频率f对应于每个归一化频率仓，请计算f = (fs*k/N) [Hz]，其中fs是采样频率。

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

1）缩放：其他库（FFTW，IMKL，KISSFFT）不执行缩放，因此在正变换和逆变换后会产生恒定的增益，因此IFFT（FFT（x））= Kx；这样做是为了避免向量乘以值。缺点是在 Matlab/octave 中正常工作的算法在 C++ 中实现后的行为方式并不相同。Eigen/FFT 有何不同：执行可逆缩放，因此 IFFT( FFT(x) ) = x。

2）真实FFT半频谱：其他库仅使用一半频谱（加上奈奎斯特仓的一个额外样本）进行真实FFT，另一半是前半部分的共轭对称。这可以为他们节省一份副本和一些内存。缺点是调用者需要对复杂与真实的 bin 数量有特殊的逻辑。Eigen/FFT 有何不同： 全频谱由正向变换返回。这通过消除真实与复杂的单独专业化来促进通用模板编程。在逆变换中，如果输出类型为实数，则实际仅使用频谱的一半。

因此，您应该期望获得增益，只需进行测试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 赫兹处的两个尖峰，正如此代码所预期的那样。我更改了这些值以更好地了解发生的情况，只需尝试其他值即可。

在您显示的图的样式中，x 轴范围[0,16)而不是 [0,32)如图所示，但是，由于您的信号是真实的，因此频谱是对称的，您可以删除一半。