使用 C 中的 fftw.h 计算 fft 和 ifft
大家好 我正在使用 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 文档。 它提到了一些缩放,但是我使用“r2c”和“c2r”变换而不是 FFT_FORWARD 和 FFT_BACKWARD。任何见解将不胜感激。
Hi all
I am using the fftw C libraries to compute the frequency spectrum for some signal processing applications on embedded systems. However, in my project I have run into a slight hinderence.
Below is a simple program I wrote to ensure I am implementing the fftw functions correctly. Basically I want to calculate the fft of a sequence of 12 numbers, then do the ifft and obtain the same sequence of numbers again. If you have fftw3 and gcc installed this program should work if you compile with:
gcc -g -lfftw3 -lm fftw_test.c -o fftw_test
Currently my fft length is the same size as the input array.
#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;
}
Which Produces the following output:
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
So obviously the ifft is producing some scaled up result. In the fftw docs found here:
fftw docs about scaling.
It mentions about some scaling, however I am using the "r2c" and "c2r" transforms rather than the FFT_FORWARD and FFT_BACKWARD. Any insight would be appreciated.
如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱,如果其他用户或者作者回复了您的评论,将不能在第一时间通知您!
发布评论
评论(3)
查看 优秀文档对于您使用的函数,您将看到您正在使用FFT_FORWARD和FFT_BACKWARD,而且正是它的预期用途。因此,您之前找到的缩放信息也适用于此处。
Looking at the great documentation for the functions you use, you will see you are using FFT_FORWARD and FFT_BACKWARD, and exactly where it is intended. Therefore, the scaling information you found previously also applies here.
很抱歉显得迂腐,但您的 out_cpx 大小不正确。它不应该是 size long,而应该是 size/2 + 1。这是因为真实信号的 FFT 是埃尔米特式的。您可以通过将 out_cpx 初始化为某个随机数(全部为 3.14159)来验证我所说的。向前和向后运行,然后打印 out_cpx 从 size/2 + 1 到 size。它不会改变。
http://www.fftw.org /fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format
Sorry to be pedantic, but your size for out_cpx is incorrect. instead of being size long, it should be size/2 + 1. This is because FFT of a real signal is Hermitian. You can verify what I say by initializing out_cpx to some random number (all 3.14159). Run both the forward and backward and then print out out_cpx from size/2 + 1 to size. It will not have changed.
http://www.fftw.org/fftw3_doc/Real_002ddata-DFT-Array-Format.html#Real_002ddata-DFT-Array-Format
r2c 和 c2r 的作用本质上与常规傅里叶变换相同。唯一的区别是输入和输出数组都需要保存一半的数字。请看一下FFTW r2c和c2r手册的最后一段。因此,归一化因子恰好是实际数组的元素数量,或者您的情况下的变量
size
(== 12)。r2c and c2r do essentially the same as the regular Fourier transform. The only difference is that both the input and output array need to hold half of the numbers. Please take a look at the last paragraph of the manual of FFTW r2c and c2r. So the normalization factor is precisely the number of elements of the real array, or the variable
size
(== 12) in your case.