在 c 中使用 fftw 进行谐波分析

Question

所以我有一个周期数列，我需要在没有一次谐波的情况下获得该数列的最高数。

系列大小为 360

所以我有

f(0)=value1
f(1)=value2
...
f(359)=value360

有人建议要做到这一点，我必须得到 360 点离散傅里叶变换 DFT (f(0), f(1)... f(359)) --> (F(0), F(1 )... F(359)) 然后将 F(0) 和 F(359) 设置为 0，这将消除一次谐波。在此之后，我应该进行 360 点逆离散傅立叶变换 iDFT，然后在结果中搜索最大值。

为此，我在 c 中使用 fft 库，但在尝试弄清楚如何正确使用它们来执行此操作时遇到了一些麻烦。我的周期数列由不复数的实数组成，所以我这样做：

#include "complex.h"
#include "fftw.h"
#include "rfftw.h"

...

fftw_real in_r[360]; //input 1d array of real numbers
fftw_complex out_c[360]; //output 1d array of complex numbers
rfftwnd_plan p_DFT; //plan to calculate the DFT
rfftwnd_plan p_iDFT; //plan to calculate the iDFT

p_DFT = rfftwnd_create_plan(1, 360, FFTW_REAL_TO_COMPLEX, FFTW_MEASURE);
p_iDFT = rfftwnd_create_plan(1, 360, FFTW_COMPLEX_TO_REAL, FFTW_MEASURE);

rffftwnd_one_real_to_complex(p_DFT,in_r,out_c);

//I GET AN ERROR IN BOTH THIS CALLS: incompatible types when assigning to type
//‘fftw_complex’ from type ‘complex double’
out_c[0]=0.0 + 0.0*_Complex_I;
out_c[359]=0.0 + 0.0*_Complex_I;

rfftwnd_one_complex_to_real(p_iDFT,out_c,in_r);

for(i=0;i<360;i++)
  max=fmaxf(in_r[i],max);

所以我有几个问题。

首先，我如何将输出数组的第一个和最后一个元素设置为 0，考虑到它很复杂但不允许我为其分配一个复数？

其次，我这样做是否正确？还是我错过了什么？

第三（这是我接下来需要做的后续工作）。我可以使用 FFTW 库来获取该系列的第一到第四谐波的幅度和相位吗？如果有，怎么做？

感谢您的帮助。

更新：

我变了

#include "fftw.h"

为

#include "fftw3.h"

所以在我的收录中我有

#include "complex.h"
#include "fftw3.h"
#include "rfftw.h"

但我遇到了同样的错误

更新 2：

我也遇到这样的错误

/usr/include/fftw.h:307:13: error: conflicting types for ‘fftw_destroy_plan’

因为在 fftw3.h 之后包含了 rfftw.h。但是如果我删除 rfftw.h 就会出现这样的错误

error: unknown type name ‘fftw_real’

fftw3.h 和 rfftw.h 之间似乎有冲突，但我无法删除 rfftw.h，因为那样我就无法使用我需要的功能。

Answer 1

参考这个：

fftw_complex out_c[360]; //output 1d array of complex numbers

...

//I GET AN ERROR IN BOTH THIS CALLS: incompatible types when assigning to type
//‘fftw_complex’ from type ‘complex double’
out_c[0]=0.0 + 0.0*_Complex_I;
out_c[359]=0.0 + 0.0*_Complex_I;

要在 C99 模式（选项 -std=c99）下在 gcc 上进行这项工作，您需要以正确的顺序包含标头：

#include <complex.h>
#include <fftw3.h>

这样fftw_complex应该被定义为原生复杂类型。

---

否则是typedef double fftw_complex[2]，所以分配是这样的：

out_c[0][0] = 0.0; // real part
out_c[0][1] = 0.0;  // imaginary part
...

---

请参阅here for details 了解如何实现复数。

---

rfftw.h在 fftw3 中已过时。来自the update-to-fftw3 page：

FFTW 2 具有单独的数据类型 fftw_plan、fftwnd_plan、rfftw_plan 和 rfftwnd_plan，用于复杂和真实的一维和多维变换，并且每种类型都有自己的“销毁”功能。在 FFTW 3 中，所有计划都是 fftw_plan 类型，并且都被 fftw_destroy_plan(plan) 销毁。