实验2_FFT的实现

电 子 科 技 大 学

实 验 报 告

学生姓名： 学 号： 指导教师：

一、实验室名称：数字信号处理实验室 二、实验项目名称：FFT 的实现 三、实验原理：

一．FFT 算法思想：

1．DFT 的定义：

对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{x[k]}可以由离散付氏变换（DFT ）求得。DFT 的定义为：

21

[][]N j

nk N

n X k x n e

π--==∑，k=0,1,…N-1

通常令2j

N

N e

W π-=，称为旋转因子。

2．直接计算DFT 的问题及FFT 的基本思想：

由DFT 的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N 点DFT 需要（N-1）2次复数乘法和N （N-1）次加法。因此，对于一些相当大的N 值（如1024）来说，直接计算它的DFT 所作的计算量是很大的。

FFT 的基本思想在于，将原有的N 点序列分成两个较短的序列，这些序列的DFT 可以很简单的组合起来得到原序列的DFT 。例如，若N 为偶数，将原有的N 点序列分成两个（N/2）点序列，那么计算N 点DFT 将只需要约[(N/2)2 ·2]=N 2/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2代表必须完成两个DFT 。上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点的DFT （假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT 运算的情况。

3．基2按时间抽取（DIT ）的FFT 算法思想：

设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

将长度为2L N =的序列[](0,1,...,1)x n n N =-，先按n 的奇偶分成两组：

12[2][][21][]

x r x r x r x r =+=，r=0,1,…,N/2-1

DFT 化为：

1/21

/21

2(21)0

/21

/21

22120

0/21

/21

1/2

2/2

[]{[]}[][2][21][][][][]N N N nk rk r k

N

N

N

n r r N N rk k rk

N

N

N

r r N N rk k rk

N N

N r r X k DFT x n x n W

x r W

x r W x r W W

x r W x r W

W

x r W ---+===--==--=====

+

+=+=

+∑∑

∑

∑

∑

∑

∑

上式中利用了旋转因子的可约性，即：2/2rk

rk

N

N W W =。又令

/2

1

/21

11

/22

2/2

[][],[][]N N

rk

rk

N

N r r X k x r W

X k x r W

--===

=∑∑，则上式可以写成： 12[][][]k

N X k X k W X k =+(k=0,1,…,N/2-1)

可以看出，12[],[]X k X k 分别为从[]X k 中取出的N/2点偶数点和奇数点序列的N/2点DFT 值，所以，一个N 点序列的DFT 可以用两个N/2点序列的DFT 组合而成。但是，从上式可以看出，这样的组合仅表示出了[]X k 前N/2点的DFT 值，还需要继续利用12[],[]X k X k 表示[]X k 的后半段本算法推导才完整。利用旋转因子的周期性，有：(/2)

/2/2

rk

r k N N N W W +=，则后半段的DFT 值表达式：

/21/21()211/21/21

00

[][][][]2N N N r k rk N N r r N X k x r W x r W X k --+==+===∑∑，同样，22[][]2N X k X k += （k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT 值可以用前半段k 值表达式获得，中间还利用到()2

2N

N

k k k N

N

W W W W +==-，

得到后半段的[]X k 值表达式为：12[][][]k

N X k X k W X k =-（k=0,1,…,N/2-1）

。 这样，通过计算两个N/2点序列12[],[]x n x n 的N/2点DFT 12[],[]X k X k ，可以组合得到N 点序列的DFT 值[]X k ，其组合过程如下图所示：

1[]X k 12[][]k

N X k W X k +

2[]X k nk

N W -1 12[][]k N X k W X k -

比如，一个N = 8点的FFT 运算按照这种方法来计算FFT 可以用下面的流程图来表示：

x(0)

x(1)

x(2)

x(3)x(4)

x(5)

x(6)

x(7)

X(7)

X(6)X(5)X(4)X(3)X(2)

X(1)X(0)

4．基2按频率抽取（DIF ）的FFT 算法思想：

设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把[]X k 按k 的奇偶分组之前，把输入按n 的顺序分成前后两半：

1

/21

1

/2

/21

/21

()2

00

/21

2

[]{[]}[][][][][]2

[[][]],0,1,...,1

2

N N N nk

nk nk

N

N

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

n n N N k nk N N n X k DFT x n x n W

x n W

x n W N x n W x n W N x n x n W W k N ---===--+==-====

+

=+

+=

++=-∑∑

∑

∑

∑

∑

因为21N N

W

=-，则有2(1)N k k N

W

=-，所以：

/21

[][[](1)[]],0,1,...,12

N k

nk

N n N X k x n x n W k N -===

+-+=-∑

按k 的奇偶来讨论，k 为偶数时：

/21

20

[2][[][]],0,1,...,12

N rn N n N X r x n x n W k N -==

++

=-∑

k 为奇数时：/21

(21)0

[21][[][]],0,1,...,12

N r n

N n N X r x n x n W k N -+=+=-+=-∑

前面已经推导过2/2rk

rk N

N W W =，所以上面的两个等式可以写为：

/21

/20

[2][[][]],0,1,...,/212

N rn

N n N X r x n x n W r N -==

++

=-∑

/21

/20

[21]{[[][]]},0,1,...,/212

N n nr N N n N X r x n x n W W r N -=+=

-+

=-∑ 通过上面的推导，[]X k 的偶数点值[2]X r 和奇数点值[21]X r +分别可以由组合而成的N/2点的序列来求得，其中偶数点值[2]X r 为输入x[n]的前半段和后半段之和序列的N/2点DFT 值，奇数点值[21]X r +为输入x[n]的前半段和后半段之差再与n

