快速傅里叶变换

时间:2020-05-06 23:37:15

一个五项余弦级数的时域信号。频率为
用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)或其逆变换的方法。傅里叶分析将信号从原始域(通常是时间或空间)转换到频域的表示或者逆过来转换。FFT会通过把DFT矩阵分解为稀疏(大多为零)因子之积来快速计算此类变换。 因此,它能够将计算DFT的复杂度从只用DFT定义计算需要的 ,其中

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

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

一般的简化理论

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

where 个乘法。

Step 2:

的乘法量上限为,其中点DFT的乘法量为点DFT的乘法量为:

,P是一个素数。

个值为的倍数,但不为

库利-图基算法

库利-图基算法是最常见的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次单位根

具有周期

  • 对称性
  • 为了简单起见,我们下面设待变换序列长度时,可以将求和区间分为两部分:

    奇数号和偶数号序列N/2点变换。由此式只能计算出

  • 算法实现

    • 蝶形结网络和位反转(Bit Reversal):
    首先将逐层加入到蝶形网络中。

    算法复杂度

    由于按蝶形结网络计算n点变换要进行log n层计算,每层计算n个点的变换,故算法的时间复杂度为互质的序列,可以采用基于剩余定理的互素因子算法将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看成是分解成cyclotomic多项式,而这些多项式的系数通常为1,0,-1。这样只需要很少的乘法量(如果有需要的话),所以winograd是可以得到最少乘法量的快速傅立叶算法,对于较小的数字,可以找出有效率的算方式。更精确地说,winograd算法让DFT可以用=点的DFT,也还没能够严谨的证明出FFT至少需要。当次无理实数的乘法即可以计算出来。更详尽,且也能趋近此下限的算法也一一被提出(Heideman & Burrus, 1986; Duhamel, 1990)。很可惜的是,这些算法,都需要很大量的加法计算,目前的硬件无法克服这个问题。

    对于所需加法量的数目,虽然我们可以在某些受限制的假设下,推得其下限,但目前并没有一个精确的下限被推导出来。1973年,Morgenstern在乘法常量趋近巨大的情形下(对大部分的FFT算法为真,但不是全部)推导出加法量的下限:,但一般来说,此假设相当的不明确。长度为是最少的。到目前为止,在长度为还少。

    还有一个问题是如何把乘法量与加法量的总和最小化,有时候称作"演算复杂度"(在这里考虑的是实际的运算量,而不是渐近复杂度)。同样的,没有一个严谨下限被证明出来。从1968年开始,的情形下,其需要。(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乘法数
    1001146242839182
    2041282514840100
    321213522610442124
    401614322711444160
    510341540286445170
    6436162030804892
    716721832327252208
    845220403316054228
    9167221623515056156
    1020882280366460160

    N < 1000

    N乘法数N乘法数N乘法数N乘法数
    63256962801927523601540
    642041044682049764202080
    6628410845621610204802360
    7030011239622410165042300
    721641203802409405123180
    8026012856025210245603100
    8148014443625613086723496
    8424816068028811607203620
    8841216858031213247844412
    9034018068033614128404580

    N > 1000

    N乘法数N乘法数N乘法数N乘法数
    1008535614408680252016540403229488
    10247436168010420268819108409637516
    11527088201612728288020060436835828
    12607640204816836336924200460836812
    13448252230415868392029900504036860

    参阅

    • 离散傅里叶变换
    • 并行快速傅里叶变换
  • 与本文近似的文章: