C 中的简单有效的FFT实现，第一部分外文翻译资料

C 中的简单有效的FFT实现，第一部分

Vlodymyr在德国科特布兰的勃兰登堡理工大学任教。 他可以myrnyy@math.tu-cottbus.de。 本文由Dobb博士的杂志提供。

本文介绍了使用C 模板元编程的Cooley-Tukey快速傅里叶变换（FFT）算法的新的有效实现。 由于FFT的递归性质，源代码比传统的实现更加可读和更快变换的效率通过不同平台上的性能基准来证明。

介绍

快速傅里叶变换（FFT）不仅是一种快速的计算数字傅里叶变换（DFT）的方法，它的复杂度O（Nlog（N））（N必须是2的幂，N = 2^P），它还是是一种线性使用性质“分化与征服”的解决非线性复杂性的许多真实数学问题的方法。

将N点信号xn的离散傅立叶变换fn定义为和：

示例一：

那是个复杂的集合。 简单来说，根据公式用于计算转换的算法将需要O（N₂）操作。 但是Danielson-Lanczos引理（1942）使用了统一g复杂的属性，给出了一个奇妙的构想，即递归构造傅里叶变换（例1）。

示例二：

其中f_k^e表示由原始xj的偶数分量形成的长度为N / 2的傅里叶变换的第k个分量，而f_k^o是由奇数分量形成的相应变换。 虽然实施例2的最后一行的k从0变化到N-1，但是fk^e和fk^o的变换在长度为N / 2的k中是周期性的。 应用于变换的相同公式fk^e和fk^e将问题减少到长度N / 4的变换等等。 这意味着，如果N是2的幂，则将以线性复杂度O（Nlog（N））计算变换。 关于FFT和高级算法的数学背景的更多信息，不限于N = 2^P，可以在许多书中找到，例如，[3,4]。

我想从算法的最简单的递归形式开始，直接来自于示例2中的关系：

FFT(x) {

n=length(x);

if (n==1) return x;

m = n/2;

X = (x_{2j})_{j=0}^{m-1};

Y = (x_{2j 1})_{j=0}^{m-1};

X = FFT(X);

Y = FFT(Y);

U = (X_{k mod m})_{k=0}^{n-1};

V = (g^{-k}Y_{k mod m})_{k=0}^{n-1};

return U V;

}

该算法只能给出FFT构造的第一印象。在源数据的偶数和奇数元素上，FFT（x）函数被递归地调用两次。之后，对数据进行一些转换。如果数据长度变为1，则递归结束。这种递归形式是有启发性的，但是绝大多数FFT实现使用了由Cooley和Tukey [2]在1965年首先实现的循环结构。Cooley-Tukey算法使用如下事实：原始长度N信号x的元素被赋予一定的“位加扰”排列，则可以使用方便的嵌套循环来执行FFT。加扰是反二进制重新索引，这意味着x_j被x_k替换，其中k是j的反二进制表示。例如，对于数据长度N = 2⁵，元素x5必须与x _{20}交换，因为5 = 001012的二进制反转为101002 = 20。稍后将考虑该数据排列的实现，因为它已经是整个FFT的一小部分。最重要的观察是，Cooley-Tukey方案实际上允许在现场执行FFT，也就是说，原始数据x被逐个元素地替换为FFT值。这是一种非常高效的记忆效率的方式来处理大量数据。清单一显示了C [5]，第53页中数值配方的Cooley-Tukey算法的经典实现。

void four1(double* data, unsigned long nn)

{

unsigned long n, mmax, m, j, istep, i;

double wtemp, wr, wpr, wpi, wi, theta;

double tempr, tempi;

// reverse-binary reindexing

n = nnlt;lt;1;

j=1;

for (i=1; ilt;n; i =2) {

if (jgt;i) {

swap(data[j-1], data[i-1]);

swap(data[j], data[i]);

}

m = nn;

while (mgt;=2 amp;amp; jgt;m) {

j -= m;

m gt;gt;= 1;

}

j = m;

};

// here begins the Danielson-Lanczos section

mmax=2;

while (ngt;mmax) {

istep = mmaxlt;lt;1;

theta = -(2*M_PI/mmax);

wtemp = sin(0.5*theta);

wpr = -2.0*wtemp*wtemp;

wpi = sin(theta);

wr = 1.0;

wi = 0.0;

for (m=1; m lt; mmax; m = 2) {

for (i=m; i lt;= n; i = istep) {

j=i mmax;

tempr = wr*data[j-1] - wi*data[j];

tempi = wr * data[j] wi*data[j-1];

data[j-1] = data[i-1] - tempr;

data[j] = data[i] - tempi;

data[i-1] = tempr;

data[i] = tempi;

}

wtemp=wr;

wr = wr*wpr - wi*wpi;

wi = wi*wpr wtemp*wpi;

}

mmax=istep;

}

}

清单一

初始信号存储在长度为2 * nn的阵列数据中，其中每个偶数元素对应于实数部分，每个奇数元素对应于复数的虚部。

递归是很好的性质

