【翻译】Fast Fourier Transform

最新推荐文章于 2024-07-05 14:02:35 发布

erutan_pku

最新推荐文章于 2024-07-05 14:02:35 发布

阅读量2.7k

点赞数 2

分类专栏：课程笔记文章标签：算法分析论文翻译快速傅里叶变换

课程笔记专栏收录该内容

13 篇文章 0 订阅

订阅专栏

这是Stefan Wörner的一篇介绍快速傅里叶变换的论文，原文下载在这里有：http://download.csdn.net/detail/u013012544/7045283，是北京大学本科生算法设计与分析2014年春季课程的选读论文之一。

【这篇论文中的公式图表太多了...所以有很多部分都是截图...】

以下是它的全文翻译：

瑞士联邦理工学院

快速傅里叶变换

数值分析研讨会

Stefan Wörner

1、历史介绍

2、连续和离散的傅里叶变换

2.1连续傅里叶变换

2.2离散傅里叶变换

2.2.1三角插值

3、FFT的推导

3.1对合数的FFT

3.2对素数的FFT

4、实现

4.1递归实现

4.2迭代实现

参考书目

1、历史介绍

快速傅里叶变换的历史是非常有趣的。它可以追溯到1805年，卡尔·弗里德里希·高斯(Carl Friedrich Gauss)试图通过采样点位置来确定小行星轨道的时候[3]，他提出了离散傅里叶变换（DFT 见定义2.2），甚至比1822年发表相关论述的傅里叶还要早。为了计算DFT，他发明的算法同库里和图基(Cooley and Tukey) [2][3]的算法等价。不过，在有生之年，高斯并没有发表他的方法，就像有其他的什么更好的方法可以处理这个问题似的。这可能也就是为什么在1866年高斯文集发表的时候并没有人留意这篇手稿。直到160年之后，库里和图基才重新发现FFT。当时，美国军方正致力于研究方法来检测苏联的核试验，其中一种方法就是去分析地震的时间序列数据。正好图基所在的肯尼迪总统科学顾问委员会[6]负责处理该问题。在1805年到1965年之间，有很多科学家都给出了计算DFT的有效算法。但是它们的一般性都不如高斯或者库里和图基的算法。详见图1

2、连续和离散的傅里叶变换

本章给出了连续傅里叶变换(CFT)和离散傅里叶变换(DFT)的定义并以三角插值问题导出了离散傅里叶变换。

2.1连续傅里叶变换

定义2.1 (CFT) 记黎曼可积函数f从[0,L]映射到复数域中，满足f(0)= f(L) .定义f的第k个复傅里叶系数：

我们称矩阵W_N为傅里叶矩阵。这表明我们需要N²次复数乘法和N(N-1)次加法来协助完成这个变换。即需要O(n²)的算术复杂度。

2.2.1三角插值

本节以对一个函数的三角插值导出离散傅里叶变换。记F = (f₀ , f₁, … , f_N-1)∈C^N为函数f : [0,L] à C的一组等距采样点。CFT是一种利用无穷级数e^{2πi k/L x}, k∈Z来近似计算f的坐标的方式。现在，我们希望仅利用e^{2πik/L x}, k=0,1,…,N-1这N个方程来近似地表达f，即：

而接下来，我们要证明上式被定义为D(j,l)的函数等价于克罗内克函数δ_jl【注：克罗内克函数返回布尔量，当两个参数相等时返回1，否则返回0。】

定理2.1还表明，(2.5)式所定义的同时也是DFT的逆变换。由于DFT与其逆变换的式子几乎一样，我们的诸如FFT的那些高效计算DFT的算法同样也适用于计算DFT的逆。

3、FFT的推导

本章介绍了FFT的理论背景并讨论了几个应用广泛的特例。

3.1对合数的FFT

记N = N₁* N₂, 其中N₁, N₂∈ N。向量 (x₀ , x₁ , … , x_N-1)∈ C^N的DFT定义如下：

