用特征进行FFT时的频率：：FFT

Question

我目前正在努力弄清楚如何准确地使用Eigen的快速傅立叶变换算法。

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

从你计算的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桶中。例如，如果索引为k，则相关的归一化频率为f=2*pi*k/N [radians]。

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

关于Eigen FFT doc的标度和全光谱特性

1)缩放:其他库(FFTW、IMKL、KISSFFT)不执行缩放，因此在正变换和逆变换之后会产生恒定的增益，因此FFT(FFT(X))= Kx；这样做是为了避免向量逐值乘。缺点是，在Matlab/octave中正确工作的算法一旦在C++中实现，其行为方式就不一样了。特征/FFT的不同之处:反向缩放，因此IFFT( FFT( x ) )=x。

2)实FFT半谱:其他库只使用一半的频谱(加一个额外的样本用于实FFT )，另一半是前半部分的共轭对称。这为他们节省了一份拷贝和一些内存。缺点是调用者需要对复杂的和真实的回收箱的数量有特殊的逻辑。特征/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)，但是，由于信号是真实的，所以频谱是对称的，你可以把它降一半。

Answer 2

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

Xk = sum_n(xn * exp(-2*pi *i*k* n/N)

哪里

• X -是傅里叶域数组.这个值代表了相应的正弦/余弦函数的振幅。
• k-是傅里叶域的指数，它也唯一地定义了频率
• i-它只是数学上的i-复数( 0+1i )
• N -是数组

的大小。

所以，在索引k，你的频率有你整个信号输入的1/k长度。特别是：

• X[0]是您的平均value
• X[1]对应的正弦/余弦函数，它恰好适合于整个domain
• X[2]对应的正弦/余弦函数，而正弦/余弦函数在您的域上适合两次。等等.

在指数k>N/2时，频率很高，实际上对应于混叠引起的较低频率。

下面是N=8的一个示例：

​

​

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