基于FFT的连续信号谱分析毕业设计论文

毕业设计（论文）题目：基于FFT的连续信号谱分析
摘要
离散傅里叶变换(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 引言 (2)
1.1 数字信号处理概述 (2)
1.2 连续信号的频谱分析 (2)
1.3 谱分析的研究意义 (3)
2 离散傅里叶变换（DFT） (3)
2.1 离散傅里叶变换的性质 (4)
2.2 利用DFT计算模拟信号的傅里叶变换 (4)
2.2.1 连续信号谱分析原理 (4)
2.2.2 对连续非周期信号的傅里叶变换的DFT逼近 (5)
2.2.3 对连续时间周期信号的傅里叶级数的DFS逼近 (7)
2.2.4 利用DFT对非周期连续时间信号傅里叶变换逼近的误差分析 (8)
3 快速傅里叶变换(FFT) (9)
3.1 FFT的来源 (9)
3.2 按时间抽选（DIT）的基-2FFT算法（库利-图基算法） (10)
3.2.1 算法原理 (10)
3.2.2 运算量 (16)
3.2.3 按时间抽选FFT算法的特点 (17)
4 数字信号处理MATLAB实现的基本知识 (19)
4.1 MATLAB简介 (19)
4.2 利用Matlab计算FFT的子函数 (20)
4.3 利用MATLAB实现信号仿真 (21)
5 总结与展望 (28)
5.1 总结 (28)
5.2 展望 (28)
致谢 (28)
参考文献 (29)
1 引言
1.1 数字信号处理概述
数字信号处理是从20世纪60年代以来，随着信息学科和计算机学科的高速发展而迅速发展起来的一门新兴学科。

它的重要性日益在各个领域的应用中表现出来。

数字信号处理是把信号用数字或符号表示的序列，通过计算机等处理设备，用数字的数值计算方法处理已达到提取有用信息便于应用的目的。

信号可以从不同角度进行分类，如周期信号和非周期信号，确定信号和随机信号，能量信号和功率信号，连续信号和离散信号等。

对于本课题，我们需要了解的是连续时间信号，即时间是连续的，幅值可以使连续的也可以是离散的，特别地，我们一般所说的连续信号指的是模拟信号，即时间，幅值均连续，它是上一种信号的特例。

于是，我们讨论将模拟信号转化成数字信号的过程。

首先把模拟信号变换为数字信号，然后在进行数字化处理，最后在恢复出模拟信号。

其方框图如图1-1所示：
前置滤波器A/D变换器信号处理器D/A变换器模拟滤波器
图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 序列对应关系 x(n)[或X(k)] X(k)[或x(n)]
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 Ω=+Ω-Ω=Ω
1
00N 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 --Ω=Ω≈⨯=⨯∑
01000
()()2N jknT k x nT X jk e π-Ω=Ω≈Ω∑ 21000()N jkn N k F X jk e
π-==Ω∑ 21
00
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
00N T n dt T -=→∑⎰ 设一个周期内的样点数为N ，则式(2-4)变为
0211000
01()()()N N jkn jknT N
n n T X jk x nT e x n e T N π----Ω==Ω≈=∑∑
在0T （一个周期）时间段内抽样间隔为T ，共有0T 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 ）为
10
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,…, 21N -
则DFT 转化成
1
()[()]()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
0()()()()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） 11222222
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)值，还需利用系数的周期性质，即
()
2
2
2
N r k rk
N N W W
+=
这样可得到
1122
()211112200
()()()()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 1()X k
2()X k
12()()k
N X k W X k +
12()()k N X k W X k -
点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）
14
444
()()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） 1(0)(0)x x =
1(1)(2)
x x =1(2)(4)
x x =1(3)(6)
x x =2(1)(3)
x x =2(0)(1)x x =2(2)(5)
x x =2(3)(7)
x x =2N 点DFT (1)
X (0)
X (2)
X (3)
X (4)
X (6)
X (7)
X (5)
X 2N 点DFT
最后将系数统一为22
k k
N N
W W =，则一个N=8点DFT 就可分成四个2点DFT ，这样就可得到3-4图。

图3-3 由两个4N 点DFT 组合成一个2N 点DFT
根据上面同样的分析知道，最后剩下的是2点DFT ，对于此例N=8，就是四个2点DFT ，其输出为3()X k ，4()X k ，5()X k ，6()X k (k=0,1)，这由式（3-13）至（3-16）可以计算出来。

由此可得出一个按时间抽选运算的大致的8点DFT 流图，如图3-5所示。

31(0)(0)(0)x x x ==
41(0)(1)(2)x x x ==
31(1)(2)(4)x x x == 41(1)(3)(6)x x x ==
4N 点DFT
4N 点DFT
1(0)X
1(1)
X 1(2)X 1(3)
X
图3-4 将一个N 点DFT 分成四个4N 点DFT
3.2.2运算量
由图3-5可知，当2L
N =时(L 表示级数)，每级都由2N 个蝶形运算组成，
每个蝶形运算有一次复乘，二次复加，因而每级运算都需2N 次复乘和N 次复加，因此L 级运算共需
复乘数 2log 22
f N N
m L N =
= 复加数 2log f a NL N N ==
因为0
1N
W =，4
N N
W j =-，这几个系数都不用乘法运算。

