用matlab实现DFT FFT

k X1 k X1 k + WN X2 (k)
k X 2 (k)X1 k − WN X2 (k) 图 2 DIT-FFT 蝴蝶因子
k WN
下图为 DIT-FFT 与 DIF-FFT 蝴蝶因子的三种形式：
8
图 3 DIT-FFT 和 DIF-FFT 的蝴蝶因子 (1)原始形式； （2）简化形式； （3）优化形式 下面以 8 点的 DIT-FFT 为例：
用 matlab 实现 DFT FFT
目录
实验目的........................................................................................................................................... 2 实验内容........................................................................................................................................... 2 1.用 MATLAB 实现 DFT ................................................................................................................. 2 2.用 MATLAB 实现 FFT，分析有限离散序列的 FFT .................................................................. 3 3.通过分别计算时间，得出 DFT 与 FFT 的算法差异 .............................................................. 7 实验原理........................................................................................................................................... 7 1. 离散傅里叶变换的快速算法 FFT ......................................................................................... 7 2. FFT 提高运算速度的原理 ..................................................................................................... 8 3. 理论分析 DFT 与 FFT 算法差异 .......................................................................................... 10 实验步骤......................................................................................................................................... 11 实验结果......................................................................................................................................... 11 实验分析......................................................................................................................................... 21 实验结论......................................................................................................................................... 25 实验体会......................................................................................................................................... 26
N
(4)
将x1 (r)和x2 (r)分别进行 N/2 点的 DFT 得X1 (k)和X 2 (k)， 且：
k X k = X1 k + WN X2 (k)
X
N 2
k + k = X1 k − WN X2 (k)
r = 0,1, … ， 2 − 1 (5)
N
重复这一过程，可得到 x（n）的 FFT。
2. FFT 提高运算速度的原理
实验Fra Baidu bibliotek容
1.用 MATLAB 实现 DFT
N 点序列 x(n) 的 DFT 为：
������−1
X k =
������ =0
nk x n WN
0≤k≤N−1
DFT 的矩阵为：
根据 DFT 公式与矩阵展开，通过 MATLAB 实现 DFT：
2
2.用 Matlab 实现 FFT
编程思想及程序框图： 原位计算 因为 DIT-FFT 与 DIF-FFT 的算法类似，这里我们以 DIT-FFT 为例。N=2M 点的 FFT 共 进行 M 级运算，且每一级都由 N/2 个蝶形运算组成，后一级的节点数据由前一级同 处一条水平线位置的节点数据产生，所以我们同样可以将后一级的节点数据储存到 前一级的节点中， 这样的方法叫做原位计算， 它大大节省了内存资源， 降低了成本， 简化了运算。 序列的倒序 无论是进行 DIT-FFT 还是 DIF-FFT 都需要进行倒序，包括输入倒序与输出倒序，以 一定的方式将数组进行重新排列。 M 倒序的方法： 首先由于 N=2 ,我们就可以用 M 位二进制数来表示节点的顺序，并且 按照奇偶时域抽取。然后，如图 1 所示，第一次按最低位 n0 的 0、1 值分解为奇偶 组，第二次按次低位 n1 的 0、1 值分解为奇偶组，以此类推。最后，所得二进制数 所对应的十进制数即为序列倒序后产生的序列。
图 1 序列倒序过程
倒序的 MATLAB 方法： 用雷德算法可以对输入信号序列进行倒序重排，流程图如下所示：
3

