如何用FFTW进行频谱分析

Question

我想更好地理解FFTW的API。FFTW是一个用于计算一维或多维离散傅里叶变换(DFT)的库。

现在，假设我有一个窦形波形x=30*sin(2*M_PI*f* I* T)，其中f是频率(例如f=1000Hz)。如果我使用FFTW的函数来分析我的波形，我希望得到一个频率f=1000Hz。

我的问题是如何使用FFTW库在c++中做到这一点?

Answer 1

您可以在FFTW的文档中找到更多详细信息。

然而，对于一维实值信号的相对简单的情况，下面是你必须做的一般步骤的总结。

通常您需要分配输入/输出缓冲区，以及FFTW用于自己记账的数据结构，库将其称为plan。这可以通过多种方式实现(详见FFTW的文档)，例如:

  #include "fftw3.h"
  // First choose a buffer size:
  //   Typically best performance with a power of 2
  //   but could be a product of small primes
  int           input_size    = 1024; 
  //   Compute corresponding number of complex output samples
  int           output_size   = (input_size/2 + 1);
  // Allocate input and output buffers
  double*       input_buffer  = static_cast<double*      >(fftw_malloc(input_size  * sizeof(double)));
  fftw_complex* output_buffer = static_cast<fftw_complex*>(fftw_malloc(output_size * sizeof(fftw_complex)));
  // Create plan
  //   Select plan creation flags
  //   see http://www.fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags
  int           flags = FFTW_ESTIMATE;
  fftw_plan     plan = fftw_plan_dft_r2c_1d(input_size, 
                                            input_buffer, 
                                            output_buffer, 
                                            flags);

一旦完成，您可以用实值数据样本填充input_buffer进行分析，并使用

执行FFT:

 fftw_execute(plan);

复值结果将存储在output_buffer中，其中output_buffer[0]对应频率为0,output_buffer[output_size-1]对应采样率的一半。该计划可以执行多次(input_buffer中的值更新，导致output_buffer中的值相应更新)。

请注意，通常fftw_complex(这是本例中用于输出的数据类型)是作为一个包含2个值的数组实现的:索引0对应于实数部分，索引1对应于虚数部分(例如output_buffer[i][0]对应于i频率分量的实数部分)。

完成后，可以使用以下命令释放已分配的资源:

  fftw_free(input_buffer);
  fftw_free(output_buffer);
  fftw_destroy_plan(plan);

请注意，如果您可以使用这些函数的float, double或long double版本。链接到相应的libfftw3f-3.lib, libfftw3-3.lib或libfftw3l-3.lib。

更新:如果你想和fftw_plan_dft_1d一起使用复值输入样本，那么你必须像这样设置实部和虚部:

for (i = 0; i < N-1; ++i) {
  t[i]=i*T;
  signal[i][0] = 0.7 * sin(2*M_PI*f*t[i]); // real-part
  signal[i][1] = 0.0; // imaginary-part
}

或者将输入样本类型更改为float, double或long double(同时使用fftw_plan_dft_r2c_1d)。

Answer 2

你的问题属于真实数据的一维dft范畴

step1:通过您提供的x = 30(2*pi*1000*id函数生成频率为1000hz的波形数据(确保采样频率为2乘以2*f，即2000，我建议您通过范围(0,1,1/2000)迭代id)，从而为您提供2000个样本数据。

step2:使用函数得到DFT输出

fftw_plan fftw_plan_dft_r2c_1d(int n, double *in, fftw_complex *out,
                                    unsigned flags);