FFTW中文参考

转自CSDN博主:FFTW中文参考

一、 简介

先看一下使用FFTW编程的方法:

  #include <fftw3.h>

...

{

fftw_complex*in, *out;

fftw_planp;

...

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 前加入#include ,这样一来fftw_complex就被定义为本机复数类型,而且与上述typedef二进制兼容(指内存排列),经 测试不能用在Windows下。C++有一个复数模板类complex,在头文件下定义。C++标准委 员会最近同意该类的存储方式与C99二进制兼容,即顺序存储,实部在前,虚部在后(见报告WG21/N1388),该解决方案在所有主流标准库实现中都能正确工作。所以实际上可以用complex 来代替fftw_complex,比如有一个复数数组complex *x,则可以将其类型转换后作为参数传递给fftw:reinterpret_cast(x)。测试如下:开 两个数组fftw_complex x1[2]和complex x2[2],然后赋相同值,在调试模式下可以看到它们的内存排列是相同的。complex类数据赋值的方式不是很直接,必须采用无名对象方式x2[i] = complex (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 *类型使用(比如赋值或读取等),变换时再将其转换为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模式自动选最优方案,速度较慢;然后使用该模式变换数据就会较快。示例代码为:


intlength = 50000;

fftw_complex* din  = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2);

fftw_complex* dout = (fftw_complex *)fftw_malloc(sizeof(double)*length * 2);

fftw_planp  = 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的内容(大小不变)以用同一个方案执行新的变换。

三、 多维复数据的DFT

    fftw_plan fftw_plan_dft_2d(int n0, intn1,

fftw_complex *in, fftw_complex*out,

int sign, unsignedflags);

fftw_plan fftw_plan_dft_3d(int n0, int n1, intn2,

fftw_complex *in, fftw_complex*out,

int sign, unsignedflags);

fftw_plan fftw_plan_dft(int rank, const int*n,

fftw_complex *in, fftw_complex *out,

int sign, unsigned flags);

多维复数据的DFT同一维复数据的DFT用法类似,数组in/out为行优先方式 存储。fftw_plan_dft是一个通用的复DFT函数,可以执行一维、二维或多维复DFT。比如对于图像(2维数据),其变换为 fftw_plan_dft_2d(height,width,85),因为原始图像数据为height×width的矩阵,即第一维(n0)为行数 height。

四、 一维实数据的DFT

函数接口

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

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

r2c版本:实输入数据,复Hermitian输出,正变换。

c2r版本:复Hermitian输入数据,实输出数据,逆变换。

n:逻辑长度,不必为物理长度。由于实数据的DFT具有 Hermitian对称性,所以只需要计算n/2+1(向下取整)个输出就可以了。比如对于r2c,输入in有n个数据,输出out有floor(n /2)+1个数据。对于原位运算,in和out为同一数组(out须强制类型转换),所以其必须足够大以容纳所有数据,长度为2*(n/2+1),in数 组的前n个数据为输入数据,后面的数据用来保存输出。

c2r逆变换在任何情况下(不管是否为原位运算)都破坏输入数组in,如果有必要可以通过在flags中加入FFTW_PRESERVE_INPUT来阻止,但这会损失一些性能,而其这个标志位目前在多维实DFT中已不被支持。

比如对于n=4,in=[1 2 3 4],out=[10 -2+2i -2 -2-2i],out具有共轭对称性,out的实际内存为10 0 -2 2 -2 0,共3个复数数据。对于n=5,in=[1 2 3 4 5],out=[15 -2.5+3.44i -2.5+0.81i -2.5-0.81i -2.5-3.44i] ,out的实际内存为15 0 -2.5 3.44 -2.5 0.81,共3个复数数据。

实数据DFT中,首个变换数据为直流分量,为实数;如果n为偶数,由 Nyquist采样定理,第n/2个变换数据也为实数;所以可以把这两个数据组合在一起,即将第n/2个变换数据作为第0个变换数据的虚部,这样一来输入 数组就和输出数组相等(n=n/2*2)。一些FFT的实现就是这么做的,但FFTW没有这么做,因为它并不能很好地推广到多维DFT中,而且存储空间的 节省也是非常小以至于可以忽略不计。

一个一维c2r和r2c DFT的替代接口可以在r2r接口中找到,它利用了半复数输出类型(即实部和虚部分开放在不通的区域里),使输出数组具有和输入数组同样的长度和类型。该接口在多维变换中用处不大,但有时可能会有一些性能的提升。

五、 多维实数据的DFT

    fftw_plan fftw_plan_dft_r2c_2d(int n0, intn1,

double *in, fftw_complex*out,

unsignedflags);

fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, intn2,

double *in, fftw_complex*out,

unsignedflags);

fftw_plan fftw_plan_dft_r2c(int rank, const int*n,

double *in, fftw_complex*out,

unsigned flags);

这是r2c接口(正变换),c2r接口(逆变换)只是简单的将输入/输出类型交换一下。用法大致同1d情况一样,需要特别注意的是复数据的存放方式。对于nnn1×…×nd-1的实输入数组(行优先),经过r2c正变换后,输出为一个nnn1×…×(nd-1/2+1)的复数(fftw_complex)数组(行优先),其中除法向下取整。由于复数数据的总长度要大于实数据,所以如果需要进行原位运算则输入实数组必须扩展以能够容纳所有复数据,即实数组的最后一维必须包含2(floor(nd-1/2)+1)个double元素。比如对于一个2维实正DFT,输入数组为nn1大小,输出复数组大小为n0×floor(n1/2+1)(由对称性),其长度大于实输入数组。所以对于原位运算,输入数组要扩展到n0×2floor(n1/2+1)。如果n1为偶数,扩展为n0×(n1+2);如果n1为奇数,扩展为n0×(n1+1);这些扩展的内存不需要赋初值,它们只用来存放输出数据。

比如对于3×3的实正DFT,in=[0 2 4;6 1 3;5 7 4],out=[32 0.5+0.86i 0.5-0.86i;-7+5.2i -1-1.73i -8.5-6.06i;-7-5.2i -8.5+6.06i -1+1.73i];out的实际内存为32,0,0.5,0.86,-7,5.2,-1,-1.73,-7,-5.19,-8.5,6.06;即为 3×2的复数组,换算成double类型为3×4,所以如果要进行原位运算,in数组大小必须为3×4,即最后一维(每行)扩展一个double元素。另 外,如果用这个out数组进行3×3的c2r逆变换,将得到实数据[0 18 36;54 9 27;45 63 36],即原始数据的9(nn1)倍,这是因为FFTW的逆变换没有归一化。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 220,367评论 6 512
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,959评论 3 396
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 166,750评论 0 357
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 59,226评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 68,252评论 6 397
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,975评论 1 308
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,592评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,497评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 46,027评论 1 319
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 38,147评论 3 340
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,274评论 1 352
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,953评论 5 347
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,623评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,143评论 0 23
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,260评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,607评论 3 375
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,271评论 2 358

推荐阅读更多精彩内容