如何用 FFT 进行频谱分析？ [关闭]

Posted 2023-02-25

【中文标题】如何用 FFT 进行频谱分析？ [关闭]

【英文标题】：how to do spectrum analyse with FFTW? [closed]

【发布时间】：2014-05-27 00:54:17

【问题描述】：

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

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

我的问题是如何使用 FFTW 库在 c++ 中执行此操作？ 任何帮助将不胜感激。

【参考方案1】：

您可以在FFTW's documentation找到更多详细信息。

但是，对于一维实值信号的相对简单的情况，以下是您必须执行的一般步骤的总结。

通常，您需要分配输入/输出缓冲区，以及 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^th 频率分量的实部）。

完成后，您可以通过以下方式释放分配的资源：

  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）。

【参考方案2】：

您的问题属于One-Dimensional DFTs of Real Data的类别

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

step2：使用函数获取dft输出

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