FFT实现的大多数已知方法都是基于避免自然FFT递归，将其替换为循环。 但是，如果在编译时解决了递归，那么递归不再是问题，就像模板类递归一样。 此外，这种递归可以提供性能优势，因为代码比普通循环更好地展开。 这个想法似乎与Todd Veldhuizen [6]的方法非常相似，他在模板元程序中重写了相同的Cooley-Tukey算法（清单1）。 嵌套循环成为具有非线性复杂度的递归模板，可以在现代工作站中花费大量时间和内存来最多编译N =

2 ^{12}。 效率相当高，这一技术实现并不适用于实际的问题，因为它们经常需要处理更多的数据。 从这两个角度来看，我试图用模板的效率找到一个“黄金分割”。

他在这里介绍的方法利用了使用模板类递归来实现Danielson-Lanczos关系的FFT的原始递归性质（示例2）。 必要的假设是信号N = 2^P的长度是一个静态常数，并作为模板参数P传递。我将算法从清单1中分为两部分，即高频抽象级别：加扰和Danielson-Lanczos部分。 清单2表示包含成员函数fft（T * data）的初始模板类GFFT，包括变换的两个部分。

templatelt;unsigned P,

typename T=doublegt;

class GFFT {

enum { N = 1lt;lt;P };

DanielsonLanczoslt;N,Tgt; recursion;

public:

void fft(T* data) {

scramble(data,N);

recursion.apply(data);

}

};

清单二现在的主要目的是使用递归模板实现DanielsonLanczos模板类，其中P是2定义N的幂。类型T是数据元素的默认类型。稍后将简要介绍函数争用的执行情况，现在可以从清单1中得到。清单3中的DanielsonLanczos模板类取决于整数N，定义递归和相同类型T中数据的当前长度。为了避免实例化模板的非线性数，我只定义一个模板classDanielsonLanczos lt;N / 2，Tgt;每递归水平因此，要实例化的模板类的总数为P 1。由于物理内存限制，常数P不能很大。例如，如果您的数据具有双重精度的复杂元素（每个元素2个8字节），那么P在32位平台上可能在1到27之间变化。情况P = 28对应于4GB的数据，并且没有其他程序变量的存储器。 P可以在64位处理器上更大，但再次受到可用物理内存的限制。这样一些实例化的模板类不应该提供任何编译问题。 Danielson-Lanczos关系的递归思想是通过成员函数应用的两个递归调用实现的：第一次使用原始信号数据，第二次移位N，每个下一个递归级别将N除以2.最后一个是专门用于N = 1并包括空成员函数适用。

templatelt;unsigned N, typename T=doublegt;

class DanielsonLanczos {

DanielsonLanczoslt;N/2,Tgt; next;

public:

void apply(T* data) {

next.apply(data);

next.apply(data N);

T tempr,tempi,c,s;

for (unsigned i=0; ilt;N; i =2) {

c = cos(i*M_PI/N);

s = -sin(i*M_PI/N);

tempr = data[i N]*c - data[i N 1]*s;

tempi = data[i N]*s data[i N 1]*c;

data[i N] = data[i]-tempr;

data[i N 1] = data[i 1]-tempi;

data[i] = tempr;

data[i 1] = tempi;

}

}

};

templatelt;typename Tgt;

class DanielsonLanczoslt;1,Tgt; {

public:

void apply(T* data) { }

};

清单三

递归完成后，数据在循环中被修改，其中cos和sin函数用于计算单位（c，s）的复数根。 所得到的（tempr，tempi）是要修改的临时复数（data[i N]，data[i N 1]）和（data[i]，data[i 1]）。 清单3中的这种简单的实现由于大量的三角函数的计算而具有差的性能。

清单1中的原始实现包含一个不同的统一的属性：它们的递归计算。 清单四是具有相同递归公式的下一个实现步骤。 只计算两个正弦函数，从1开始，下一个根（wr，wi）从（wpr，wpi）：（wr，wi） =（wr，wi）*（wpr，wpi）反复计算。 Danielson Lanczostemplate课程的最终专业化水平保持不变。消除三角函数调用

清单1中的原始实现包含一个不同的统一的属性：它们的递归计算。 清单四是具有相同递归公式的下一个实现步骤。 只计算两个正弦函数，从1开始，下一个根（wr，wi）从（wpr，wpi）：（wr，wi） =（wr，wi）*（wpr，wpi）反复计算。 Danielson Lanczostemplate的阶数最终保持不变。

templatelt;unsigned N, typename T=doublegt;

class DanielsonLanczos {

DanielsonLanczoslt;N/2,Tgt; next;

public:

void apply(T* data) {

next.apply(data);

next.apply(data N);

T wtemp,tempr,tempi,wr,wi,wpr,wpi;

wtemp = sin(M_PI/N);

wpr = -2.0*wtemp*wtemp;

wpi = -sin(2*M_PI/N);

wr = 1.0;

wi = 0.0;

for (unsigned i=0; ilt;N; i =2) {

lt;

剩余内容已隐藏，支付完成后下载完整资料

---

资料编号：[141618]，资料为PDF文档或Word文档，PDF文档可免费转换为Word