c++ - 使用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 ...)。

以下是将上述代码段放入单个文件中:

#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)范围一律划分为N箱。例如，如果index为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.

因此，您应该期望获得 yield ，只需测试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 Hz和27 Hz的两个尖峰。我更改了值以更好地了解发生了什么，只需尝试其他任何值即可。



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