FFTW3.3.10库与QT结合的使用
FFTW(Fastest Fourier Transform in the West)是世界上最快的FFT, 实测计算长度为10000的double数组, 单次运行时间在2ms左右。为了详细了解FFTW以及为编程方便,特将用户手册看了一下,并结合手册制作了以下FFTW中文参考。其中大部分是原文重点内容的翻译,并加入了一些注解。
先看一下使用FFTW编程的方法:
#include <fftw3.h> ... { fftw_complex *in, *out; fftw_plan p; ... in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); // 输入数据in赋值 p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(p); // 执行变换 ... fftw_destroy_plan(p); fftw_free(in); fftw_free(out); } |
大致是先用fftw_malloc分配输入输出内存,然后输入数据赋值,然后创建变换方案(fftw_plan),然后执行变换(fftw_execute),最后释放资源,还是比较简单的。
一维复数据的DFT
1. 数据类型
fftw_complex默认由两个double组成,在内存中顺序排列,实部在 前,虚部在后,即typedef double fftw_complex[2]。FFTW文档指出如果有一个支持C99标准的C编译器(如gcc),可以在#include <fftw3.h>前加入#include <complex.h>,这样一来fftw_complex就被定义为本机复数类型,而且与上述typedef二进制兼容(指内存排列),经 测试不能用在Windows下。C++有一个复数模板类complex<T>,在头文件<complex>下定义。C++标准委 员会最近同意该类的存储方式与C99二进制兼容,即顺序存储,实部在前,虚部在后(见报告WG21/N1388),该解决方案在所有主流标准库实现中都能正确工作。所以实际上可以用complex <double> 来代替fftw_complex,比如有一个复数数组complex<double> *x,则可以将其类型转换后作为参数传递给fftw:reinterpret_cast<fftw_complex*>(x)。测试如下:开 两个数组fftw_complex x1[2]和complex<double> x2[2],然后赋相同值,在调试模式下可以看到它们的内存排列是相同的。complex<T>类数据赋值的方式不是很直接,必须采用无名对象方式x2[i] = complex <double>(1,2) 或成员函数方式x2[i].real(1);x2[i].imag(2);不能直接写x2[i].real=1;x2[i].imag=2。 fftw_complex赋值方式比较直接:x1[i][0]=1;x1[i][1]=2。最后,考虑到数据对齐(见后),最好使用 fftw_malloc分配内存,所以可以将其返回的指针转换为complex <double> *类型使用(比如赋值或读取等),变换时再将其转换为fftw_complex*。
2. 函数接口
fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); |
n为数据个数,可以为任意正整数,但如果为一些小因子的乘积计算起来可以更有效,不过即使n为素数算法仍然能够达到O(nlogn)的复杂度。FFTW对N=2a 3b 5c 7d 11e 13f的变换处理得最好,其中e+f=0/1,其它幂指数可以为任意值。
如果in和out指针相同为原位运算,否则为非原位运算。
sign可以为正变换FFTW_FORWARD(-1),也可以为逆变换FFTW_BACKWORD(+1),实际上就是变换公式中指数项的符号。需注意FFTW的逆变换没有除以N,即数据正变换再反变换后是原始数据的N倍。
flags参数一般情况下为FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW会先计算一些FFT并测量所用的时间,以便为大小为n的变换寻找最优的计算方法。依据 机器配置和变换的大小(n),这个过程耗费约数秒(时钟clock精度)。FFTW_ESTIMATE则相反,它直接构造一个合理的但可能是次最优的方案。总体来说,如果你的程序需要进行大量相同大小的FFT,并且初始化时间不重要,可以使用FFTW_MEASURE,否则应使用 FFTW_ESTIMATE。FFTW_MEASURE模式下in和out数组中的值会被覆盖,所以该方式应该在用户初始化输入数据in之前完成。
不知道上述说法是不是这个意思:先用FFTW_MEASURE模式自动选最优方案,速度较慢;然后使用该模式变换数据就会较快。示例代码为:
int length = 50000; fftw_complex* din = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); fftw_complex* dout = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2); fftw_plan p = fftw_plan_dft_1d(length, din, din, FFTW_FORWARD, FFTW_MEASURE); fftw_execute(p); // 输入数据din赋值 // ... fftw_execute(p); // 读取变换结果 // ... fftw_destroy_plan(p); fftw_free(din); fftw_free(dout); |
实验发现第一个fftw_execute耗费了数秒,而第二个fftw_execute则瞬间完成,说明上述猜想可能是对的。
创建完方案(fftw_plan)后,就可以用fftw_execute对指定的 数据in/out做任意次变换。如果想变换一个相同大小(N相等)但数据不同的另外一个数组in,可以创建一个新方案,FFTW会自动重用上次方案的信息。这一点其实是非常好的,比如你首先用FFTW_MEASURE模式创建了一个最优的变换方案,只要变换数据的大小不变,你可以用 fftw_plan_dft_1d创建新的方案以对新数据执行变换,同时新变换仍然是最优的。
一个fftw_plan只能对固定的in/out进行变换, 但可以在变换后改变in的内容(大小不变)以用同一个方案执行新的变换。
软件整体界面如下:
1、可设置正弦信号的频率、采样频率、采样点数、直流信号、是否加噪声。
2、点击“计算FFT”按钮可生成频谱。
3、还可以进行自动计算。信号频谱按10Hz依次增加,再依次计算这个频率对应的频谱。运行效果如下:
将原始信号修改为 sin(f0) + sin(f1) + sin(f2) + sin(f3) 模拟4个不同频率的信号合成信号。
再看它的FFT变换结果。
源码出售,完整代码,包教会