基于FFT的连续信号谱分析设计
毕业设计（论文）题目：基于FFT的连续信号谱分析
学院：电气与电子信息工程学院
专业名称：电子信息工程
学号：
学生姓名：XX
指导教师：
2013年05月25日
摘要
离散傅里叶变换(DFT)的快速算法FFT的出现，使DFT在数字通信、语音信号处理、图像处理、功率谱估计、系统分析与仿真等各个领域中都得到了广泛的应用。各种应用一般都以卷积和相关运算的具体计算为依据，或者以DFT为连续傅里叶变换的近似为基础。
本文主要涉及用FFT对连续信号的频谱分析，概述了信号的频谱分析，介绍了谱分析的重要性，连续信号谱分析的过程，FFT算法的思想及性质；利用matlab 软件编制信号产生子程序，对典型信号进行谱分析并用仿真实现，绘制不同采样下的时域波形和频谱特性；根据谱分析的结果验证DFT的共轭对称性；了解可能出现的分析误差及其原因。通过matlab软件，我们演示了部分基本信号的波形和变换，使我们可以直观的了解和掌握信号与系统，数字信号处理的一些基本知识。
关键词：谱分析；DFT；FFT；matlab；连续信号
ABSTRACT
Digital signal processing course is a basic course of the telecommunications, in which the signal spectrum analysis is very common in practical applications. The emergence of the fast algorithm FFT of DFT, makes DFT have been widely used in various fields of system analysis and simulation such as in digital communications, speech signal processing, image processing, power spectrum estimation and so on. The various applications are generally based on the specific calculation of the convolution and correlation calculation, or the approximation of continuous Fourier Transform .
This paper introduce the signal spectrum analysis,and summarizes the importance of spectrum analysis ,the processes of spectral, the ideas and nature of FFT (Fast Fourier Transform); use the matlab software to develop signal generation subroutine to achieve the typical signal spectral analysis and simulation, draw time-domain waveform and spectrum characteristics under different sampling; verify the conjugate symmetry of DFT Based on the results of spectral analysis; understand errors and their causes of the possible analysis. In that experiment,we demonstrate some waveforms and transforms of the basic signals, so that we can intuitively understand and grasp the Basic knowledge of signals and systems and digital signal processing.
Key words:spectrum analysis;Discrete Fourier Transform;Fast Fourier Transform; matlab;Continuous signal
目录
1 引言 (1)
1.1 数字信号处理概述 (1)
1.2 连续信号的频谱分析 (1)
1.3 谱分析的研究意义 (2)
2 离散傅里叶变换（DFT） (3)
2.1 离散傅里叶变换的性质 (3)
2.2 利用DFT计算模拟信号的傅里叶变换 (3)
2.2.1 连续信号谱分析原理 (3)
2.2.2 对连续非周期信号的傅里叶变换的DFT逼近 (4)
2.2.3 对连续时间周期信号的傅里叶级数的DFS逼近 (6)
2.2.4 利用DFT对非周期连续时间信号傅里叶变换逼近的误差分析 (7)
3 快速傅里叶变换(FFT) (9)
3.1 FFT的来源 (9)
3.2 按时间抽选（DIT）的基-2FFT算法（库利-图基算法） (10)
3.2.1 算法原理 (10)
3.2.2 运算量 (15)
3.2.3 按时间抽选FFT算法的特点 (16)
4 数字信号处理MATLAB实现的基本知识 (19)
4.1 MATLAB简介 (19)
4.2 利用Matlab计算FFT的子函数 (19)
4.3 利用MATLAB实现信号仿真 (20)
5 总结与展望 (29)
5.1 总结 (29)
5.2 展望 (29)
致谢 (30)
参考文献 (31)
1 引言
1.1 数字信号处理概述
数字信号处理是从20世纪60年代以来，随着信息学科和计算机学科的高速发展而迅速发展起来的一门新兴学科。它的重要性日益在各个领域的应用中表现出来。数字信号处理是把信号用数字或符号表示的序列，通过计算机等处理设备，用数字的数值计算方法处理已达到提取有用信息便于应用的目的。
信号可以从不同角度进行分类，如周期信号和非周期信号，确定信号和随机信号，能量信号和功率信号，连续信号和离散信号等。对于本课题，我们需要了解的是连续时间信号，即时间是连续的，幅值可以使连续的也可以是离散的，特别地，我们一般所说的连续信号指的是模拟信号，即时间，幅值均连续，它是上一种信号的特例。
于是，我们讨论将模拟信号转化成数字信号的过程。首先把模拟信号变换为数字信号，然后在进行数字化处理，最后在恢复出模拟信号。其方框图如图1-1所示：
图1-1 数字信号处理简单方框图
图1-1表示的是模拟信号数字处理系统的方框图，实际的系统并不一定要包括它的所有框图。例如，有些系统只需数字输出，因而就不需要D/A变换器。对于数字信号处理器，可以是数字计算机或微处理机，通过软件编程对输入信号进行预期处理，这是一种软件实现方法，本课题，我们主要运用matlab软件进行处理。
自从1965年库利（Cooley）和图基（Tukey）在《计算数学》上发表了“用机器计算复序列傅里叶级数的一种算法”即“快速傅里叶变换算法”以来，数字信号处理这样一门学科蓬勃发展，逐渐形成了一整套较为完整的学科领域和理论体系，其中就包括谱分析与快速傅里叶变换（FFT），快速卷积与相关算法。
1.2 连续信号的频谱分析
信号的谱分析，就是计算信号的傅里叶变换，即将信号源发出的信号强度按频率顺序展开，使其成为频率的函数，并考察其变化规律。傅里叶变换是声学、语音、电信等领域中的一种重要方法。我们知道，为增大了数字信号处理的机动性，使数字信号处理可以在频域采用数值计算的方式进行，对有限长序列采用离散傅里叶变换(DFT)实现频域的离散化。但直接计算DFT，其计算量与区间长度的平方成正比，当长度较大时，计算量很大。而快速傅里叶变换的出现使得DFT的运算量下降明显，从而极大提高了DFT的计算效率。为数字信号处理技术应用于各种信号实时处理创造了良好的条件，大大推动了数字信号处理技术的发展。
1.3 谱分析的研究意义
现代社会，通信与传感，计算及仿真技术紧密结合，信息成为社会焦点。随着我国科学技术的发展和国内外竞争加剧，对通信水平的要求也日益突出，如果通信水平跟不上，社会成员的合作程度就会受到限制，生产力也会受到影响。而频谱分析是处理通信信息的关键一环。生活中对信号进行频谱分析具有十分重要的意义。通过对信号进行频谱分析，我们可以得到信号的频谱结构，了解信号的频率成分或系统的特征。在此基础上，可实现对信号的跟踪控制，从而实现对系统状态的早期预测，发现潜在的危险并诊断可能发生故障的原因，对系统参数进行识别及校正。因此，频谱分析是揭示信号特征的重要方法，也是处理信号的重要手段。而且这些方法和手段已经广泛应用于通信，雷达，地震，声纳，生物医学，物理，化学，音乐等领域。
信号的频谱分析不仅在现实生活中具有重要意义，同样在教学过程中也是不可或缺的。由于频谱分析仪价格昂贵，高等院校只有少数实验室配有频谱仪。但电子信息类教学，如果没有频谱仪辅助观察，学生只能从书中抽象理解信号特性，可能对教学产生影响。
2 离散傅里叶变换（DFT ）
我们知道，有限长序列在数字信号处理中很重要，而反映它的“有限长”特点的一种工具是离散傅里叶变换。理论上离散傅里叶变换除作为有限长序列的一种傅里叶表示及其重要之外，且在实际应用中有着各种高效快速计算DFT ，如快速傅里叶变换（FFT ），因而离散傅里叶变换在数字信号处理领域中起着核心作用。
2.1 离散傅里叶变换的性质
本节讨论DFT 的一些性质，它们本质上是和周期序列的DFS 概念有关，而且是由有限长序列及其DFT 表示式隐含的周期性得出的。
1、线性性质
设()X k 和()y n 是长度均为N 的两个有限长序列，它们的离散傅里叶变换分别为()X k 和()Y k ，则线性组合的离散傅里叶变换为
2、共轭对称性
设*()x n 为()x n 的复共轭序列，则有
**[()]()DFT x n X N k =-
()()jt X j x t e dw ∞
-Ω-∞
Ω=?
1
()2jt x t e d π
∞
Ω-∞
=
Ω?
3、()x n 与()X k 的奇，偶，虚，实关系如表2-1所示。
2.2 利用DFT 计算模拟信号的傅里叶变换 2.2.1 连续信号谱分析原理
我们知道，对连续信号进行频谱分析时，要对()x t 进行时域采样，得到()x n ,在对()x n 进行DFT ，得到()X k ，其中，()X k 是()x n 的傅里叶变换()jw X e 在频率区
上的等间隔采样。又因若信号为有限长，则其频谱为无线宽。因此，为防间[0,2]
止时域采样后产生频谱失真，对频谱范围较宽的信号，常采用预滤波法滤除幅度较小的高频成分，使连续信号的带宽小于折叠频率。而对于无线长信号，通常截取有限个点进行傅里叶变换。从而可看出，对信号谱分析只是近似的。下面我们做具体介绍。
表2-1 序列对应关系
2.2.2 对连续非周期信号的傅里叶变换的DFT 逼近
连续时间非周期信号()x t 的傅里叶变换对为：
()()j t X j x t e dw ∞
-Ω-∞Ω=? (2-1)
1()2j t x t e d π
∞
Ω-∞
=
Ω?
(2-2)
计算方法如下：
(1)将x(t)在x 轴上等间隔分段，每一段用一个矩形脉冲代替，脉冲的幅度为其起点的抽样值()()()x t x nT x n ==,然后把所有矩形脉冲的面积相加，而
t nT → dt T →
n dt T ∞
∞
-∞
=-∞
→
∑?
从而可得频谱密度()()j t
X j x t dw e
∞
-Ω-∞Ω=?的近似值为：
()()j nT
n X j x nT T
e
∞
-Ω=-∞
Ω≈
?∑ (2-3)
(2)将序列()()x n x nT =截断成从0t =开始长度为0T 的有限长序列，包含有N 个抽样，则上式(2-3)变为：
1
0()()N jnT n X j T x nT e --Ω=Ω≈?∑
由于时域抽样的抽样频率为1/s f T =，则频域产生以s f 为周期的周期延拓。如果频域是限带信号，则有可能不产生混叠，成为连续周期频谱序列，频域周期为1/s f T =（即时域的抽样频率）。
(3)为了数值计算，在频域上也要离散化（抽样），即在频域的一个周期(s f )中也分成N 段，即取N 个样点0s f NF =，每个样点间的间隔为0F 。频域抽样，则频域的积分公式(2-2)就变成求和式，而时域就得到原已截断的离散时间序列的周期延拓序列，其时域周期为，这时0k Ω=Ω，则
000(1)d k k Ω=+Ω-Ω=Ω
100
N k d -∞
-∞
=Ω→Ω∑?
参量关系为
001//s T F N f NT
===
这样，进过上面三个步骤，时域，频域都是离散序列。推导如下：
经过前两步：时域抽样，截断，即
()()j nT
n X j x nT T
e
∞
-Ω=-∞
Ω≈
?∑
01()()2s j nT
x nT X j d e πΩΩ≈
ΩΩ?
再由第三步：频域抽样，则得
01
00()()[()]N jknT n X jk T x nT e T DFT x n --Ω=Ω≈?=?∑
1
()()2N jknT k x nT X jk e
π
-Ω=Ω
≈Ω∑
21
000
()N jkn
N
k F X jk e
π-==Ω∑
21
1
()N jkn N
s k f X jk e
N
π-==?
Ω∑
0[()]s f IDFT x jk =?Ω 即
0()[()]X jk T DFT x n Ω≈?
01
()[()]x n IDFT X jk T
=
Ω 此为从离散傅里叶变换法求连续非周期信号的傅里叶变换的抽样值的方法。其模型如图2-1。
N N *()()x()()()()*()()()*D()()()X a N
jw jw jw
a
N
t x n n d n n n j X X k k x x x
e e e X
X
????→???→???→????→←????↓
????→Ω→???→→←????取一个周期抽样截断周期延拓周期延拓
取一个周期
卷积周期延拓
图2-1 连续非周期信号的傅里叶变换模型
2.2.3 对连续时间周期信号的傅里叶级数的DFS 逼近
连续时间周期信号x(t)的傅里叶级数为
00001()()jk t
T X jk X t dt T e -Ω=
Ω? (2-4)
00
()()jk
t
k x t X jk e
∞
=-∞
Ω=
∑Ω (2-5)
这里T 0为连续时间周期信号的周期。
由于满足：时域周期→频域离散 时域连续→频域非周期
所以，需将连续周期信号与DFS 联系起来，就需要进行如下操作：
先对时域抽样()()()x n x nT x t ==，当t nT =时，则
(1)dt n T nT T =+-=
1
N T n dt T -=→∑?
设一个周期内的样点数为N ，则式(2-4)变为
211000
01()()()N N jkn jknT N n n T X jk x nT e x n e T N π
----Ω==Ω≈=∑∑
在0T （一个周期）时间段内抽样间隔为T ，共有0
T N T
=
个抽样点。 再将频域离散序列加以截断，使它成为有限长序列，如果这个截断长度正好等于一个周期（即时域抽样造成的频域周期延拓的一个周期），则式(2-5)变成 01
00()()N jknT n x nT X jk e -Ω=≈Ω∑
1
200
()N jk
n N
n X jk e
π
-==∑Ω (2-6)
按照DFT 的定义，可得
01
()[()]X jk DFS x n N
Ω=
? 这就是用DFS 来逼近连续时间周期信号傅里叶级数对的公式。
2.2.4 利用DFT 对非周期连续时间信号傅里叶变换逼近的误差分析
我们已经讨论了用数学表达式分析了逼近过程，需关注的是最后得到的()n X n 是否是()a X t 的抽样？以及()n X n 的DFT 的系数()N X k 是否是()a X t 的频谱()a X j Ω的准确抽样？如果不是的话，误差产生在哪些运算过程中，如何减少这些误差？ 前面已经谈到，从理论上说，信号的时域长度和频域带宽不可能同时是有限长的，也就是说，时域是有限长的，则频域带宽是无限的，反之亦对。而DFT 是有限长离散时间序列与有限长离散频域序列之间的 傅里叶变换对，并隐含周期性，因而用DFT 来逼近连续时间傅里叶变换对肯定是有误差的。 从模型中可看出，最主要的是三个处理：
(1)时域抽样成x(n)则频域就会产生周期性延拓。若频域是限带宽的，只要抽样频率s f 大于信号最高频率h f 的两倍，就不会产生频域的频域的频域响应混叠现象。由于频域限带，故时域()a X t 一定是无线长的。
(2)时域截短成有限长序列()()x n d n ，此序列只在01n N ≤≤-的范围内和原()x n 相同，在此范围外为0（当用矩形序列()N R n 截短时）
。时域截断就是相乘，则在频域表现为周期性卷积，即()*()jw jw X e D e ，卷积的结果造成频谱的“扩散”，也就是频谱的“泄露”，使得形成的信号频谱与原来()jw X e 的频谱并不相同，也就是说，时域截断造成频域泄露，使频域产生误差。
(3)由于频域仍是连续的，要对频域抽样，即频域也要离散化。频域抽样造成时域()()x n d n 的周期延拓，只要频域抽样间隔0F 满足0s
f F N
≤
（即频域一个周期抽样点数M ≥N ），则这个周期延拓是不会产生时域混叠失真的。
因而，（1）中如果频域不是有限带宽的，则时域抽样会造成频域混叠失真，此时则需尽量选取时域抽样频域足够大，以减小频域混叠失真。（2）中时域截断造成频域泄露，泄露也会造成频域的混叠，因而要选择截断后的N 足够大，使得频域卷积结果更接近()jw X e （即当N →∞时，等于不产生截断）。此外，如果不采用突变形式的矩形窗截断，而是采用其他缓变窗（如哈明窗等），则可减小泄露，但这类
截断会使()()x n d n 在01n N ≤≤-范围内不等于原来的x(n)，使时域产生误差。（3）中如果频域抽样的抽样间隔0s
f F N
≥
，则造成时域()()x n d n 的周期延拓的混叠失真，使得在01n N ≤≤-范围内，时域序列也不等于x(n)。
总之，时域采用矩形窗d(n)截断，则时域只是造成截断区域01n N ≤≤-外的误差，在区域内()()()x n d n x n =就是()a X t 的抽样。频域若是限带信号，则只有时域截断时形成频域的泄露而造成频域的误差；如果频域不是限带信号，则时域抽样和时域截断都会造成频域的误差。
3 快速傅里叶变换(FFT)
快速傅里叶变换(FFT)只是离散傅里叶变换（DFT ）的一种快速计算方法，其并不是一种新的变换。如前面所讲，有限长序列的特点是其频域可离散化成有限长序列。DFT 的计算在信号处理中非常有用，再有，信号的谱分析对通信，图像传输，声纳等都是很重要的。此外，在系统的分析，设计和实现中都会用到DFT 计算。但是，在很长一段时间里，由于DFT 的计算量极大，即使采用计算机也很难对实际问题进行处理，所以没有得到真正的运用。直到库利和图基提出机器计算傅里叶级数的一种算法，并经过后来人们的改进发展，完善了一套高速有效的运算方法，使得DFT 的计算量得到简化，运算时间一般可以缩短一，二个数量级，从而使DFT 的运算在实践中得到真正广泛应用。
3.1 FFT 的来源
若x(n)是有限长序列，点数为N ，其DFT 为 1
0()()N nk N
k X k x n W
-==∑，其中，2j
N
N W e
π-=，k=0,1,…N-1 （3-1）
反变换（IDFT ）为
1
1
()()N nk N
k x n X k W
N
--==
∑，其中，n=0,1,…,N-1 （3-2）
二者的差别只在于N W 的指数符号不同，以及差一个常数因子。因而，直接计算DFT ，乘法次数和加法次数都是和2N 成正比的，当N 很大时运算量显而易见，例如，当N ＝8，DFT 需64次复乘，而当N ＝1024时，DFT 需一百多万次复运算，这对实用性很强的信号处理来说，对计算速度要求是太高了。
下面讨论减少运算量的方法。观察DFT 的运算形式可看出，利用系数nk N W 的以下特性，就可以减少DFT 的运算量：
1、nk N W 的共轭对称性
*()nk nk
N N W W -=
2、nk N W 的周期性
()()
nk n N k n k N N N N W W W ++==
3、nk N W 的可约性
nk nkm N
Nm
W
W
=，nk
nk m
N
N
m
W
W =
由此可得出
()()N k n N n k nk N N N W W W ---==，21N
N W =-，()
2
W N k k
N
N
W +=-
这样，通过这些性质，使DFT 运算中有些项合并；利用nk N W 的对称性，周期性及可约性,将较长序列的DFT 分解成较短序列的DFT ，从而使DFT 计算过程转化成许多迭代运算过程。而前面已经说到，DFT 的运算量是与2N 成正比的，所以N 越小，运算量越低。
快速傅里叶变换算法正是基于此而发展起来的。它的运算基本上可以分成两大类，即按时间抽选法和按频率抽选法。根据论文要求，我们只考虑前一种方法。
3.2 按时间抽选（DIT ）的基-2FFT 算法（库利-图基算法） 3.2.1 算法原理
取序列点数2L N =，L 为正整数。如果不满足这个条件，可以加上若干零值点，使达到这一要求，即为基-2FFT 。
对N 点序列，将()x n 按n 的奇偶性分组：
1(2)()x r x r =，2(21)()x r x r += ，其中r=0,1,…,
21
N -
则DFT 转化成
1
0()[()]()N nk
N
n X k DFT x n x n W -===∑
11
2
2(21)0
(2)(21)N N rk r k
N N r r x r W x r W --+===++∑∑ 12
222120
()()()()
N
N rk rk
N
N r r x r W
x r W -===+∑∑
（3-3）
利用系数nk N
W 的可约性，即422
j
N
N
N W e
W π-==，则
1122
12122
2
()()()()()N
N
rk k rk k
N N N N r r X k x r W W x r W X k W X k --===+=+∑∑ （3-4）
式中2()X k 与2()X k 分别是1()X r 及2()X r 的N/2点DFT ：
112
2
1112
2
()()(2)N
N rk rk
N N
r r X k x r W x r W --====∑∑ （3-5） 1122
2222
2
()()(21)N
N rk rk
N N
r r X k x r W x r W --====+∑∑ （3-6） 由式（3-4）可以看出，一个N 点DFT 分解成两个N/2点的DFT ，它们按（3-4）式又可组合成一个N 点DFT 。但是，1()x r ，2()x r 以及1()X k ，2()X k 都是N/2点的序列，即r,k 满足r,k=0,1,…,21N -。而X(k)却有N 点，而用式（3-4）只得X(k)前一半，要用1()X k ，2()X k 来表式全部X(k)值，还需利用系数的周期性质，即
()
22
2
N r k rk
N N
W W +=
这样可得到
112
2
()2
11112200
()()()()2N
N N r k rk N N r r N X k x r W x r W X k --+==+===∑∑ （3-7）
同理可得
22()()2
N
X k X k +
= （3-8）
式（3-7）、式（3-8）表明后半部段k 值所对应的1()X k ，2()X k 与前半部段k 值所对应的1()X k ，2()X k 分别相等。
在考虑nk N W 的性质：
()
2
2N N
k k k
N
N N N
W W W W +==- （3-9）
这样，将式（3-7）、式（3-8）、式（3-9）带入式（3-4）中，就可将X(k)表示成前后两段：
前半部分X(k)，(k=0,1,…,21N -)
12()()()k
N X k k k W X X =
+
（3-10）
后半部分X(k)，(k=0,1,…,21N -)
212()()()222
N
k N N N N
X k X k W X k ++=+++
12()()k N X k W X k =- ，k=0,1,…,21N - （3-11）
这样，只要求出0到(21N -)区间的所有1()X k ，2()X k 值，即可求出0到(1N -）区间内的所有()X k 值，这就大大节省了运算。
式（3-10）、（3-11）的运算可以用如下的蝶形信号流图3-1表示。
图3-1蝶形运算流
图
采用这种表示法，其过程可图3-2说明。图3-2表示N=8的情况，输出值X(0)到X(3)是由式（3-10）给出的，输出值X （4）到X （7）是由式（3-11）给出。
据此，一个N 点DFT 分成两个2N 点DFT 时，若直接计算2N 点DFT ，则每个
2N 点DFT 只需要24N 次复数乘法，
(1)22
N N
-次复数加法。同时，把两个2N 点DFT 合成为一个DFT 时，有2N 个蝶形运算，还需要2N 次复数乘法及N 次复数加法。因而通过这第一步分解后，工作量总体得到大幅降低。
既然如此，由于2L N =，因而2N 仍是偶数，可以更进一步把每个2N 点子序列按其奇偶性再分成两个4N 点子序列。先将1()x r 进行分解：
1314(2)()
(21)()
x l x l x l x l =??
+=? （3-12）
其中，l=0,1,…,41N -
114
4
2(21)
1112
2
()(2)(21)
N N lk k l N N l r X k x l W W x l --+===++∑∑
1144
34
4
24
()()N
N lk k
k
N N
N
l l x l W W x l W
--===+∑∑
342
()()k
N X k W X k =+，k=0,1,…,41N -
且 1342()()()4
k
N N X k X k W X k +
=-，k=0,1,…,41N -
其中 14334
()()N lk
N l X k x l W -==∑ （3-13）
1
4444
()()N
lk
N l X k x l W == （3-14）
图3-2 将N 点DFT 分为
两个2N 点的DFT 图
图3-3给出N=8时，将一个2N 点的DFT 分成两个4N 点DFT ，由这两个4N 点DFT 组合成一个2N 点DFT 流图。
同时，2()x r 也可以进行同样的分解，得到
2562
()()()
k
N X k X k W X k =+ 2362()()()
4
k N N X k X k W X k +
=-
其中 11445254
4
()(2)()N N
lk lk
N N
l l X k x l W x l W --====∑∑ （3-15） 11446264
4
()(21)()N
N lk lk
N N
l l X k x l W x l W --===+=∑∑ （3-16） 最后将系数统一为22
k k N N
W W =，则一个N=8点DFT 就可分成四个2点DFT ，这
样就可得到3-4图。
图3-3 由两个4N 点DFT 组合成一个2N 点DFT