毕 业 设 计（论文）
题 目：基于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 )
X (j Ω) =⎰x (t ) e -jt Ωdw -∞∞
x (t ) =1
2π⎰∞
-∞e jt Ω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 ) 的傅里叶变换X (e jw ) 在频率区间[0,2π]上的等间隔采样。又因若信号为有限长，则其频谱为无线宽。因此，为防止时域采样后产生频谱失真，对频谱范围较宽的信号，常采用预滤波法滤除幅度较小的高频成分，使连续信号的带宽小于折叠频率。而对于无线长信号，通常截取有
限个点进行傅里叶变换。从而可看出，对信号谱分析只是近似的。下面我们做具体介绍。
表2-1 序列对应关系
2.2.2 对连续非周期信号的傅里叶变换的DFT 逼近
连续时间非周期信号x (t ) 的傅里叶变换对为：
X (j Ω) =⎰x (t ) e -j Ωt dw (2-1) -∞∞
1x (t ) =2π⎰∞
-∞e j Ωt d Ω (2-2)
计算方法如下：
(1)将x(t)在x 轴上等间隔分段，每一段用一个矩形脉冲代替，脉冲的幅度为其起点的抽样值x (t ) =x (nT ) =x (n ) , 然后把所有矩形脉冲的面积相加，而
t →nT
dt →T
⎰
∞-∞dt →n =-∞∑T
∞
∞从而可得频谱密度X (j Ω) =⎰x (t ) e -∞-j Ωt dw 的近似值为：
x (nT ) e -j ΩnT X (j Ω) ≈
n =-∞∑∞⨯T (2-3)
(2)将序列x (n ) =x (nT ) 截断成从t =0开始长度为T 0的有限长序列，包含有N 个抽样，则上式(2-3)变为：
X (j Ω) ≈T ⨯∑x (nT ) e -jnT Ω
n =0N -1
由于时域抽样的抽样频率为f s =1/T ，则频域产生以f s 为周期的周期延拓。如果频域是限带信号，则有可能不产生混叠，成为连续周期频谱序列，频域周期为f s =1/T （即时域的抽样频率）。
(3)为了数值计算，在频域上也要离散化（抽样），即在频域的一个周期(f s ) 中也分成N 段，即取N 个样点f s =NF 0，每个样点间的间隔为F 0。频域抽样，则频
域的积分公式(2-2)就变成求和式，而时域就得到原已截断的离散时间序列的周期延拓序列，其时域周期为，这时Ω=k Ω0，则
d Ω=(k +1) Ω0-k Ω0=Ω0
⎰
参量关系为 ∞-∞d Ω→∑Ω0 k =0N -1
T 0=1/F 0=N /f s =NT
这样，进过上面三个步骤，时域，频域都是离散序列。推导如下：
经过前两步：时域抽样，截断，即
X (j Ω) ≈
n =-∞∑∞x (nT ) e -j ΩnT ⨯T
x (nT ) ≈1Ωs j ΩnT X (j Ω) d Ωe ⎰02π
再由第三步：频域抽样，则得
X (jk Ω0) ≈T ⨯∑x (nT ) e -jknT Ω0=T ⨯DFT [x (n )]
n =0N -1
Ω0N -1x (nT ) ≈X (jk Ω0) e jknT Ω0 ∑2πk =0
=F 0∑X (jk Ω0) e
k =0N -1jkn 2πN 2πjkn 1N -1
=f s ⨯∑X (jk Ω0) e N N k =0
=f s ⨯IDFT [x (jk Ω0)]
即
X (jk Ω0) ≈T ⨯DFT [x (n )]
x (n ) =1IDFT [X (jk Ω0)] T
此为从离散傅里叶变换法求连续非周期信号的傅里叶变换的抽样值的方法。其模型如图2-1。
取一个周期抽样截断周期延拓*(n ) ←−−−−−−−→(n ) (t ) −−−→x (n ) −−−→x(n ) d (n ) −−−−→−x a x N x N 周期延拓
↓
−→X (e )*D(e ) →X a (j Ω) →X (e ) −−卷积jw jw jw X 取一个周期*(k ) ←−−−−−−−→−X N (k ) N 周期延拓
图2-1 连续非周期信号的傅里叶变换模型
2.2.3 对连续时间周期信号的傅里叶级数的DFS 逼近
连续时间周期信号x(t)的傅里叶级数为
X (jk Ω0) =1T 0-jk Ω0t X (t ) dt e ⎰0T (2-4)
x (t ) =
k =-∞∑∞jk t X (jk Ω0) e Ω0
(2-5)
这里T 0为连续时间周期信号的周期。
由于满足：时域周期→频域离散
时域连续→频域非周期
所以，需将连续周期信号与DFS 联系起来，就需要进行如下操作：
先对时域抽样x (n ) =x (nT ) =x (t ) ，当t =nT 时，则
dt =(n +1) T -nT =T
⎰
T 0
dt →∑T
n =0
N -1
设一个周期内的样点数为N ，则式(2-4)变为
2π-jkn T N -11N -1-jknT Ω0
X (jk Ω0) ≈∑x (nT ) e =∑x (n ) e N
T 0n =0N n =0
在T 0（一个周期）时间段内抽样间隔为T ，共有N =
T 0
个抽样点。 T
再将频域离散序列加以截断，使它成为有限长序列，如果这个截断长度正好等于一个周期（即时域抽样造成的频域周期延拓的一个周期），则式(2-5)变成 x (nT ) ≈∑X (jk Ω0) e jknT Ω0
n =0N -1
=∑X (jk Ω0) e
n =0
N -1
jk
2π
n N
(2-6)
按照DFT 的定义，可得 X (jk Ω0) =
1
⨯DFS [x (n )] N
这就是用DFS 来逼近连续时间周期信号傅里叶级数对的公式。
2.2.4 利用DFT 对非周期连续时间信号傅里叶变换逼近的误差分析
我们已经讨论了用数学表达式分析了逼近过程，需关注的是最后得到的X n (n ) 是否是X a (t ) 的抽样？以及X n (n ) 的DFT 的系数X N (k ) 是否是X a (t ) 的频谱X a (j Ω) 的准确抽样？如果不是的话，误差产生在哪些运算过程中，如何减少这些误差？ 前面已经谈到，从理论上说，信号的时域长度和频域带宽不可能同时是有限长的，也就是说，时域是有限长的，则频域带宽是无限的，反之亦对。而DFT 是有限长离散时间序列与有限长离散频域序列之间的 傅里叶变换对，并隐含周期性，因而用DFT 来逼近连续时间傅里叶变换对肯定是有误差的。 从模型中可看出，最主要的是三个处理：
(1)时域抽样成x(n)则频域就会产生周期性延拓。若频域是限带宽的，只要抽样频率f s 大于信号最高频率f h 的两倍，就不会产生频域的频域的频域响应混叠现
象。由于频域限带，故时域X a (t ) 一定是无线长的。
(2)时域截短成有限长序列x (n ) d (n ) ，此序列只在0≤n ≤N -1的范围内和原。时域截断就是相乘，则x (n ) 相同，在此范围外为0（当用矩形序列R N (n ) 截短时）
在频域表现为周期性卷积，即X (e jw )*D (e jw ) ，卷积的结果造成频谱的“扩散”，也就是频谱的“泄露”，使得形成的信号频谱与原来X (e jw ) 的频谱并不相同，也就是说，时域截断造成频域泄露，使频域产生误差。
(3)由于频域仍是连续的，要对频域抽样，即频域也要离散化。频域抽样造成时域x (n ) d (n ) 的周期延拓，只要频域抽样间隔F 0满足F 0≤
f s
（即频域一个周期抽N
样点数M ≥N ），则这个周期延拓是不会产生时域混叠失真的。
因而，（1）中如果频域不是有限带宽的，则时域抽样会造成频域混叠失真，此时则需尽量选取时域抽样频域足够大，以减小频域混叠失真。（2）中时域截断造成频域泄露，泄露也会造成频域的混叠，因而要选择截断后的N 足够大，使得频域卷积结果更接近X (e jw ) （即当N →∞时，等于不产生截断）。此外，如果不采用突变形式的矩形窗截断，而是采用其他缓变窗（如哈明窗等），则可减小泄露，但这类截断会使x (n ) d (n ) 在0≤n ≤N -1范围内不等于原来的x(n)，使时域产生误差。（3）中如果频域抽样的抽样间隔F 0≥
f s
，则造成时域x (n ) d (n ) 的周期延拓的混叠失真，N
使得在0≤n ≤N -1范围内，时域序列也不等于x(n)。
总之，时域采用矩形窗d(n)截断，则时域只是造成截断区域0≤n ≤N -1外的误差，在区域内x (n ) d (n ) =x (n ) 就是X a (t ) 的抽样。频域若是限带信号，则只有时域截断时形成频域的泄露而造成频域的误差；如果频域不是限带信号，则时域抽样和时域截断都会造成频域的误差。
3 快速傅里叶变换(FFT)
快速傅里叶变换(FFT)只是离散傅里叶变换（DFT ）的一种快速计算方法，其并不是一种新的变换。如前面所讲，有限长序列的特点是其频域可离散化成有限长序列。DFT 的计算在信号处理中非常有用，再有，信号的谱分析对通信，图像传输，声纳等都是很重要的。此外，在系统的分析，设计和实现中都会用到DFT 计算。但是，在很长一段时间里，由于DFT 的计算量极大，即使采用计算机也很难对实际问题进行处理，所以没有得到真正的运用。直到库利和图基提出机器计算傅里叶级数的一种算法，并经过后来人们的改进发展，完善了一套高速有效的运算方法，使得DFT 的计算量得到简化，运算时间一般可以缩短一，二个数量级，从而使DFT 的运算在实践中得到真正广泛应用。
3.1 FFT的来源
若x(n)是有限长序列，点数为N ，其DFT 为
X (k ) =∑x (n ) W
k =0N -1
nk N
，其中，W N =e
-j
2πN
，k=0,1,„N-1 （3-1）
反变换（IDFT ）为
1N -1-nk
x (n ) =∑X (k ) W N ，其中，n=0,1,„,N-1 （3-2）
N k =0
二者的差别只在于W N 的指数符号不同，以及差一个常数因子。因而，直接计算DFT ，乘法次数和加法次数都是和N 2成正比的，当N 很大时运算量显而易见，例如，当N ＝8，DFT 需64次复乘，而当N ＝1024时，DFT 需一百多万次复运算，这对实用性很强的信号处理来说，对计算速度要求是太高了。
nk 下面讨论减少运算量的方法。观察DFT 的运算形式可看出，利用系数W N 的以
下特性，就可以减少DFT 的运算量：
nk
1、W N 的共轭对称性
nk *-nk
(W N ) =W N
nk
2、W N 的周期性
nk (n +N ) k n (k +N )
W N =W N =W N
nk
3、W N 的可约性
nk nkm nk
，W N W N =W Nm =W N
由此可得出
W
(N -k ) n
N
=W
(N -n ) k N
=W
-nk N ，W N =-1，W
N
(k +N )
N
k
=-W N
nk
这样，通过这些性质，使DFT 运算中有些项合并；利用W N 的对称性，周期性
及可约性, 将较长序列的DFT 分解成较短序列的DFT ，从而使DFT 计算过程转化成许多迭代运算过程。而前面已经说到，DFT 的运算量是与N 2成正比的，所以N 越小，运算量越低。
快速傅里叶变换算法正是基于此而发展起来的。它的运算基本上可以分成两大类，即按时间抽选法和按频率抽选法。根据论文要求，我们只考虑前一种方法。
3.2 按时间抽选（DIT ）的基-2FFT 算法（库利-图基算法） 3.2.1 算法原理
取序列点数N =2L ，L 为正整数。如果不满足这个条件，可以加上若干零值点，使达到这一要求，即为基-2FFT 。
对N 点序列，将x (n ) 按n 的奇偶性分组：
x (2r ) =x 1(r ) ，x (2r +1) =x 2(r ) ，其中r=0,1,„, N 2-1
则DFT 转化成
nk
X (k ) =DFT [x (n )]=∑x (n ) W N
n =0N -1
2rk (2r +1) k
=∑x (2r ) W N +∑x (2r +1) W N r =0N
2
r =0N -12r =0
N -1
N -12
=
∑x (r )(W
1r =0
2N
2rk N
2rk
) +∑x 2(r )(W N ) （3-3） 4π
N
利用系数W 的可约性，即W =e
N -12r =0
nk N
-j
=W N ，则
rk k
X (k ) =∑x 1(r ) W N +W N ∑x 2(r ) W N rk =X 1(k ) +W N k X 2(k ) （3-4）
r =0
N
-12
式中X 2(k ) 与X 2(k ) 分别是X 1(r ) 及X 2(r ) 的N/2点DFT ：
rk rk
（3-5） X 1(k ) =∑x 1(r ) W N =∑x 1(2r ) W N
r =0
r =0
N
-12
N -12
rk rk
（3-6） X 2(k ) =∑x 2(r ) W N =∑x 2(2r +1) W N
r =0
r =0
N
-12N -12
由式（3-4）可以看出，一个N 点DFT 分解成两个N/2点的DFT ，它们按（3-4）式又可组合成一个N 点DFT 。但是，x 1(r ) ，x 2(r ) 以及X 1(k ) ，X 2(k ) 都是N/2点的序列，即r,k 满足r,k=0,1,„, N 2-1。而X(k)却有N 点，而用式（3-4）只得X(k)前一半，要用X 1(k ) ，X 2(k ) 来表式全部X(k)值，还需利用系数的周期性质，即
rk W N =W N
r (k +)
这样可得到
N r (k +N ) rk X 1(k +) =∑x 1(r ) W N =∑x 1(r ) W N =X 1(k ) （3-7）
2r =0r =0
同理可得
X 2(k +
N
-12N -12
N
) =X 2(k ) （3-8） 2
式（3-7）、式（3-8）表明后半部段k 值所对应的X 1(k ) ，X 2(k ) 与前半部段k 值所对应的X 1(k ) ，X 2(k ) 分别相等。
nk 在考虑W N 的性质：
(k +N )
N
W N k k =W N W N =-W N （3-9）
这样，将式（3-7）、式（3-8）、式（3-9）带入式（3-4）中，就可将X(k)表示成前后两段：
前半部分X(k)，(k=0,1,„, N 2-1)
X (k ) =
k
X 1(k ) +X 2(k ) W N （3-10）
后半部分X(k)，(k=0,1,„, N 2-1)
N
k +N N N
X (k +) =X 1(k +) +W N 2X 2(k +)
222
k
=X 1(k ) -W N k=0,1,„, N 2-1 （3-11） X 2(k ) ，
这样，只要求出0到(N 2-1) 区间的所有X 1(k ) ，X 2(k ) 值，即可求出0到(N -1）区间内的所有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 分成两个N 2点DFT 时，若直接计算N 2点DFT ，则每个
N 2点DFT 只需要N 2次复数乘法，
N N
(-1) 次复数加法。同时，把两个N 2点22
DFT 合成为一个DFT 时，有N 2个蝶形运算，还需要N 2次复数乘法及N 次复数加法。因而通过这第一步分解后，工作量总体得到大幅降低。
既然如此，由于N =2L ，因而N 2仍是偶数，可以更进一步把每个N 2点子序列按其奇偶性再分成两个N 4点子序列。先将x 1(r ) 进行分解：
⎧x 1(2l ) =x 3(l )
（3-12） ⎨
⎩x 1(2l +1) =x 4(l )
其中，l=0,1,„, N 4-1
2lk k (2l +1)
X 1(k ) =∑x 1(2l ) W N +W N ∑x 1(2l +1)
l =0
r =0
N
-14l =0
N -14
N -14N -14
lk k
=∑x 3(l ) W +W N
l =0
∑x (l ) W
4
k
N
k
X 4(k ) ，k=0,1,„, N 4-1
=X 3(k ) +W
且 X 1(k +
N k
) =X 3(k ) -W N X 4(k ) ，k=0,1,„, N 4-1
4
N
-14l =0
lk
其中 X 3(k ) =∑x 3(l ) W N （3-13）
1
l =0
lk
X 4(k ) =x 4(l ) W （3-14）
图3-2 将N 点DFT 分为两个N 2点的DFT 图
图3-3给出N=8时，将一个N 2点的DFT 分成两个N 4点DFT ，由这两个N 4点DFT 组合成一个N 2点DFT 流图。
同时，x 2(r ) 也可以进行同样的分解，得到
k
X 2(k ) =X 5(k ) +W N X 6(k )
X 2(k +
N k
) =X 3(k ) -W N X 6(k )
4
N
-14l =0
N -14l =0
lk lk
其中 X 5(k ) =∑x 2(2l ) W N （3-15） =∑x 5(l ) W lk lk X 6(k ) =∑x 2(2l +1) W N （3-16） =∑x 6(l ) W N
l =0
l =0
N
-14N -14
k 2k
=W N 最后将系数统一为W N ，则一个N=8点DFT 就可分成四个2点DFT ，这
样就可得到3-4图。
图3-3 由两个N 4点DFT 组合成一个N 2点DFT
根据上面同样的分析知道，最后剩下的是2点DFT ，对于此例N=8，就是四个2点DFT ，其输出为X 3(k ) ，X 4(k ) ，X 5(k ) ，X 6(k ) (k=0,1)，这由式（3-13）至（3-16）可以计算出来。由此可得出一个按时间抽选运算的大致的8点DFT 流图，如图3-5所示。
图3-4 将一个N 点DFT 分成四个N 4点DFT
3.2.2运算量
由图3-5可知，当N =2时(L表示级数) ，每级都由N 2个蝶形运算组成，每
L
个蝶形运算有一次复乘，二次复加，因而每级运算都需N 2次复乘和N 次复加，因此L 级运算共需
复乘数 m f =
N N
L =log 2N 22
复加数 a f =NL =N log 2N
因为W N =1，W N
N =-j ，这几个系数都不用乘法运算。此外，当N 较大时，
这些特例相对而言就很少。直接DFT 算法运算量，复数乘法：N 2，复数加法：
N (N -1) 。直接计算DFT 与FFT 算法的计算量之比为M =
2N
。 log 2N
图3-5 N=8按时间抽选FFT 运算流图
3.2.3 按时间抽选FFT 算法的特点
(1)原位运算（同址运算）
从上图可以看出该运算是有规律可循的，其每级运算都是由N 2个蝶形运算构成，下述迭代运算由每一个蝶形结构实现：
r
⎧⎪X m (k ) =X m -1(k ) +X m -1(j ) W N
（3-17） ⎨r
⎪⎩X m (j ) =X m -1(k ) -X m -1(j ) W N
式（3-17）的蝶形运算如图3-6所示。
从图3-5中看出，某一列的任何两点k 和j 的运算后，得到结果为下一列k,j 两节点的变量，而与其他节点变量无关，通过原位运算，经过蝶形运算，其结果为另一列数据，它们以蝶形为单位仍存储在这一组存储器中，中间无需其他存储器。这样只需N 个存储结构单元，虽然进入蝶形的组合关系有差异，但下一级的运算也可采用这种运算方式。这种结构只需少量存储单元，从而降低设备成本。
图3-6 蝶形运算结构
(2)倒位序规律
按照原位计算，FFT 的输出X(k)在存储单元中按顺序排列，即X(0),X(1)，... ，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，按x(0),x(4)，... ，x(7)的顺序依次存入存储单元，看起来好像是“杂乱无序”的，实际是有规律可循的，即倒位序。
造成该原因是由输入x(n)按标号n 的偶奇的不断分组引起。这种不断分成奇偶子序列的过程可用二进制树状图3-7来描述。这就是DIT 的FFT 算法输入序列的 序数成为倒位序的原因。
(3)倒位序的实现
一般实际运算中，总是先向存储单元中按顺序输入序列，为形成倒位序的排列，这通过变址来完成。若序列序号n 用二进制数表示，顺序为(n 2, n 1, n 0) ，那么其倒位序二进制数应当是(n 0, n 1, n 2) 。表3-1列出了N=8时自然顺序以及相应的倒位序的二进制数。
(4)蝶形运算节点的“距离”
仍以N=8时FFT 运算流图为例，其输入是倒位序的，输出是自然顺序的，第一级每个蝶形两节点间“距离”是4，由此类推得，对N =2L 点FFT ，当输出为正常顺序，输入为倒位序时，其第ｍ级运算，每个蝶形两节点“距离”则为2m -1。
nk (5)W N 的确定
对第m 级，一个蝶形运算两节点的“距离”为2m -1，于是第m 级的一个蝶形计算可写成
r
X m (k ) =X m -1(k ) +X m -1(k +2m -1) W N
r
X m (k +2m -1) =X m -1(k ) -X m -1(k +2m -1) W N
图3-7 描述倒位序的树状图
表3-1 自然顺序及相应的倒位序
4 数字信号处理 MATLAB实现的基本知识
4.1 MATLAB简介
MATLAB软件最初是由美国新墨西哥大学计算机系Cleve Moler 教授用FORTRAN
语言编写的矩阵计算软件，目前该软件是Math Works公司的产品。MATLAB 全称是Matrix Laboratory，该软件以矩阵的形式来描述和处理所有的计算机数据，主要由MATLAB 语言、工作环境、图形系统、数学函数库、API 函数和工具箱等六大部分组成、前面五部分是其核心内容工作箱子是其辅助内容。第一部分是MATLAB 语言，与C 语言、FORTRAN 等高级语言近似，用于实现各种算法，MATLAB 包含的编译器还可以将MATLAB 的算法和应用程序文件转变成独立可执行的应用程序。第二部分是MATLAB 工作环境，用户可以完成程序设计、数值计算、图形绘制、输入输出、文件管理等多项功能。第三部分是图形系统，可以完成2D 和3D 数据显示、动画生成、视频处理、图形显示等功能。第四部分是基础与通用的MATLAB 数据函数库，其功已经大大超越了常用高级语言的基本数学函数库。第五部分是MATLAB API函数，这使得C 语言、FORTRAN 语言可以与MATLAB 混合应用，这种交互可以取长补短，提高程序的运行效率，丰富程序开发的手段。第六部分是各种功能极其强大的专业工具箱，涉及到数据获取、数字信号处理、数字图像处理、控制系统分析和设计、最优化方法、系统辨识、智能计算、金融财务分析以及生物信息处理等内容。这些工具箱的算法是开放可拓展的。第三方和用户可根据需要进行拓展。
目前，MATLAB 已经在自动控制原理、数字信号处理、数字图像处理、时间序列分析、数理统计、动态系统仿真等课程教学中得到广泛应用，是理工类本科学生和研究生必须掌握的基本工具。不仅如此，MATLAB 在科学研究中被公认为准确、可靠、方便的科学计算软件；在工业部门，MATLAB 被认为是高效仿真、分析和开发的软件工具。基于此，本论文以MATLAB 为基本，做频谱分析。
4.2利用Matlab 计算FFT 的子函数
用Matlab 方法计算信号的DFT 时，主要是用函数fft (x , N ) 和ifft (x , N ) 。对于这两个函数，如果N 为2的正整数幂，则可以得到本章中介绍的基2 FFT 快速算法；如果N 既不是2的正整数幂，也不是质数，则函数将N 分解成质数，得到较慢的混
合基 FFT算法；如果 N 为质数，则fft 函数采用原来的 DFT 算法。 实验涉及的MATLAB 子函数
1、fft
功能：一维快速傅里叶变换
调用格式：y =fft (x ) ；利用FFT 算法计算矢量x 的离散傅里叶变换，当x 为矩阵时，y 为矩阵x 每一列的FFT 。当x 的长度为2的幂次方时，则fft 函数采用基2的FFT 算法，否则采用稍慢的混合基算法。
y =fft (x , n ) ；采用n 点FFT 。当x 的长度小于n 时，fft 函数在x 的尾部补零，
以构成n 点数据；当x 的长度大于n ，fft 函数会截断序列x 。当x 为矩阵时，fft 函数按类似的方式处理列长度。
2、ifft
功能：一维快速傅里叶变换
调用格式：y =ifft (x ) ; 用于计算矢量x 的IFFT 。当x 为矩阵时，计算所得的y 为矩阵x 中每一列的IFFT 。
y =ifft (x , n ) ; 采用n 点IFFT 。当length (x ) n 时，
将x 截断，使length (x ) =n 。
用MATLAB 提供的子函数进行快速傅里叶变换时，从理论学习可知，DFT 是唯一在时域和频域均为离散序列的变换方法，它适用于有限长序列。尽管这种变换方法是可以数值计算的，但如果只是简单的按照定义进行数据处理，当序列长度很大时，则将占用很大的内存空间，运算时间将很长。
快速傅里叶变换是用于DFT 运算的高效运算方法的统称，FFT 只是其中的一种，FFT 主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N 的序列分解成多个短序列，如基2算法，大大缩短了运算时间。
4.3 利用MATLAB 实现信号仿真
用FFT 进行无限长序列的频谱计算，首先要将无限长序列截断成一个有限长序列。序列长度的取值对频谱有较大的影响，带来的问题是引入频谱的泄漏和波动。
例1 给出信号为x (t ) =cos(20πt ) +cos(8πt ) +cos(16πt ) ，是一个连续信号，所以我们根据其最高频率确定采样速率f s =64Hz ，以及由频率分辨率选择采样点数N=64、32、16三种情况进行分析。
解答 程序如下:
N=input('选择FFT 变换区间长度N ：16 or 32 or 64:'); %**输入FFT 变换区间长度N**%
fs=64; n=0:N-1; t=n/fs;
x=cos(20*pi*t)+cos(8*pi*t)+cos(16*pi*t); subplot(2,2,1);
plot(t,x); %**作信号时域波形**% xlabel('t'); ylabel('x(t)'); title('x(t)时域波形');
f=fft(x,N); %**计算FFT**% subplot(2,2,2); %**绘制x(n)时域图像**% M=N*0.4; n=0:N-1; stem(n,x,'.'); n=0:N-1; m=zeros(N); hold on; plot(n,m); t=max(x); xlabel('n') ylabel('x(n)') title('x(n)时域波形');
subplot(2,2,3); %**绘制N 点FFT 图**% n=0:N-1; stem(n,f,'.'); t=max(f); xlabel('K')
string2=['x(n)的N=',sprintf('%2.0f',N),'点FFT']; title(string2); ylabel('|X(k)|'); 运行如下：当输入16时，
图4-1 N=16时x(t)的DFT
当输入32时，
图4-2 N=32时x(t)的DFT
当输入64时，
图4-3 N=32时x(t)的DFT
分析：
1、T s 相同，当N 越大，在所给范围内等间隔抽样点数越多，且时域信号的长度L NT s 保留得越长，则分辨率越高，频谱特性误差越小。反之，当N 越小时，等间隔采样点数越少，则有可能漏掉某些重要的信号分量，称为栅栏效应。
2、图4-1中，在点数n=1,2,14,15外，还有其他值，而在N=32，N=64时只有几根谱线有值，其它处为0，原因是图4-1中信号发生泄漏，信号的频率成分泄漏到周围一些离散频率点上。
例2 求正，余弦信号谱。假设采样频率为64hz ，采样点N=32，正弦信号频率为10HZ ，余弦信号频率为20HZ ，并验证有关性质。
解答 程序如下： clear all; N=32;n=0:N-1; Fs=64;T=1/Fs; x1=cos(20*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N));
subplot(4,2,1);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('余弦序列'); subplot(4,2,2);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs; x1=sin(10*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,3);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('正弦序列'); subplot(4,2,4);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs;
x1=sin(10*2*pi*n*T)+cos(20*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,5);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('合序列'); subplot(4,2,6);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs;
x1=cos(20*2*pi*n*T)+j*sin(10*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,7);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('复数序列'); subplot(4,2,8);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); 结果如下：
余弦序列
n
正弦序列
32点FFT 幅频曲线
X 1(k )
k
32点FFT 幅频曲线
x 1(n )
n 合序列
X 1(k )
x 1(n )
k
32点FFT 幅频曲线
n
复数序列
X 1(k )
x 1(n )
k
32点FFT 幅频曲线
n
X 1(k )
x 1(n )
k
图4-4 各种序列频谱
分析：
1、对单一连续信号进行傅里叶变换，首先，将时域信号进行采样，显然，f s f h ，满足题意。由x(t)到x(n)，可得最小采样点数为10，又因使用基-2FFT 算法，故最小采样点数为16，这里N 取32符合要求。
2、图4-4中，正弦或余弦信号频谱只有两个点不等于0，这正说明它们自身的频率。
3、由合序列可知，两信号时域叠加，其频谱亦为线性叠加。
4、复数序列实数部分的离散傅立叶变换是原来序列离散傅立叶变换的共轭对称分量；复数序列虚数部分的离散傅立叶变换是原来序列离散傅里叶变换的共轭反对称分量。这里讨论的x(n)与X(k)均为有限长序列，所以这里的对称性是指关于
N 2点的对称性。
例3 绘制一个正弦信号频谱图并作相关分析，参数自定。 解答 程序如下： clear;
fs=100;N=1024; %采样频率和数据点数 A=20;B=30;; n=0:N-1;t=n/fs; %时间序列 x=A*sin(2*pi*B*t); %信号
y=fft(x,N); %对信号进行傅里叶变换 yy=abs(y); %求得傅里叶变换后的振幅 yy=yy*2/N; %幅值处理 f=n*fs/N; %频率序列
subplot(2,2,1),plot(f,yy); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图1 fs=100 N=1024'); grid on; %两种信号叠加
x=A*sin(2*pi*B*t)+2*A*sin(2*pi*1.5*B*t); %信号 y=fft(x,N); %对信号进行傅里叶变换 yy=abs(y); %求得傅里叶变换后的振幅 yy=yy*2/N; %幅值处理 f=n*fs/N; %频率序列
subplot(2,2,2),plot(f,yy); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图2 fs=100,N=1024 两种信号叠加'); grid on; %加噪声之后的图像
x=A*sin(2*pi*B*t)+28*randn(size(t)); y=fft(x,N); yy=abs(y); yy=yy*2/N; %幅值处理
subplot(2,2,3),plot(f(1:N/2.56),yy(1:N/2.56)); xlabel('频率/\itHz'); ylabel('振幅'); title('图3 fs=100,N=1024混入噪声');
grid on; %改变采样点数
N=128; n=0:N-1;t=n/fs; %时间序列
x=A*sin(2*pi*B*t); %信号
y=fft(x,N); %对信号进行傅里叶变换
yy=abs(y); %求得傅里叶变换后的振幅
yy=yy*2/N; %幅值处理
f=n*fs/N; %频率序列
subplot(2,2,4),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图4 fs=100 N=128');
结果如下：
图4-5 信号频谱
分析：
1、图1和图2取相同的采样频率fs=100和数据点数N=1024，不同的是图2采用两种不同幅值和频率的正弦信号叠加。从图中可以看出图2可以明显的看出含有两种不同的频率成分的信号，幅值也不相同。由此可以看出，不同频率的正弦信号叠加，在频域当中互相分离互不影响。
2、从图1和图2可以看出:整个频谱图是以fs/2频率为对称轴的。由此可以知道傅里叶变换数据的对称性。因此在用傅里叶变换做频谱分析时，我们只需做出前一半频谱图即可。
3、图3为混入噪声之后的频谱。可以看出噪声分布在整个频率轴上，并且由于噪声中含有与原信号频率相同的成分，叠加之后导致幅值增加。加大噪声的幅值之后将分辨不出原信号的频率。
4、图4减少了数据点数，N=128。与图1相比较，采用128点和1024点的相同频率的振幅是有不同的表现值。因此振幅的大小与所用采样点数有关。一定范围内采样点数越多。信号的幅值越接近真实值。
5 总结与展望
5.1 总结
通过这次毕业设计，使我充分理解了信号处理的重要性，及MATLAB 软件编程方面的知识。在学习过程中，虽然遇到了许多问题，但在老师的指导，自己翻阅大量资料，一步步找出问题所在，同时也暴露出之前在学习数字信号处理等科目的不足。实践是检验真理的唯一标准，通过动手操作，理论与实际结果相结合，才能更好的理解和掌握知识，为我所用。
对于我们通信专业，我认为信号处理是基本专业知识，这个专业以后就业的方向也很多，就业面很广。我们毕业以后工作，可以进入设备制造商、 运营商、专有服务提供商以及银行等领域工作。当然, 就业形势每年都会变化，所以关键还是要看自己。可以从事硬件方面，比如说PCB ，别小看这门技术，平时我们在试验时制作的简单，这一技术难点就在于板的层数越多，要做的越稳定就越难，这可是非常有难度的，如果学好了学精了，也是非常好找工作的。也可以从事软件方面，这实际上也要求我们具备比较好的模电，数电等基础知识。
在今后社会发展和学习过程中，我相信信号处理这方面知识一定大有用处。只有不懈努力，遇到问题不要退缩，然后一一进行解决，才能取得成功。
5.2 展望
目前，信号处理课群体系正逐步成熟，并得到国内高校的认可，其体现了理论与实践的有机结合，体现了原理、方法和技术的有机结合。在介绍数字信号处理的理论和方法的基础上，进行MATLAB 仿真实验，再进行基于DSP 系统的开发应用实验。
信息处理前景可谓前途无可限量, 随着社会的发展, 信息交流将越发频繁, 在信息充斥的世界里, 信息处理将越发被需要, 发展前景广阔。
致 谢
本文是在桂静宜老师的谆谆教诲和悉心指导下完成，桂老师在学术上有着严谨的教学作风，实事求是的治学态度，让我受益匪浅。生活上桂老师平易近人，和蔼可亲，是我学习和生活中的榜样。在毕业设计研究之中，桂老师给了我最及时和最有效的指导，这使得我最终克服各种困难，顺利地完成了论文。在此，谨向我的老师表示感谢。
同时，也感谢我的同学一直对我的鼓励和支持，在这次论文中也给了我不少建议，使我在学习上不断进步。最后，《基于FFT 的连续信号谱分析》得以基本完成。
参考文献
[1] 高西全, 丁玉美著. 数字信号处理（第三版). 西安电子科技大学出版社
[2]刘小东著. MATLAB7.X的系统分析与设计——信号处理. 西安电子科技大学出版社
[3]郭文彬, 桑林编著. 庞兴华审订. 通信原理—基于MATLAB 的计算机仿真. 第一版. 北京：北京邮电大学出版社2006年6月
[4]郑南宁. 数字信号处理（第二版）. 西安交通大学出版社. 1995
[5]邓华编著. MATLAB通信仿真及应用实例详解. 第一版. 北京：人民邮电出版社. 2003年9月
[6]邵玉斌编著. Matlab/Simulink通信系统建模与仿真实例分析. 第1版. 北京：清华大学出版社2008年6月
[7]唐向宏.MATLAB 及在电子信息类课程中的应用. 第四版. 西安：电子科技大学出版社，2006年
[8]樊昌信, 詹道庸, 徐炳祥, 吴成柯. 通信原理. 第四版. 北京:国防工业出版社. 2000年
[9]Schilling R J and Harris S L. Fundamentals of digital signal processing using MATLAB .西安：西安交通大学出版社. 2005
[10]楼顺天, 李博菡. 基于MATLAB 的系统分析与设计——信号处理. 西安电子科技大学出版社. 2000
[11]Paulo S R.Diniz ,Eduardo A.B.da Silva . Digital Signal Processing and System Analysis and Design. 电子工业出版社. 2002
[12]刘树棠译. Oppenheim A V,Willsky S and Nawab S H.信号与系统（第二版). 西安交通大学出版社. 2001
[13]胡广书. 数字信号处理-理论、算法、与实现(第二版). 北京：清华大学出版社.2003
[14]党宏社. 信号与系统实验(MATLAB版). 西安：西安电子科技大学出版社.2007
[15]陈怀琛. 数字信号处理教程—MATLAB 释义与实现. 北京：电子工业出版社.2005
毕 业 设 计（论文）
题 目：基于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 )
X (j Ω) =⎰x (t ) e -jt Ωdw -∞∞
x (t ) =1
2π⎰∞
-∞e jt Ω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 ) 的傅里叶变换X (e jw ) 在频率区间[0,2π]上的等间隔采样。又因若信号为有限长，则其频谱为无线宽。因此，为防止时域采样后产生频谱失真，对频谱范围较宽的信号，常采用预滤波法滤除幅度较小的高频成分，使连续信号的带宽小于折叠频率。而对于无线长信号，通常截取有
限个点进行傅里叶变换。从而可看出，对信号谱分析只是近似的。下面我们做具体介绍。
表2-1 序列对应关系
2.2.2 对连续非周期信号的傅里叶变换的DFT 逼近
连续时间非周期信号x (t ) 的傅里叶变换对为：
X (j Ω) =⎰x (t ) e -j Ωt dw (2-1) -∞∞
1x (t ) =2π⎰∞
-∞e j Ωt d Ω (2-2)
计算方法如下：
(1)将x(t)在x 轴上等间隔分段，每一段用一个矩形脉冲代替，脉冲的幅度为其起点的抽样值x (t ) =x (nT ) =x (n ) , 然后把所有矩形脉冲的面积相加，而
t →nT
dt →T
⎰
∞-∞dt →n =-∞∑T
∞
∞从而可得频谱密度X (j Ω) =⎰x (t ) e -∞-j Ωt dw 的近似值为：
x (nT ) e -j ΩnT X (j Ω) ≈
n =-∞∑∞⨯T (2-3)
(2)将序列x (n ) =x (nT ) 截断成从t =0开始长度为T 0的有限长序列，包含有N 个抽样，则上式(2-3)变为：
X (j Ω) ≈T ⨯∑x (nT ) e -jnT Ω
n =0N -1
由于时域抽样的抽样频率为f s =1/T ，则频域产生以f s 为周期的周期延拓。如果频域是限带信号，则有可能不产生混叠，成为连续周期频谱序列，频域周期为f s =1/T （即时域的抽样频率）。
(3)为了数值计算，在频域上也要离散化（抽样），即在频域的一个周期(f s ) 中也分成N 段，即取N 个样点f s =NF 0，每个样点间的间隔为F 0。频域抽样，则频
域的积分公式(2-2)就变成求和式，而时域就得到原已截断的离散时间序列的周期延拓序列，其时域周期为，这时Ω=k Ω0，则
d Ω=(k +1) Ω0-k Ω0=Ω0
⎰
参量关系为 ∞-∞d Ω→∑Ω0 k =0N -1
T 0=1/F 0=N /f s =NT
这样，进过上面三个步骤，时域，频域都是离散序列。推导如下：
经过前两步：时域抽样，截断，即
X (j Ω) ≈
n =-∞∑∞x (nT ) e -j ΩnT ⨯T
x (nT ) ≈1Ωs j ΩnT X (j Ω) d Ωe ⎰02π
再由第三步：频域抽样，则得
X (jk Ω0) ≈T ⨯∑x (nT ) e -jknT Ω0=T ⨯DFT [x (n )]
n =0N -1
Ω0N -1x (nT ) ≈X (jk Ω0) e jknT Ω0 ∑2πk =0
=F 0∑X (jk Ω0) e
k =0N -1jkn 2πN 2πjkn 1N -1
=f s ⨯∑X (jk Ω0) e N N k =0
=f s ⨯IDFT [x (jk Ω0)]
即
X (jk Ω0) ≈T ⨯DFT [x (n )]
x (n ) =1IDFT [X (jk Ω0)] T
此为从离散傅里叶变换法求连续非周期信号的傅里叶变换的抽样值的方法。其模型如图2-1。
取一个周期抽样截断周期延拓*(n ) ←−−−−−−−→(n ) (t ) −−−→x (n ) −−−→x(n ) d (n ) −−−−→−x a x N x N 周期延拓
↓
−→X (e )*D(e ) →X a (j Ω) →X (e ) −−卷积jw jw jw X 取一个周期*(k ) ←−−−−−−−→−X N (k ) N 周期延拓
图2-1 连续非周期信号的傅里叶变换模型
2.2.3 对连续时间周期信号的傅里叶级数的DFS 逼近
连续时间周期信号x(t)的傅里叶级数为
X (jk Ω0) =1T 0-jk Ω0t X (t ) dt e ⎰0T (2-4)
x (t ) =
k =-∞∑∞jk t X (jk Ω0) e Ω0
(2-5)
这里T 0为连续时间周期信号的周期。
由于满足：时域周期→频域离散
时域连续→频域非周期
所以，需将连续周期信号与DFS 联系起来，就需要进行如下操作：
先对时域抽样x (n ) =x (nT ) =x (t ) ，当t =nT 时，则
dt =(n +1) T -nT =T
⎰
T 0
dt →∑T
n =0
N -1
设一个周期内的样点数为N ，则式(2-4)变为
2π-jkn T N -11N -1-jknT Ω0
X (jk Ω0) ≈∑x (nT ) e =∑x (n ) e N
T 0n =0N n =0
在T 0（一个周期）时间段内抽样间隔为T ，共有N =
T 0
个抽样点。 T
再将频域离散序列加以截断，使它成为有限长序列，如果这个截断长度正好等于一个周期（即时域抽样造成的频域周期延拓的一个周期），则式(2-5)变成 x (nT ) ≈∑X (jk Ω0) e jknT Ω0
n =0N -1
=∑X (jk Ω0) e
n =0
N -1
jk
2π
n N
(2-6)
按照DFT 的定义，可得 X (jk Ω0) =
1
⨯DFS [x (n )] N
这就是用DFS 来逼近连续时间周期信号傅里叶级数对的公式。
2.2.4 利用DFT 对非周期连续时间信号傅里叶变换逼近的误差分析
我们已经讨论了用数学表达式分析了逼近过程，需关注的是最后得到的X n (n ) 是否是X a (t ) 的抽样？以及X n (n ) 的DFT 的系数X N (k ) 是否是X a (t ) 的频谱X a (j Ω) 的准确抽样？如果不是的话，误差产生在哪些运算过程中，如何减少这些误差？ 前面已经谈到，从理论上说，信号的时域长度和频域带宽不可能同时是有限长的，也就是说，时域是有限长的，则频域带宽是无限的，反之亦对。而DFT 是有限长离散时间序列与有限长离散频域序列之间的 傅里叶变换对，并隐含周期性，因而用DFT 来逼近连续时间傅里叶变换对肯定是有误差的。 从模型中可看出，最主要的是三个处理：
(1)时域抽样成x(n)则频域就会产生周期性延拓。若频域是限带宽的，只要抽样频率f s 大于信号最高频率f h 的两倍，就不会产生频域的频域的频域响应混叠现
象。由于频域限带，故时域X a (t ) 一定是无线长的。
(2)时域截短成有限长序列x (n ) d (n ) ，此序列只在0≤n ≤N -1的范围内和原。时域截断就是相乘，则x (n ) 相同，在此范围外为0（当用矩形序列R N (n ) 截短时）
在频域表现为周期性卷积，即X (e jw )*D (e jw ) ，卷积的结果造成频谱的“扩散”，也就是频谱的“泄露”，使得形成的信号频谱与原来X (e jw ) 的频谱并不相同，也就是说，时域截断造成频域泄露，使频域产生误差。
(3)由于频域仍是连续的，要对频域抽样，即频域也要离散化。频域抽样造成时域x (n ) d (n ) 的周期延拓，只要频域抽样间隔F 0满足F 0≤
f s
（即频域一个周期抽N
样点数M ≥N ），则这个周期延拓是不会产生时域混叠失真的。
因而，（1）中如果频域不是有限带宽的，则时域抽样会造成频域混叠失真，此时则需尽量选取时域抽样频域足够大，以减小频域混叠失真。（2）中时域截断造成频域泄露，泄露也会造成频域的混叠，因而要选择截断后的N 足够大，使得频域卷积结果更接近X (e jw ) （即当N →∞时，等于不产生截断）。此外，如果不采用突变形式的矩形窗截断，而是采用其他缓变窗（如哈明窗等），则可减小泄露，但这类截断会使x (n ) d (n ) 在0≤n ≤N -1范围内不等于原来的x(n)，使时域产生误差。（3）中如果频域抽样的抽样间隔F 0≥
f s
，则造成时域x (n ) d (n ) 的周期延拓的混叠失真，N
使得在0≤n ≤N -1范围内，时域序列也不等于x(n)。
总之，时域采用矩形窗d(n)截断，则时域只是造成截断区域0≤n ≤N -1外的误差，在区域内x (n ) d (n ) =x (n ) 就是X a (t ) 的抽样。频域若是限带信号，则只有时域截断时形成频域的泄露而造成频域的误差；如果频域不是限带信号，则时域抽样和时域截断都会造成频域的误差。
3 快速傅里叶变换(FFT)
快速傅里叶变换(FFT)只是离散傅里叶变换（DFT ）的一种快速计算方法，其并不是一种新的变换。如前面所讲，有限长序列的特点是其频域可离散化成有限长序列。DFT 的计算在信号处理中非常有用，再有，信号的谱分析对通信，图像传输，声纳等都是很重要的。此外，在系统的分析，设计和实现中都会用到DFT 计算。但是，在很长一段时间里，由于DFT 的计算量极大，即使采用计算机也很难对实际问题进行处理，所以没有得到真正的运用。直到库利和图基提出机器计算傅里叶级数的一种算法，并经过后来人们的改进发展，完善了一套高速有效的运算方法，使得DFT 的计算量得到简化，运算时间一般可以缩短一，二个数量级，从而使DFT 的运算在实践中得到真正广泛应用。
3.1 FFT的来源
若x(n)是有限长序列，点数为N ，其DFT 为
X (k ) =∑x (n ) W
k =0N -1
nk N
，其中，W N =e
-j
2πN
，k=0,1,„N-1 （3-1）
反变换（IDFT ）为
1N -1-nk
x (n ) =∑X (k ) W N ，其中，n=0,1,„,N-1 （3-2）
N k =0
二者的差别只在于W N 的指数符号不同，以及差一个常数因子。因而，直接计算DFT ，乘法次数和加法次数都是和N 2成正比的，当N 很大时运算量显而易见，例如，当N ＝8，DFT 需64次复乘，而当N ＝1024时，DFT 需一百多万次复运算，这对实用性很强的信号处理来说，对计算速度要求是太高了。
nk 下面讨论减少运算量的方法。观察DFT 的运算形式可看出，利用系数W N 的以
下特性，就可以减少DFT 的运算量：
nk
1、W N 的共轭对称性
nk *-nk
(W N ) =W N
nk
2、W N 的周期性
nk (n +N ) k n (k +N )
W N =W N =W N
nk
3、W N 的可约性
nk nkm nk
，W N W N =W Nm =W N
由此可得出
W
(N -k ) n
N
=W
(N -n ) k N
=W
-nk N ，W N =-1，W
N
(k +N )
N
k
=-W N
nk
这样，通过这些性质，使DFT 运算中有些项合并；利用W N 的对称性，周期性
及可约性, 将较长序列的DFT 分解成较短序列的DFT ，从而使DFT 计算过程转化成许多迭代运算过程。而前面已经说到，DFT 的运算量是与N 2成正比的，所以N 越小，运算量越低。
快速傅里叶变换算法正是基于此而发展起来的。它的运算基本上可以分成两大类，即按时间抽选法和按频率抽选法。根据论文要求，我们只考虑前一种方法。
3.2 按时间抽选（DIT ）的基-2FFT 算法（库利-图基算法） 3.2.1 算法原理
取序列点数N =2L ，L 为正整数。如果不满足这个条件，可以加上若干零值点，使达到这一要求，即为基-2FFT 。
对N 点序列，将x (n ) 按n 的奇偶性分组：
x (2r ) =x 1(r ) ，x (2r +1) =x 2(r ) ，其中r=0,1,„, N 2-1
则DFT 转化成
nk
X (k ) =DFT [x (n )]=∑x (n ) W N
n =0N -1
2rk (2r +1) k
=∑x (2r ) W N +∑x (2r +1) W N r =0N
2
r =0N -12r =0
N -1
N -12
=
∑x (r )(W
1r =0
2N
2rk N
2rk
) +∑x 2(r )(W N ) （3-3） 4π
N
利用系数W 的可约性，即W =e
N -12r =0
nk N
-j
=W N ，则
rk k
X (k ) =∑x 1(r ) W N +W N ∑x 2(r ) W N rk =X 1(k ) +W N k X 2(k ) （3-4）
r =0
N
-12
式中X 2(k ) 与X 2(k ) 分别是X 1(r ) 及X 2(r ) 的N/2点DFT ：
rk rk
（3-5） X 1(k ) =∑x 1(r ) W N =∑x 1(2r ) W N
r =0
r =0
N
-12
N -12
rk rk
（3-6） X 2(k ) =∑x 2(r ) W N =∑x 2(2r +1) W N
r =0
r =0
N
-12N -12
由式（3-4）可以看出，一个N 点DFT 分解成两个N/2点的DFT ，它们按（3-4）式又可组合成一个N 点DFT 。但是，x 1(r ) ，x 2(r ) 以及X 1(k ) ，X 2(k ) 都是N/2点的序列，即r,k 满足r,k=0,1,„, N 2-1。而X(k)却有N 点，而用式（3-4）只得X(k)前一半，要用X 1(k ) ，X 2(k ) 来表式全部X(k)值，还需利用系数的周期性质，即
rk W N =W N
r (k +)
这样可得到
N r (k +N ) rk X 1(k +) =∑x 1(r ) W N =∑x 1(r ) W N =X 1(k ) （3-7）
2r =0r =0
同理可得
X 2(k +
N
-12N -12
N
) =X 2(k ) （3-8） 2
式（3-7）、式（3-8）表明后半部段k 值所对应的X 1(k ) ，X 2(k ) 与前半部段k 值所对应的X 1(k ) ，X 2(k ) 分别相等。
nk 在考虑W N 的性质：
(k +N )
N
W N k k =W N W N =-W N （3-9）
这样，将式（3-7）、式（3-8）、式（3-9）带入式（3-4）中，就可将X(k)表示成前后两段：
前半部分X(k)，(k=0,1,„, N 2-1)
X (k ) =
k
X 1(k ) +X 2(k ) W N （3-10）
后半部分X(k)，(k=0,1,„, N 2-1)
N
k +N N N
X (k +) =X 1(k +) +W N 2X 2(k +)
222
k
=X 1(k ) -W N k=0,1,„, N 2-1 （3-11） X 2(k ) ，
这样，只要求出0到(N 2-1) 区间的所有X 1(k ) ，X 2(k ) 值，即可求出0到(N -1）区间内的所有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 分成两个N 2点DFT 时，若直接计算N 2点DFT ，则每个
N 2点DFT 只需要N 2次复数乘法，
N N
(-1) 次复数加法。同时，把两个N 2点22
DFT 合成为一个DFT 时，有N 2个蝶形运算，还需要N 2次复数乘法及N 次复数加法。因而通过这第一步分解后，工作量总体得到大幅降低。
既然如此，由于N =2L ，因而N 2仍是偶数，可以更进一步把每个N 2点子序列按其奇偶性再分成两个N 4点子序列。先将x 1(r ) 进行分解：
⎧x 1(2l ) =x 3(l )
（3-12） ⎨
⎩x 1(2l +1) =x 4(l )
其中，l=0,1,„, N 4-1
2lk k (2l +1)
X 1(k ) =∑x 1(2l ) W N +W N ∑x 1(2l +1)
l =0
r =0
N
-14l =0
N -14
N -14N -14
lk k
=∑x 3(l ) W +W N
l =0
∑x (l ) W
4
k
N
k
X 4(k ) ，k=0,1,„, N 4-1
=X 3(k ) +W
且 X 1(k +
N k
) =X 3(k ) -W N X 4(k ) ，k=0,1,„, N 4-1
4
N
-14l =0
lk
其中 X 3(k ) =∑x 3(l ) W N （3-13）
1
l =0
lk
X 4(k ) =x 4(l ) W （3-14）
图3-2 将N 点DFT 分为两个N 2点的DFT 图
图3-3给出N=8时，将一个N 2点的DFT 分成两个N 4点DFT ，由这两个N 4点DFT 组合成一个N 2点DFT 流图。
同时，x 2(r ) 也可以进行同样的分解，得到
k
X 2(k ) =X 5(k ) +W N X 6(k )
X 2(k +
N k
) =X 3(k ) -W N X 6(k )
4
N
-14l =0
N -14l =0
lk lk
其中 X 5(k ) =∑x 2(2l ) W N （3-15） =∑x 5(l ) W lk lk X 6(k ) =∑x 2(2l +1) W N （3-16） =∑x 6(l ) W N
l =0
l =0
N
-14N -14
k 2k
=W N 最后将系数统一为W N ，则一个N=8点DFT 就可分成四个2点DFT ，这
样就可得到3-4图。
图3-3 由两个N 4点DFT 组合成一个N 2点DFT
根据上面同样的分析知道，最后剩下的是2点DFT ，对于此例N=8，就是四个2点DFT ，其输出为X 3(k ) ，X 4(k ) ，X 5(k ) ，X 6(k ) (k=0,1)，这由式（3-13）至（3-16）可以计算出来。由此可得出一个按时间抽选运算的大致的8点DFT 流图，如图3-5所示。
图3-4 将一个N 点DFT 分成四个N 4点DFT
3.2.2运算量
由图3-5可知，当N =2时(L表示级数) ，每级都由N 2个蝶形运算组成，每
L
个蝶形运算有一次复乘，二次复加，因而每级运算都需N 2次复乘和N 次复加，因此L 级运算共需
复乘数 m f =
N N
L =log 2N 22
复加数 a f =NL =N log 2N
因为W N =1，W N
N =-j ，这几个系数都不用乘法运算。此外，当N 较大时，
这些特例相对而言就很少。直接DFT 算法运算量，复数乘法：N 2，复数加法：
N (N -1) 。直接计算DFT 与FFT 算法的计算量之比为M =
2N
。 log 2N
图3-5 N=8按时间抽选FFT 运算流图
3.2.3 按时间抽选FFT 算法的特点
(1)原位运算（同址运算）
从上图可以看出该运算是有规律可循的，其每级运算都是由N 2个蝶形运算构成，下述迭代运算由每一个蝶形结构实现：
r
⎧⎪X m (k ) =X m -1(k ) +X m -1(j ) W N
（3-17） ⎨r
⎪⎩X m (j ) =X m -1(k ) -X m -1(j ) W N
式（3-17）的蝶形运算如图3-6所示。
从图3-5中看出，某一列的任何两点k 和j 的运算后，得到结果为下一列k,j 两节点的变量，而与其他节点变量无关，通过原位运算，经过蝶形运算，其结果为另一列数据，它们以蝶形为单位仍存储在这一组存储器中，中间无需其他存储器。这样只需N 个存储结构单元，虽然进入蝶形的组合关系有差异，但下一级的运算也可采用这种运算方式。这种结构只需少量存储单元，从而降低设备成本。
图3-6 蝶形运算结构
(2)倒位序规律
按照原位计算，FFT 的输出X(k)在存储单元中按顺序排列，即X(0),X(1)，... ，X(7)的顺序排列，但是这时输入x(n)却不是按自然顺序存储的，按x(0),x(4)，... ，x(7)的顺序依次存入存储单元，看起来好像是“杂乱无序”的，实际是有规律可循的，即倒位序。
造成该原因是由输入x(n)按标号n 的偶奇的不断分组引起。这种不断分成奇偶子序列的过程可用二进制树状图3-7来描述。这就是DIT 的FFT 算法输入序列的 序数成为倒位序的原因。
(3)倒位序的实现
一般实际运算中，总是先向存储单元中按顺序输入序列，为形成倒位序的排列，这通过变址来完成。若序列序号n 用二进制数表示，顺序为(n 2, n 1, n 0) ，那么其倒位序二进制数应当是(n 0, n 1, n 2) 。表3-1列出了N=8时自然顺序以及相应的倒位序的二进制数。
(4)蝶形运算节点的“距离”
仍以N=8时FFT 运算流图为例，其输入是倒位序的，输出是自然顺序的，第一级每个蝶形两节点间“距离”是4，由此类推得，对N =2L 点FFT ，当输出为正常顺序，输入为倒位序时，其第ｍ级运算，每个蝶形两节点“距离”则为2m -1。
nk (5)W N 的确定
对第m 级，一个蝶形运算两节点的“距离”为2m -1，于是第m 级的一个蝶形计算可写成
r
X m (k ) =X m -1(k ) +X m -1(k +2m -1) W N
r
X m (k +2m -1) =X m -1(k ) -X m -1(k +2m -1) W N
图3-7 描述倒位序的树状图
表3-1 自然顺序及相应的倒位序
4 数字信号处理 MATLAB实现的基本知识
4.1 MATLAB简介
MATLAB软件最初是由美国新墨西哥大学计算机系Cleve Moler 教授用FORTRAN
语言编写的矩阵计算软件，目前该软件是Math Works公司的产品。MATLAB 全称是Matrix Laboratory，该软件以矩阵的形式来描述和处理所有的计算机数据，主要由MATLAB 语言、工作环境、图形系统、数学函数库、API 函数和工具箱等六大部分组成、前面五部分是其核心内容工作箱子是其辅助内容。第一部分是MATLAB 语言，与C 语言、FORTRAN 等高级语言近似，用于实现各种算法，MATLAB 包含的编译器还可以将MATLAB 的算法和应用程序文件转变成独立可执行的应用程序。第二部分是MATLAB 工作环境，用户可以完成程序设计、数值计算、图形绘制、输入输出、文件管理等多项功能。第三部分是图形系统，可以完成2D 和3D 数据显示、动画生成、视频处理、图形显示等功能。第四部分是基础与通用的MATLAB 数据函数库，其功已经大大超越了常用高级语言的基本数学函数库。第五部分是MATLAB API函数，这使得C 语言、FORTRAN 语言可以与MATLAB 混合应用，这种交互可以取长补短，提高程序的运行效率，丰富程序开发的手段。第六部分是各种功能极其强大的专业工具箱，涉及到数据获取、数字信号处理、数字图像处理、控制系统分析和设计、最优化方法、系统辨识、智能计算、金融财务分析以及生物信息处理等内容。这些工具箱的算法是开放可拓展的。第三方和用户可根据需要进行拓展。
目前，MATLAB 已经在自动控制原理、数字信号处理、数字图像处理、时间序列分析、数理统计、动态系统仿真等课程教学中得到广泛应用，是理工类本科学生和研究生必须掌握的基本工具。不仅如此，MATLAB 在科学研究中被公认为准确、可靠、方便的科学计算软件；在工业部门，MATLAB 被认为是高效仿真、分析和开发的软件工具。基于此，本论文以MATLAB 为基本，做频谱分析。
4.2利用Matlab 计算FFT 的子函数
用Matlab 方法计算信号的DFT 时，主要是用函数fft (x , N ) 和ifft (x , N ) 。对于这两个函数，如果N 为2的正整数幂，则可以得到本章中介绍的基2 FFT 快速算法；如果N 既不是2的正整数幂，也不是质数，则函数将N 分解成质数，得到较慢的混
合基 FFT算法；如果 N 为质数，则fft 函数采用原来的 DFT 算法。 实验涉及的MATLAB 子函数
1、fft
功能：一维快速傅里叶变换
调用格式：y =fft (x ) ；利用FFT 算法计算矢量x 的离散傅里叶变换，当x 为矩阵时，y 为矩阵x 每一列的FFT 。当x 的长度为2的幂次方时，则fft 函数采用基2的FFT 算法，否则采用稍慢的混合基算法。
y =fft (x , n ) ；采用n 点FFT 。当x 的长度小于n 时，fft 函数在x 的尾部补零，
以构成n 点数据；当x 的长度大于n ，fft 函数会截断序列x 。当x 为矩阵时，fft 函数按类似的方式处理列长度。
2、ifft
功能：一维快速傅里叶变换
调用格式：y =ifft (x ) ; 用于计算矢量x 的IFFT 。当x 为矩阵时，计算所得的y 为矩阵x 中每一列的IFFT 。
y =ifft (x , n ) ; 采用n 点IFFT 。当length (x ) n 时，
将x 截断，使length (x ) =n 。
用MATLAB 提供的子函数进行快速傅里叶变换时，从理论学习可知，DFT 是唯一在时域和频域均为离散序列的变换方法，它适用于有限长序列。尽管这种变换方法是可以数值计算的，但如果只是简单的按照定义进行数据处理，当序列长度很大时，则将占用很大的内存空间，运算时间将很长。
快速傅里叶变换是用于DFT 运算的高效运算方法的统称，FFT 只是其中的一种，FFT 主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N 的序列分解成多个短序列，如基2算法，大大缩短了运算时间。
4.3 利用MATLAB 实现信号仿真
用FFT 进行无限长序列的频谱计算，首先要将无限长序列截断成一个有限长序列。序列长度的取值对频谱有较大的影响，带来的问题是引入频谱的泄漏和波动。
例1 给出信号为x (t ) =cos(20πt ) +cos(8πt ) +cos(16πt ) ，是一个连续信号，所以我们根据其最高频率确定采样速率f s =64Hz ，以及由频率分辨率选择采样点数N=64、32、16三种情况进行分析。
解答 程序如下:
N=input('选择FFT 变换区间长度N ：16 or 32 or 64:'); %**输入FFT 变换区间长度N**%
fs=64; n=0:N-1; t=n/fs;
x=cos(20*pi*t)+cos(8*pi*t)+cos(16*pi*t); subplot(2,2,1);
plot(t,x); %**作信号时域波形**% xlabel('t'); ylabel('x(t)'); title('x(t)时域波形');
f=fft(x,N); %**计算FFT**% subplot(2,2,2); %**绘制x(n)时域图像**% M=N*0.4; n=0:N-1; stem(n,x,'.'); n=0:N-1; m=zeros(N); hold on; plot(n,m); t=max(x); xlabel('n') ylabel('x(n)') title('x(n)时域波形');
subplot(2,2,3); %**绘制N 点FFT 图**% n=0:N-1; stem(n,f,'.'); t=max(f); xlabel('K')
string2=['x(n)的N=',sprintf('%2.0f',N),'点FFT']; title(string2); ylabel('|X(k)|'); 运行如下：当输入16时，
图4-1 N=16时x(t)的DFT
当输入32时，
图4-2 N=32时x(t)的DFT
当输入64时，
图4-3 N=32时x(t)的DFT
分析：
1、T s 相同，当N 越大，在所给范围内等间隔抽样点数越多，且时域信号的长度L NT s 保留得越长，则分辨率越高，频谱特性误差越小。反之，当N 越小时，等间隔采样点数越少，则有可能漏掉某些重要的信号分量，称为栅栏效应。
2、图4-1中，在点数n=1,2,14,15外，还有其他值，而在N=32，N=64时只有几根谱线有值，其它处为0，原因是图4-1中信号发生泄漏，信号的频率成分泄漏到周围一些离散频率点上。
例2 求正，余弦信号谱。假设采样频率为64hz ，采样点N=32，正弦信号频率为10HZ ，余弦信号频率为20HZ ，并验证有关性质。
解答 程序如下： clear all; N=32;n=0:N-1; Fs=64;T=1/Fs; x1=cos(20*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N));
subplot(4,2,1);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('余弦序列'); subplot(4,2,2);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs; x1=sin(10*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,3);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('正弦序列'); subplot(4,2,4);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs;
x1=sin(10*2*pi*n*T)+cos(20*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,5);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('合序列'); subplot(4,2,6);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); grid on; N=32;n=0:N-1; Fs=64;T=1/Fs;
x1=cos(20*2*pi*n*T)+j*sin(10*2*pi*n*T); k=0:N-1; X1=abs(fft(x1,N)); subplot(4,2,7);stem(n,x1);
xlabel('n');ylabel('x1(n)');title('复数序列'); subplot(4,2,8);stem(k,X1); xlabel('k');ylabel('X1(k)'); title('32点FFT 幅频曲线'); 结果如下：
余弦序列
n
正弦序列
32点FFT 幅频曲线
X 1(k )
k
32点FFT 幅频曲线
x 1(n )
n 合序列
X 1(k )
x 1(n )
k
32点FFT 幅频曲线
n
复数序列
X 1(k )
x 1(n )
k
32点FFT 幅频曲线
n
X 1(k )
x 1(n )
k
图4-4 各种序列频谱
分析：
1、对单一连续信号进行傅里叶变换，首先，将时域信号进行采样，显然，f s f h ，满足题意。由x(t)到x(n)，可得最小采样点数为10，又因使用基-2FFT 算法，故最小采样点数为16，这里N 取32符合要求。
2、图4-4中，正弦或余弦信号频谱只有两个点不等于0，这正说明它们自身的频率。
3、由合序列可知，两信号时域叠加，其频谱亦为线性叠加。
4、复数序列实数部分的离散傅立叶变换是原来序列离散傅立叶变换的共轭对称分量；复数序列虚数部分的离散傅立叶变换是原来序列离散傅里叶变换的共轭反对称分量。这里讨论的x(n)与X(k)均为有限长序列，所以这里的对称性是指关于
N 2点的对称性。
例3 绘制一个正弦信号频谱图并作相关分析，参数自定。 解答 程序如下： clear;
fs=100;N=1024; %采样频率和数据点数 A=20;B=30;; n=0:N-1;t=n/fs; %时间序列 x=A*sin(2*pi*B*t); %信号
y=fft(x,N); %对信号进行傅里叶变换 yy=abs(y); %求得傅里叶变换后的振幅 yy=yy*2/N; %幅值处理 f=n*fs/N; %频率序列
subplot(2,2,1),plot(f,yy); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图1 fs=100 N=1024'); grid on; %两种信号叠加
x=A*sin(2*pi*B*t)+2*A*sin(2*pi*1.5*B*t); %信号 y=fft(x,N); %对信号进行傅里叶变换 yy=abs(y); %求得傅里叶变换后的振幅 yy=yy*2/N; %幅值处理 f=n*fs/N; %频率序列
subplot(2,2,2),plot(f,yy); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图2 fs=100,N=1024 两种信号叠加'); grid on; %加噪声之后的图像
x=A*sin(2*pi*B*t)+28*randn(size(t)); y=fft(x,N); yy=abs(y); yy=yy*2/N; %幅值处理
subplot(2,2,3),plot(f(1:N/2.56),yy(1:N/2.56)); xlabel('频率/\itHz'); ylabel('振幅'); title('图3 fs=100,N=1024混入噪声');
grid on; %改变采样点数
N=128; n=0:N-1;t=n/fs; %时间序列
x=A*sin(2*pi*B*t); %信号
y=fft(x,N); %对信号进行傅里叶变换
yy=abs(y); %求得傅里叶变换后的振幅
yy=yy*2/N; %幅值处理
f=n*fs/N; %频率序列
subplot(2,2,4),plot(f(1:N/2.56),yy(1:N/2.56)); %绘出随频率变化的振幅
xlabel('频率/\itHz'); ylabel('振幅'); title('图4 fs=100 N=128');
结果如下：
图4-5 信号频谱
分析：
1、图1和图2取相同的采样频率fs=100和数据点数N=1024，不同的是图2采用两种不同幅值和频率的正弦信号叠加。从图中可以看出图2可以明显的看出含有两种不同的频率成分的信号，幅值也不相同。由此可以看出，不同频率的正弦信号叠加，在频域当中互相分离互不影响。
2、从图1和图2可以看出:整个频谱图是以fs/2频率为对称轴的。由此可以知道傅里叶变换数据的对称性。因此在用傅里叶变换做频谱分析时，我们只需做出前一半频谱图即可。
3、图3为混入噪声之后的频谱。可以看出噪声分布在整个频率轴上，并且由于噪声中含有与原信号频率相同的成分，叠加之后导致幅值增加。加大噪声的幅值之后将分辨不出原信号的频率。
4、图4减少了数据点数，N=128。与图1相比较，采用128点和1024点的相同频率的振幅是有不同的表现值。因此振幅的大小与所用采样点数有关。一定范围内采样点数越多。信号的幅值越接近真实值。
5 总结与展望
5.1 总结
通过这次毕业设计，使我充分理解了信号处理的重要性，及MATLAB 软件编程方面的知识。在学习过程中，虽然遇到了许多问题，但在老师的指导，自己翻阅大量资料，一步步找出问题所在，同时也暴露出之前在学习数字信号处理等科目的不足。实践是检验真理的唯一标准，通过动手操作，理论与实际结果相结合，才能更好的理解和掌握知识，为我所用。
对于我们通信专业，我认为信号处理是基本专业知识，这个专业以后就业的方向也很多，就业面很广。我们毕业以后工作，可以进入设备制造商、 运营商、专有服务提供商以及银行等领域工作。当然, 就业形势每年都会变化，所以关键还是要看自己。可以从事硬件方面，比如说PCB ，别小看这门技术，平时我们在试验时制作的简单，这一技术难点就在于板的层数越多，要做的越稳定就越难，这可是非常有难度的，如果学好了学精了，也是非常好找工作的。也可以从事软件方面，这实际上也要求我们具备比较好的模电，数电等基础知识。
在今后社会发展和学习过程中，我相信信号处理这方面知识一定大有用处。只有不懈努力，遇到问题不要退缩，然后一一进行解决，才能取得成功。
5.2 展望
目前，信号处理课群体系正逐步成熟，并得到国内高校的认可，其体现了理论与实践的有机结合，体现了原理、方法和技术的有机结合。在介绍数字信号处理的理论和方法的基础上，进行MATLAB 仿真实验，再进行基于DSP 系统的开发应用实验。
信息处理前景可谓前途无可限量, 随着社会的发展, 信息交流将越发频繁, 在信息充斥的世界里, 信息处理将越发被需要, 发展前景广阔。
致 谢
本文是在桂静宜老师的谆谆教诲和悉心指导下完成，桂老师在学术上有着严谨的教学作风，实事求是的治学态度，让我受益匪浅。生活上桂老师平易近人，和蔼可亲，是我学习和生活中的榜样。在毕业设计研究之中，桂老师给了我最及时和最有效的指导，这使得我最终克服各种困难，顺利地完成了论文。在此，谨向我的老师表示感谢。
同时，也感谢我的同学一直对我的鼓励和支持，在这次论文中也给了我不少建议，使我在学习上不断进步。最后，《基于FFT 的连续信号谱分析》得以基本完成。
参考文献
[1] 高西全, 丁玉美著. 数字信号处理（第三版). 西安电子科技大学出版社
[2]刘小东著. MATLAB7.X的系统分析与设计——信号处理. 西安电子科技大学出版社
[3]郭文彬, 桑林编著. 庞兴华审订. 通信原理—基于MATLAB 的计算机仿真. 第一版. 北京：北京邮电大学出版社2006年6月
[4]郑南宁. 数字信号处理（第二版）. 西安交通大学出版社. 1995
[5]邓华编著. MATLAB通信仿真及应用实例详解. 第一版. 北京：人民邮电出版社. 2003年9月
[6]邵玉斌编著. Matlab/Simulink通信系统建模与仿真实例分析. 第1版. 北京：清华大学出版社2008年6月
[7]唐向宏.MATLAB 及在电子信息类课程中的应用. 第四版. 西安：电子科技大学出版社，2006年
[8]樊昌信, 詹道庸, 徐炳祥, 吴成柯. 通信原理. 第四版. 北京:国防工业出版社. 2000年
[9]Schilling R J and Harris S L. Fundamentals of digital signal processing using MATLAB .西安：西安交通大学出版社. 2005
[10]楼顺天, 李博菡. 基于MATLAB 的系统分析与设计——信号处理. 西安电子科技大学出版社. 2000
[11]Paulo S R.Diniz ,Eduardo A.B.da Silva . Digital Signal Processing and System Analysis and Design. 电子工业出版社. 2002
[12]刘树棠译. Oppenheim A V,Willsky S and Nawab S H.信号与系统（第二版). 西安交通大学出版社. 2001
[13]胡广书. 数字信号处理-理论、算法、与实现(第二版). 北京：清华大学出版社.2003
[14]党宏社. 信号与系统实验(MATLAB版). 西安：西安电子科技大学出版社.2007
[15]陈怀琛. 数字信号处理教程—MATLAB 释义与实现. 北京：电子工业出版社.2005