热门问题
时间线
聊天
视角
快速傅里叶变换
快速計算序列的離散傅立葉變換或其逆變換的方法 来自维基百科,自由的百科全书
Remove ads
Remove ads
快速傅里叶变换(英语:Fast Fourier Transform, FFT),是快速计算序列的离散傅里叶变换(DFT)或其逆变换的方法[1]。傅里叶分析将信号从原始域(通常是时间或空间)转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏(大多为零)因子之积来快速计算此类变换。[2] 因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的 ,降低到 ,其中 为数据大小。



快速傅里叶变换广泛的应用于工程、科学和数学领域。这里的基本思想在1965年才得到普及,但早在1805年就已推导出来。[3] 1994年美国数学家吉尔伯特·斯特朗把FFT描述为“我们一生中最重要的数值算法”[4],它还被IEEE科学与工程计算期刊列入20世纪十大算法。[5]
Remove ads
定义和速度
用FFT计算DFT会得到与直接用DFT定义计算相同的结果;最重要的区别是FFT更快。(由于舍入误差的存在,许多FFT算法还会比直接运用定义求值精确很多,后面会讨论到这一点。)
令 x0, ...., xN-1 为复数。DFT由下式定义
直接按这个定义求值需要次运算:Xk 共有 N 个输出,每个输出需要 N 项求和。直接使用DFT运算需使用N个复数乘法(4N个实数乘法)与N-1个复数加法(2N-2个实数加法),因此,计算使用DFT所有N点的值需要N2复数乘法与N2-N 个复数加法。FFT则是能够在次操作计算出相同结果的任何方法。更准确的说,所有已知的FFT算法都需要次运算(技术上O只标记上界),虽然还没有已知的证据证明更低的复杂度是不可能的。[6]
要说明FFT节省时间的方式,就得考虑复数相乘和相加的次数。直接计算DFT的值涉及到 N2 次复数相乘和 N(N−1) 次复数相加(可以通过削去琐碎运算(如乘以1)来节省次运算)。众所周知的基2库利-图基算法,N 为2的幂,可以只用 (N/2)log2(N) 次复数乘法(再次忽略乘以1的简化)和 Nlog2(N) 次加法就可以得到相同结果。在实际中,现代计算机通常的实际性能通常不受算术运算的速度和对复杂主体的分析主导[7],但是从到的总体改进仍然能够体现出来。
Remove ads
一般的简化理论
假设一个M*N型矩阵S可分解成列向量以及行向量相乘:
若有个相异的非平凡值( where )
有个相异的非平凡值
则S共需要个乘法。
Step 1:
Step 2:
简化理论的变型:
也是一个M*N的矩阵。
若有个值不等于0,则的乘法量上限为。
Remove ads
快速傅里叶变换乘法量的计算
假设,其中彼此互素
点DFT的乘法量为,则点DFT的乘法量为:
假设,P是一个素数。
若点的DFT需要的乘法量为
且当中()
有个值不为及的倍数,
有个值为及的倍数,但不为的倍数,
则N点DFT的乘法量为:
Remove ads
库利-图基算法
库利-图基算法是最常见的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算法联合使用。
Remove ads
下面,我们用N次单位根来表示。
的性质:
- 周期性,具有周期,即
- 对称性:。
- 若是的约数,
为了简单起见,我们下面设待变换序列长度。根据上面单位根的对称性,求级数时,可以将求和区间分为两部分:
和是两个分别关于序列奇数号和偶数号序列N/2点变换。由此式只需计算出的前N/2个点,对于后N/2个点可以通过复指数函数的对称性快速求解。即,注意和都是周期为N/2的函数,由单位根的对称性,有以下变换公式:
- 。
这样,一个N点变换就分解成了两个N/2点变换。照这样可继续分解下去。这就是库利-图基快速傅里叶变换算法的基本原理。根据主定理不难分析出此时算法的时间复杂度为
Remove ads
- 蝶形结网络和位反转(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个点列作为输入传送到蝶形结网络中,注意将因子逐层加入到蝶形网络中。
Remove ads
由于按蝶形结网络计算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)。
Remove ads
实数或对称资料专用的算法
在许多的运用当中,要进行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可用下列的算法代替:
- DCT
- DST
- DHT
- 直接做折积(direct computing)
- 分段式DFT折积(sectioned DFT convolution)
- 威诺格拉德快速傅里叶变换算法
- 沃尔什(Walsh、Hadamard)转换
- 数论转换
单一路径延迟回授系统
单一路径延迟回授系统(英语: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的运算和中间结果的存储位置。

首先考虑系统输入为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)和寄存器使用量。
Remove ads
- Radix-2 SDF: 2(log4N-1) 个乘法器、4log4N 个加法器、存储器大小为 N-1。
- Radix-4 SDF: log4N-1 个乘法器、8log4N 个加法器、存储器大小为 N-1。
- Radix-22 SDF: log4N-1 个乘法器、4log4N 个加法器、存储器大小为 N-1。
- Radix-2 MDC: 2(log4N-1) 个乘法器、4log4N 个加法器、存储器大小为 1.5N-2。
- Radix-4 MDC: 3(log4N-1) 个乘法器、8log4N 个加法器、存储器大小为 2.5N-4。
- 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)。
快速算法查表
For a N-point DFT:
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 |
其他N的算法:
(1)若N可被两个以上互素的数相乘得到
设,且与互素(无须与为素数):
则当输入信号长度为N所需之的乘法器=
其中为-point DFT乘法量,为-point DFT乘法量。
例如:可以拆成,其中5与13互素,
借由查表可得知5-points DFT乘法量为10,13-points DFT乘法量为52,
所以13-points DFT乘法量为 。
扩展:
设,且互素,
-point DFT乘法量为,
则N-point DFT可分解为个-point DFT+个-point DFT++个-point DFT,
总乘法量为
(2)若N不可被两个互素的数相乘得到
设,且与不互素:
且当中(,),
有个值不为及的倍数,
有个值为或的倍数,但不为的倍数。
则N-point DFT所需之的乘法量=
其中为-point DFT乘法量,为-point DFT乘法量。
例如:16-point DFT,
乘法量=
or
乘法量=
FFT的应用
快速傅里叶变换可将卷积化乘法并以对数级数降低运算量,已成数字信号与科学计算核心工具。
其应用涵盖:
MP3、AAC 与 JPEG 均把信号分块后先作 FFT/MDCT,以把时域样本转到频域;高能量系数集中于低频,小幅子载波便能重现重要消息。接着结合人耳/视觉掩蔽模型及 Zig-zag 扫描,将低感知值系数量化为零,大幅降低比特率。最后用霍夫曼或算术编码存储系数与量化步阶,达到有损压缩与可调画质/音质的目的。若不经 FFT,逐点卷积将耗费 O(N²) 运算,无法在即时设备上完成。
现代蜂窝网络与 Wi-Fi 系统皆采正交分频多工 (OFDM)。基站先将并发调变符号经 IFFT 合成等间隔子载波;终端则以 FFT 还原频域资料。此架构可在多径通道中保持子载波互不干扰,同时透过可变子载波间距支撑 5G NR 的 eMBB 与 mMTC 场景。运用循环前缀可把时域卷积化为频域乘法,等效频率选择通道简化为每载波一次复数乘法,计算量对比传统单载波大幅减少。
FFT 让示波器、功率分析仪与机房电力质量监控能在毫秒级显示功率频谱。将长度 N 的窗函数样本取 FFT 后,平方其复数幅度并适当标准化,即得单/双边功率谱密度 (PSD)。工程师可由此快速识别谐波、噪声底或机械转速边带;科学领域亦用来分析地震波与星体光谱。无须昂贵的扫频式分析器,即可精确估测频率分量与相位。
空间滤波、去模糊及模板匹配常利用“乘频域、回时域”策略:对影像与滤波器 PSF 取 FFT,相乘后作 IFFT,即得卷积结果,计算量由 O(N²) 降为 O(N log N)。在盲去模糊或维纳解卷积中,亦可用频域除法抑制模糊并保留高频细节。此方法支持任意遮罩大小,远优于滑动视窗直接卷积,适用医学影像与天文摄影等高清晰度场景。
雷达为兼顾高清晰与大探测距离,常发射长码调变脉冲并于接收端做匹配滤波。将回波与已知码序列分别取 FFT,相乘后 IFFT,即可在频域完成相关;此“脉冲压缩”技术能在固定峰值功率下获得窄等效脉冲宽,提升距离分辨率与信杂比,并大幅降低即时硬件乘加数。
在多项式乘法或卷积为瓶颈的应用(例如格基同态加密、FFT-based 大数计算)中,先将系数以 NTT/FFT 转到频域,再做逐点乘法,最后逆变换还原,即可把 O(N²) 乘法降至 O(N log N)。加密函数库 SEAL、HElib 均内置此优化,使百万阶多项式运算在一般 CPU 上于毫秒内完成,直接左右现代隐私计算性能。
由此可见,FFT已由数学算法演变为支撑信息社会不可或缺的基础技术。
参阅
参考资料
延伸阅读
Wikiwand - on
Seamless Wikipedia browsing. On steroids.
Remove ads