在 C 中使用 fftw.h 计算 fft 和 ifft答案

【问题标题】：computing fft and ifft with fftw.h in C在 C 中使用 fftw.h 计算 fft 和 ifft
【发布时间】：2011-08-14 16:24:14
【问题描述】：

大家好我正在使用 fftw C 库来计算嵌入式系统上某些信号处理应用程序的频谱。但是，在我的项目中，我遇到了一点阻碍。

以下是我编写的一个简单程序，用于确保正确实现 fftw 功能。基本上我想计算12个数字序列的fft，然后做ifft并再次获得相同的数字序列。如果您安装了 fftw3 和 gcc，则该程序应该可以工作，如果您使用以下代码进行编译：

gcc -g -lfftw3 -lm fftw_test.c -o fftw_test

目前我的 fft 长度与输入数组的大小相同。

#include <stdio.h>
#include <stdlib.h>
#include <sndfile.h>
#include <stdint.h>
#include <math.h>
#include <fftw3.h>

int main(void)
{
double array[] = {0.1, 0.6, 0.1, 0.4, 0.5, 0, 0.8, 0.7, 0.8, 0.6, 0.1,0};
//double array2[] = {1, 6, 1, 4, 5, 0, 8, 7, 8, 6, 1,0};
double *out;
double *err;
int i,size = 12;

fftw_complex *out_cpx;

fftw_plan fft;
fftw_plan ifft;
out_cpx = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*size);
out = (double *) malloc(size*sizeof(double));
err = (double *) malloc(size*sizeof(double));

fft = fftw_plan_dft_r2c_1d(size, array, out_cpx, FFTW_ESTIMATE);  //Setup fftw plan for fft
ifft = fftw_plan_dft_c2r_1d(size, out_cpx, out, FFTW_ESTIMATE);   //Setup fftw plan for ifft

fftw_execute(fft);
fftw_execute(ifft);

//printf("Input:    \tOutput:    \tError:\n");
printf("Input:    \tOutput:\n");
for(i=0;i<size;i++)
{
err[i] = abs(array[i] - out[i]);    
printf("%f\t%f\n",(array[i]),out[i]);
//printf("%f\t%f\t%f\n",(array[i]),out[i],err[i]);
}

fftw_destroy_plan(fft);
fftw_destroy_plan(ifft);
fftw_free(out_cpx);
free(err);
free(out);
return 0;
}

产生以下输出：

Input:      Output:
0.100000    1.200000
0.600000    7.200000
0.100000    1.200000
0.400000    4.800000
0.500000    6.000000
0.000000    0.000000
0.800000    9.600000
0.700000    8.400000
0.800000    9.600000
0.600000    7.200000
0.100000    1.200000
0.000000    0.000000

显然，ifft 正在产生一些放大的结果。在此处找到的 fftw 文档中： fftw docs about scaling。它提到了一些缩放，但是我使用的是“r2c”和“c2r”转换而不是 FFT_FORWARD 和 FFT_BACKWARD。任何见解将不胜感激。

【问题讨论】：

标签： c signal-processing fft fftw

【解决方案1】：

查看the great documentation for the functions you use，您将看到您正在使用 FFT_FORWARD 和 FFT_BACKWARD，并且准确地知道它的用途。因此，您之前找到的缩放信息也适用于此。

【讨论】：

好的，谢谢，是的，我确实在某个时候读过。这是漫长的一天......是的，解决了：out[i] = out[i]/size;

【解决方案2】：

抱歉，您的 out_cpx 大小不正确。它应该是 size/2 + 1，而不是 size/2 + 1。这是因为真实信号的 FFT 是 Hermitian。您可以通过将 out_cpx 初始化为某个随机数（全部为 3.14159）来验证我所说的内容。向前和向后运行，然后从 size/2 + 1 打印出 out_cpx 到 size。它不会改变。

http://www.fftw.org/fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format

【讨论】：

非常感谢您的意见：) 现在阅读文档。

【解决方案3】：

r2c 和 c2r 与常规傅里叶变换基本相同。唯一的区别是输入和输出数组都需要保存一半的数字。请看FFTW r2c and c2r手册的最后一段。因此，归一化因子正是实际数组的元素数，或者在您的情况下是变量 size (== 12)。

【讨论】：