出于惯例，上式前面的因子1/n被省略。【傅里叶变换对前面的系数没有要求，只要与其逆变换的系数乘起来得1/n就是对的。】我们可以通过如下换元方法将一维的DFT求和转化为二维的：

j = j(a , b) = aN₁ + b; 0 ≤ a < N2, 0 ≤ b < N1 (3.2)

k = k(c , d) = cN₂ + d; 0 ≤ c < N1, 0 ≤ d < N2 (3.3)

在引用[6][7]中指出，上面的算法已经达到同类算法的下限，即不存在更优的算法可以处理这类问题。

3.2对素数的FFT

最后一节中的算法并不对于N是素数的情况适用。但是在引用[5]中给出了两种高效计算DFT的方法。其中第一种可以在N-1是一个因子很多的合数时候运用，在这种情况下，我们可以通过处理三个FFT算法来解决。第二种方法通过创建一个较大的数组N’=2ⁿ > N,然后通过一般的FFT算法来实现。我不打算提供更多细节因为这些算法使用了一些技巧通过FFT来处理一些循环相关的函数，而这已经超出了我们讨论的范围。如果想要了解更多的话可以去看引用[5]。总而言之，这些算法在N为素数的时候依然可以很快地计算出DFT的值。比如说上面提到的第二种算法的算术复杂度为O(N’log₂(N’))

4、实现

这一章中给出了两种FFT在N=2ⁿ时候的实现（这就是所谓的基2算法），并且分别讨论了这两种实现的利弊。为了更好地了解这两种算法是怎么实现的，我们可以看一下下面的式子。这是式(3.4)和(3.5)的特例。同样的，前面的1/N也被省略了。

算法1 FFT的递归实现

  function y = fft_rec(x)
  n = length(x);
  if n == 1
    y = x;
  else
    m = n/2;
    y_top = fft_rec(x(1:2:(n-1)));
    y_bottom = fft_rec(x(2:2:n));
    d = exp(-2 * pi * i / n) .^ (0:m-1);
    z = d .* y_bottom;
    y = [ y_top + z , y_top - z ];
  end

4.2 迭代实现

FFT的迭代实现并不像递归实现那样直白。(4.1)式中，数据按奇偶分成了两部分，如果继续运用该式，数据还会被继续分割下去。如图4.1展示的是N=8的时候，进行log₂N=3次分割的过程：

图4.1 迭代FFT时的数据排列[1]

这种方法也可以推广到N=2ⁿ的情况。这表明，如果数据如下排序，DFT可以如图4.2所示进行计算。其中，将两个元素进行归并的方式与算法1中的相同，F(2,0)表示(x₀,x₁,…,x₇)的DFT.

图4.2 FFT N = 8时的方案[4]

接下来需要解决的问题就是如何重排数据了，而这可以通过所谓的位反转轻易搞定。这也就是说k的索引序就是k的二进制表示然后按位进行一下反转后新得到的二进制值。如表4.1所示：

表4.1 位反转 [4]

定理4.1：按图4.1中所描述那种分割方式进行排列等价于位反转。

证明：

（这个证明只涉及N = 8的情况，不过这个证明是有启发性的并且显然地对N=2ⁿ也成立。）

我们如下表示数k的二进制表示及其位反转：

即为k的按位反转。

证毕。

Matlab提供了一个计算位反转的指令：” bitrevorder(x)”,于是我们得到如算法二所示的matlab代码：

算法2 FFT的迭代实现[4]

  function y = fft_it(x)
  n = length(x);
  x = x(bitrevorder(1:n));
  q = round(log(n)/log(2));
  for j = 1:q
    m = 2^(j-1);
    d = exp(-2 *pi * i /m).^(0:m-1);
    for k = 1:2^(q-j)
      s = (k-1)*2*m+1;               % 开始位置索引
      e = k*2*m;                     % 结束位置索引
      r = s + (e-s+1)/2;             % 中间位置索引
      y_top = x(s:(r-1));
      y_bottom = x(r:e);
      z = d .* y_bottom;
      y = [y_top + z, y_top - z];
      x(s:e) = y;
    end
  end

引用：

【略】