蝴蝶因子的变化规律 在 DIT-FFT 中， 每一级都由 N/2 个蝶形运算构成， 每个蝶形运算包含一个蝴蝶因子， 每一级的蝶形因子又有一定的变化规律： 设 L 表示自左而右的运算级次（L=1,2,3,…,M） ，序数 R，次数 K。 每个蝶形运算的两个输入量相距 B=2^(L-1)个点。 假设 N=8， 则 M=3，这时有： L=1 时， B=1， S=N/2， R=0， K=1: N/2 则有 P=(K-1)*1，所以
图 4 8 点的 DIT-FFT 对于 8 点的 FFT，我们需要 M=log 2 8=3 阶运算，每一阶有四个蝴蝶因子，在每个蝴蝶因
9
子中需要 1 次复数乘法与两次复数加法，因为每个蝴蝶因子都一样，所以 FFT 通过重复 运算大大简化了算法与复杂度。
3. 理论分析 DFT 与 FFT 算法差异 DFT
7
WN
nk +N/2
nk = −WN
（3）
快速傅里叶变换算法正是基于这样的基本思想而发展起来的， 它的算法基本可以分称两 大类：时间抽取法（DIT-FFT）和频率抽取（DIF-FFT） 。由于 DIF-FFT 算法思想基本一 致，只是划分方式略有差异，所以这里以 DIT-FFT 算法为例进行说明。 当 N 是 2 的整数次方时，称为基 2 的 FFT 算法。 首先将序列 x(n)分解为两组，偶数项为一组，奇数项为一组： x 2r = x1 (r) x 2r + 1 = x2 (r) r = 0,1, … ， 2 − 1
实验原理
1. 离散傅里叶变换的快速算法 FFT
N 点序列 x(n) 的 DFT 为： X k =
������−1 ������ =0 nk x n WN
0≤k≤N−1
nk
（1）
nk 由于系数WN = e−j N n(N −k)
2π
是一个周期函数：
−nk = WN
WN
= WN
k(N −n)
（2）
且是对称的：
FFT 算法将长序列的 DFT 分解为短序列的 DFT。N 点的 DFT 先分解为 2 个 N/2 点的 DFT， 每个 N/2 点的 DFT 又分解为 N/4 点的 DFT， 以此类推。 最小变换的点数即所谓的 “基数” ， 这个基数是 2 点的 DFT（又叫蝴蝶因子） 。 如公式（5）所示，对于 X（k）的计算可以分为两部分，因为在上下两个式子中都出现 了X1 k 与X2 (k)，因此只需要一次复数乘法与两次复数加法即可分别得到 N/2 点的 DFT， 从而得到总的 N 点的 DFT。如图 2 所示。

FFT
k X k = X1 k + WN X2 (k) N k X + k = X1 k − WN X2 (k) 2
r = 0,1, … ，
N −1 2
FFT
算法复杂度的计算： 对于 N = 2M 点的 FFT,需要M = log 2 N阶的重复运算； 每一阶有 N/2 个蝴蝶因子； 每一个蝴蝶因子需要 1 此复数乘与 2 次复数加。 复数乘法 N 点 FFT N log 2 N 2 复数加法 N log 2 N
1
实验目的
1. 通过研究 DFT,FFT 性质，用语言实现 DFT, FFT。不使用 MATLAB 现有的 FFT 函数， 自己编写具体算法。 2. 掌握 FFT 基 2 时间抽选法，理解其提高减少乘法运算次数提高运算速度的原理。 3. 设计实验， 得出 DFT 和 FFT 算法差异的证明， 如复杂度等(精度、 不同长度的序列等） 。
P WN =WN/4 ， J=0 J
4
所以对于一般情况 N=2 ，第 L 级的旋转因子为 P=(K-1)*B。 从而得到 DIT-FFT 的 Matlab 流程图：
M
5
DIT-FFT 的 MATLAB 实现为（以对 x(n)=cos(
nπ ) 进行 6
16 点的变化为例） ：
6
3.通过分别计算时间，得出 DFT 与 FFT 的算法差异：
P WN =WN ， J
J=0,1,2,3
L=2 时， B=2， S=N/4， R=0: N/2:N-1， K=1: N/4 则有 P=(K-1)*2，所以
P WN =WN/2 ， J=0,2 J
L=3 时， B=4， S=N/8， R=0: N/4:N-1， K=1: N/8 则有 P=(K-1)*4，所以
������−1
X k =
������ =0
nk x n WN
0≤k≤ N−1
DFT 算法复杂度的计算： 复数乘法 一次 X(k) N 次 X(k) (N 点 DFT) N 复数加法 N-1
N
2
N(N-1)
对于 N= 512、1024 和 8192 的 DFT，分别需要复数乘法（次） 2 2 18 N = 512 = 2 = 262144 2 2 20 N = 1024 = 2 = 1048576 2 2 26 N = 8192 = 2 = 67108864 这是一个巨大的数字，很难在实际中应用。
法一： 对 N=512，1024，2048 和 4096 点的离散时间信号 x(n)，用 Matlab 语言编程分别以 DFT 和 FFT 计算 N 个频率样值 X(k)，比较两者所用时间的大小。设 x(n)=cos(
nπ )。 6
利用 TIC,TOC 函数，分别跟踪变换时间。 法二： 利用 clock 返回当前日期向量的时间， 通过 etime （） 函数返回走过的日期向量的时间。 用 MATLAB 实现如下：