DSP 2006-2007
快速傅立叶算法分析
The College of Computer Science
Beijing University of Technology S200607097 杨涛
Email : [email protected]
数字信号处理
S200607097
杨涛
目录
一、傅立叶分析理论的提出 ................................................................................................... 3 二、傅立叶分析的类别及其联系 ........................................................................................... 3 （一）连续周期信号的傅立叶级数 ................................................................................... 3 （二）连续非周期信号的傅立叶变换 ............................................................................... 4 （三）离散周期信号的傅立叶级数 ................................................................................... 4 （四）离散非周期信号的傅立叶变换 ............................................................................... 4 （五）离散傅立叶变换 ....................................................................................................... 5 三、快速傅立叶变换的算法分析 ........................................................................................... 5 （一）快速傅立叶变换提出的原因 ................................................................................... 5 （二）快速傅立叶变换的原理 ........................................................................................... 6 （三）基 2 快速傅立叶算法分析 ...................................................................................... 6 四、快速傅立叶算法的实现 ................................................................................................... 6 （一）算法设计 ................................................................................................................... 6 （1）计算相位因子Wn .................................................................................................. 7 （2）蝶形运算 ................................................................................................................. 7 （3）重排信号序列 ......................................................................................................... 8 （二）结果分析 ................................................................................................................... 9 五、参考文献 ........................................................................................................................... 9
2
数字信号处理
S200607097
杨涛
一、傅立叶分析理论的提出
在物理学研究中， 科学家发现当白光通过玻璃棱镜可被分解为不同颜色的光， 每一种颜 色对应一种可见光谱的特定频率， 因此由光到颜色的分析就是一种频率分析。 将这种思想应 用到对信号性质的研究中，就引出了傅立叶分析理论。 将信号展开成傅立叶级数的基本目的就是要将时间变量 t 的函数分解为不同的频率分 量，这些基本的构造分量正弦和余弦函数。数学家已经证明，形如（1）式的三角函数的和 式就足以表示绝大多数曲线。 cos sin （1）
对周期函数进行这种分解，被称为傅立叶级数；对有限能量信号函数，这种分解被称为 傅立叶变换。三角函数展开起源于 18 世纪关于简谐振动以及类似的物理现象的研究。1808 年傅立叶在著作《热的分析理论》中详细地研究了三角函数，并利用三角函数解决了很多热 传导问题。在随后的二百年间，傅立叶分析理论得到不断的发展，并被应用于众多的领域， 例如用来消除信号中的高频噪声或进行数据压缩等等。
二、傅立叶分析的类别及其联系
根据信号的各种分类和相应特点， 我们可以将施加其上的傅立叶分析分为诸如连续信号 傅立叶级数和傅立叶变换、 离散信号傅立叶级数和傅立叶变换； 根据对信号的处理策略上的 改进，又由离散傅立叶变换（DFT）衍生了快速傅立叶变换（FFT）等等。在刚刚接触这些术 语和缩写时常常会感到混乱和迷惑。 为了更加清晰地认识其理论脉络， 下面将简要介绍它们 各自的应用范围，以及它们之间的联系与区别。
（一）连续周期信号的傅立叶级数
傅立叶级数是周期信号的基本数学表示， 它是相关正弦信号或负指数信号的线性加权和。 出于数学上分析和表示的方便，常根据欧拉公式（2）将傅立叶级数表示为复指数信号和的 信号。 e cos sin （2） Dirichlet 条件保证了除不连续时刻外，级数式（3）一定和初始信号 x(t)相吻合。 x t c c e x t e
F
3
F
TP TP
dt 4
其中系数 c 确定波形的形状，基频F 是给定周期TP 的倒数，确定基本周期，选取适当 的基频F 和系数 c 就可以构造不同类型的周期信号。
3
数字信号处理
S200607097
杨涛
（二）连续非周期信号的傅立叶变换
由上面的公式（3）我们可以知道，周期信号的频率表示是离散的图形，谱线之间有相 等的间隔，这些间隔等于基本频率。如果信号的周期趋近于无限，信号变为非周期的，谱线 间隔将趋近 0，因此信号的频率表示变成连续函数。非周期信号通过周期为TP 的重复得到周 期信号， 我们可以看到傅立叶级数和傅立叶变换之间的本质区别是后者的谱是连续的， 因此 在公式（5）中合成非周期信号的过程使用积分运算代替了（3）式中的求和运算。 x t X F X F e x t e
F
dF 5
F
dt 6
（三）离散周期信号的傅立叶级数
从上面的介绍我们可以知道，由于连续周期信号的频率范围是从 ∞延伸到∞。然而离 散时间信号的频率范围在（ π，π）或（0，2π），以 N 为周期，频率分量间隔2π/N。因 此离散周期信号的傅立叶级数最多包含 N 个频率分量，这是连续周期信号与离散周期信号 的傅立叶级数之间的主要差别。
N
x t
N
c e
/N
7
c
1 N
x t e
/N
8
公式（7）称为离散时间傅立叶级数（DTFS） 。
（四）离散非周期信号的傅立叶变换
与连续信号的傅立叶分析相似， 我们也可以同样从对周期信号的频谱分析推广到非周期 信号的情况。 1 X ω e dω 9 x n 2π X ω x n e 10
公 式 （ 10 ） 中 的 X ω 以 2π 为 周 期 ， 因 为 任 何 离 散 时 间 信 号 的 频 率 范 围 限 制 在 （ π，π）或（0，2π）。信号在时间域是离散的，因此公式（10）用求和号代替了公式 （6）中的积分号，被称为离散时间傅立叶变换（DTFT） 。
4
数字信号处理
S200607097
杨涛
（五）离散傅立叶变换
从上文对离散时间傅立叶变换的分析可以看出X ω 是序列 x n 的频域等效表示，然而 X ω 是频率的连续函数。由于在计算机中无法表示连续频域，考虑使用X ω 的等间隔频率 样本X（2πk/N）, k 0,1, … , N 1,来表示频域信息,就引出了离散傅立叶变换（DFT） 。需要 注意的是，DFT 和上小节介绍的 DTFT 之间的区别，后者的频率域是连续的。离散傅立叶变 换在频域分析中有着重要的地位，因为它将时域的离散信号转换为离散频率域的表示形式， 使得在计算机以及其他数字系统中应用傅立叶变换成为可能。公式（11）和（12）分别表示 了傅立叶变换和反变换。
N
X k
N
x n e
/N
k
0,1,2 … , N
1 11
x n
1 N
X k e
/N
k
0,1,2 … , N
1 12
三、快速傅立叶变换的算法分析
快速傅立叶变换是离散傅立叶变换的改进和优化，利用 DFT 的相关数学性质，通过使 用各种巧妙的算法提高计算性能。 由于在计算机上实现傅立叶变换需要巨大的计算量， 正是 快速傅立叶变换的出现才使得这些理论和方法能够在计算机上实现和普及。
（一）快速傅立叶变换提出的原因
由第二章的介绍我们知道离散傅立叶变换能够在计算机上进行实现，在诸如线性滤波、 相关分析和谱分析等众多领域有着广泛的应用。然而在实际解决问题时，DFT 的巨大数据运 算量却使得它难以最大地发挥傅立叶变换的优势。 给定长度为 N 的数据序列x n ，将它的 DFT 定义为：
N
X k
x n WN 0
k
N
1 13
/N （14） 其中 WN e 通过对公式（13）的计算过程进行分析可以发现，公式首先进行了 N 次的复数乘法
x n WN ，每次复数乘法可以分解为 4 次实数乘法；随后公式进行了 N‐1 次的复数加法，共 进行了 4（N‐1）次实数加法。以此对 N 个数值总共进行 4N 次实数乘法和 4N（N‐1）次实 数加法运算。根据欧拉公式（2）计算WN 还要进行 2N 次三角函数运算。 ，当信号序列的 N 数值比 从上面的分析可知，离散傅立叶变换的时间复杂度是 O（N ） 较大，DFT 算法的计算量是相当巨大的。1965 年，Cooley 和 Turkey 发表论文，提出了 FFT 的概念和算法。由此，众多的数学家、工程师和科学家们都不断地努力寻找各种方法来降低 DFT 的计算量，提出了许多有效的快速傅立叶变换算法，正是有赖于这些高效的算法，傅立 叶变换才得以在各领域获得广泛应用。
5
数字 字信号处理
S S200607097
杨涛 杨
（二 二）快速 速傅立叶 叶变换的原 原理
要为 DFT 提 提出高效的算 算法需要抓住 住问题的关键 键。细致分析 DFT 的公式 式，可以发现 现算法 效率 率低的主要原 原因是没有 有充分利用相 相位因子WN 的对称性 （WN = WN ）和周 期性 N （WN N =WN ）这 W 这两个基本性 性质。基于对 对问题的分析 析，现在已经 经有很多种快 快速算法被提出， 其中 中最重要的一 一类是分解征 征服方法。 类方法的主要 此类 要思路是将被 被计算序列分 分解为较短的子序 列进 DFT 计算 进行 算，使用最为普遍的主要有基 2、基 4 和劈分基 F 算法。本 FFT 本文将以基 2 的快 速算 算法为例进行 行分析，并在 在下文中详述按照时间和频 频率抽取的实 实现方法。
N/
（三 三）基 2 快速傅 2 傅立叶算 算法分析
分解征服方 方法的原理是 是将长度为 N 的信号序列 列逐次地分解 解为计算长度 度较短的序列进行 离散 散傅立叶变换 换。将 N 可 可因子分解为 N=r r r …. r 其中r 为素 为 素数时，算法 法非常有效 。当 r r r …. r r时,称 称这样具有整齐计算模式的 DFT 为 r 以基的 FFT。下面设定 r=2， 的 r 因此 此信号序列长 N=2 。 长度 首先如果将 将信号序列 x(n)劈分为长度 度各位 N/2 的数据子序列 n 和f n , f n 和f n 分 的 列f 别对 对应偶序号样 样本和奇序列 列样本。这种按照因子 2 对 x(n)抽取的 的方法称为按 按时间抽取算 算法， 区别于上 因为 为这种方法是 是对时域信号 号序列进行划分。 上面的时域信 信号抽取， 我 我们也可以对 对频域 信号 X(k)划分为 X(2k)和 X(2k 号 k+1)，并把这 这样的方法称 称为按频率抽 抽取算法。 经过 v=log N次重复上述 述过程， 到每个子 直到 序列 2 个信号 列为 号为止，并进 DFT 计算 进行 算。由此 整个 个信号序列的 的计算过程形 形如一棵满二 二叉树结 构，如图 1 所示 示。从图中可 可以看出，8 点的基 2FFT 需要 3 层计 T 计算，每层分 分别由2 个 个子序列 组成 成，如图 1 中 中第一层有 4 个子块。先 先计算 4 个 2 点 DFT， 后计算 2 个 4 点 DFT， 然后 最后计算 最 图表 1 一个 8 点 DFT。 个 信号序列 由于经过快 快速傅立叶变 变换的计算之后， 列的顺序会出 出现变化， 因 因此在最后一步要 使用 用一定的算法 法将结果进行 行重新排序。相 相反，也可以 以在计算前对 对时域信号重 重新排列，输 输出的 结果 果即为自然顺 顺序，图表 1 所使用的就是这种方法。 。 使用快速算 算法，最终将 将乘法运算量 量减为 log N 次，加法运算 次 算量为Nlog N 次，可以有效地 提高 DFT 的计 高了 计算速度，下 下面一章将详细分析基 2 算法的程序实 算 实现并分析试 试验结果。
N
四、快速 立叶算法 实现 四 速傅立 法的实
（一 一）算法 法设计
基 2 快速傅 傅立叶算法的 的设计是根据相位因子WN 的对称性和 和周期性，最 最大限度地减 减小不 N 必要 要的计算量，将信号按照 照时间或频率抽取后进行蝶 蝶形运算。
6
数字 字信号处理
S S200607097
杨涛 杨
（1 1）计算 算相位因子 子
从公式（13）我们可以 以看到，在 DFT 的计算 算过程中有大 大量的工作用 用来计算WN ， WN k=0,1,…,N‐1。然 然而从第三章第二小节的论 论述可知WN = WN 且 N N =WN ，以此我们只需 且W N 计算 N/2 个相位 算 位因子就可以 以完全得到其它因子的值，没有必要进 进行重复的计 计算。 已知WN e
/N N/
， 欧拉公式 ） 由欧 （2） 可得WN
cos 2πi/N c
sin 2πi/N j ， i=0,1,…,
N
1，
通过 过循环语句可 可以计算这 N/ 个相位因子并将其保存 /2 存在复数数组 组中。
（2 2）蝶形 形运算
计算 基 2 算法中 中最关键的一 一步就是蝶形运算， 无论是 是按照时域抽 抽取还是按照 照频域抽取， 过程 程的本质是一 一样的，都是 是提取两个复数，与相位因 因子WN 乘积 积之后再作加 加减运算，最后生 成两 两个新的复数 数。 图表 3 图表 2 图表 2 显示 示的是按照时 时间抽取的蝶 蝶形计算方式 式，首先提取复数（a,b） 用WN 乘以 b， ， 以 然后 a 进行加 后与 加减运算，得到新的复数（A,B） 表 3 显示的是 。图表 是按照频率抽 抽取的蝶形计算方 式，首先提取复 复数（a,b） 然后计算它们 ，然 们的和生成新 新的复数 A，用WN 乘以 以它们的差得到新 的复 B。 复数
外层循环：共进 进行v个层次 次的计算过 过程 中层循环：每个层次 次过程分为 为k块进行 行计算 内层循 循环：每个 个块内蝶形 形计算N/k k个元素
图表 4
蝶形运算的 的整体过程正 正如图表 4 所显示，运算分 v=log N 分为 N层，在每一 一层中要对信号序 列进 进行分块，由于是基 2 的快 快速算法，每 每次分块均是 是依据 2 的幂 幂次，同时也 也决定了每一个块 内的 的元素数量。 按照时间抽取 按 取和按频率抽 抽取的方法的 的区别就在每 每层中分块顺 顺序和块内计 计算的 差异 异。 图示 5 显示 示了按照时间 间抽取的基 2 算法 蝶形 形运算的整体 体过程，每个 个子块内的元 元素个 数从 2 开始，随 从 随后每级计算 算时块内元素 素数增 加一 一倍，块数减 减少一倍。按 按照频率抽取 取的算 图表 5 表 法恰 恰恰相反，子 子块内的元素 N 开始，随后 素从 每级 级计算时块内 内元素数减少 少一半，块数 数增加 一倍 倍。在块内就 就分别使用图表 2 和图表 3 所 表 示的 的蝶形运算过 过程进行计算 算，每一级计 计算的 输出 出结果又会成 成为下一级的 的输入数据。
7
数字信号处理
S200607097
杨涛
在每一级的块内蝶形计算中， 每个复数会与距离半个块长度的另一个复数进行计算， 因 此每一级计算的结果都与初始序列不一致。 但这其中是有规律可循的， 下面就将详细分析信 号序号的变化，并根据变化的特点设计算法将序号恢复到初始状态。
（3）重排信号序列
从上小节对蝶形运算的分析和实现可知， 经过快速计算后， 信号序列原来的顺序已经发 生了变化。根据计算过程我们看出，每经过一次计算，信号序号的二进制表示的后 v‐i+1 位 将进行循环右移操作，所以经过 v‐1 次数据抽取之后，序号的二进制表示将会反转（其中 v=log N，i 为抽取次数） 。如下图所示，红黄蓝三色的矩形分表表示二进制序号的高中低三 位，没有进过抽取时信号的序号是 3。当进行一次抽取后，二进制表示的后三位循环右移， 序号变为 5；当进行第二次抽取后，二进制表示的后两位循环右移，序号变为 6。观察最后 结果的二进制表示， 我们发现它正是初始序号二进制表示的倒序， 所以只要通过反转结果的 二进制表示就可以恢复原来的顺序输出。 图表 6 0 1 1 数据抽取 1 1 0 1 数据抽取 2 1 1 0 下面就来分析对二进制表示的反转方法。 思路一： 可以通过将十进制数转换为二进制数表示， 将保存在数组中的每个二进制位交 换顺序，最后将变换后的二进制数重新转换为十进制数。这种思路简单，然而其中涉及到除 法、取模等运算，此外还要利用附加空间，计算效率不高，将会使 FFT 的效果打折扣。 思路二： 使用十进制数的移位和按位与操作可以避免使用乘法和除法操作， 巧妙地提取 和累加，最后完成二进制表示的反转。具体步骤如下： （a） 把 1 向左位移（
& 第一位是 1， 将 1 左移 v‐i‐1=2 位 ， & 第二位是 1， & 第三位是 0，
结果变量加 结果变量加 因此最后序 i=0 True i=2 False i=1 True 4，暂时为 4 2，暂时为 6 号结果为 6 由于信号的排序变化是在进行蝶形运算中产生的， 所以在蝶形运算之前和之后进行位序 重排都是可以的。 本文所用的设计在时域抽取的算法中先对初始序号进行了位序重排， 在频 0 0 1 0 1 0
v‐i‐1=1 位，
将 1 左移
1
0
0
不需要改变 结果变量，
8
数字信号处理
S200607097
杨涛
域抽取的算法中对结果序号进行了位序重排，输出的结果是一致的。
（二）结果分析
图示 8 显示的分别由时间抽取算法 TD_FFT2 和频率抽取算法 FD_FFT2 计算的 FFT 结果， 初始 的输入数据为 A=[8，7，6，5，4，3，2，1]。 为了验证程序计算的正确性和准确性，使用 MATLAB 7 中的基 2 函数 fft2(A)计算的结果如下： ans= 36.0000 4.0000 ‐ 9.6569i 4.0000 ‐ 4.0000i 4.0000 ‐ 1.6569i 4.0000 4.0000 + 1.6569i 4.0000 + 4.0000i 4.0000 + 9.6569i 从而确认所编写的程序实现了傅立叶变换的功 能并达到了一定精度。
图表 8
五、参考文献
[1] 数字信号处理：原理、算法与应用 作者：John G.Proakis Dimitris G.Manolakis 张晓林 译 肖创柏 译审 [2] An algorithm for the machine calculation of complex Fourier series 作者：J. W. Cooley and J. W. Tukey [3] Trigonometric Delights chapter15: Fourier’s Theorem 作者：Eli Maor 普林斯顿大学 [4] A First Course in Wavelets with Fourier analysis 作者：Albert Boggess Francis J.Narcowich [5] An Introduction to Fourier Theory 作者：Forrest Moffman [6] An Intuitive Explanation of Fourier Theory 作者：Steven Lehar [7] Fourier Transforms in Crystallography 作者：Kevin Cowtan University of York, UK [8] FFT Tutorial University of Rhode Island Department of Electrical and Computer Engineering [9] Fourier series & Fourier Transform-- from Wolfram MathWorld [10] 相关参考网站 http://www.fftw.org/links.html M. Frigo and S. G. Johnson http://www.earlevel.com/Digital%20Audio/FFT.html 快速傅立叶的原理分析 http://www.relisoft.com/Science/Physics/fft.html 快速傅立叶算法分析 http://www-history.mcs.st-andrews.ac.uk/Printonly/Fourier.html 傅立叶简介
9
DSP 2006-2007
快速傅立叶算法分析
The College of Computer Science
Beijing University of Technology S200607097 杨涛
Email : [email protected]
数字信号处理
S200607097
杨涛
目录
一、傅立叶分析理论的提出 ................................................................................................... 3 二、傅立叶分析的类别及其联系 ........................................................................................... 3 （一）连续周期信号的傅立叶级数 ................................................................................... 3 （二）连续非周期信号的傅立叶变换 ............................................................................... 4 （三）离散周期信号的傅立叶级数 ................................................................................... 4 （四）离散非周期信号的傅立叶变换 ............................................................................... 4 （五）离散傅立叶变换 ....................................................................................................... 5 三、快速傅立叶变换的算法分析 ........................................................................................... 5 （一）快速傅立叶变换提出的原因 ................................................................................... 5 （二）快速傅立叶变换的原理 ........................................................................................... 6 （三）基 2 快速傅立叶算法分析 ...................................................................................... 6 四、快速傅立叶算法的实现 ................................................................................................... 6 （一）算法设计 ................................................................................................................... 6 （1）计算相位因子Wn .................................................................................................. 7 （2）蝶形运算 ................................................................................................................. 7 （3）重排信号序列 ......................................................................................................... 8 （二）结果分析 ................................................................................................................... 9 五、参考文献 ........................................................................................................................... 9
2
数字信号处理
S200607097
杨涛
一、傅立叶分析理论的提出
在物理学研究中， 科学家发现当白光通过玻璃棱镜可被分解为不同颜色的光， 每一种颜 色对应一种可见光谱的特定频率， 因此由光到颜色的分析就是一种频率分析。 将这种思想应 用到对信号性质的研究中，就引出了傅立叶分析理论。 将信号展开成傅立叶级数的基本目的就是要将时间变量 t 的函数分解为不同的频率分 量，这些基本的构造分量正弦和余弦函数。数学家已经证明，形如（1）式的三角函数的和 式就足以表示绝大多数曲线。 cos sin （1）
对周期函数进行这种分解，被称为傅立叶级数；对有限能量信号函数，这种分解被称为 傅立叶变换。三角函数展开起源于 18 世纪关于简谐振动以及类似的物理现象的研究。1808 年傅立叶在著作《热的分析理论》中详细地研究了三角函数，并利用三角函数解决了很多热 传导问题。在随后的二百年间，傅立叶分析理论得到不断的发展，并被应用于众多的领域， 例如用来消除信号中的高频噪声或进行数据压缩等等。
二、傅立叶分析的类别及其联系
根据信号的各种分类和相应特点， 我们可以将施加其上的傅立叶分析分为诸如连续信号 傅立叶级数和傅立叶变换、 离散信号傅立叶级数和傅立叶变换； 根据对信号的处理策略上的 改进，又由离散傅立叶变换（DFT）衍生了快速傅立叶变换（FFT）等等。在刚刚接触这些术 语和缩写时常常会感到混乱和迷惑。 为了更加清晰地认识其理论脉络， 下面将简要介绍它们 各自的应用范围，以及它们之间的联系与区别。
（一）连续周期信号的傅立叶级数
傅立叶级数是周期信号的基本数学表示， 它是相关正弦信号或负指数信号的线性加权和。 出于数学上分析和表示的方便，常根据欧拉公式（2）将傅立叶级数表示为复指数信号和的 信号。 e cos sin （2） Dirichlet 条件保证了除不连续时刻外，级数式（3）一定和初始信号 x(t)相吻合。 x t c c e x t e
F
3
F
TP TP
dt 4
其中系数 c 确定波形的形状，基频F 是给定周期TP 的倒数，确定基本周期，选取适当 的基频F 和系数 c 就可以构造不同类型的周期信号。
3
数字信号处理
S200607097
杨涛
（二）连续非周期信号的傅立叶变换
由上面的公式（3）我们可以知道，周期信号的频率表示是离散的图形，谱线之间有相 等的间隔，这些间隔等于基本频率。如果信号的周期趋近于无限，信号变为非周期的，谱线 间隔将趋近 0，因此信号的频率表示变成连续函数。非周期信号通过周期为TP 的重复得到周 期信号， 我们可以看到傅立叶级数和傅立叶变换之间的本质区别是后者的谱是连续的， 因此 在公式（5）中合成非周期信号的过程使用积分运算代替了（3）式中的求和运算。 x t X F X F e x t e
F
dF 5
F
dt 6
（三）离散周期信号的傅立叶级数
从上面的介绍我们可以知道，由于连续周期信号的频率范围是从 ∞延伸到∞。然而离 散时间信号的频率范围在（ π，π）或（0，2π），以 N 为周期，频率分量间隔2π/N。因 此离散周期信号的傅立叶级数最多包含 N 个频率分量，这是连续周期信号与离散周期信号 的傅立叶级数之间的主要差别。
N
x t
N
c e
/N
7
c
1 N
x t e
/N
8
公式（7）称为离散时间傅立叶级数（DTFS） 。
（四）离散非周期信号的傅立叶变换
与连续信号的傅立叶分析相似， 我们也可以同样从对周期信号的频谱分析推广到非周期 信号的情况。 1 X ω e dω 9 x n 2π X ω x n e 10
公 式 （ 10 ） 中 的 X ω 以 2π 为 周 期 ， 因 为 任 何 离 散 时 间 信 号 的 频 率 范 围 限 制 在 （ π，π）或（0，2π）。信号在时间域是离散的，因此公式（10）用求和号代替了公式 （6）中的积分号，被称为离散时间傅立叶变换（DTFT） 。
4
数字信号处理
S200607097
杨涛
（五）离散傅立叶变换
从上文对离散时间傅立叶变换的分析可以看出X ω 是序列 x n 的频域等效表示，然而 X ω 是频率的连续函数。由于在计算机中无法表示连续频域，考虑使用X ω 的等间隔频率 样本X（2πk/N）, k 0,1, … , N 1,来表示频域信息,就引出了离散傅立叶变换（DFT） 。需要 注意的是，DFT 和上小节介绍的 DTFT 之间的区别，后者的频率域是连续的。离散傅立叶变 换在频域分析中有着重要的地位，因为它将时域的离散信号转换为离散频率域的表示形式， 使得在计算机以及其他数字系统中应用傅立叶变换成为可能。公式（11）和（12）分别表示 了傅立叶变换和反变换。
N
X k
N
x n e
/N
k
0,1,2 … , N
1 11
x n
1 N
X k e
/N
k
0,1,2 … , N
1 12
三、快速傅立叶变换的算法分析
快速傅立叶变换是离散傅立叶变换的改进和优化，利用 DFT 的相关数学性质，通过使 用各种巧妙的算法提高计算性能。 由于在计算机上实现傅立叶变换需要巨大的计算量， 正是 快速傅立叶变换的出现才使得这些理论和方法能够在计算机上实现和普及。
（一）快速傅立叶变换提出的原因
由第二章的介绍我们知道离散傅立叶变换能够在计算机上进行实现，在诸如线性滤波、 相关分析和谱分析等众多领域有着广泛的应用。然而在实际解决问题时，DFT 的巨大数据运 算量却使得它难以最大地发挥傅立叶变换的优势。 给定长度为 N 的数据序列x n ，将它的 DFT 定义为：
N
X k
x n WN 0
k
N
1 13
/N （14） 其中 WN e 通过对公式（13）的计算过程进行分析可以发现，公式首先进行了 N 次的复数乘法
x n WN ，每次复数乘法可以分解为 4 次实数乘法；随后公式进行了 N‐1 次的复数加法，共 进行了 4（N‐1）次实数加法。以此对 N 个数值总共进行 4N 次实数乘法和 4N（N‐1）次实 数加法运算。根据欧拉公式（2）计算WN 还要进行 2N 次三角函数运算。 ，当信号序列的 N 数值比 从上面的分析可知，离散傅立叶变换的时间复杂度是 O（N ） 较大，DFT 算法的计算量是相当巨大的。1965 年，Cooley 和 Turkey 发表论文，提出了 FFT 的概念和算法。由此，众多的数学家、工程师和科学家们都不断地努力寻找各种方法来降低 DFT 的计算量，提出了许多有效的快速傅立叶变换算法，正是有赖于这些高效的算法，傅立 叶变换才得以在各领域获得广泛应用。
5
数字 字信号处理
S S200607097
杨涛 杨
（二 二）快速 速傅立叶 叶变换的原 原理
要为 DFT 提 提出高效的算 算法需要抓住 住问题的关键 键。细致分析 DFT 的公式 式，可以发现 现算法 效率 率低的主要原 原因是没有 有充分利用相 相位因子WN 的对称性 （WN = WN ）和周 期性 N （WN N =WN ）这 W 这两个基本性 性质。基于对 对问题的分析 析，现在已经 经有很多种快 快速算法被提出， 其中 中最重要的一 一类是分解征 征服方法。 类方法的主要 此类 要思路是将被 被计算序列分 分解为较短的子序 列进 DFT 计算 进行 算，使用最为普遍的主要有基 2、基 4 和劈分基 F 算法。本 FFT 本文将以基 2 的快 速算 算法为例进行 行分析，并在 在下文中详述按照时间和频 频率抽取的实 实现方法。
N/
（三 三）基 2 快速傅 2 傅立叶算 算法分析
分解征服方 方法的原理是 是将长度为 N 的信号序列 列逐次地分解 解为计算长度 度较短的序列进行 离散 散傅立叶变换 换。将 N 可 可因子分解为 N=r r r …. r 其中r 为素 为 素数时，算法 法非常有效 。当 r r r …. r r时,称 称这样具有整齐计算模式的 DFT 为 r 以基的 FFT。下面设定 r=2， 的 r 因此 此信号序列长 N=2 。 长度 首先如果将 将信号序列 x(n)劈分为长度 度各位 N/2 的数据子序列 n 和f n , f n 和f n 分 的 列f 别对 对应偶序号样 样本和奇序列 列样本。这种按照因子 2 对 x(n)抽取的 的方法称为按 按时间抽取算 算法， 区别于上 因为 为这种方法是 是对时域信号 号序列进行划分。 上面的时域信 信号抽取， 我 我们也可以对 对频域 信号 X(k)划分为 X(2k)和 X(2k 号 k+1)，并把这 这样的方法称 称为按频率抽 抽取算法。 经过 v=log N次重复上述 述过程， 到每个子 直到 序列 2 个信号 列为 号为止，并进 DFT 计算 进行 算。由此 整个 个信号序列的 的计算过程形 形如一棵满二 二叉树结 构，如图 1 所示 示。从图中可 可以看出，8 点的基 2FFT 需要 3 层计 T 计算，每层分 分别由2 个 个子序列 组成 成，如图 1 中 中第一层有 4 个子块。先 先计算 4 个 2 点 DFT， 后计算 2 个 4 点 DFT， 然后 最后计算 最 图表 1 一个 8 点 DFT。 个 信号序列 由于经过快 快速傅立叶变 变换的计算之后， 列的顺序会出 出现变化， 因 因此在最后一步要 使用 用一定的算法 法将结果进行 行重新排序。相 相反，也可以 以在计算前对 对时域信号重 重新排列，输 输出的 结果 果即为自然顺 顺序，图表 1 所使用的就是这种方法。 。 使用快速算 算法，最终将 将乘法运算量 量减为 log N 次，加法运算 次 算量为Nlog N 次，可以有效地 提高 DFT 的计 高了 计算速度，下 下面一章将详细分析基 2 算法的程序实 算 实现并分析试 试验结果。
N
四、快速 立叶算法 实现 四 速傅立 法的实
（一 一）算法 法设计
基 2 快速傅 傅立叶算法的 的设计是根据相位因子WN 的对称性和 和周期性，最 最大限度地减 减小不 N 必要 要的计算量，将信号按照 照时间或频率抽取后进行蝶 蝶形运算。
6
数字 字信号处理
S S200607097
杨涛 杨
（1 1）计算 算相位因子 子
从公式（13）我们可以 以看到，在 DFT 的计算 算过程中有大 大量的工作用 用来计算WN ， WN k=0,1,…,N‐1。然 然而从第三章第二小节的论 论述可知WN = WN 且 N N =WN ，以此我们只需 且W N 计算 N/2 个相位 算 位因子就可以 以完全得到其它因子的值，没有必要进 进行重复的计 计算。 已知WN e
/N N/
， 欧拉公式 ） 由欧 （2） 可得WN
cos 2πi/N c
sin 2πi/N j ， i=0,1,…,
N
1，
通过 过循环语句可 可以计算这 N/ 个相位因子并将其保存 /2 存在复数数组 组中。
（2 2）蝶形 形运算
计算 基 2 算法中 中最关键的一 一步就是蝶形运算， 无论是 是按照时域抽 抽取还是按照 照频域抽取， 过程 程的本质是一 一样的，都是 是提取两个复数，与相位因 因子WN 乘积 积之后再作加 加减运算，最后生 成两 两个新的复数 数。 图表 3 图表 2 图表 2 显示 示的是按照时 时间抽取的蝶 蝶形计算方式 式，首先提取复数（a,b） 用WN 乘以 b， ， 以 然后 a 进行加 后与 加减运算，得到新的复数（A,B） 表 3 显示的是 。图表 是按照频率抽 抽取的蝶形计算方 式，首先提取复 复数（a,b） 然后计算它们 ，然 们的和生成新 新的复数 A，用WN 乘以 以它们的差得到新 的复 B。 复数
外层循环：共进 进行v个层次 次的计算过 过程 中层循环：每个层次 次过程分为 为k块进行 行计算 内层循 循环：每个 个块内蝶形 形计算N/k k个元素
图表 4
蝶形运算的 的整体过程正 正如图表 4 所显示，运算分 v=log N 分为 N层，在每一 一层中要对信号序 列进 进行分块，由于是基 2 的快 快速算法，每 每次分块均是 是依据 2 的幂 幂次，同时也 也决定了每一个块 内的 的元素数量。 按照时间抽取 按 取和按频率抽 抽取的方法的 的区别就在每 每层中分块顺 顺序和块内计 计算的 差异 异。 图示 5 显示 示了按照时间 间抽取的基 2 算法 蝶形 形运算的整体 体过程，每个 个子块内的元 元素个 数从 2 开始，随 从 随后每级计算 算时块内元素 素数增 加一 一倍，块数减 减少一倍。按 按照频率抽取 取的算 图表 5 表 法恰 恰恰相反，子 子块内的元素 N 开始，随后 素从 每级 级计算时块内 内元素数减少 少一半，块数 数增加 一倍 倍。在块内就 就分别使用图表 2 和图表 3 所 表 示的 的蝶形运算过 过程进行计算 算，每一级计 计算的 输出 出结果又会成 成为下一级的 的输入数据。
7
数字信号处理
S200607097
杨涛
在每一级的块内蝶形计算中， 每个复数会与距离半个块长度的另一个复数进行计算， 因 此每一级计算的结果都与初始序列不一致。 但这其中是有规律可循的， 下面就将详细分析信 号序号的变化，并根据变化的特点设计算法将序号恢复到初始状态。
（3）重排信号序列
从上小节对蝶形运算的分析和实现可知， 经过快速计算后， 信号序列原来的顺序已经发 生了变化。根据计算过程我们看出，每经过一次计算，信号序号的二进制表示的后 v‐i+1 位 将进行循环右移操作，所以经过 v‐1 次数据抽取之后，序号的二进制表示将会反转（其中 v=log N，i 为抽取次数） 。如下图所示，红黄蓝三色的矩形分表表示二进制序号的高中低三 位，没有进过抽取时信号的序号是 3。当进行一次抽取后，二进制表示的后三位循环右移， 序号变为 5；当进行第二次抽取后，二进制表示的后两位循环右移，序号变为 6。观察最后 结果的二进制表示， 我们发现它正是初始序号二进制表示的倒序， 所以只要通过反转结果的 二进制表示就可以恢复原来的顺序输出。 图表 6 0 1 1 数据抽取 1 1 0 1 数据抽取 2 1 1 0 下面就来分析对二进制表示的反转方法。 思路一： 可以通过将十进制数转换为二进制数表示， 将保存在数组中的每个二进制位交 换顺序，最后将变换后的二进制数重新转换为十进制数。这种思路简单，然而其中涉及到除 法、取模等运算，此外还要利用附加空间，计算效率不高，将会使 FFT 的效果打折扣。 思路二： 使用十进制数的移位和按位与操作可以避免使用乘法和除法操作， 巧妙地提取 和累加，最后完成二进制表示的反转。具体步骤如下： （a） 把 1 向左位移（
& 第一位是 1， 将 1 左移 v‐i‐1=2 位 ， & 第二位是 1， & 第三位是 0，
结果变量加 结果变量加 因此最后序 i=0 True i=2 False i=1 True 4，暂时为 4 2，暂时为 6 号结果为 6 由于信号的排序变化是在进行蝶形运算中产生的， 所以在蝶形运算之前和之后进行位序 重排都是可以的。 本文所用的设计在时域抽取的算法中先对初始序号进行了位序重排， 在频 0 0 1 0 1 0
v‐i‐1=1 位，
将 1 左移
1
0
0
不需要改变 结果变量，
8
数字信号处理
S200607097
杨涛
域抽取的算法中对结果序号进行了位序重排，输出的结果是一致的。
（二）结果分析
图示 8 显示的分别由时间抽取算法 TD_FFT2 和频率抽取算法 FD_FFT2 计算的 FFT 结果， 初始 的输入数据为 A=[8，7，6，5，4，3，2，1]。 为了验证程序计算的正确性和准确性，使用 MATLAB 7 中的基 2 函数 fft2(A)计算的结果如下： ans= 36.0000 4.0000 ‐ 9.6569i 4.0000 ‐ 4.0000i 4.0000 ‐ 1.6569i 4.0000 4.0000 + 1.6569i 4.0000 + 4.0000i 4.0000 + 9.6569i 从而确认所编写的程序实现了傅立叶变换的功 能并达到了一定精度。
图表 8
五、参考文献
[1] 数字信号处理：原理、算法与应用 作者：John G.Proakis Dimitris G.Manolakis 张晓林 译 肖创柏 译审 [2] An algorithm for the machine calculation of complex Fourier series 作者：J. W. Cooley and J. W. Tukey [3] Trigonometric Delights chapter15: Fourier’s Theorem 作者：Eli Maor 普林斯顿大学 [4] A First Course in Wavelets with Fourier analysis 作者：Albert Boggess Francis J.Narcowich [5] An Introduction to Fourier Theory 作者：Forrest Moffman [6] An Intuitive Explanation of Fourier Theory 作者：Steven Lehar [7] Fourier Transforms in Crystallography 作者：Kevin Cowtan University of York, UK [8] FFT Tutorial University of Rhode Island Department of Electrical and Computer Engineering [9] Fourier series & Fourier Transform-- from Wolfram MathWorld [10] 相关参考网站 http://www.fftw.org/links.html M. Frigo and S. G. Johnson http://www.earlevel.com/Digital%20Audio/FFT.html 快速傅立叶的原理分析 http://www.relisoft.com/Science/Physics/fft.html 快速傅立叶算法分析 http://www-history.mcs.st-andrews.ac.uk/Printonly/Fourier.html 傅立叶简介
9