利用FFTW实现快速傅里叶变换

时间:2022-03-25 07:00:51

摘要: 本文对FFTW的工作原理及机制作了详细的分析,并通过实例讲解了FFTW实现FFT的具体方法及性能的测试。

关键词: FFTW 工作原理 使用方法 性能评测

引言

快速傅里叶变换(Fast Fourier transform简称FFT)作为时域和频域转换的分析工具,广泛应用在数字通信、音频信号处理、图像处理、 谱估计、 雷达信号处理、地震勘探等各领域中。在FFT的软件实现中,比较有名的库有fftpack,mkl及FFTW等。

一、FFT及FFTW

FFT是离散傅里叶变换的快速算法,也可用于计算离散傅里叶变换的逆变换。快速傅里叶变换有广泛的应用,如数字信号处理、计算大整数乘法、求解偏微分方程,等等。

对于复数序列{x(n)}(n=0,1,2,...,N-1),则该复数序列的离散傅里叶变换为[1]:

X(K)= x(n)W,K=0,1,2,...N-1

其中W =e 。

直接按照上述变换公式的计算复杂度是O(n );快速傅里叶变换(FFT)可以计算出与直接计算相同的结果,但只需要计算复杂度O(nlog(n))。而实现FFT的算法有AMD Core Math Library (ACML),intel mkl,cwplib,dfftpack,kissfft,mixfft,monnier,jmffte,numutils,FFTW3等,FFTW以其在各个平台上的优秀表现,应用相当的广泛。

FFTW是Fastest Fourier Transform in the West(西方快速Fourier变换)的缩写。它是计算离散Fourier变换(DFT)的快速C程序的一个完整集合,它可计算一维或多维、实和复数据以及任意规模的DFT。FFTW还包含对共享和分布式存储系统的并行变换。FFTW由麻省理工学院计算机科学实验室超级计算技术组开发。

二、FFTW的工作原理

FFTW在各种平台上都能实现相当高的性能,主是由于以下两个FFTW工作特点[2]。

1.将长度为N的序列分解成长度为N 和N 的短序列进行计算。这些短序列的计算由一系列经过高度优化的C代码片段来完成,这些代码片段称为solvers。solvers由planner来管理与组合起来完成计算,其中一种组合称为一个plan。

2.对于同一个长度为N,其分解方式有多种,不同的组合路径对于运算的速度影响不同,在不同的机器与平台下其表现也会有所不同,因此需要确定哪种组合的效果最佳。FFTW通过对各种组合的评估,选择出其中执行速度最快的一个plan。将此plan记录下来,以后每次做同样长度的变换时,只要使用相当的plan即可,这样使FFTW能适应各种运行环境下的最优。

FFTW的核心工作模块主要由planner,solver和plan来组成。其中solver是解决变换的一种独立的方法;plan是由solver生成的包含计算变换的路径和所有输入输出数组的指针的一个数据结构,可以供后面的算法使用;planner则是管理solver生成plan,并通过一定的算法找出执行速度最快的plan,并将其返回给用户。具体的执行情况如下:

首先planner 对一系列的solver进行初始化,然后对给出一个指定长度变换,planner依次调用每个solver来生成plan,当此solver可以解决的话就返回一个plan指针,不能解决时则生成一个空指针。Planner对生成的各个plan进行执行并测量,选择执行速度最快的plan返回给用户。而用户则调用返回的plan来执行傅里叶变换。可以将此包含执行路径及组合算法的plan保存下来,重复使用。

由于FFTW的 planner是一些独立的解决特定问题的plan的模块。故同样的planner可以用于复数的DFT,实数的DFT及其它的变换。

三、FFTW的使用方法

在使用FFTW之前应该配置好其开发环境,本文以VC6.0为例来讲述其环境的配置方法。其步骤如下:

1.下载FFTW for windwos版本(ftp://ftp.省略/pub/fftw/fftw-3.2-dll.zip)。

2.解压fftw-3.2-dll.zip到一指定目录,如c:\fftw3。在此目录下有编译所需要的动态链接库,库的定义文件等,但缺少静态链接的lib文件,可利用def文件和VC的命令来生成。

3.使用VC下的命令lib来生成VC静态链接时需要的lib文件。如下命令将生成libfftw3-3.lib文件。

lib/def:libfftw3-3.def

4.在VC的包含文件路径中加入fftw3头文件的路径,如c:\fftw3。

5.在新建的工程中添加上生成的libfftw3-3.lib。

6.然后在工程就可以使用FFTW库。

FFTW的使用过程要使用到的数据结构为fftw_complex 。它的定义如下:

typedef double fftw_complex;

其是一个复数,其数据的存储格式为:部分存储实数部分,部分存储虚数部分。此数据结构由fftw_malloc和fftw_free函数分配与释放。其函数原型为:

#include

void *fftw_malloc(size_t n);

void fftw_free(void *p);

在使用过程中使用到的函数有fftw_plan_dft_*(*号代表各种方式的变换)和fftw_execute函数。其中fftw_plan_dft_*可实现生成一维及多维的及时变换和复杂变换的plan。fftw_execute则根据生成的plan来执行变换。其函数原型为:

void fftw_execute(const fftw_plan plan);

FFTW函数的使用步骤如下:

首先,使用fttw_malloc分配变换中要使用的输入输出数组。

第二步,利用它生成一个plan,用来包含FFTW计算FFT是要包含的所有的数据。

当输入输出数组使用同一个数组时,将执行原位运算。

第三步,利用先前生成的plan作为参数,调用fftw_execute执行傅里叶变换。当需要多次执行时,可以重复调用fftw_execute来完成任务。

最后,当不使用时,可以调用fftw_destroy_plan,fftw_free函数将相应的plan和输入输出数组销毁。

具体的使用实例如下所示:

/* 定义输入输出数组及fftw_plan*/

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);

/*调用planner生成最优plan */

p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);

/*给输入数组in添加数据*/

...

/*执行变换*/

fftw_execute(p); /* 根据需要可反复调用 */

/*执行完后在输出数组out中取出变换后的输出数据 */

...

/*不使用该plan,将其和输入输出数组销毁 */

fftw_destroy_plan(p);

fftw_free(in); fftw_free(out);

四、FFTW的性能评测

在2.2 GHz 32位Dual Core AMD Opteron Processor 275处理器上,分别对AMD Core Math Library (ACML),intel mkl,cwplib,dfftpack,kissfft,mixfft,monnier,jmffte,numutils,FFTW3对双精度一维的不同长度的DFT进行测试。其测试的速度如图1所示[3]。

从上图可以看出其性能明显要优于其它的FFT算法的执行速度。在实践中使用FFTW将会显著缩短算法的执行时间,提高整个系统的实时性能。

本人在做地震实时相关处理项目中使用FFTW做相关运算的处理,取得了显著的性能提升。由此可见在实践中使用FFTW实现FFT是切实可行的,并且能够高效地执行。

参考文献:

[1]刘益成,孙祥娥.数字信号处理.北京:电子工业出版社,2004.

[2]王刚.遥感图像快速傅立叶变换新算法的研究与实现.中科院遥感应用研究所,2003.

[3]FFTW User Manual.http://www.省略.

注:“本文中所涉及到的图表、注解、公式等内容请以PDF格式阅读原文。”

本文为全文原貌 未安装PDF浏览器用户请先下载安装 原版全文

上一篇:信息技术与初中物理教学的有效整合 下一篇:Protel设计电路中元件属性的相关性研究