快速傅里叶变换FFT算法及其应用 [当文网提供]

快速傅里叶变换FFT算法及其应用

摘要

本文较为系统地阐述了快速傅里叶变换的算法原理及其在数字信号处理等工程技术中的应用。根据抽取方法的不同，一维基2 FFT算法分为两种：频域抽取的FFT算法和时频域抽取的FFT算法。第1节阐述了这两种FFT算法的原理。第2节给出了两种算法的编程思想和步骤。第3节阐述了一维非基2 FFT的两种算法：Cooley-tukey FFT算法和素因子算法(Prime Factor Algorithm)的思想原理，给出了在把一维非基2 DFT的多层分解式转化为二层分解的过程中，如何综合运用这两种算法以达到总运算次数最少的方案；并以20点DFT为例描述了非基2 FFT算法实现的一般步骤。第4节介绍了一维FFT算法在计算连续时间信号的傅里叶变换、离散信号的线性卷积、离散信号压缩和滤波等数字信号处理中的典型应用。第5节把一维FFT变换推广到二维FFT变换，并在一维FFT算法的基础上，给出了二维FFT算法的原理和实现过程。最后在附录中给出了一维DFT 的基2 FFT 算法（包括频域抽取的FFT和IFFT算法、时域抽取的FFT和IFFT 算法），一维任意非基2 FFT算法，二维DFT的基2 FFT 算法以及二维DFT的任意非基2 FFT 算法的详细的Visual C++程序。

本文通过各种流程图和表格，较为深入系统地阐述了FFT的算法原理；运用Matlab编程，通过大量生动的实例，图文并茂地列举出了FFT算法的各种应用，并在每个实例中都附上了完整的Matlab程序，可供读者参考。由于篇幅所限，本文未涉及FFT变换以及其应用的数学理论背景知识。

关键词：FFT算法的应用，一维基2 FFT算法，频域抽取，时域抽取，非基2 FFT算法，Cooley-Tukey算法，素因子算法，线形卷积，信号压缩和滤波，二维FFT算法

快速傅里叶变换FFT算法及其应用

摘要

目录

摘要 0

目录 (2)

1 一维DFT的快速算法—FFT (3)

1.1频域抽取的基2算法 (3)

1.1.1 正变换的计算 (3)

1.1.2 逆变换的计算 (6)

1.2时域抽取的基2算法 (7)

2 一维基2 FFT算法编程 (8)

3 一维任意非基2 FFT算法 (12)

3.1C OOLEY-T UKEY FFT算法 (12)

3.2素因子算法（PFA） (13)

3.3一维任意非基2FFT算法 (15)

4 一维FFT算法的应用 (18)

4.1利用FFT计算连续时间信号的傅里叶变换 (18)

4.2利用FFT计算离散信号的线性卷积 (21)

4.3利用FFT进行离散信号压缩 (23)

4.4利用FFT对离散信号进行滤波 (26)

4.5利用FFT提取离散信号中的最强正弦分量 (29)

5 二维DFT的快速变换算法及应用简介 (34)

5.1二维FFT变换及其算法介绍 (34)

5.2二维FFT变换算法的应用 (35)

参考文献 (35)

附录 (36)

1．一维DFT的基2FFT算法V ISUAL C++程序 (36)

（1）频域抽取的FFT和IFFT算法 (36)

（2）时域抽取的FFT和IFFT算法 (41)

2．一维任意非基2FFT算法V ISUAL C++程序 (46)

快速傅里叶变换FFT 算法及其应用

3．二维DFT 的基2 FFT 算法V ISUAL C++程序 ............................................... 51 4．二维DFT 的任意非基2 FFT 算法V ISUAL C++程序 .. (59)

1 一维DFT 的快速算法—FFT

当序列[]f n 的点数不超过N 时，它的N 点D FT 定义为

21

[][]01N i

k n

N

n F k f n e k

N π

--==

≤≤-∑

（1）

反变换ID FT 定义为 21

1[][]01N i

k n

N

k f n F k e

n

N N

π-==

≤≤-∑

（2）

二者形式相似，快速算法的原理一样，这里先就其正变换进行讨论。令

2/i N

N W e

π-=，当k 依次取为0,1,2,,1N - 时，可表示为如下的方程组：

