课程设计任务书
学生姓名：专业班级：
指导教师：工作单位：
题目：8点基于DIT的FFT的实现
初始条件：
具备Matlab编程能力；
熟悉基于DIT的FFT的实现原理；
提供编程所需要的计算机一台。
要求完成的主要任务：（包括课程设计工作量及其技术要求，以及说明书撰写等具体要求）
1、编写一个8点的基于DIT的FFT函数，不能使用matlab自带的FFT实现函数；
2、并调用该函数实现16点的FFT运算，用matlab自带函数对运行结果结果进行验证；
3、完成符合学校要求的设计说明书。
时间安排：
一周，其中3天程序设计，2天程序调试
指导教师签名：年月日
系主任（或责任教师）签名：年月日
目录
摘要................................................................................................................................................ I 1 概述 (1)
1.1 快速傅立叶变换(FFT)简介 (1)
1.2 MATLAB简介 (2)
2 直接计算DFT的问题及改进 (3)
2.1直接计算DFT的运算量 (3)
2.2 改进措施 (4)
3 按时间抽选的基-2FFT算法（DIT-FFT) (5)
3.1 DIT-FFT算法原理 (5)
3.2 DIT-FFT的运算量 (11)
3.3 DIT-FFT算法的特点 (12)
3.4 N=16时的DIT-FFT算法 (14)
4 MATLAB程序代码 (15)
4.1 N=8点DIT-FFT代码 (15)
4.2 N=16点DIT-FFT代码 (16)
5 MATLAB仿真结果及验证 (17)
5.1 DIT-FFT函数调试 (17)
5.2 DIT-FFT函数运行结果 (18)
5.3调用系统函数验证 (19)
6 心得体会 (21)
参考文献 (22)
摘要
此次课设目的是利用MATLAB实现8点基于DIT的FFT的仿真，不使用MATLAB 自带的FFT实现函数。本文先就直接计算傅立叶变换（DFT）存在的问题进行讨论，之后详细介绍了快速傅立叶变换（FFT）的原理以及推导过程，给出了8点FFT的蝶形流图以及MATLAB仿真的程序代码，并通过调用该函数代码计算16点的FFT。最后给出了仿真调试结果和此次课设的总结。
关键词：FFT；MATLAB；仿真
Abstract
The aim of this Course Design is to use MATLAB to achieve 8-point DIT-FFT simulation, and can not use the built-in MATLAB FFT function to realize. The beginning of this article discuss the problems of direct calculation of the Fourier transform (DFT) , and then introduces the principle of Fast Fourier Transform (FFT) and the process of derivation. Then there is given butterfly flow diagram of 8-point FFT and the MATLAB simulation program code, and realize 16-point FFT calculation by calling the function code. Finally, enumerate the simulation results and make the summary of this curriculum design.
Keywords: FFT; MATLAB; Simulation
1 概述
1.1 快速傅立叶变换(FFT)简介
傅立叶变换，表示能将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。傅立叶变换是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。傅立叶变换是声学、语音、电信和信号处理等领域中一种重要的分析工具。在不同的研究领域，傅立叶变换具有多种不同的变体形式，如连续傅立叶变换和离散傅立叶变换。
离散傅立叶变换（DFT），是傅立叶变换在时域和频域上都呈现离散的形式，将时域信号的采样变换为在离散时间傅立叶变换（DTFT）频域的采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作经过周期延拓成为周期信号再作变换。有限长序列可以通过离散傅立叶变换(DFT）将其频域也离散化成有限长序列。但其计算量太大，很难实时地处理问题，因此引出了快速傅立叶变换(FFT)。
1965年，Cooley和Tukey提出了计算离散傅立叶变换（DFT）的快速算法，将DFT的运算量减少了几个数量级。从此，对快速傅立叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。FFT在离散傅立叶反变换、线性卷积和线性相关等方面也有重要应用。
计算离散傅立叶变换的快速方法，有按时间抽取的FFT算法（DIT-FFT）和按频率抽取的FFT算法（DIF-FFT)。前者是将时域信号序列按偶奇分排，后者是将频域信号序列按偶奇分排。它们都借助于的两个特点：一是周期性；二是对称性。这样，便可以把离散傅立叶变换的计算分成若干步进行，计算效率大为提高。快速傅立叶变换（FFT），是离散傅立叶变换的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
1.2 MATLAB简介
MATLAB是美国Math Works公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink 两大部分。
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国Math Works公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。
MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JA V A的支持。
优势特点:
1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来；
2) 具有完备的图形处理功能，实现计算结果和编程的可视化；
3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握；
4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。
2 直接计算DFT 的问题及改进
DFT 在数字信号处理中有着重要的作用。然而直接计算DFT 的运算量非常大，它与序列长度的平方成正比，因此制约了DFT 的应用。快速傅立叶变换（Fast Fourier Transform,FFT ）是实现DFT 的一种快速算法。
2.1 直接计算DFT 的运算量
离散傅立叶变换对为：
正变换： ∑-===1
0)()]([)(N n nk
N W n x n x DFT k X , 1...,2,1,0-=N n （2.1.1）
反变换： W nk N N k k X N k X IDFT n x --=∑==10
)(1)]([)( , 1,...,2,1,0-=N k （2.1.2）
其中，N W 为旋转因子且N
j N e
W π2-=，x(n)表示输入的离散数字信号序列，X(k)为输入序
列x(n)对应的N 个离散频率点的相对幅度。
比较两式，正变换与反变换的差别就在于旋转因子的指数差个负号及少一个比例因子。因此DFT 与IDFT 的计算量极为相似，所以只需以正变换为例来考虑直接计算DFT 时所存在的问题。
一般情况下，序列N W 、x(n)及其离散频谱X(k)都是复序列，因此,要计算离散频谱X(k)的一个值就需要N 次复数乘法与N-1次复数加法运算，而计算一个完整的N 点的X(k)（对应1,...,2,1,0-=N k ），就需要2
N 次复数乘法与N(N-1)次复数加法运算，当N 很大时，
2≈1)-N(N N ，因此直接计算DFT 的运算量就几乎与2N 成正比，随着N 的增加，运算量将急剧增大，即使采用计算机，也很难实时处理。
当然，以上的分析的DFT 计算量与实际的运算量稍有出入，例如，计算10
=W N 时就不需要做乘法运算，但是当N 很大时，这种情况对整个DFT 的计算量影响很小，一般不做特别统计。
2.2 改进措施
FFT 主要利用DFT 旋转因子的周期性与对称性来减少运算量：
周期性： W W W n
k N N
k
N n N
nk
N )()(++== （2.2.1）
对称性： W W n N k N nk nk N
)
(-N W *)(-== （2.2.2） 利用周期性与对称性，一方面可以在DFT 的运算过程中把有些项进行合并，另一方面可以把长序列的DFT 分解成若干短序列的DFT 。因为DFT 的运算量与变换长度的平方成正比，如果可以把一个长序列DFT 分解成若干短序列DFT 再进行计算，就可以大大减少运算量。
3 按时间抽选的基-2FFT 算法（DIT-FFT)
常用的FFT 算法有两大类，一类是按时间抽取的FFT 算法（简称DIT-FFT ），另一类是按频率抽取的FFT 算法（简称DIF-FFT ）。
最早提出的基-2FFT 算法，使DFT 的运算效率提高了1~2个数量级，从而为DFT 由理论研究到实际应用创造了条件。
3.1 DIT-FFT 算法原理
按时间抽取的FFT 算法基本思想是：时域下的序列x(n)按序列n 的奇偶分组，频域下序列X(k)按序列k 前后分组。有限长序列x(n)，设其长度2L
N =，L 为整数，若不满足该条件，加零补足。显然N 为偶数，可以按序列序号分为奇、偶两组序列，长度分别为N/2，如：
)}1(),...,1(),0({)(-=N x x x n x
按n 的奇偶分组，对x(n)重新排列，得：
2)}-x(N x(2),...,｜x(0),)1(),...,3(),1({-N x x x
令
12/,...,2,1,0)},1(),...,2(),0({)2(-=-=N r N x x x r x （3.1.1） 12/,...,2,1,0)},1(),...,3(),1({)12(-=-=+N r N x x x r x （3.1.2）
再令
)2()(1r x n x = （3.1.3） )12()(2+=r x n x （3.1.4）
N 点序列x(n)的DFT 为：
W W W W
W
W W W rk
N
N r k N
rk
N N r k r N N r rk N
N r nk
N
n nk
N n nk
N
N n r x r x r x r x n x n x n x n x DFT k X 212/0
221
2/0
1)12(12/0
212/0
1
0)()()12()2()()()()]
([)(∑∑∑∑∑∑∑-=-=+-=-=-=+=
++
=+
===为奇数
为偶数
因为
W W rk
N rk N j
rk N
j rk
N e
e 2/2
/2222===--π
π （3.1.5）
所以
)
()()()()(211
2/0
2
/21
2/0
2/1k X k X r x r x k X W W W W k
N N r rk
N rk
N
N r rk N +=+=
∑
∑
-=-= （3.1.6）
其中)(1k X 与)(2k X 分别是)(1n x 与)(2n x 的N/2点的DFT ，即：
1
2
,...2,1,0,)()]([)(12
,...
2,1,0,)()]([)(1
2/0
2/2221
2/0
2/111-==
=-==
=∑
∑
-=-=N k r x n x DFT k X N k r x n x DFT k X N r rk
N N r rk
N W W （3.1.7）
式（3.1.6）把一个N 点的DFT 分解成了两个N/2点的DFT 的组合，但是)(1k X 与)(2k X 分别是)(1n x 与)(2n x 的N/2点的DFT ，r 与k 的取值满足r,k=0,1,2,...,N/2-1，而X(k)是一个N 点的DFT ，因此式（3.1.7）只计算X （k)前N/2点的值。利用旋转因子的周期性来计算 X(k)后N/2的值。由式（2.1.1）得： W W N
k r N nk
N )2
(2/2/+
= （3.1.8）
这样X(k)后N/2点的值为：
)2
()2()2(2
22/1N k X N k X N k X W N
k N +++=++
考虑周期性：
)
()()()2(11
2/0
2
/112/0
)
2
(2/11k X r x r x N k X N r rk
N N r N
k r N W W ==
=+∑
∑-=-=+ （3.1.9）
同理可得：
)()2
(22k X N
k X =+
（3.1.10） 式（3.1.9）与式（3.1.10）表明)(1k X 与)(2k X 前N/2点的值与后N/2的值相同，实际上就是DFT 隐含的周期性。 再考虑对称性：
W W k
N N
k
k N
N -==+W W 2N
N )2( （3.1.11） 所以：
)
()()
2()2()2(212)
2(2/1k X k X N k X N k X N k X W W k
N N
k N -=+++=++ （3.1.12） 这样利用式（3.1.11）与式（3.1.12）就可以把一个完整的N 点的DFT 分解成两个N/2点的DFT 来运算。上述讨论的运算过程可以用图3.1所示的信号流图来表示。 )(1k X )()(21k X k X W k
N + )(2k X )()(21k X k X W k N -
图3.1 DIT-FFT 蝶形运算符号
如图3.1所示，任一支路上的值等于支路起始节点的值乘以支路传输系数，任一节点上的值等于所有输入支路值之和。从图中可以看出，每个蝶形的运算需要一次复数乘法运算和两次复数加法运算。
在图3.1中，输入为两个N/2点的DFT 输出为一个N 点的DFT 结果，输入输出点数一致。运用这种表示方法，8点的DFT 可以用图3.2来表示。
当N=8为例，采用蝶形图表示，DFT 的分解运算如图3.2所示。其中)3()0(X X →由式（3.1.7）给出，)7()4(X X →由式（3.1.12）给出。
根据式（3.1.11）与式（3.1.12），一个N 点的DFT 可以由两个N/2点的DFT 运算构成，再结合图3.1的蝶形信号流图可以得到图3.2的8点DFT 的第一次分解。该分解可以用以下几个步骤来描述：
1.将N 点的输入序列按奇偶分为2组分别为N/2点的序列；
2.分别对1中的每组序列进行DFT 变换得到两组点数为N/2的DFT 变换值X1和X2；
3.按照蝶形信号流图将2的结果组合为一个N 点的DFT 变换结果。
经过一次分解，计算一个完整的N 点DFT 需要2/2/)2/(222N N N ≈+次复数乘法，以及2/)12/(2N N N N =+-次复数加法，运算量减少了将近一半，由于N/2依然是偶数，故可将N/2点的DFT 按同样方法分解成两个N/4点的DFT 。
与第一次分解相同，把序列)(1r x 按r 的奇偶次序分解成两个N/4点的序列)(3l x 与)(4l x ，即：
)
0(x )2(x )4(x )
6(x )1(x )3(x )5(x )
7(x DFT
N 点2/DFT
N 点2/)
0(1X )1(1X )
2(1X )
3(1X )
0(2X )1(2X )2(2X )
3(2X 0N
W 1N
W 2N
W
3N
W 1
-1-1-1
-)0(X )1(X )
2(X )
3(X )4(X )5(X )6(X )
7(X 图3.2 DIT-FFT 一次分解
14
,...,2,1,0,)()12()
()2(4131-=??
?=+=N l l x l x l x l x （3.1.13） 代入)(1r x 的N/2点的DFT 表达式中有：
)
()()()()12()2()()]
([)(42/31
4/0
4
/42
/1
4/0
4/31
4/0)12(2/11
4/0
22/11
2/0
2/111k X k X l x l x l x l x r x r x DFT k X W W W W W W W k
N N l lk
N k
N N l lk N N l k
l N N l lk
N N r rk
N +=+=
++
==
=∑
∑
∑
∑
∑
-=-=-=+-=-= （3.1.14）
且由对称性和周期性有：
)()()4
(42/31k X k X N k X W k
N -=+ （3.1.15） 其中：
14
,...
2,1,0,)()]([)(14,...
2,1,0,)()]([)(1
4/0
4/4441
4/04/333-==
=-==
=∑∑
-=-=N k l x l x DFT k X N k l x l x DFT k X N l rk
N N l rk
N W W （3.1.16）
所以，当N=8时，上一步分解出的一个N/2点的DFT 可以分解成两个N/4点的DFT ，运算流图如图3.3所示。
x (0) N/4点 )0(3X )0(1X DFT
x(4) )1(3X )1(1X
x(2) N/4点 )0(4X W N 1
2/ -1 )2(1X DFT
x(6) )1(4X W N 12/ -1 )3(1X
图3.3 由两个N/4点DFT 组合成1个N/2点DFT
同理，)k (2X 也做同样的分解，得：
)()()(62/52k X k X k X W k
N += （3.1.17） )()()4
(62/52k X k X N k X W k
N -=+ （3.1.18） 其中：
14
,...
2,1,0,)()]([)(14,...2,1,0,)()]([)(1
4/0
4/6661
4/04/555-==
=-==
=∑
∑
-=-=N k l x l x DFT k X N k l x l x DFT k X N l rk
N N l rk
N W W （3.1.19）
这样就可以把两个N/2点的DFT 分解成四个N/4点的DFT 。当N=8时，经过两次分解后的运算流图如图3.4所示。
图3.4 DIT-FFT 二次分解
由于满足L N 2=，所以经过二次分解后N/4仍然是偶数，可以继续分解。根据上述讨论，经过L-1次分解后，就可以把一个N 点的DFT 分解成N/2个两点的DFT ，而每一个两点的DFT 可以根据如图3.1所示的蝶形图来进行计算。当N=8时，于是就可以把一个完整的8点的DFT 分解成4个两点的DFT 进行运算。计算一个完整的N=8点的DIT-FFT 蝶形图
如图3.5所示。
图3.5 N=8 DIT-DFT 的蝶形图
3.2 DIT-FFT 的运算量
观察图3.5，可见当L N 2=时，经过L-1次分解，整个DIT-FFT 运算有L 级蝶形，每一级蝶形有N/2个蝶形运算，每一个蝶形运算有一次复数乘法好人两次复数加法，所以整个运算流图的运算量为：
复数乘法：
bN N
N L m F 12
12=??
= (3.2.1) 复数加法：
bN N N
L a F 122
=??
= (3.2.2) 直接计算DFT 需要2N 次复数乘法与N(N-1)次复数加法运算，直接计算DFT 与DIT-FFT 的复数乘法运算量之比为：
)
0(x )
3(x )5(x )1(x )
4(x )
2(x )
6(x )7(x 0N
W
1N
W 2N
W
3N
W
N
W
0N
W
0N
W
0N
W
1
-1
-1-1
-1-1
-1-1-1
-1
-1
-1
-0N
W
0N W
2N
W
2N
W )
0(X )
1(X )2(X )3(X )4(X )
7(X )
6(X )5(X
bN N
bN N N L N N 1212
222=
= (3.2.3) 3.3 DIT-FFT 算法的特点
1. 原位运算
在DIT-FFT 的蝶形图中，取第m 级且两输入节点分别在第k 、j 行的蝶形为例，讨论DIT-FFT 的原位运算规律。如图3.6所示，蝶形运算的关系式可表示为：
??
??
?-=+=----)()()()()()(1111j X k X j X j X k X k X m r N m m m r
N m m W W (3.2.4)
从式中可以看出，第m 级蝶形第k 行与第j 行的输出，只与第m-1级的第k 行与第j 行的输出有关，换言之，第m-1级的第k 行与第j 行的输出)(1k X m -与)(1j X m -在运算流图中的作用就是用来计算第m 级的第k 行与第j 行的输出)(k X m 与)(j X m 。这样当计算完)(k X m 与)(j X m 后，)(1k X m -与)(1j X m -在运算流图不再起作用，因此就可以把)(k X m 与)(j X m 直接存放在原来存放)(1k X m -与)(1j X m -的存储单元中。同理，可以把第m 级蝶形的N 个输出值直接存放在第m-1级蝶形输出的N 个存储单元中，这样从第一级的输入x （n ）开始，到最后一级输出X(k)，只需要N 个存储单元。
)(1k X )()(21k X k X W k
N +
)(2k X )()(21k X k X W k N - W r
N -1
图3.6 按时间抽取蝶形运算结构
2. 倒序规律
从图3.5可以看出，按原位计算时，蝶形图的输出正好是自然顺序X(0),X(1),...,X(7)，但是输入却不是自然顺序，而是x(0),x(4),x(6)...，表面看起来好像是“混乱无序”的，实际上是有规律的，即倒序的排列方法。
倒序的形成原因是FFT 不断对序列进行奇偶分组造成的，重新排列了序列的存放顺序，
因此它是按码位倒置顺序排放的。由于N=2M ，所以倒序数可用M 位二进制数表示。第一次分组，按0n 的0和1分成奇偶两组，0n =0相当是偶序列，0n =1相当是奇序列。第二次分组按1n 的0和1分成两组……依次类推，直到M 次分组。
0 0
1 0 0 1 1
0 0 ()2210n n n 1 1 0 1 1
图3.7 倒序形成的树状结构
3. 蝶形运算两节点之间的“距离”
从图3.5可以看出，第一级蝶形每个蝶形运算两节点之间的“距离”为1，第二级每个蝶形运算两节点之间的“距离”为2，第三级蝶形每个蝶形运算两节点之间的“距离”为4。依次类推，对于L N 2=的DIT-FFT ，可以得到第M 级蝶形每一个蝶形运算两节点之间的“距离”为12-M 。
4. 旋转因子的变化规律
以图3.5的8点FFT 为例，每一级的每一个蝶形运算都要乘以一个旋转因子W r
N ，r 是旋转因子的指数。在第一级蝶形，r=0；在第二级蝶形，r=0,1；在第三级的蝶形，r=0,1,2,3;依次类推，对于第M 级蝶形，旋转因子的指数为 ：
12,...,2,1,0,21-=?=--M M L J J r
这样就可以算出每一级的旋转因子。对于M 级的任一蝶形运算所对应的旋转因子的指数，可以 如下方法得到：
(1) 将待求的蝶形输入节点中上面节点的行标号值k 写成L 位二进制数；
(2) 将此二进制数乘以M L -2，即将L 位二进制数左移L-M 位，右边的空位补零，然后
从低位到高位截取L 位，即得到要求的指数r 所对应的二进制数。
3.4 N=16时的DIT-FFT 算法
先将序列x(n)奇偶分组得：
)2()(1r x r x =
)12()(2+=r x r x r=0,1,2,...,7 (3.4.1) 将DFT 运算也分为两组：
)
()()()()12()2()()()()]([)(21618
7
21687
1)12(12/021
2/01
0k X k X r x r x r x r x n x n x n x n x DFT k X W W W W W
W
W W W k
rk
r k
rk r k r N
N r rk N
N r nk
N
n nk
N n nk
N
N n +=+=++
=
+
=
==∑∑∑∑∑∑∑==+-=-=-=为奇数
为偶数 (3.4.2)
其中)(1k X 与)(2k X 分别是)(1n x 与)(2n x 的8点DFT ，即：
7
,...,2,1,0,)()]([)(7
,...,2,1,0,)()]([)(7
82227
08111======∑∑==k r x n x DFT k X k r x n x DFT k X r rk
r rk
W W (3.4.3)
这样，一个16点的DFT 就被分解为两个8点的DFT ，即：
)
()()8()()()(21612161k X k X k X k X k X k X W W k
k
-=++=，k=0,1,2,...,7 (3.4.4)
对两个8点的DFT 再分别做进一步分解，每个8点的DFT 分解成两个4点的DFT ，即：
)
()()4()()()(421631421631k X k X k X k X k X k X W W k
k
-=++=，k=0,1,2,...,7 (3.4.5)
对四个4点的DFT 再分别做进一步分解，每个4点的DFT 分解成两个2点的DFT ，即：
)
()()2()()()(641653641653k X k X k X k X k X k X W W k
k
-=++=，k=0,1,2,...,7 (3.4.6)
所以，N=16点的DFT 最终可以分解成8个2点的DFT 。
4 MATLAB 程序代码
4.1 N=8点DIT-FFT 代码
根据蝶形图一级一级依次分解，可以把一个完整的8点的DFT 分解成4个两点的DFT 进行运算。
function y=fft8(x) %定义根据蝶形图计算8点DIT-FFT 的函数 Wn=exp(-j*2*pi/8); %定义旋转因子计算公式N
j
N e W π2-=
x1(1)=x(1)+x(5); %计算第一级蝶形图输出x(5)x(1)=(1)x 0
1W N +
x1(2)=x(1)-x(5); %x(5)-x(1)=(2)x 01W N x1(3)=x(3)+x(7); %x(7)x(3)=(3)x 0
1W N + x1(4)=x(3)-x(7); %x(7)-x(3)=(4)x 0
1W N x1(5)=x(2)+x(6); %x(6)x(2)=(5)x 0
1W N + x1(6)=x(2)-x(6); %x(6)-x(2)=(6)x 0
1W N x1(7)=x(4)+x(8); %x(8)x(4)=(7)x 0
1W N + x1(8)=x(4)-x(8); %x(8)-x(4)=(8)x 01W N
x2(1)=x1(1)+x1(3); %计算第二级蝶形图输出(3)x (1)x =(1)x 10
12W N + x2(2)=x1(2)+x1(4)*(Wn.^2); %(4)x (2)x =(2)x 1212W N + x2(3)=x1(1)-x1(3); %(3)x -(1)x =(3)x 1012W N x2(4)=x1(2)-x1(4)*(Wn.^2); %(4)x -(2)x =(4)x 1212W N x2(5)=x1(5)+x1(7); %(7)x (5)x =(5)x 10
12W N + x2(6)=x1(6)+x1(8)*(Wn.^2); %(8)x (6)x =(6)x 1212W N + x2(7)=x1(5)-x1(7); %(7)x -(5)x =(7)x 1012W N x2(8)=x1(6)-x1(8)*(Wn.^2); %(8)x -(6)x =(8)x 1212W N
y(1)=x2(1)+x2(5); %计算第三级蝶形图输出(5)x (1)x =y(1)20
2W N +
y(2)=x2(2)+x2(6)*(Wn.^1); %(6)x (2)x =y(2)21
2W N + y(3)=x2(3)+x2(7)*(Wn.^2); %(7)x (3)x =y(3)22
2W N + y(4)=x2(4)+x2(8)*(Wn.^3); %(8)x (4)x =y(4)232W N + y(5)=x2(1)-x2(5); %(5)x -(1)x =y(5)202W N y(6)=x2(2)-x2(6)*(Wn.^1); %(6)x -(2)x =y(6)212W N y(7)=x2(3)-x2(7)*(Wn.^2); %(7)x -(3)x =y(7)222W N y(8)=x2(4)-x2(8)*(Wn.^3); %(8)x -(4)x =y(8)232W N
4.2 N=16点DIT-FFT 代码
调用计算N=8点时编写的函数fft8来实现N=16点的DIT-FFT 运算。可以先将N=16点的DFT 根据奇偶分解成两个8点的DFT ，然后再分别用函数fft8对这两组8点的DFT 进行计算。
function y=fft16(x) %定义计算16点DIT-FFT 的函数 Wn=exp(-j*2*pi/16); %定义旋转因子计算公式N
j N e W π2-=
y1=fft8(x(1:2:16)); %计算偶数组的8点FFT y2=fft8(x(2:2:16)); %计算奇数组的8点FFT y(1:8)=y1+y2.*(Wn.^[0:7]); %计算前八个点 y(9:16)=y1-y2.*(Wn.^[0:7]); %计算后八个点
DSP 课程作业 用C 语言编写FFT 程序 1，快速傅里叶变换FFT 简介 快速傅氏变换（FFT ），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 我们假设 x(n)为N 项的复数序列，由DFT 变换，任一X （m ）的计算都需要N 次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的X （m ）,即N 点DFT 变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT 中，利用WN 的周期性和对称性，把一个N 项序列（设N=2k,k 为正整数），分为两个N/2项的子序列，每个N/2点DFT 变换需要（N/2）2次运算，再用N 次运算把两个N/2点的DFT 变换组合成一个N 点的DFT 变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要Nlog2N 次的运算，N 在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT 的优越性。 2，FFT 算法的基本原理 FFT 算法的基本思想：利用DFT 系数的特性，合并DFT 运算中的某些项，吧长序列的DFT —>短序列的DFT ，从而减少其运算量。 FFT 算法分类:时间抽选法DIT: Decimation-In-Time ；频率抽选法DIF: Decimation-In-Frequency 按时间抽选的基-2FFT 算法 1、算法原理 设序列点数 N = 2L ，L 为整数。 若不满足，则补零。N 为2的整数幂的FFT 算法称基-2FFT 算法。将序列x(n)按n 的奇偶分成两组： 则x(n)的DFT: ()()()()12221x r x r x r x r =+=0,1,...,12N r =-()()()() 111000N N N nk nk nk N N N n n n X k x n W x n W x n W ---=====+∑∑∑
创作编号： GB8878185555334563BT9125XW 创作者： 凤呜大王* 《测试信号分析及处理》课程作业 快速傅里叶变换 一、程序设计思路 快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M 级，其中N M 2log =；在输入序列()i x 中是按码位倒序排列的，输出序列()k X 是按顺序排列；每级包含 2N 个蝶形单元，第i 级有i N 2 个群，每个群有12-i 个蝶形单元； 每个蝶形单元都包含乘r N W 和r N W -系数 的运算，每个蝶形单元数据的间隔为12-i ，i 为第i 级； 同一级中各个群的系数W 分布规律完全相同。 将输入序列()i x 按码位倒序排列时，用到的是倒序算法——雷德算法。 自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。 若已知某数的倒序数是J ，求下一个倒序数，应先判断J 的最高位是 否为0，与2 N k =进行比较即可得到结果。如果J k >，说明最高位为0， 应把其变成1，即2 N J +，这样就得到倒序数了。如果J k ≤，即J 的最高 位为1，将最高位化为0，即2N J -，再判断次高位；与4N k =进行比较， 若为0，将其变位1，即4 N J +，即得到倒序数，如果次高位为1，将其化 为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。
注：因为0的倒序数为0，所以可从1开始进行求解。 二、程序设计框图 （1）倒序算法——雷德算法流程图 (2)FFT算法流程
实验三 用FFT 对信号进行频谱分析 一 实验目的 1 能够熟练掌握快速离散傅立叶变换的原理及应用FFT 进行频谱分析的基本方法； 2了解用FFT 进行频谱分析可能出现的分析误差及其原因； 二 实验原理 1.用DFT 对非周期序列进行谱分析 单位圆上的Z 变换就是序列的傅里叶变换，即 ()() j j z e X e X z ω ω== （3-1） ()j X e ω是ω的连续周期函数。对序列()x n 进行N 点DFT 得到()X k ，则()X k 是在区间 []0,2π上对()j X e ω的N 点等间隔采样，频谱分辨率就是采样间隔 2N π。因此序列的傅里叶变换可利用DFT （即FFT ）来计算。 用FFT 对序列进行谱分析的误差主要来自于用FFT 作频谱分析时，得到的是离散谱，而非周期序列的频谱是连续谱，只有当N 较大时，离散谱的包络才能逼近连续谱，因此N 要适当选择大一些。 2.用DFT 对周期序列进行谱分析 已知周期为N 的离散序列)(n x ，它的离散傅里叶级数DFS 分别由式（3-2）和（3-3） 给出： DFS ： ∑-=-=1 2)(1N n kn N j k e n x N a π ， n =0,1,2,…,N -1 （3-2） IDFS ： ∑-== 10 2)(N k kn N j k e a n x π ， n =0,1,2,…,N -1 （3-3） 对于长度为N 的有限长序列x (n )的DFT 对表达式分别由式（3-4）和（3-5）给出： DFT ： ∑-=-= 10 2)()(N n kn N j e n x k X π ， n =0,1,2,…,N -1 （3-4） IDFT ： ∑-==1 2)(1)(N k kn N j e k X N n x π ， n =0,1,2,…,N -1 （3-5） FFT 为离散傅里叶变换DFT 的快速算法，对于周期为N 的离散序列x (n )的频谱分析便可
按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。 按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 [ 对N=M 2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p W N ，p 称为旋转因子的指数。例如N ＝8 ＝3 2 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：p W N =J /4W N =J 2L W J=0
第二级：L=2，有2个旋转因子：p W N =J /2W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N ＝M 2的一般情况，第L 级共有1-L 2个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1-L 2 －1 L 2=M 2×M -L 2= N ·M -L 2 故： 按照上面两式可以确定第L 级运算的旋转因子 \ 3、同一级中，同一旋转因子对应蝶形数目 第L 级FFT 运算中，同一旋转因子用在L -M 2 个蝶形中； 4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L 级中，蝶距：D=L 2； 5、同一蝶形运算两输入数据的距离 在输入倒序，输出原序的FFT 变换中，第L 级的每一个蝶形的2个输入数据相距：B=1 -L 2。 6、码位颠倒 输入序列x(n)经过M 级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。 '
FFT相关原理及使用注意事项 在信号分析与处理中，频谱分析是重要的工具。FFT（Fast Fourier Transform，快速傅立叶变换）可以将时域信号转换至频域，以获得信号的频率结构、幅度、相位等信息。该算法在理工科课程中都有介绍，众多的仪器或软件亦集成此功能。FFT实用且高效，相关原理与使用注意事项也值得好好学习。 一、何为FFT 对于模拟信号的频谱分析，首先得使用ADC（模拟数字转换器）进行采样，转换为有限序列，其非零值长度为N，经DFT（离散傅立叶变换）即可转化为频域。DFT变换式为： 在上式中，N点序列的DFT需要进行N2次复数乘法和次复数加法，运算量大。 FFT是DFT的快速算法，利用DFT运算中的对称性与周期性，将长序列DFT分解为短序列DFT 之和。最终运算量明显减少，使得FFT应用更加广泛。 FFT基于一个基本理论：任何连续的波形，都可以分解为不同频率的正弦波形的叠加。FFT将采样得到的原始信号，转化此信号所包含的正弦波信号的频率、幅度、相位，为信号分析提供一个创新视觉。 例如在日常生活中有使用到的AM（Amplitude Modulation，幅度调制）广播，其原理是将人的声音（频率约20Hz至20kHz，称为调制波）调制到500kHz～1500kHz正弦波上（称为载波）中，载波的幅度随调制波的幅度变化。声音经这样调制后，可以传播得更远。在AM的时域波形（波形电压随时间的变化曲线），载波与调制波特征不易体现，而在FFT后的幅频曲线中则一目了然。如下图为1000kHz载波、10kHz调制波的AM调制信号，时域信号经FFT后其频率能量出现在990kHz、1.01MHz频率处，符合理论计算。 图 1 调制波10kHz、载波1000kHz的AM时域与频域曲线 二、FFT相关知识 现实生活中的模拟信号，大多都是连续复杂的，其频谱分量十分丰富。正如在数学中常量π，其真实值是个无理数。当用3.14来替代π时，计算值与真实值就会有偏差。在使用FFT这个工具时，受限于采样时的频率Fs、采样点长度N、ADC的分辨率n bit等因素的制约，所得到的信息会有所缺失与混淆。 1.奈奎斯特区与波形混叠 FFT分析结果中，存在一个那奈奎斯特区的概念，其宽度为采样率的一半Fs/2，信号频
matlab的FFT函数 相关语法： Y = fft(X) Y = fft(X,n) Y = fft(X,[],dim) Y = fft(X,n,dim) 定义如下： 相关的一个例子： Fs = 1000; % 采样频率 T = 1/Fs; % 采样时间 L = 1000; % 总的采样点数 t = (0:L-1)*T; % 时间序列（时间轴） %产生一个幅值为0.7频率为50HZ正弦+另外一个信号的幅值为1频率为120Hz 的正弦信号 x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); y = x + 2*randn(size(t)); % 混入噪声信号 plot(Fs*t(1:50),y(1:50)) %画出前50个点 title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % 求得最接近总采样点的2^n,这里应该是2^10=1024 Y = fft(y,NFFT)/L; %进行fft变换（除以总采样点数，是为了后面精确看出原始信号幅值） f = Fs/2*linspace(0,1,NFFT/2+1);%频率轴(只画到Fs/2即可，由于y为实数，后面一半是对称的) % 画出频率幅度图形，可以看出50Hz幅值大概0.7,120Hz幅值大概为1. plot(f,2*abs(Y(1:NFFT/2+1))) title('Single-Sided Amplitude Spectrum of y(t)') xlabel('Frequency (Hz)') ylabel('|Y(f)|')
快速傅立叶变换（FFT）的C++实现收藏 标准的离散傅立叶DFT 变换形式如： y k=Σj=0n-1a jωn-kj = A (ωn-k). （ωn k为复数1 的第k 个n 次方根，且定义多项式A (x)=Σj=0n-1a j x j） 而离散傅立叶逆变换IDFT (Inverse DFT)形式如：a j=(Σk=0n-1y kωn kj)/n . yk=Σj=0n-1 ajωn-kj = A (ωn-k). （ωnk 为复数1 的第k 个n 次方根，且定义多项式 A (x) = Σj=0n-1 ajxj ） 而离散傅立叶逆变换IDFT (Inverse DFT)形式如：aj=(Σk=0n-1 ykωnkj)/n . 以下不同颜色内容为引用并加以修正： 快速傅立叶变换（Fast Fourier Transform，FFT）是离散傅立叶变换（Discrete Fourier transform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 设Xn 为N 项的复数序列，由DFT 变换，任一Xi 的计算都需要N 次复数乘法和N -1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的Xi ，即N 点DFT 变换大约就需要N2 次运算。当N =1024 点甚至更多的时候，需要N2 = 1048576 次运算，在FFT 中，利用ωn 的周期性和对称性，把一个N 项序列（设N 为偶数），分为两个N / 2 项的子序列，每个N / 2点DFT 变换需要(N / 2)2 次运算，再用N 次运算把两个N / 2点的DFT 变换组合成一个N 点的DFT 变换。这样变换以后，总的运算次数就变成N + 2 * (N / 2)2 = N + N2 / 2。继续上面的例子，N =1024 时，总的运算次数就变成了525312 次，节省了大约50% 的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要N * log2N 次的运算，N = 1024 点时，运算量仅有10240 次，是先前的直接算法的1% ，点数越多，运算量的节约就越大，这就是FFT 的优越性。 FFT 的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把Xi 放到合适的位置，设i 和j 互为s = log2N 位二进制的回文数，假设s = 3, i = (110)2 = 6, j = (011)2 = 3, 如果i ≠j , 那么Xi 和Xj 应该互换位置。（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学C++ 的时候就有这样的习题。）当“位反转置换”完成后，先将每一个Xi 看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有2 项），合并实际上是“蝶形运算”（Butterfly Operation, 参考《算法导论》吧^_^），继续合并（第二次的每个多项式有4 项），直到只剩
实验七 快速傅立叶变换（FFT ）实验 一 实验目的 1． 熟悉CCS 集成开发环境； 2． 了解FFT 的算法原理和基本性质； 3． 熟悉DSP 中cmd 文件的作用及对它的修改； 4． 学习用FFT 对连续信号和时域信号进行频谱分析的方法； 5． 利用DSPLIB 中现有的库函数； 6． 了解DSP 处理FFT 算法的特殊寻址方式； 7． 熟悉对FFT 的调试方法。 二 实验内容 本实验要求使用FFT 变换对一个时域信号进行频谱分析，同时进行IFFT 。这里用到时域信号可以是来源于信号发生器输入到CODEC 输入端，也可以是通过其他工具计算获取的数据表。本实验使用Matlab 语言实现对FFT 算法的仿真，然后将结果和DSP 分析的结果进行比较，其中原始数据也直接来自Matlab 。 三 实验原理 一个N 点序列][k x 的DFT ][m X ,以及IDFT 分别定义为： 1,,1,0,][][1 0-==∑-=N m W k x m X km N N k 1,,1,0,][1 ][1 -== --=∑ N k W m X N k x km N N m 如果利用上式直接计算DFT,对于每一个固定的m,需要计算N 次复数乘法，N-1次加法,对于N 个不同的m,共需计算N 的2次方复数乘法,N*(N-1)次复数加法.显然,随着N 的增加,运算量将急剧增加, 快速傅里叶算法有效提高计算速度（本例使用基2 FFT 快速算法）,利用FFT 算法只需(N/2)logN 次运算。 四 知识要点 . 1、 CMD 文件的功能及编写 2、 一种特殊的寻址方式：间接寻址 间接寻址是按照存放在某个辅助寄存器的16位地址寻址的。C54x 的8个辅助寄存器 （AR0—AR7）都可以用来寻址64K 字数据存储空间中的任何一个存储单元。 3、 TMS320C54x DSPLIB 中关于FFT 变换的一些函数的调用（SPRA480B.pdf ） 利用DSPLIB 库时，在主程序中要包含头文件：54xdsp.lib 4、 FFT 在CCS 集成开发环境下的相关头文件 #include
MATLAB仿真实验 傅里叶变换与信号频谱图 本实验将简要介绍如何利用FFT函数描绘指定信号的频谱图像。 一、相关函数 1、FFT函数 离散傅里叶（Fourier）变换函数。 【语法】 Y = fft(X) Y = fft(X,n) Y = fft(X,[],dim) Y = fft(X,n,dim) 相关函数：IFFT(x)逆傅里叶变换。 【例1】画出函数y(t)的图像。 t = 0:0.001:0.6; x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(size(t)); plot(1000*t(1:50),y(1:50)) title('Signal Corrupted with Zero-Mean Random Noise') xlabel('time (milliseconds)')
Signal Corrupted with Zero-Mean Random Noise Frequency content of y
【图像】
Power spectral density Frequency (Hz)（a）时域图f(t)（b）频域图F(ω) 图3Sin(100πt)＋2Sin(280πt)的频谱图
4、()(100)(280)f t Sin t Cos t ππ=g 的频谱图 调制信号sin(100)t π，载波cos(280)t π。 【程序】 Frequency (Hz) Frequency (Hz) 图（c ）Cos (280πt )频谱图 图4 Sin (100πt )Cos (280πt )的频谱图
FFT算法的DSP实现 对于离散傅里叶变换(DFT)的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT非常有效。 一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列，这2 N点的复数序列与原始的2N点的实数输入序列的DFT输出一致。 使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。这样利用实数FFT 算法来计算实输入序列的DFT的速度几乎是一般FFT算法的两倍。下面用这种方法来实现一个256点实数FFT（2N=256）运算。 1.实数FFT运算序列的存储分配 如何利用有限的DSP系统资源，合理的安排好算法使用的存储器是一个比较重要的问题。本文中，程序代码安排在0x3000开始的存储器中，其中0x3000~0x3080存放中断向量表，FFT程序使用的正弦表和余弦表数据(.data段)安排在0xc00开始的地方，变量(.bss段定义) 存放在0x80开始的地址中。另外，本文中256点实数FFT程序的数据缓冲位0x2300~0x23ff , FFT变换后功率谱的计算结果存放在0x2200~0x22ff中。连续定位.cmd文件程序如下： MEMORY { PAGE 0: IPROG: origin = 0x3080,len=0x1F80 VECT: lorigin=0x3000,len=0x80 EPROG: origin=0x38000,len=0x8000 PAGE 1: USERREGS: origin=0x60,len=0x1c BIOSREGS: origin=0x7c,len=0x4 IDA TA: origin=0x80,len=0xB80 EDA TA: origin=0xC00,len=0x1400 } SECTIONS { .vectors: { } > VECT PAGE 0 .sysregs: { } > BIOSREGS PAGE 1 .trcinit: { } > IPROG PAGE 0 .gblinit: { } > IPROG PAGE 0 .bios: { } > IPROG PAGE 0 frt: { } > IPROG PAGE 0 .text: { } > IPROG PAGE 0 .cinit: { } > IPROG PAGE 0 .pinit: { } > IPROG PAGE 0 .sysinit: { } > IPROG PAGE 0 .data: { } > EDATA PAGE 1 .bss: { } > IDATA PAGE 1 .far: { } > IDATA PAGE 1 .const: { } > IDATA PAGE 1 .switch: { } > IDATA PAGE 1
MATLAB关于FFT频谱分析的程序 %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值');
title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel('频率(Hz)'); ylabel('均方根谱'); title('正弦信号y=2*pi*10t均方根谱'); grid; %求功率谱 power=sq.^2; figure(1); subplot(234); plot(f,power); xlabel('频率(Hz)'); ylabel('功率谱'); title('正弦信号y=2*pi*10t功率谱'); grid; %求对数谱 ln=log(sq); figure(1); subplot(235); plot(f,ln);
FFT算法 FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。 假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。好了，说了半天，看着公式也晕，下面圈圈以一个实际的信号来做说明。 假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为的交流信号。用数学表达式就是如下： S=2+3*cos(2*pi*50*t-pi*30/180)+*cos(2*pi*75*t+pi*90/180) 式中cos参数为弧度，所以-30度和90度要分别换算成弧度。 我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每
FFT实现自相关函数 N=38; noise=(randn(1,N)+1i*randn(1,N))/sqrt(2); f1=0.11; f2=0.15; f3=0.23; SNR1=20; SNR2=18; SNR3=17; A1=10^(SNR1/20); A2=10^(SNR2/20); A3=10^(SNR3/20); signal1=A1*exp(1i*2*pi*f1*(0:N-1)); signal2=A2*exp(1i*2*pi*f2*(0:N-1)); signal3=A3*exp(1i*2*pi*f3*(0:N-1)); un=signal1+signal2+signal3+noise; Uk=fft(un,2*N); Sk=(1/N)*abs(Uk).^2;r0=ifft(Sk); r1=[r0(N+2:2*N),r0(1:N)]; r=xcorr(un,N-1,'biased'); r11=real(r1); r12=imag(r1); r1=real(r); r2=imag(r); m=1-N:N-1; subplot(2,2,1); stem(m,r11,'o'); xlabel('m'); ylabel('实部'); title('基于FFT的自相关函数'); subplot(2,2,2); stem(m,r12,'o'); xlabel('m'); ylabel('虚部'); subplot(2,2,3); stem(m,r1); xlabel('m'); ylabel('实部');
title('基于直接计算的自相关函数'); subplot(2,2,4); stem(m,r2); xlabel('m'); ylabel('虚部');
盛年不重来，一日难再晨。及时宜自勉，岁月不待人。 DSP 课程作业 用C 语言编写FFT 程序 1，快速傅里叶变换FFT 简介 快速傅氏变换（FFT ），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 我们假设 x(n)为N 项的复数序列，由DFT 变换，任一X （m ）的计算都需要N 次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的X （m ）,即N 点DFT 变换大约就需要N^2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT 中，利用WN 的周期性和对称性，把一个N 项序列（设N=2k,k 为正整数），分为两个N/2项的子序列，每个N/2点DFT 变换需要（N/2）2次运算，再用N 次运算把两个N/2点的DFT 变换组合成一个N 点的DFT 变换。这样变换以后，总的运算次数就变成N+2（N/2）2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要Nlog2N 次的运算，N 在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT 的优越性。 2，FFT 算法的基本原理 FFT 算法的基本思想：利用DFT 系数的特性，合并DFT 运算中的某些项，吧长序列的DFT —>短序列的DFT ，从而减少其运算量。 FFT 算法分类:时间抽选法DIT: Decimation-In-Time ；频率抽选法DIF: Decimation-In-Frequency 按时间抽选的基-2FFT 算法 1、算法原理 设序列点数 N = 2L ，L 为整数。 若不满足，则补零。N 为2的整数幂的FFT 算法称基-2FFT 算法。将序列x(n)按n 的奇偶分成两组： 则x(n)的DFT: ()() ()() 12221x r x r x r x r =+=0,1,...,12N r =-()()()() 111 000 N N N nk nk nk N N N n n n X k x n W x n W x n W ---=====+∑∑∑
《测试信号分析及处理》课程作业 快速傅里叶变换 一、程序设计思路 快速傅里叶变换的目的是减少运算量，其用到的方法是分级进行运算。全部计算分解为M 级，其中N M 2log =；在输入序列()i x 中是按码位倒序排列的，输出序列()k X 是按顺序排列；每级包含2N 个蝶形单元，第i 级有i N 2 个群，每个群有12-i 个蝶形单元； 每个蝶形单元都包含乘r N W 和r N W -系数的运算，每个蝶形 单元数据的间隔为12-i ，i 为第i 级； 同一级中各个群的系数W 分布规律完全相同。 将输入序列()i x 按码位倒序排列时，用到的是倒序算法——雷德算法。 自然序排列的二进制数，其下面一个数总比上面的数大1，而倒序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位仅为而得到的。 若已知某数的倒序数是J ，求下一个倒序数，应先判断J 的最高位是否为0，与2 N k =进行比较即可得到结果。如果J k >，说明最高位为0，应把其变成1，即2 N J +，这样就得到倒序数了。如果J k ≤，即J 的最高位为1，将最高位化为0，即2N J -，再判断次高位；与4N k =进行比较，若为0，将其变位1，即4 N J +，即得到倒序数，如果次高位为1，将其化为0，再判断下一位……即从高位到低位依次判断其是否为1，为1将其变位0，若这一位为0，将其变位1，即可得到倒序数。若倒序数小于顺序数，进行换位，否则不变，防治重复交换，变回原数。 注：因为0的倒序数为0，所以可从1开始进行求解。 二、程序设计框图 （1）倒序算法——雷德算法流程图
(2)FFT算法流程
快速傅里叶变换 FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 现在圈圈就根据实际经验来说说FFT结果的具体物理意义。一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示 采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高
班级：学号：姓名 实验二FFT算法的MATLAB实现 （一）实验目的： （1）掌握用matlab进行FFT在数字信号处理中的高效率应用。 （2）学习用FFT对连续信号和时域离散信号进行谱分析。 （二）实验内容及运行结果： 题1：若x(n)=cos(nπ/6)是一个N=12的有限序列，利用MATLAB计算它的DFT 并进行IDFT变换同时将原图与IDFT变换后的图形进行对比。当求解IFFT变换中，采样点数少于12时，会产生什么问题。 程序代码： N=12; n=0:11; Xn=cos(n*pi/6); k=0:11; nk=n'*k; WN=exp(-j*2*pi/N) WNnk=WN.^nk XK=Xn*WNnk; figure(1) stem(Xn) figure(2) stem(abs(XK)) 运行结果：
IFFT变换中，当采样点数少于12时图像如下图显示：
分析：由图像可以看出，当采样点数小于12时，x(n)的频谱不变，周期为6，而XK 的频谱图发生改变。 题2：对以下序列进行谱分析 132()()103()8470x n R n n n x n n n =+≤≤?? =-≤≤??? 其他n 选择FFT 的变换区间N 为8和16点两种情况进行频谱分析，分别打印其幅频特 性曲线并进行对比、分析和讨论。 ㈠ 程序代码： x=ones(1,3);nx=0:2; x1k8=fft(x,8); F=(0:length(x1k8)-1)'*2/length(x1k8); %进行对应的频率转换 stem(f,abs(x1k8));%8点FFT title('8点FFTx_1(n)'); xlabel('w/pi'); ylabel('幅度'); N=8时：
Matlab编程实现FFT变换及频谱分析的程序代码 内容 1．用Matlab产生正弦波,矩形波,以及白噪声信号，并显示各自时域波形图 2．进行FFT变换，显示各自频谱图，其中采样率，频率、数据长度自选 3．做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱 4．用IFFT傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图 源程序 %*************************************************************** **********% % FFT实践及频谱分析% %*************************************************************** **********% %*************************************************************** **********% %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值'); title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y);
FFT原理与实现 在数字信号处理中常常需要用到离散傅立叶变换(DFT)，以获取信号的频域特征。尽管传统的DFT算法能够获取信号频域特征，但是算法计算量大，耗时长，不利于计算机实时对信号进行处理。因此至DFT被发现以来，在很长的一段时间内都不能被应用到实际的工程项目中，直到一种快速的离散傅立叶计算方法——FFT，被发现，离散是傅立叶变换才在实际的工程中得到广泛应用。需要强调的是，FFT并不是一种新的频域特征获取方式，而是DFT的一种快速实现算法。本文就FFT的原理以及具体实现过程进行详尽讲解。 DFT计算公式 其中x(n)表示输入的离散数字信号序列，WN为旋转因子，X(k)一组N点组成的频率成分的相对幅度。一般情况下，假设x(n)来自于低通采样，采样频率为fs，那么X(k)表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。因为DFT计算得到的一组离散频率幅度只实际上是在频率轴上从成周期变化的，即X(k+N)=X(k)。因此任意取N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X(k)的周期特性，于是我们又可以认为X(k)表示了从零频率开始，频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。 N点DFT的计算量 根据(1)式给出的DFT计算公式，我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法，计算N各点的X(k)共需要N^2次复数乘法和N*(N-1)次复数加法。当x(n)为实数的情况下，计算N点的DFT需要2*N^2次实数乘法，2*N*(N-1)次实数加法。 旋转因子WN的特性 1. W的对称性 N W的周期性 2. N W的可约性 3. N