此外，当N 较大时，
这些特例相对而言就很少。

直接DFT 算法运算量，复数乘法：2N ，复数加法：
31(0)(0)(0)x x x ==
31(1)(2)(4)x x x ==
41(0)(1)(2)x x x ==
41(1)(3)(6)x x x ==
52(0)(0)(1)x x x == 52(1)(2)(5)x x x == 62(0)(1)(3)x x x == 62(1)(3)(7)x x x ==
4N 点DFT
4N 点DFT
4N 点DFT
4N 点DFT
(0)
X (1)X
(2)X (3)
X (4)X (5)
X (6)X (7)
X
(1)N N -。

直接计算DFT 与FFT 算法的计算量之比为22log N
M N
=。

图3-5 N=8按时间抽选FFT 运算流图
3.2.3 按时间抽选FFT 算法的特点
(1)原位运算（同址运算）
从上图可以看出该运算是有规律可循的，其每级运算都是由2N 个蝶形运算构成，下述迭代运算由每一个蝶形结构实现：
1111()()()()()()r
m m m N
r
m m m N
X k X k X j W X j X k X j W ----⎧=+⎪⎨=-⎪⎩ （3-17） 式（3-17）的蝶形运算如图3-6所示。

从图3-5中看出，某一列的任何两点k 和j 的运算后，得到结果为下一列k,j 两节点的变量，而与其他节点变量无关，通过原位运算，经过蝶形运算，其结果为另一列数据，它们以蝶形为单位仍存储在这一组存储器中，中间无需其他存储器。

这样只需N 个存储结构单元，虽然进入蝶形的组合关系有差异，但下一级的运算也可采用这种运算方式。

这种结构只需少量存储单元，从而降低设备成本。

x(0) x(4) x(6) x(1) x(5) x(2) x(3) x(7)
X(0) X(1) X(2) X(3) X(4) X(6) X(5) X(7)
图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 用二进制数表示，顺序为210(,,)n n n ，那么其倒位序二进制数应当是012(,,)n n n 。

表3-1列出了N=8时自然顺序以及相应的倒位序的二进制数。

(4)蝶形运算节点的“距离”
仍以N=8时FFT 运算流图为例，其输入是倒位序的，输出是自然顺序的，第一级每个蝶形两节点间“距离”是4，由此类推得，对2L N =点FFT ，当输出为正常顺序，输入为倒位序时，其第ｍ级运算，每个蝶形两节点“距离”则为12m -。

(5)nk
N W 的确定
对第m 级，一个蝶形运算两节点的“距离”为12m -，于是第m 级的一个蝶形计算可写成
111()()(2)m r
m m m N
X k X k X k W ---=++ 1111(2)()(2)m m r
m m m N
X k X k X k W ----+=-+ 1()m X k -
1()m X j -
11()()()r
m m m N
X k X k X j W --=+11()()()r
m m m N
X j X k X j W --=-
图3-7 描述倒位序的树状图
表3-1 自然顺序及相应的倒位序
4 数字信号处理 MATLAB 实现的基本知识
4.1 MATLAB 简介
MATLAB 软件最初是由美国新墨西哥大学计算机系Cleve Moler 教授用
自然顺序
二进制数 倒位序二进制数 倒位序顺序
0 000 000 0 1 001 100 4 2 010 010 2 3 011 110 6 4 100 001 1 5 101 101 5 6 110 011 3 7
111
111
7
x(000)
x(100) x(010) x(110)
x(111)
x(011) x(101) x(001)
1
1
1
1
1
1
1
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时，主要是用函数(,)
ifft x N。

对
fft x N和(,)
于这两个函数，如果N为2的正整数幂，则可以得到本章中介绍的基2 FFT快速算法；如果N既不是2的正整数幂，也不是质数，则函数将N分解成质数，得到较慢的混合基 FFT算法；如果 N 为质数，则fft函数采用原来的 DFT 算法。

实验涉及的MATLAB子函数
1、fft
功能：一维快速傅里叶变换
调用格式：()
=；利用FFT算法计算矢量x的离散傅里叶变换，当x
y fft x
为矩阵时，y为矩阵x每一列的FFT。

当x的长度为2的幂次方时，则fft函数采用基2的FFT算法，否则采用稍慢的混合基算法。

=；采用n点FFT。

当x的长度小于n时，fft函数在x的尾部补y fft x n
(,)
零，以构成n点数据；当x的长度大于n，fft函数会截断序列x。

当x为矩阵时，fft函数按类似的方式处理列长度。

2、ifft
功能：一维快速傅里叶变换
调用格式：()
=;用于计算矢量x的IFFT。

当x为矩阵时，计算所得的
y ifft x。