N W 相乘序列的N/2点DFT 值。

令1[][][]2

N x n x n x n =++

，2[][[][]]2n

N N x n x n x n W =-+，则有：

/21

/21

1/2

2/20

[2][],[21][],0,1,...,

12

N N rn

rn

N N n n N

X r x n W

X r x n W r --===

+=

=-∑

∑

这样，也可以用两个N/2点DFT 来组合成一个N 点DFT ，组合过程如下图所示：

[]x n [][]2

N

x n x n ++

[]2N x n + -1 n N W [[][]]2

n

N N x n x n W -+

二．在FFT 计算中使用到的MATLAB 命令：

函数fft(x)可以计算R 点序列的R 点DFT 值；而fft(x,N)则计算R 点序列的N 点DFT ，若R>N ，则直接截取R 点DFT 的前N 点，若R

四、实验目的：

离散傅氏变换（DFT ）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT 变换到时域。FFT 是DFT 的一种快速算法。在数字信号处理系统中，FFT 作为一个非常重要的工具经常使用，甚至成为DSP 运算能力的一个考核因素。

本实验通过直接计算DFT ，利用FFT 算法思想计算DFT ，以及使用MATLAB 函数中的FFT 命令计算离散时间信号的频谱，以加深对离散信号的DFT 变换及FFT 算法的理解。

五、实验内容：

a) 计算实数序列5()cos

,025616

x n n n π

=≤≤的256点DFT 。 b) 计算周期为1kHz 的方波序列（占空比为50％，幅度取为＋/-512，采样

频率为25kHz ，取256点长度） 256点DFT 。

六、实验器材（设备、元器件）：

安装MATLAB 软件的PC 机一台，DSP 实验演示系统一套。

七、实验步骤：

（1） 先利用DFT 定义式，编程直接计算2个要求序列的DFT 值。

（2） 利用MATLAB 中提供的FFT 函数，计算2个要求序列的DFT 值。 （3） (拓展要求)不改变序列的点数，仅改变DFT 计算点数（如变为计

算1024点DFT 值），观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。

（4） 利用FFT 的基本思想（基2－DIT 或基2－DIF ），自己编写FFT 计

算函数，并用该函数计算要求序列的DFT 值。并对前面3个结果进行对比。

（5） (拓展要求)尝试对其他快速傅立叶变换算法（如Goertzel 算法）

进行MATLAB 编程实现，并用它来计算要求的序列的DFT 值。并与前面的结果进行对比。

（6） （拓展要求）在提供的DSP 实验板上演示要求的2种序列的FFT

算法（基2－DIT ），用示波器观察实际计算出来的频谱结果，并与理论结果对比。

八、实验数据及结果分析：

注：本次实验在寝室电脑上完成，所用MATLAB版本为MATLAB R2010b 程序：

（2）对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序序列生成的代码同（1）。这里给出FFT算法的程序：

（3）对要求的2种序列用MATLAB中提供的FFT函数进行计算的程

结果：

（1）对2种要求的序列直接进行DFT计算的频域波形

50

100

150200250

300

-1-0.500.5

1Time index n

A m p l i t u d e Sinusoidal, time-domain sequence

Frequency index k

M a g n i t u d e

Magnitude of DFT samples(direct method)

050

100150200250300

-1000

-50005001000

Time index n

A m p l i t u d e

Periodic square wave, time-domain sequence

4

Frequency index k

M a g n i t u d e

Magnitude of DFT samples(direct method)

（2） 对2种要求的序列进行基2－DIT 和基2－DIF FFT 算法频域波形

（3） 对2种要求的序列用MATLAB 中提供的FFT 函数计算的频域波形。

Frequency index k

M a g n i t u d e

Magnitude of DFT samples-Sinusoidal wave(fft method)

4

Frequency index k

M a g n i t u d e

Magnitude of DFT samples-Periodic square wave(fft method)

（4）（拓展要求）分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。

（5）（拓展要求）保持序列点数不变，改变DFT 计算点数（变为1024点），观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。

M a g n i t u d e

Frequency index k

1024 DFT samples of sinusoidal sequence

4

M a g n i t u d e

Frequency index k

1024 DFT samples of square wave sequence

可见，将DFT 计算点数由256点增加到1024点后，频谱变得更密了。这是由于“栅栏效应”——通过栅栏观察频谱，可能造成一些频谱丢失。增加DFT 计算点数相当于时域补零，使得谱线更密。但是补零并不能提高频谱分辨率，提高频谱分辨率的有效方法是增加采样信号的实际长度。

九、实验结论：

1. 直接用DFT 定义和用FFT 算法计算的DFT 结果一致，但是FFT 算法的效率更高。

2. 时域补零并不能提高频谱分辨率，有效的方法是增加采样信号的实际长度。

十、总结及心得体会：

通过本实验，加强了MATLAB 使用技能，练习了用定义和FFT 算法计算序列的DFT ；对于FFT 算法，用基2-DIT 实现，由于时间原因，没有做基2-DIF ，但基本思想一致，只是变成频域抽样；实验加深了我对FFT 算法思想的理解，学到了不少新知识、新技巧，将许多知识运用到了实践中。

十一、对本实验过程及方法、手段的改进建议：

1.建议在实验指导书中增加一些内容，如编程思想、码位倒序等，因为具体的编程实现和仅了解算法思想之间还是有一定差别的。

2.建议适当调整实验时间段，不要和考试时间段重合，让同学们有机会认真研究实验内容，并有时间和老师交流。

报告评分：

指导教师签字：