00

01

02

0(

11011

12

1(1)2021222(1)

(1)0

(

1)1(1)(1)

[0][0][1][2][1]

[1][0][1][2][

1]

[2][0][1][2][

1][1][0][1]

[1]

N N N N N N N N N N

N N N N N

N N N N N

N N

F f W f W f W f N W F f W f W f W f

N W F f W f

W f

W f N W F N f W f W f N W -------?=

++++-=++++-=

+++

+

-?-=+

++-

?

??

??

??

（3）

由上式可见，直接按照定义计算N 点序列的N 点D FT 时，每行含N 个复乘和N 个加，从而直接按定义计算点的总计算量为2N 个复乘和2N 个加。当N 较大时，2N 很大，计算量过大不仅耗时长，还会因字长有限而产生较大的误差，甚至造成计算结果不收敛。所谓快速傅里叶变换就是能大大减少计算量而完成全部点计算的算法。下面介绍两种经典的D FT 的快速算法：频域抽取的F F T 算法和时域抽取的F F T 算法。

1.1 频域抽取的基2算法

1.1.1 正变换的计算

这里仅介绍基2算法，即是2的整次幂的情况。由定义

1

[][]01N kn

N

n F k f n W k N -==

≤≤-

∑

（4）

1 一维DFT 的快速算法—FFT

把[]f n 分成两半，即[]f n 和[/2]f n N +(0/21)n N ≤≤-，代入（4）式得 /21/21

(/2)

[][][/2]01N N

kn

k n N N

N

n n F k f n W f n N W k N --+===

+

+≤≤-∑

∑

（5）

由于

(/2)

/2

(1)k n N kn kN k kn

N

N N

N

W W W W +==-

（5）式两项又可合并为

/21

[]{[](

1)[/2]}

01N k

k n

N n F k f n f n N

W k

N -

==

+-+≤≤-

∑

（6） 当2k r =为偶数时，注意到(1)1k -=，222/kn rn i rn N N N W W e π-==/2rn N W =，（6）式变为

/21

/20/21

/2

[2]([][/2])()()

0/21

N rn

N n N rn

N n F r f n f n N W g n W G r r N -

=-==

++=

=≤≤-∑

∑

（7）

当21k r =+为奇数时， (21)2(21)//2kn r n i r n N n rn N N N N W W e W W π+-+===，（6）式变为

/21

/2

0/21

/20

[21]{([][/2])}()()

0/21

N n rn

N N n N rn

N n F r f n f n N W W p n W P r r N -

=-=+=

-+=

=≤≤-∑

∑

（8）

这样就把一个N 点序列（[]f n ）的N 点D FT （[]F k ）计算化成了两个/2N 点序列（[]g n 和[]p n ）的/2N 点D FT （[]G r 和[]P r ）计算。由[]f n 划分成[]g n 和[]p n 的计算量为N 个加，即

[][/2]f n f n N ++和[][/2],0/21f n f n N n N -+≤≤-

和/2N 个乘，即

([][/2]),0/21n

N f n f n N W n N -+≤≤-

由于[]g n 算出的/2N 点D FT ，是[]f n 的N 点D FT （[]F k ）中k 为偶数的那一半，由[]p n 算出的则是k 为偶数的那一半，故需要把偶数k 的[]F k 抽出来放在一起作为[]g n 的D FT （()G r ）输出，同时把奇数k 的[]F k 抽出来放在一起作为[]p n 的D FT （()P r ）输出。由于k 是频域变量，故这种算法称为频域抽取的

F F T

算法。

快速傅里叶变换FFT 算法及其应用

接着，两个/2N 点D FT 仍可用上述方法各经/4N 个乘/2N 个加划分成两个

/4

N 点D FT （同时还要做相应的频域抽取），从而共划分成4个/2N 点D FT ，

总划分计算量仍是N 个加和/2N 个乘。这样的划分可一步步做下去，不难看出，每步的总划分计算量都是N 个加和/2N 个乘。

经过1M -步的划分后就划成了/2N 个如下两点D FT 的计算问题

00

01

22

10

110

2

2

2()A aW bW a b

B aW

bW

a b W ?????=+=+??=+=-??? （9）

上式所需计算量是2个加和1个乘，于是完成/2N 个两点D FT 的总计算量仍是N 个加和/2N 个乘。从而完成全部N 点D FT 的总计算量2log M N N N ?=个加和2/2(/2)log M N N N ?=个乘，这比直接按定义计算所需的2N 个乘和加要少得多。例如，1021024N ==，10M =，用F F T 算法计算所需的乘法个数为

/2M N ?51024

=?，而直接按定义计算所需的乘法个数为210241024N =?，二

者相差1024/5200≈倍。若直接计算需半小时，而用F F T 计算只需9s 即可完成，可见其效率之高，而且N 越大，F F T 的效率提高越明显。

图1 频域抽取的8点FFT 计算流图

一般情况下，由于做了1M -次分奇偶的抽取，此算法最后的/2N 个两点

D FT

计算出的[]F k 不是顺序抽取的。次序的变化可用二进码来说明：第一次抽

取所分的奇偶是由二进码第1位是1或0来区分的，该位为0时为偶数，该位为1时为奇数，第二次抽奇偶是由二进码第2位是1或0来区分的……，每次抽取

1 一维DFT 的快速算法—FFT

都是把偶数项放在前（左）边，把奇数项放在后（右）边，从而抽取以后数的二进码是按照二进制位从左向右依次排列的，和普通二进制数从右向左依次排列的的规律正好相反，所以称为倒位序。在计算出[]F k 之后要把倒位序变成顺序。 1.1.2 逆变换的计算

所谓逆变换是指由[]F k 求[]f n 的计算，若直接按定义

1

1[][]01N kn

N

k f n F k W

n N N

--==

≤≤-∑

做计算，则除了求和号和正变换相同的计算量外，每算一个[]f n 都还需再多做一个乘1/N 的乘法运算。故按定义完成全部N 点D FT 的总计算量是2N 个加和

(1)N N +个乘。下面从图导出它的快速算法，先讨论第

3列的2点D FT 的逆运

算如何完成。由式（7）得，

2()a b A

a b W B

+=??-?=? 由上式不难解出

20

21()

21()

2

a A B W

b A B W --=

+=

- （10）

图2 频域抽取的8点IFFT 计算流图

快速傅里叶变换FFT 算法及其应用

此计算过程如图2所示，可以看出：左边各列的划分计算也都是由/2N 个碟形运算来完成的，只是各碟形运算所乘的相移因子W 不同。把每个碟形运算都用图的办法变成对应的逆运算，并把它们按输入在左、输出在右重新排列，就得到了全部N 点ID FT 的计算流图。给出了8N =的示例，图中先对顺序输入的[]F k 做1M -次的频域抽取，

并把3个乘1/2的运算合成了一个乘1/8的运算放在

了最前边，然后就开始做求逆的碟形运算。

1.2 时域抽取的基2算法

比较正变换D FT 和反变换ID FT 的定义式

1

0[][]01N kn

N

n F k f n W k N -==

≤≤-∑

1

1[][]01N kn

N

k f n F k W

n N N

--==

≤≤-∑

可见，去掉乘1/N 的运算，把1W -换成W ，交换[]F k 、[]f n 和k 、n ，反变换定义式就变成了正变换的定义式。对图2做这些变换，则得到图3的流程图。对图1做这些变换，则得到图4的流程图。这就是时域抽取的算法流图。进行碟形运算之前，先要对顺序的时域输入序列进行1M -次的奇偶抽取，故称为时域抽取的F F T 算法。

图3 时域抽取的8点FFT 计算流图

2 一维基2 FFT 算法编程

比较图2和图3不难看出，两种算法的计算量是完全一样的。这里先算出

/2

N 个两点的D FT

00

01

2

2

()(/2)()(/2)

A f n W f n N W f n f n N ??=++=++

1011

22()(/2)()(/2)B f n W f n N W f n f n N ??=++=-+ （11）

图4 时域抽取的8点IFFT 计算流图

然后把/2N 个两点的D FT 组合成/4N 个4点的D FT ，再把/4N 个4点的

D FT 组合成/8N 个8点的D FT ，经过1M -次的组合之后，就得到了顺序N 点D FT

计算结果。算法原理参考文献【4】。

2 一维基2 FFT 算法编程

以频域抽取的基2 FFT 正变换为例，对FFT 的信号流图进行讨论，以找到FFT 算法的规律。

1）分级

在进行D FT 变换的过程中，从N 点D FT 到两点D FT 共分了M 级，如图1所示，从左到右依次为1m =级，2m =级，…，m M =级。

2）倒位序

快速傅里叶变换FFT 算法及其应用

在频域抽取的基2 FFT 算法中，输出数据不是按照序列的先后顺序排列的，这是由于变换过程中，输出按奇、偶抽取的缘故。如果将序列[]x n 中标号n 用二进制值0212()M M n n n -- 表示，那么在FFT 信号流图输入端，[]x n 位于1202()M M n n n -- 处，

称为倒序。以8点FFT 为例，顺序和倒序的关系如表1所示。

从表1可以看出，一个自然顺序二进制数，是在最低位加1，逢2向左移位；而倒序数的顺序是在最高位加1，逢2向右移位。用i 表示顺序数，j 表示倒序数，k 表示位权重。对于一个倒序数j 来说，下一个倒序数可以按下面的方法求：先对最高位加1，相当于十进制运算/2j N +。如果/2j N <，说明二进制最高位为0，则直接由/2j N +得到下一个倒序值；如果/2j N ≥，说明二进制最高 位为1，则将j 的最高位变为0，通过/2j j N ?-实现，同时令/2k k ?，接着判断次高位是否为0，直到位为0时，令j j k =+ 。归结起来算法流程图如图5所示。

2 一维基2 FFT 算法编程

图5 倒位序程序流程图 图6 频域抽取FFT 程序流程图 3）蝶形运算单元和同址计算

频域抽取的信号流图中，基本的运算结构如图7所示，该运算结构的形状像蝴蝶，故称为“蝶形运算单元”。

在这一结构中，DFT 和IDFT 运算关系分别为

1111[][][]

[]([][])m m m r m m m N

F p F p F q F p F p F q W ----=+???=-??

1111[]([][])/2[]([][])/2

r

m m m N r m m m N f p f p f q W f p f p f q W ------?=+??=-?? （12）

快速傅里叶变换FFT 算法及其应用

(a) 两点DFT 的计算

(b) 两点IDFT 的计算

图7 频域抽取FFT 算法流图中第m 级碟形单元

而在时域抽取的信号流图中，基本的运算结构如图8所示。在这一结构中，DFT 和IDFT 运算关系分别为

1111[][][][][][]r

m m m N r m m m N

F p F p F q W F p F p F

q W ----?=+??=-??

1111[]([][])/2

[]([][])/2

m m m r m m m N f p f p f q f p f p f q W -----=+???=-?? （13）

(a) 两点DFT 的计算

(b) 两点IDFT 的计算

图8 时域抽取FFT 算法流图中第m 级碟形单元

其中，p 、q 分别表示该蝶形运算单元的上下节点的序号。可以看出参与运算的输入序号p ，输出序号仍为q ，并且该运算不涉及到其它的点，因此我们可以将输出的结果仍然放在数组中，称这样的操作为同址计算。也就是说，共同占有同一个存储单元。

4）寻址和相移因子r N W 的计算

时域抽取基2 FFT 信号流图中，每一级有个蝶形单元。每一级的个蝶形单元又可以分为若干组，每一组具有相同的结构和因子的分布。

如图1所示，第1级分为1组，第2级分为2组，…，第m 级分为12m -组。

3 一维任意非基2 FFT 算法

在第m 级中，相邻组之间的间距（也即每个分组所含节点数）为12M m +-，每个蝶形单元的上下节点之间的距离（也即每个分组所含碟形单元数）为2M m -。每组的相移因子

22cos(

)sin(

)r

N W r i r N

N

ππ=-，其中m-1

r=(l-1)2

,

? 0,1,l =,21M m --

综合以上各步骤，得到频域抽取FFT 程序流程图如图6所示。采用类似的步骤可得到频域抽取IFFT 流程图、时域抽取FFT 与IFFT 流程图（略）。

频域抽取FFT 算法、时域抽取FFT 算法的Visual C++源程序分别见附录1.（1），1.（2）。在Matlab 中，傅里叶变换及其逆变换分别用dwt()和idwt()函数实现。

3 一维任意非基2 FFT 算法

3.1 Cooley-Tukey FFT 算法

F F T

的核心是将一层运算映射为二层运算。作N 点F F T 时，若N 不是素数，

则N 可分解为12N N N =，那么由[]f n 的D FT

1

[][]01N nk

N

n F k f n W k N -==

≤≤-

∑

（14）

通过映射：

21211221121

12201,0101,01n N n n n N n N k k N k k N k N =+≤≤-≤≤-=+≤≤-≤≤-

（15）

可得到

212112211121221122()()

()

N n n k N k N n k N N n k n k N n k nk

N N

N

W W W +++++==

而12N N N =，1

2

N N N W W =，2

1

N N N W W =，可化简为

11

21

22

1

2

n k n k n k nk N N N N W W W W = （16）

从而式（14）转化为

2122

21112

1211

1

12120

[,]([,])N N n k n k n k

N N

N n n F k k W

W

f n n W --===∑∑

（17）

其中

112201,01k N k N ≤≤-≤≤-。

以20点F F T 为例，1220,5,4N N N ===，映射方式为：124n n n =+，

快速傅里叶变换FFT 算法及其应用

125k k k =+，则计算流图如图9

所示。

3.2 素因子算法(Prime Factor Algorithm, PFA)

当Cooley-Tukey FFT 算法中的1N 、2N 互素的话，相移因子21

n k N W 可通过选择

12,,n n 12,k k 前的特殊系数而消去，FFT 的算法进一步的简化。

121

1

22

121122

01,0101,01n A n B n n N n N k C k D k k N k N =+≤≤

-≤≤

-

=+≤

≤-≤≤- （18）

当A B C D 、、、满足以下条件

0m od 0m od AD N BC N

≡≡ （19）

图9 Cooley-T ukey 20 点FFT 算法

3 一维任意非基2 FFT 算法

2

1

m o d m o d A C N N B D N N ≡≡

（19）

时，有

1212

11122122

211122

11

2

2

12

()()

()

An Bn C k D k nk

N N

AC n k AD n k BC n k BD n k N N n k N n k N n k

n k N N W W W W

W W ++++++==== （20）

于是式（14）化简为 2122

1121211

1

12120

[,]([,])N N n k n k

N N n n F k k W

f n n W --===

∑∑ （21）

从而消去了相移因子21

n k N W 。同样以20点F F T 为例，修改映射方式为：

20,N = 125,4N N ==

121245,04,03n n n n n =+≤≤≤≤ （22） 12121625,04,03k k k k k =+≤≤≤≤ （23） 则由式（22）得到的映射关系如表2，由式（23）得到的映射关系如表3，计算流图如图10所示。

快速傅里叶变换FFT 算法及其应用

3.3 一维任意非基2 FFT 算法

对于非素数N 点D FT ，对N 做因式分解

12l N N N N = （24） 令22l l N N N = ，则12l N N N =。于是式（24）中多层F F T 转化为二层F F T 。如果1N 与2l N 互素，那么采用PFA 算法进行映射，否则采用Cooley-Tukey FFT 算法（此时需乘以相移因子）。22l l N N N = 可采用同样的方法分解成2N 个3l N 点

D FT

，其中33l l N N N = ，依此类推，直到D FT 长度为l N 。 可以证明，复乘的总次数不大于

12()l N N N N l +++- （25）

例如，若63337N ==??，则复乘总次数至多为63(3373)630???-=。而用基2的F F T 算法计算64点D FT ，需要的复乘总次数为192。这说明，将N 分解得越细，运算量越少。实际中，常常将输入序列补零，使得N 成为2的幂次数，这样能够减少复乘运算的次数。

3 一维任意非基2 FFT 算法

图10 素因子(PFA)20 点FFT 算法

再以20点F F T 为例，进行如下三层因式分解

123522

N N N N ==??

即15N =，23224N =?=，由于5与4互素，所以第一层采用PFA 算法，相应的二层映射为

12312345,04,03

n n n n n =+≤≤≤≤

1231231625,04,03k k k k k =+≤≤≤≤

由于2与2不互素，所以第二层采用Cooley-Tukey FFT 算法，相应的二层映射为

2323232,01,01n n n n n =+≤≤≤≤

快速傅里叶变换FFT 算法及其应用

2323232,01,01k k k k k =+≤≤≤≤

总的F F T 变换公式如式（26），计算流图如图11所示。

图11 多层分解时20 点FFT 算法

321333222113

23

2

13211

1

1

1231230

[,,]{([,,])}N N N n k n k n k n k

N N N N n n n F k k k W

W

W

f n n n W ---====

∑∑∑

4 一维FFT 算法的应用

3332221132111

4

2

4

2

12350

{([,,])}n k n k n k n k

n n n W

W

W

f n n n W ====

∑∑∑

（26）

比较正变换D FT 和反变换ID FT 的定义式可知，在正变换前加上乘1/N 的运算，把W 换成1W -，并交换[]f n 、[]F k 和n 、k ，就得到反变换的算法。

一维任意非基2 FFT 算法的Visual C++源程序见附录2。在Matlab 中，傅里叶变换及其逆变换也分别用dwt()和idwt()函数实现。

4 一维FFT 算法的应用

4.1 利用FFT 计算连续时间信号的傅里叶变换

设()x t 是连续时间信号，并假设0t <时()0x t =，则其傅里叶变换由下式给出

()()i t

X x t e

dt ωω∞

-=

?

令Γ是一个固定的正实数，N 是一个固定的正整数。当

,0,1,2,,1k k N ω=Γ=- 时，利用FFT 算法可计算()X ω。

已知一个固定的时间间隔T ，选择T 足够小，使得每一个T 秒的间隔

(1)nT t n T

≤<+内，()x t 的变化很小，则式中积分可近似为

(1)0()()()

n T iw t

nT

n X e

dt x nT ω∞

+-==

∑?

(

1)0

1[

]()

i t t n T

t n T

n e x n T i ωω

∞

-=+==-=

∑ 0

1()i T

i nT

n e

e

x nT i ωωω

-∞

-=-=

∑ （27）

假设N 足够大，对于所有n N ≥的整数，幅值()x nT 很小，则式（27）变为 1

1()()i T N i nT

n e

X e

x nT i ωωωω

---=-=

∑ （28）

当2/k N T ωπ=时，式（28）两边的值为 2/2/1

2/0

211(

)()[]2/2/i k N i k N

N i nk N

n k e

e

X e

x nT X k N T

i k N T

i k N T

ππππππ----=--=

=

∑ （29）

其中[]X k 代表抽样信号[]()x n x nT =的N 点D FT 。最后令2/N T πΓ=，则上式

快速傅里叶变换FFT 算法及其应用

变为

2/1()[]

0,1,2,,12/i k N

e

X k X k k N i k NT

ππ--Γ=

=- （30）

首先用FFT 算法求出[]X k ，然后可用上式求出0,1,2,,1k N =- 时的

()X k Γ。

应该强调的是，式（28）只是一个近似表示，计算得到的()X ω只是一个近似值。通过取更小的抽样间隔T ，或者增加点数N ，可以得到更精确的值。如果

B ω>时，幅度谱()X ω很小，对应于奈奎斯特抽样频率2s B ω=，抽样间隔T 选

择/B π比较合适。如果已知信号只在时间区间10t t ≤≤内存在，可以通过对

1n T t >时的抽样信号[]()x n x nT =补零，使N

足够大。

例1 利用FFT 计算傅里叶变换 如图12所示的信号

1

02()0

t t x t -≤

?其它

其傅里叶变换为：

2

cos()sin()

()2i X e

i

ω

ωωωωω

--=

利用下面的命令，可得到()x t 的近

似值和准确值。 图12 连续时间信号x(t)

N=input('Input N:'); T=input('Input T:'); %计算X(w)近似值 t=0:T:2;

x=[t-1 zeros(1,N-length(t))]; X=fft(x);

gamma=2*pi/(N*T); k=0:10/gamma;

Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X; %计算真实值X(w)

w=0.05:0.05:10;

Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w); plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact)); legend('近似值','真实值'); xlabel('频率(rad/s)');ylabel('|X|')

运行程序后输入N=128，T=0.1，此时0.4909Γ=，得到实际的和近似的傅里