使用 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我要找的是图片中 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();
}

结果如下图这似乎有点不对劲.. 缩放当然没有多大意义，但事实上有两个尖峰值让我有点怀疑..

Answer 1

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

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

哪里

X - 是傅里叶域数组。这些值表示相应 sine/cosine 函数的振幅。
k - 是傅立叶域中的索引，也唯一定义了频率
i - 这只是数学上的 i - 复数 ( 0+1i )
N - 是数组的大小

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

X[0]是你的平均值
X[1] 对应于 sine/cosine 函数，在您的整个域中恰好适合一次
X[2] 对应于适合您域两次的 sine/cosine 函数 ...等等...

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

这里是 N=8 的例子：

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

Answer 2

根据你对 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 return 后者）。因此，您需要将范围 [0,2*pi) 始终如一地划分为 N 个箱子。例如，如果索引为 k，则关联的归一化频率为 f=2*pi*k/N [radians]。

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

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

1) Scaling: Other libraries (FFTW,IMKL,KISSFFT) do not perform scaling, so there is a constant gain incurred after the forward&inverse transforms , so IFFT(FFT(x)) = Kx; this is done to avoid a vector-by-value multiply. The downside is that algorithms that worked correctly in Matlab/octave don't behave the same way once implemented in C++. How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.

2) Real FFT half-spectrum: Other libraries use only half the frequency spectrum (plus one extra sample for the Nyquist bin) for a real FFT, the other half is the conjugate-symmetric of the first half. This saves them a copy and some memory. The downside is the caller needs to have special logic for the number of bins in complex vs real. How Eigen/FFT differs: The full spectrum is returned from the forward transform. This facilitates generic template programming by obviating separate specializations for real vs complex. On the inverse transform, only half the spectrum is actually used if the output type is real.

所以，你应该期望有所收获，只需进行测试ifft(fft(x)) == x（测试为"error power" << "signal power"）。您可以除以 N 以获得标准化版本。

另一方面，您看到的两个尖峰是因为第 2 点。您上面 post 的图只是变换的一侧，如果信号是真实的，另一侧是对称的。您可以删除输出的高半部分。

此代码：

#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)，但是，由于您的信号是真实的，因此频谱是对称的，您可以掉一半。

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

Frequencies when performing FFT with Eigen::FFT

c++

math

fft

eigen

eigen3