跳转到内容

快速傅里叶变换

本页使用了标题或全文手工转换
维基百科,自由的百科全书
(差异) ←上一版本 | 最后版本 (差异) | 下一版本→ (差异)
一个五项余弦级数的时域信号。频率为 的倍数。振幅为 1, 2, 3, 4, 5。时间步长为 0.001 s
用FFT(~40,000次运算)五项余弦级数及用显式积分(~1亿次运算)得出的DFT。时间窗口是10秒。FFT是用FFTW3( http://www.fftw.org/页面存档备份,存于互联网档案馆) )计算的。显式积分DFT是用 https://sourceforge.net/projects/amoreaccuratefouriertransform/页面存档备份,存于互联网档案馆) 计算的。
缩放后的五项余弦级数的DFT。注意到显式积分更细的步长大小比FFT更精确地再现了峰值(4)和频率(56.569 Hz),代价是速度慢了上千倍。

快速傅里叶变换(英语:Fast Fourier Transform, FFT),是快速计算序列的离散傅里叶变换(DFT)或其逆变换的方法[1]傅里叶分析将信号从原始域(通常是时间或空间)转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解稀疏(大多为零)因子之积来快速计算此类变换。[2] 因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的 ,降低到 ,其中 为数据大小。

快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及,但早在1805年就已推导出来。[3] 1994年美国数学家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法[4],它还被IEEE科学与工程计算期刊列入20世纪十大算法。[5]

定义和速度

[编辑]

用FFT计算DFT会得到与直接用DFT定义计算相同的结果;最重要的区别是FFT更快。(由于舍入误差的存在,许多FFT算法还会比直接运用定义求值精确很多,后面会讨论到这一点。)

x0, ...., xN-1复数。DFT由下式定义

直接按这个定义求值需要 O(N2) 次运算:Xk 共有 N 个输出,每个输出需要 N 项求和。直接使用DFT运算需使用N个复数乘法(4N 个实数乘法)与N-1个复数加法(2N-2个实数加法),因此,计算使用DFT所有N点的值需要N2复数乘法与N2-N 个复数加法。FFT则是能够在 O(N log N) 次操作计算出相同结果的任何方法。更准确的说,所有已知的FFT算法都需要 O(N log N) 次运算(技术上O只标记上界),虽然还没有已知的证据证明更低的复杂度是不可能的。[6]

要说明FFT节省时间的方式,就得考虑复数相乘和相加的次数。直接计算DFT的值涉及到 N2 次复数相乘和 N(N−1) 次复数相加(可以通过削去琐碎运算(如乘以1)来节省 O(N) 次运算)。众所周知的基2库利-图基算法N 为2的幂,可以只用 (N/2)log2(N) 次复数乘法(再次忽略乘以1的简化)和 Nlog2(N) 次加法就可以得到相同结果。在实际中,现代计算机通常的实际性能通常不受算术运算的速度和对复杂主体的分析主导[7],但是从 O(N2) 到 O(N log N) 的总体改进仍然能够体现出来。

一般的简化理论

[编辑]

假设一个M*N型矩阵S可分解成列向量以及行向量相乘:

个相异的非平凡值( where ) 

个相异的非平凡值

S共需要个乘法。

Step 1:

Step 2:

简化理论的变型:

也是一个M*N的矩阵。

个值不等于0,则的乘法量上限为

快速傅立叶变换乘法量的计算

[编辑]

假设,其中彼此互质

点DFT的乘法量为,则点DFT的乘法量为:

假设,P是一个质数。

点的DFT需要的乘法量为

当中(

个值不为的倍数,

个值为的倍数,但不为的倍数,

则N点DFT的乘法量为:

库利-图基算法

[编辑]

库利-图基算法是最常见的FFT算法。这一方法以分治法为策略递归地将长度为离散傅里叶变换分解为长度为个较短序列的离散傅里叶变换,以及与旋转因子的复数乘法。

这种方法以及FFT的基本思路在1965年J. W. Cooley和J. W. Tukey合作发表An algorithm for the machine calculation of complex Fourier series之后开始为人所知。但后来发现,实际上这两位作者只是重新发明了高斯在1805年就已经提出的算法(此算法在历史上数次以各种形式被再次提出)。

库利-图基算法最有名的应用,是将序列长为N 的DFT分割为两个长为N/2 的子序列的DFT,因此这一应用只适用于序列长度为2的幂的DFT计算,即基2-FFT。实际上,如同高斯和Cooley与Tukey都指出的那样,Cooley-Tukey算法也可以用于序列长度N 为任意因数分解形式的DFT,即混合基FFT,而且还可以应用于其他诸如分裂基FFT等变种。尽管Cooley-Tukey算法的基本思路是采用递归的方法进行计算,大多数传统的算法实现都将显式的递归算法改写为非递归的形式。另外,因为Cooley-Tukey算法是将DFT分解为较小长度的多个DFT,因此它可以同任一种其他的DFT算法联合使用。

设计思想

[编辑]

下面,我们用N次单位根来表示

的性质:

  1. 周期性具有周期,即
  2. 对称性
  3. 的约数,

为了简单起见,我们下面设待变换序列长度。根据上面单位根的对称性,求级数时,可以将求和区间分为两部分:

是两个分别关于序列奇数号和偶数号序列N/2点变换。由此式只需计算出的前N/2个点,对于后N/2个点可以通过复指数函数的对称性快速求解。即,注意都是周期为N/2的函数,由单位根的对称性,有以下变换公式:

这样,一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据主定理不难分析出此时算法的时间复杂度为

算法实现

[编辑]
  • 蝶形结网络和位反转(Bit Reversal):
首先将个输入点列按二进制进行编号,然后对各个编号按位倒置并按此重新排序。例如,对于一个8点变换,
001倒置以后变成 100
000 → 000
001 → 100
010 → 010
011 → 110
100 → 001
101 → 101
110 → 011
111 → 111
倒置后的编号为{0,4,2,6,1,5,3,7}。
然后将这n个点列作为输入传送到蝶形结网络中,注意将因子逐层加入到蝶形网络中。

算法复杂度

[编辑]

由于按蝶形结网络计算n点变换要进行log n层计算,每层计算n个点的变换,故算法的时间复杂度为

其他算法

[编辑]

Cooley-Tukey算法之外还有其他DFT的快速演算法。对于长度互质的序列,可以采用基于中国剩余定理互质因子算法将N长序列的DFT分解为两个子序列的DFT。与Cooley-Tukey算法不同的是,互质因子算法不需要旋转因子。

Rader-Brenner算法是类似于Cooley-Tukey算法,但是采用的旋转因子都是纯虚数,以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey所取代,与Rader-Brenner演算法相比,有一样多的乘法量,却有较少的加法量,且不牺牲数值的准确性。

Bruun以及QFT演算法是不断的把DFT分成许多较小的DFT运算。(Rader-Brenner以及QFT演算法是为了2的指数所设计的演算法,但依然可以适用在可分解的整数上。Bruun演算法则可以运用在可被分成偶数个运算的数字)。尤其是Bruun演算法,把FFT看成是,并把它分解成的形式。

另一个从多项式观点的快速傅立叶变换法是Winograd算法。此演算法把分解成cyclotomic多项式,而这些多项式的系数通常为1,0,-1。这样只需要很少的乘法量(如果有需要的话),所以winograd是可以得到最少乘法量的快速傅立叶演算法,对于较小的数字,可以找出有效率的算方式。更精确地说,winograd演算法让DFT可以用点的DFT来简化,但减少乘法量的同时,也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。

Rader演算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋折积来表示原本的DFT,如此就可利用折积用一对基本的FFT来计算DFT。另一个prime-size的FFT演算法为chirp-Z演算法。此法也是将DFT用折积来表示,此法与Rader演算法相比,能运用在更一般的转换上,其转换的基础为Z转换(Rabiner et al., 1969)。

实数或对称资料专用的演算法

[编辑]

在许多的运用当中,要进行DFT的资料是纯实数,如此一来经过DFT的结果会满足对称性:

=

而有一些演算法是专门为这种情形设计的(e.g. Sorensen, 1987)。另一些则是利用旧有的演算法(e.g. Cooley-Tukey),删去一些不必要的演算步骤,如此省下了记忆体的使用,也省下了时间。另一方面,也可以把一个偶数长度且纯实数的DFT,用长度为原本一半的复数型态DFT来表示(实数项为原本纯实数资料的偶数项,虚数项则为奇数项)。

在Matlab上用一次N点FFT计算两个N点实数信号x,y的FFT:

function [X,Y] = Real2FFT( x,y )

% N-point FFT is used for 2 real signals both with length N

N = length(x);

z = x+y*i ;

Z = fft(z);

X = 0.5*(Z+conj(Z(1+mod(0:-1:1-N,N))));

Y = 0.5*(Z-conj(Z(1+mod(0:-1:1-N,N))))/i;

end

一度人们认为,用离散哈特利转换(Discrete Hartley Transform)来处理纯实数的DFT会更有效率,但接着人们发现,对于同样点数的纯实数DFT,经过设计的FFT,可以比DHT省下更多的运算。Bruun演算法是第一个试着从减少实数输入的DFT运算量的演算法,但此法并没有成为人们普遍使用的方法。

对于实数输入,且输入为偶对称或奇对称的情形,可以更进一步的省下时间以及记忆体,此时DFT可以用离散余弦转换离散正弦转换来代替(Discrete cosine/sine transforms)。由于DCT/DST也可以设计出FFT的演算法,故在此种情形下,此方法取代了对DFT设计的FFT演算法。

DFT可以应用在频谱分析以及做折积的运算,而在此处,不同应用可以用不同的演算法来取代,列表如下:

用来做频谱分析的情况下,DFT可用下列的演算法代替:

  • DCT
  • DST
  • DHT
  • 正交基底的扩展(orthogonal basis expantion)包括正交多项式(orthogonal polynomials)以及CDMA.
  • Walsh(Hadamard)转换
  • Haar转换
  • 小波(wavelet)转换
  • 时频分布(time-frequency distribution)

用来做折积的情况下,DFT可用下列的演算法代替:

单一路径延迟回授系统

[编辑]

单一路径延迟回授系统(英语:Single-path delay feedback,缩写为SDF)是一种用于实现快速傅立叶转换(Fast Fourier Transform)中运算的管线式架构(pipeline architecture)。此架构的特点是资料从输入一直到最终的输出只会通过单一的路径,并使用蝶形结(butterfly diagram)做为资料的运算单元。经过蝶形结后的输出资料会先暂存在系统的移位暂存器(shift registers)或是随机存取记忆体(RAM)当中,而由于只有单一路径的关系,对于每一级的蝶形结架构我们只需要一个复数乘法器来进行和FFT当中的旋转因子(twiddle factor)的乘积运算。[8]

SDF的优点在于其架构简单并且利用硬体实现的成本较低,原因在于其单一路径的设计,使得资料会依序进行处理并输出,并不会在系统中加入多个复杂的运算单元以进行平行化的运算,事实上SDF的架构在管线式快速傅立叶转换处理器(pipelined FFT processors)中具有最高的记忆体使用效率,而由于记忆体单元数量会随快速傅立叶转换的级数而指数增长,记忆体使用效率进而成为影响电路面积和功耗的主要因素之一。因此在实作大型快速傅立叶转换处理器的情况下,SDF架构一直是被优先考量采用的选择。然而其缺点也同样来自于单一路径的设计,导致其吞吐量(throughput)和多路径延迟交换线路系统(英语:Multiple-path delay commutator,缩写为MDC)相比较低。

接著我们将以8个点的基-2单一路径延迟回授系统(8-point R2SDF)为例进行其演算法细节的介绍,并且展示在硬体实作中每一个cycle的运算和中间结果的存储位置。

8-point R2SDF Architecture

首先考虑系统输入为x0, x1, ..., x7,并在前8个cycle依序输入。

  • Cycle 0~3: x0~x3依序存入第一级BF-2上方的暂存器中。
  • Cycle 4: x0和x4输入第一级BF-2后输出x0'和x4',x0'存入第二级BF-2上方的暂存器中,x4'存入第一级BF-2上方的暂存器中。
  • Cycle 5: x1和x5输入第一级BF-2后输出x1'和x5',x1'存入第二级BF-2上方的暂存器中,x5'存入第一级BF-2上方的暂存器中。
  • Cycle 6: x2和x6输入第一级BF-2后输出x2'和x6',x6'存入第一级BF-2上方的暂存器中,x0'和x2'输入第二级BF-2后输出x0''和x2'',x0''存入第三级BF-2上方的暂存器中,x2''存入第二级BF-2上方的暂存器中。
  • Cycle 7: x3和x7输入第一级BF-2后输出x3'和x7',x7'存入第一级BF-2上方的暂存器中,x1'和x3'输入第二级BF-2后输出x1''和x3'',x3''存入第二级BF-2上方的暂存器中,x0''和x1''输入第三级BF-2后输出x0'''和x1''',x0'''直接输出,x1'''存入第三级BF-2上方的暂存器中。(此为第一笔输出,此后每一个cycle皆会输出一笔资料。)
  • Cycle 8: x4'乘上W0后存入第二级BF-2上方的暂存器中,x2''乘上W0后存入第三级BF-2上方的暂存器中,x1'''直接输出。
  • Cycle 9: x5'乘上W1后存入第二级BF-2上方的暂存器中,x3''乘上W2后和x2''输入第三级BF-2后输出x2'''和x3''', x2'''直接输出,x3'''存入第三级BF-2上方的暂存器中。
  • Cycle 10: x6'乘上W2后和x4'输入第二级BF-2后输出x4''和x6'',x6''存入第二级BF-2上方的暂存器中,x4''存入第三级BF-2上方的暂存器中,x3'''直接输出。
  • Cycle 11: x7'乘上W3后和x5'输入第二级BF-2后输出x5''和x7'',x7''存入第二级BF-2上方的暂存器中,x4''和x5''输入第三级BF-2后输出x4'''和x5''',x4'''直接输出,x5'''存入第三级BF-2上方的暂存器中。
  • Cycle 12: x6''乘上W0后存入第三级BF-2上方的暂存器中,x5'''直接输出。
  • Cycle 13: x7''乘上W2后和x6''输入第三级BF-2后输出x6'''和x7''',x6'''直接输出,x7'''存入第三级BF-2上方的暂存器中。
  • Cycle 14: x7'''直接输出。(此为最后一笔输出。)

此外,我们可以在蝶形结中的加法器和乘法器加上额外的暂存器(register)来降低系统的关键路径(critical path),来避免某些时刻不同级之间的繁重运算必须挤在同一个cycle中运算完成,进而严重影响整体电路在timing上的表现。但如此的管线式设计(pipeline design)也会增加整体电路的延迟(latency)和暂存器使用量。

相似架构之硬体实作比较

[编辑]
  1. Radix-2 SDF: 2(log4N-1) 个乘法器、4log4N 个加法器、记忆体大小为 N-1。
  2. Radix-4 SDF: log4N-1 个乘法器、8log4N 个加法器、记忆体大小为 N-1。
  3. Radix-22 SDF: log4N-1 个乘法器、4log4N 个加法器、记忆体大小为 N-1。
  4. Radix-2 MDC: 2(log4N-1) 个乘法器、4log4N 个加法器、记忆体大小为 1.5N-2。
  5. Radix-4 MDC: 3(log4N-1) 个乘法器、8log4N 个加法器、记忆体大小为 2.5N-4。
  6. Radix-4 SDC: log4N-1 个乘法器、3log4N 个加法器、记忆体大小为 2N-2。

复杂度以及运算量的极限

[编辑]

长久以来,人们对于求出快速傅立叶变换的复杂度下限以及需要多少的运算量感到很有兴趣,而实际上也还有许多问题需要解决。即使是用较简单的情形,即点的DFT,也还没能够严谨的证明出FFT至少需要(不比小)的运算量,目前也没有发现复杂度更低的演算法。通常数学运算量的多寡会是运算效率好坏最主要的因素,但在现实中,有许多因素也会有很大的影响,如快取记忆体以及CPU均有很大的影响。

在1978年,Winograd率先导出一个较严谨的FFT所需乘法量的下限:。当时,DFT只需要次无理实数的乘法即可以计算出来。更详尽,且也能趋近此下限的演算法也一一被提出(Heideman & Burrus, 1986; Duhamel, 1990)。很可惜的是,这些演算法,都需要很大量的加法计算,目前的硬体无法克服这个问题。

对于所需加法量的数目,虽然我们可以在某些受限制的假设下,推得其下限,但目前并没有一个精确的下限被推导出来。1973年,Morgenstern在乘法常数趋近巨大的情形下(对大部分的FFT演算法为真,但不是全部)推导出加法量的下限:。Pan(1986)在假设FFT演算法的不同步的情形有其极限下证明出加法量的下限,但一般来说,此假设相当的不明确。长度为的情形下,在某些假设下,Papadimitriou(1979)提出使用Cooley-Tukey演算法所需的复数加法量是最少的。到目前为止,在长度为情况,还没有任何FFT的演算法可以让复数的加法量比还少。

还有一个问题是如何把乘法量与加法量的总和最小化,有时候称作"演算复杂度"(在这里考虑的是实际的运算量,而不是渐近复杂度)。同样的,没有一个严谨下限被证明出来。从1968年开始,点DFT而言,split-radix FFT演算法需要最少的运算量,在的情形下,其需要个乘法运算以及加法运算。最近有人导出更低的运算量:。(Johnson and Frigo, 2007; Lundy and Van Buskirk, 2007)

大多数尝试要降低或者证明FFT复杂度下限的人都把焦点放在复数资料输入的情况,因其为最简单的情形。但是,复数资料输入的FFT演算法,与实数资料输入的FFT演算法,离散馀弦转换(DCT),离散哈特列转换(DHT),以及其他的演算法,均有很大的关连性。故任何一个演算法,在复杂度上有任何的改善的话,其他的演算法复杂度也会马上获得改善(Duhamel & Vetterli, 1990)。

快速演算法查表

[编辑]

当输入信号长度为N时:

N = 1~60

N 乘法数 加法数 N 乘法数 N 乘法数 N 乘法数
1 0 0 11 46 24 28 39 182
2 0 4 12 8 25 148 40 100
3 2 12 13 52 26 104 42 124
4 0 16 14 32 27 114 44 160
5 10 34 15 40 28 64 45 170
6 4 36 16 20 30 80 48 92
7 16 72 18 32 32 72 52 208
8 4 52 20 40 33 160 54 228
9 16 72 21 62 35 150 56 156
10 20 88 22 80 36 64 60 160

N < 1000

N 乘法数 N 乘法数 N 乘法数 N 乘法数
63 256 96 280 192 752 360 1540
64 204 104 468 204 976 420 2080
66 284 108 456 216 1020 480 2360
70 300 112 396 224 1016 504 2300
72 164 120 380 240 940 512 3180
80 260 128 560 252 1024 560 3100
81 480 144 436 256 1308 672 3496
84 248 160 680 288 1160 720 3620
88 412 168 580 312 1324 784 4412
90 340 180 680 336 1412 840 4580

N > 1000

N 乘法数 N 乘法数 N 乘法数 N 乘法数
1008 5356 1440 8680 2520 16540 4032 29488
1024 7436 1680 10420 2688 19108 4096 37516
1152 7088 2016 12728 2880 20060 4368 35828
1260 7640 2048 16836 3369 24200 4608 36812
1344 8252 2304 15868 3920 29900 5040 36860

参阅

[编辑]

参考资料

[编辑]
  1. ^ 杨毅明. 数字信号处理(第2版). 北京: 机械工业出版社. 2017年: 第95页. ISBN 9787111576235. 
  2. ^ Charles Van Loan, Computational Frameworks for the Fast Fourier Transform (SIAM, 1992).
  3. ^ Heideman, M. T.; Johnson, D. H.; Burrus, C. S. Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine. 1984, 1 (4): 14–21. doi:10.1109/MASSP.1984.1162257. 
  4. ^ Strang, Gilbert. Wavelets. American Scientist. May–June 1994, 82 (3): 253. JSTOR 29775194. 
  5. ^ Dongarra, J.; Sullivan, F. Guest Editors Introduction to the top 10 algorithms. Computing in Science Engineering. January 2000, 2 (1): 22–23. ISSN 1521-9615. doi:10.1109/MCISE.2000.814652. 
  6. ^ Johnson and Frigo, 2007
  7. ^ Frigo & Johnson, 2005
  8. ^ Shousheng He; Torkelson, M. A new approach to pipeline FFT processor. IEEE Comput. Soc. Press. 1996. ISBN 978-0-8186-7255-2. doi:10.1109/IPPS.1996.508145. 

延伸阅读

[编辑]