毕业论文-快速傅里叶变换算法及其在信号处理中的应用

题 专 业学 号

姓 指 导学 院

2015 届毕业设计(论文)

目 快速傅里叶变换算法及其在信号处理中的应用

班 级 2011电子信息工程02 1104030231 名 周汝耀 教 师 华夏讲师 名 称 电气信息学院

2011 年 6 月 9 日

武汉工程大学毕业设计(论文)说明书

摘 要

快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得 的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。 傅里叶变换的理论与方法在“数理方程”、“线性系统分析”、“信号处理、仿真”等很多学科领域都有着广泛应用,由于计算机只能处理有限长度的离散的序列,所以真正在计算机上运算的是一种离散傅里叶变换. 虽然傅里叶运算在各方面计算中有着重要的作用，但是它的计算过于复杂，大量的计算对于系统的运算负担过于庞大，使得一些对于耗电量少，运算速度慢的系统对其敬而远之，然而，快速傅里叶变换的产生，使得傅里叶变换大为简化，在不牺牲耗电量的条件下提高了系统的运算速度，增强了系统的综合能力，提高了运算速度，因此快速傅里叶变换在生产和生活中都有着非常重要的作用，对于学习掌握都有着非常大的意义。

关键词：快速傅氏变换；快速算法；简化；广泛应用

武汉工程大学毕业设计(论文)说明书

Abstract

Fast Fourier Transform (FFT), is a discrete fast Fourier transform algorithm, which is based on the Discrete Fourier Transform of odd and even, false, false, and other characteristics of the Discrete Fourier Transform algorithms improvements obtained. Its Fourier transform theory has not found a new, but in the computer system or the application of digital systems Discrete Fourier Transform can be said to be a big step into. Fourier transform theory and methods in the \equation\other areas have a wide range of applications, as the computer can only handle a limited length of the sequence of discrete, so true On the computer's operation is a discrete Fourier transform. Fourier Although all aspects of computing in the calculation has an important role, but its calculation was too complicated, a lot of computing system for calculating the burden is too large for some Less power consumption, the slow speed of operation of its system at arm's length, however, have the fast Fourier transform, Fourier transform greatly simplifying the making, not in power at the expense of the conditions to increase the speed of computing systems, and enhance the system The comprehensive ability to improve the speed of operation, the Fast Fourier Transform in the production and life have a very important role in learning to master all have great significance.

KeyWords：Fast Fourier Transform; fast algorithm; simplified; widely used

II

武汉工程大学毕业设计(论文)说明书

目 录

摘要……..……………………….………………………………………………….I Abstract ………….…………… ……………………………………………………………II 1.绪论

1.1 选题背景................................................................................................................ 1 1.2 课题研究的意义.................................................. 2

2.快速傅里叶变换原理及性质

2.1快速傅里叶变换原理 .............................................. 3 2.2快速傅里叶变换的优越性 .......................................... 4 2.3快速傅里叶变换的意义 ............................................ 4

3.快速傅里叶变换的算法

3.1快速傅里叶变换算法 .............................................. 6 3.2 Cooley=Tukey FFT算法 ........................................... 8 3.3 Rader-Brennr FFT算法 ........................................... 9 3.4 Goertsel 算法.................................................. 10

4.快速傅里叶变换在信号处理中的理论应用

4.1利用FFT计算连续时间信号的傅里叶变换 ........................... 13 4.2利用FFT计算离散信号的线性卷积 ................................. 17 4.3利用FFT进行离散信号压缩 ....................................... 19 4.4利用FFT对离散信号进行滤波 ..................................... 22 4.5利用FFT提取离散信号中的最强正弦分量 ........................... 24 5.快速傅里叶变换在数字信号分析与处理的实际应用 5.1快速傅里叶变换在喇曼光谱信号噪声平滑中的应用 ................... 29 5.2采用异步实现的快速傅里叶变换处理器 ............................. 31 5.3快速傅里叶算法在哈特曼夏克传感器波前重构算法中的应用 ........... 33

致谢 .................................................... 36 参考文献 ................................................ 37

武汉工程大学毕业设计(论文)说明书

1 绪论

傅立叶变换在生产生活中的重要性非常突出，它将原本难以处理的时域信号相比比较容易地转换成易于分析的频域信号，我们可以利用一些专业工具对这些频域信号进行处理、加工，使信号转化为可以对其进行各种数学变换的数学公式，对其进行处理。最后还可以根据傅立叶反变换将这些频域信号转换成原来的时域信号，这是一种特殊的积分变换。它能够将满足一定条件的某个函数表示成为正弦基函数的线性组合或者积分。然而，它在运算上过于复杂，过于宏大的运算过程，对于一些相对简单的低功耗处理器来说，难以自如应对，因此，快速傅里叶变换则显出了它的优越性。快速傅氏变换（FFT），即离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的[1]。对于计算机处理信号方面上是一大进步。系统的速度不但取决于其本身的速度，而且还在相当大的程度上取决于运用的算法，算法运算量的大小直接影响到对设备的控制质量。通过傅立叶变换（DFT），运用测试软件进行检测，我们可以看出，快速傅里叶变换大大的提高了运算速度，它为各系统的设计方案提供了简单算法，有着非常重要的意义。 1. 1 选题背景

近十多年来，数字信号处理技术同大规模集成电路、数字计算机等，都有了突飞猛进的发展，日新月异，早已成为了一门具有强大生命力的技术科学。因为它本身就具有一系列的优点，所以能够有效地促进工程技术领域的技术改造和科学发展，应用领域也更加地广泛、深入，越来越受到人们的重视。

在信号处理中，离散傅里叶变换（Discrete Fourier Transform，DFT）是比较常用的变换方法之一，它在各种数字信号处理系统中扮演着及其重要的角色。由于离散傅里叶变换（DFT）而发现了频率离散化，可以直接用它来分析信号的频谱、计算滤波器的频率响应、以及实现信号通过线系统的卷积运算等，因而在信号的谱分析等方面有着非常大的作用。

傅里叶变换已经有一百多年的历史了，我们熟知频域分析往往比时域分析更优越，不仅简单明了，而且易于分析较为复杂的信号。但需要用较为精准的数字方法，即DFT，进行谱分析，在快速傅氏变换（FFT）出现以前是不切实际的。

1

武汉工程大学毕业设计(论文)说明书

由于DFT的计算量太大，即使运用计算机也很难对问题进行实时的有效处理，所以DFT并没有得到真正的应用。直到1965年库利（J.W.Cooly）和图基（J.W.Tukey）首次发现DFT的一种快速算法，局面才发生根本性的变化。继库利和图基算法出来之后，桑德（G.Sander）等快速算法也相继出现，又经过其他学者一步步改进，很快就出现了通用型的快速傅里叶变换，简称FFT。快速傅里叶变换（Fast Fourier Transform，FFT）并非与离散傅里叶变换完全不同的另一种变换，而是为了减少DFT计算次数而诞生的一种快速、有效的算法。应当指出的是，也是因为当时电子数字计算机的“落后”条件也促成了这个算法的提出。它使得DFT的运算量大大的缩小简化，它推动了近30年来信号处理技术止步不前的前进发展，成为了数字信号处理应用领域里强有力的工具，为DFT乃至数字信号处理技术的实际应用创造了良好的条件，从而使DFT在实际使用中得以广泛的应用[2]。

数字信号处理器（DSP），是一种可编程的高性能处理器。近年来发展尤为迅速，它不仅应用于数字信号处理方面，而且在图像处理、语音处理、通信等领域得到广泛的应用。之前通用的微处理器在运算速度上已经很难适应信号实时处理的高要求。DSP处理器中集成了高速的乘法硬件，能快速、准确地进行大量数据的乘法以及加法的运算。数字信号处理区别于普通的科学计算与分析，它强调运算的实时性。除了需要普通微处理器所强调的高速运算和控制能力之外，鉴于实时数字信号处理的特点，在处理器结构、指令系统、指令流程上做了很大程度上的改进。

1. 2 课题研究的意义

如上所述，基于对DSP的快速傅里叶变换算法的研究，从而使FFT算法能够有效地在DSP芯片上实现。DSP芯片的出现，使FFT的实现更加方便。多数的DSP芯片都能够在一个指令周期内完成一次乘法和加法，并且提供了专门的FFT指令，完成一次指令的周期只需10ns，使得FFT算法在DSP芯片上实现的速度更加快速。快速傅里叶变换为频谱分析、卷积、相关数字滤波器设计与实现与功率谱计算、传递函数建模、图像处理等，提供了快速有效的运算方法。FFT技术应用DSP芯片，从而可以提供使调制、解调、压缩、解压缩的数据传输更为高效的信号处理解决方案，因而广泛应用于雷达、图像处理、通信、生物医学和声纳领域。

2

武汉工程大学毕业设计(论文)说明书

2.快速傅里叶变换原理及性质

数字信号中的傅里叶变换,通常是采用离散型傅里叶变换(DFT)。DFT 存在的

缺点就是计算量太大,不易进行实时处理。比如，计算一个N 点的DFT ,一般需要N次复数乘法和N(N-1)次复数加法运算.因此,当N较大或要求对信号进行实时处理时,往往很难实现达到所需的运算速度。1965年,J.W.Cooly和J.W.Tukey发现了DFT的一种快速算法,经过后来学者的进一步改进, 很快便形成了一套高效的运算方法,即现在通用的快速傅里叶变换, 简称FFT( The Fast Fourier

nkWNTransform)。快速傅里叶变换的实质是利用式(3-1)中的权函数的对称性和

2周期性,把N点DFT进行一系列分解和组合,使整个DFT的计算过程变成一系列叠代运算过程,使DFT的运算量大大简化,为DFT及数字信号的实时处理和应用创造了非常良好的条件[3]。 2. 1 快速傅里叶变换原理

快速傅里叶变换原理：

1. 将长序列DFT分解为短序列的DFT

2. 利用旋转因子的对称性、周期性、可约性。将时域序列逐次分解为一组子序列，依据旋转因子的特性，由子序列的DFT来实现整个序列的DFT[4]。

其中：快速傅里叶变换分为两种，分为基2时间抽取算法和基2频率抽取算法

基2-时间抽取(Decimation in time)FFT算法

?x[2r]x[k]???x[2r?1]N?1

其中:r=0,1,2? 2 （2-1）

基2-频率抽取(Decimation in frequency)FFT算法

?X[2m]X[m]???X[2m?1] （2-2）

3

武汉工程大学毕业设计(论文)说明书

2. 2 快速傅里叶变换的优越性

设xn为N项的复数序列，依据DFT变换，任一

x(m)的计算都需要有N次复

数乘法和（N?1）次复数加法，而且一次复数乘法等同于四次实数乘法和两次实数加法，同样的，一次复数加法等同于两次实数加法，即使我们把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的

x(m)2, 即N点DFT变换大约就需要N次运算。当N?1024点

甚至更多的时候，需要N2=1048576次运算。在FFT中，利用WN的对称性和周期性，把一个N项序列（设N?2k,k为正整数），分为两个N/2项的子序列，而且每个N/2点的DFT变换需要?N/2?次运算，再运用N次运算把两个N/2点的

2DFT 变换重新组合成一个N点的DFT变换。如此变换以后，总的运算次数就变成了N?2(N/2)2?N?N2/2。承接上面的例子，当N?1024时，总的运算次数就变成了525312次，这样看来，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要Nlog2次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性.当然，FFT提高了运算速度,但是,也对参与运算的样本序列作出了限制,即要求样本数为2^N点.离散傅里叶变换DFT则无上述限制[5]。 2. 3 快速傅里叶变换的意义

傅立叶变换的物理意义：傅立叶变换是数字信号处理技术领域一项很重要的算法。要知道傅立叶变换算法的意义，首先要了解傅立叶变换原理的意义。傅立叶变换原理表明：任何连续测量的时序或信号，都能够表示成为不同频率的正弦波信号的无限叠加。而利用该原理而创立的傅立叶变换算法则利用直接能测量到的原始信号，并以累加方式来计算该信号中不同正弦波信号的频率、相位和振幅。

与傅立叶变换算法对应的是反傅立叶变换算法。该反变换从本质上说也就是一种累加处理，这样便可以将单独改变的正弦波信号转换成一个信号。因此，也可以说，傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号（信号的频谱），我们可以利用一些专业工具对这些频域信号进行加工、处理。

4

n武汉工程大学毕业设计(论文)说明书

最后还可以利用傅立叶反变换将这些频域信号转换成原来的时域信号。

从现代数学的眼光来看，傅里叶变换其实就是一种特殊的积分变换。它能够将满足一定条件下的某个函数表示成为正弦基函数的线性组合或者积分形式。在不同的研究领域里，傅里叶变换具有多种形式各异的变体形式，如连续傅里叶变换和离散傅里叶变换。

在数学领域，尽管最初傅立叶分析是作为热过程的解析分析的工具，但是其思想方法仍然具有典型的还原论和分析主义的特征。\任意\的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式，而正弦函数在物理上是被充分研究而相对简单的函数类：1. 傅立叶变换是线性算子,若赋予适当的范数,它还是酉算子;2. 傅立叶变换的逆变换容易求出,而且形式与正变换非常类似;3. 正弦基函数是微分运算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;4. 离散形式的傅立叶的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;5. 著名的卷积定理指出:傅立叶变换可以化复变换可以利用数字计算机快速的算出(其算法称为快速傅立叶变换算法(FFT))。

正是由于上述的良好性质,傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用[6]。

5

武汉工程大学毕业设计(论文)说明书

3. 快速傅里叶变换的算法

3.1 快速傅里叶变换算法

快速傅里叶变换FFT 是离散傅里叶变换DFT 的一种快速算法，只有FFT 才能在现实中有实际应用的意义。

错误！未找到引用源。X(n)??x0(k)eN;n?0,1,2?N?1

K?0

(3-1)

由(3-1)式可知，对每一个n，计算X(ｎ)须作N次复数乘法及N-1次复数加法，要完成这组变换共需N2错误！未找到引用源。次乘法及N(N-1)次复数加法。但以下介绍的快速傅里叶变换的算法，可大大减少运算次数，提高工作效率。

当N?2r时，n和k可用二进制数表示：

r?1r?2n?2n?2nr?2?r?1

N?1?nk?n0?nr?1nr?2n0错误！未找到引用源。

(3-2)

r?1r?2k?2k?2kr?2?r?1

?k0?kr?1kr?2k0 （3-3）

又记 W?e X(nr?1nr?2???N,则（3-1）式可改写为

10?0n0)??k?1k1?0?1kr?1?0x0(kr?1kr?2k0)Wp （3-4）

式中：P?nk?(2r?1kr?1?2r?2kr?2? WP?W(2 W因为W2rr?1?k0)?(2r?1nr?1?2r?2nr?2??n0)

nr?1?2r?2nr?2??n0)2r?1kr?1W(2r?1nr?1?2r?2nr?2??n0)2r?2kr?2

K0(2r?1nr?1?2r?2nr?2??n0) (3-5)

???W2r?[eN]N?1所以（3-2）可改成Z

10?0X(nr?1nr?2n0)?k?11k1?0?1kr?1?0x0(kr?1kr?2k0)W(2r?1nr?1?2r?2nr?2??n0)2r?1kr?1W(2r?1nr?1?2r?2nr?2??n0)2r?2kr?2WK0(2r?1nr?1?2r?2nr?2??n0) （3-6）

x2(n0n1kr?3k0)??kr?2?0x0(n0kr?2k0)W(2n1?n0)2r?2kr?2 (3-7)

X(nr?1nr?2n0)?xr(n0n1nr?1)

6

武汉工程大学毕业设计(论文)说明书

则式（3-7）即为式（3-6）的分解形式。将初始数据代入式（3-7）的第一个等式，可得每一组计算数据，一般将痗L-1组计算数据代入式（3-7）的第L个等式，计算后可得第L组计算数据（L＝１，２，?，γ），计算公式也可表示为

x1(n0kr?2k0)??k1r?1?0x0(kr?1kr?2k0)W(2r?1nr?1?2r?2nr?2??n0)k0=

xl?1(n0n1nr?20kr?1kr?2k0)?xl?1(n0n1nr?20kr?1kr?2 k0)WP （3-8）

式中P?2r?1nl?1?2r?2nl?2?（3-9） ?2r?1n0

0根据式（3-8），第L个数组中每个xl(k)?xl(n?r1n?r2nk?r1k?r20)k

的计算只依赖于上一个数组的两个数据这两个数据的标号相差2Y?1?N/2l，即

j?i?n/2l，而且这两个数据只用于计算第L个数组中标号的数据（等号右端为二进制数）。当nl?1分别取０和１时，分别有k?i,k?j?i?n/2l。因此，用上一组的两个数据计算所得的两个新数据仍可储存在原来位置，计算过程中只需要N个存储器。将xl(i)与xl(i?n/2l)称为第L个数组中的对偶结点对。计算每个对偶结点对只需一次乘法，事实上由式（3-8）可得

xl(i)?xl?1(i)?[i?Np1]W 2l

xl(i?NNp2)?x(i)?x[i?]W l?1l?12l2l （3-10）

r?2r?lr?lr?2r?lnP?2n?...?2P?2?2n?...?2n0别为式0l?2l?2式中：1 ；2（3-9）

中

nl?1R?1pP?P?2?P211?N/2，取０，１时对应的P值。因于是对偶结点的W有

如下关系：

WP2?WP1?N2?[e???N]P1?N2??WP1,因此式（3-8）可表示为

N]Wpl2NN（3-11）xl(i?l)?xl?1(i)?xl?1[i?l]Wp22

xl(i)?xl?1(i)?xl?1[i? 7

武汉工程大学毕业设计(论文)说明书

P的求法：在xl(i)中，i写成二进制数n0n1?nl?1kr?l?1?k0右移r?l位，就成为

0?0n0n1?nl?1颠倒位序得p?nl?1?n1n00?0(l?1,2?r)式（3-7）中，前面的r个等式，每个等式均对应一组数据进行计算，每组数据都有N/２对结点，根据式（3-11），每对结点只需作１次乘法和２次加法，因此，每组数据只需N/2次乘法和N次加法，因而完成r组数据的计算共需Nr/2次乘法和Nr次加法。 3. 2 Cooley-Tukey FFT算法

作N点FFT时，若N不是素数，FFT的核心是将一层运算映射为二层运算。则N可分解为N?N1N2，那么由f[n]的DFT

nk F[k]??f[n]WNn?0N?1 0?k?N?1 （3-12）

通过映射： 可得到

(N2n1?n2)(k1?N1k2)(N2n1k1?N1N2n1k2?n2k1?N1n2k2)nk WN?WN?WNn?N2n1?n2k?k1?N1k20?n1?N1?1,0?n2?N2?10?k1?N1?1,0?k2?N2?1 （3-13）

N1而N?N1N2，WN?WN2，WNN2?WN1，可化简为

n1k1n2k1n2k2nk WN ?WNWWNN2 （3-14）1从而式（3-12）转化为 F[k1,k2]?其中

N2?1n2?0?Wn2k2N2(Wn2k1NN1?1n1?0?f[n,n]W12n1k1N1 ) （3-15）

0?k1?N1?1,0?k2?N2?1。

以20点FFT为例，N?20,N1?5,N2?4，映射方式为：n?4n1?n2，

k?k1?5k2，则计算流图如图3-1所示。

8

武汉工程大学毕业设计(论文)说明书

n1 k2 f[0] 0 W0 0 F[0] f[4] 1 W0 k 1 =0 1 F[5] 0 f[8] 2 n 2 =0 W2 F[10] f[12] 3 W0 3 F[15] f[16] 4 W0 0 F[1] f[1] 0 W1 k =1 1 F[6] 1f[5] 1 W2 2 F[11] f[9] 2 n 2 =1 W3 3 F[16] f[13] 3 f[17] 4 W0 0 F[2] k=2 W2 1 1 F[7] f[2] 0 W4 2 F[12] f[6] 1 W6 3 F[17] n2=2 f[10] 2 f[14] 3 W0 0 F[3] k1=3 f[18] 4 W3 1 F[8] W6 2 F[13] f[3] 0 W9 3 F[18] f[7] 1 n2=3 0f[11] 2 W k =4 0 F[4] 1f[15] 3 W4 1 F[9] f[19] 4 W8 2 F[14] W12 3 F[19] 图3-1 Cooley-Tukey 20 点FFT算法

3. 3 Rader-Brenner FFT算法

Rader-Brenner 算法是类似于 Cooley-Tukey 算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey所取代，与Rader-Brenner算法相比，有一样多的乘法量，却有较少的加法量，且不牺牲数值的准确性[7]。

Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。(Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字)。尤其是Bruun算法，把FFT看成是z?1 ，并把它分解成z?1 与z的形式。

另一个从多项式观点的快速傅立叶变换法是Winograd 算法。此算法把

9

NM2M?az?1

M武汉工程大学毕业设计(论文)说明书

zN?1 分 解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。 这

样只需要很少的乘法量(如果有需要的话)，所以winograd是可以得到最少乘法量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd算法让DFT可以用2K点的DFT来简化，但减少乘法量的同时，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。

Rader算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋折积来表示原本的DFT，如此就可利用折积用一对基本的FFT来计算 DFT。另一个prime-size的FFT算法为chirp-Z算法。此法也是将DFT用折积来表示，此法与Rader算法相比，能运用在更一般的转换 上，其转换的基础为Z转换[8]。 3.4 Goertsel 算法

如前所述, N点时域序列x(n)的离散付里叶变换式为

?kn， k?0,1,2,....,N?1 （3-16） X(k)??x(n)WNn?0N?1这N点频域序列是同时被算出的，不可能只计算其中某一个或几个指定点。Goetzel 算法是为了解决这个问题而提出的。这个算法把离散付里叶变换看作一组滤波器，将输入端的时域序列与其中一个滤波器的冲激响应序列进行卷积运

k算，求滤波器的输出序列，即得X(k)序列的一点。这种算法利用旋转因子WN的

周期性，使DFT运算化为线性滤波运算[9]。

由于

W故式（3-16）可化为 X(k)?W?kNNN?1m?0?kNN?e?j(2?)(?kN)N?1

?x(m)WkmN?k(N?m)??x(m)WNk?0,1,2,....,N?1 （3-17） m?0N?1定义序列yk(n)为

?k(n?m) yk(n)??x(m)WN （3-18）

m?0N?1可见yk(n)是由两个序列卷积而得到的序列。

yk(n)?x(n)?hk(n) （3-19） 其中，x(n)是输入的N点序列，另一个序列被看作滤波器的冲激响应序列

10

武汉工程大学毕业设计(论文)说明书

?knhk(n)?WNu(n)。 （3-20）

对比式(3-17）和式(3-18），可知：按式（3-19）进行卷积运算，当n?N时，滤波器的输出yk(N)就是X(k)：

X(k)?yk(n)|n?N （3-21） 对式(3-20)进行Z变换，可得滤波器的系统函数 Hk(z)?11?W?k?1z （3-22）

这是一个一阶系统。图3-2示出这个系统的信号流图，相应的差分方程为

?k yk(n)?WNyk(n?1)?x(n)， y(?1)?0 (3-23)

按照此式进行递推运算，到了 n?N 时刻，即可依据式（3-21）得到X(k)。

?k按照式(3-22)进行运算时，可先算好旋转因子WN，储存起来。每次递推包含一

次复数乘法。按式(3-16)直接计算N点离散付里叶变换,需要4N2次实数乘法和

N(4N-2)次实数加法。按照上述Goertzel算法，所需的实数乘法和实数加法都

是4N2次。所以当N不大时，上述算法的效率稍差[10]。下面介绍改进的Goertzel 算法，这种算法所需的实数乘法次数约为直接方法的一半。

图3-2 用一阶系统实现Goertzel 算法 图 3-3 用二阶系统实现Goert算法

k?1把式（3-22）的分子和分母都乘上因子(1?WNz)，就得到第k个滤波器的

系统函数为

k?11?WNz Hk(z)? ?k?1k?1(1?WNz)(1?WNz) 11

武汉工程大学毕业设计(论文)说明书

k?11?WNz ? （3-24）

1?2cos((2?/N)k)z?1?z?2与此相应的信号流图示于图3-3。由式（3-24)可见，滤波器是一个二阶系统，有一对复数共轭极点和一个复数零点。

为了便于运算，在图3-3所示的流图中，设立状态变量 u和v。按照图3-3计算X(k)时，步骤有二，即： 1.实现一对复数极点

输入点依次取 x(0),x(1),x(2),...,x(n?1)，进行递推运算。每次运算中，更新状态变量u和v。作N次迭代所需的计算量是2N次实数乘法和4N次实数加法。 2.实现复数零点。

x(n)是一个N点序列， n?0,1,2,...,N?1。在x(N)?0点上。计算状态变量u和v。

这时，按照图3-2算出滤波器的输出yk(N)，此即X(k)。所需的计算量是4次实数乘法和4N次实数加法[11]。

综上所述，计算一点X(k)需要进行2(N?2)次实数乘法和4(N?1)次实数加法。这种算法要求的乘法次数约为直接算法的一半。在这种较为有效的方案中，

k仍具有这样的优点，即必须计算和存储的系数只有cos((2?/N)k)和WN。

还要说明图3-3所示的算法的另一个优点。当输入序列为实序列时，离散付里叶变换序列X(k)是对称的，即X(k)?X*(N?k)。容易证明，图3-3的网络形式在计算X(N?k)时和计算X(k)时具有完全相同的极点，但前者的零点系数与后者的零点系数成复共轭关系。由于零点仅在最后的迭代中实现，所以诸极点要求的

2N次乘法和 4N 次加法可以用来计算离散付里叶变换的两个值。因此，若用Goertzel 算法计算离散付里叶变换的所有N个点的值，需要的乘法次数近似为

N2，加法次数近似为2N2。然而，它同直接计算离散付里叶变换一样，计算量

仍然正比于N2。

12

武汉工程大学毕业设计(论文)说明书

4.快速傅里叶变换在信号处理中的理论应用

4.1 连续时间信号的快速傅里叶变换

设x(t)是连续时间信号，并假设t?0时x(t)?0，则其傅里叶变换由下式给出：

X(?)??x(t)e?i?tdt （4-1）

0?令?是一固定的正实数，N是一个固定的正整数。当??k?,k?0,1,2,,N?1时，利用FFT算法可计算X(?)。

已知一个固定的时间间隔T，选择让T足够小，使得每一个T秒的间隔

nT?t?(n?1)T内，x(t)的变化很小，则式中积分可近似为

X(?)??(?n?0??(n?1)TnTe?iwtdt)x(nT)

??[n?0?1?i?tt?(n?1)Te]t?nTx(nT) i?1?e?i?T ?i??en?0??i?nT x(nT) （4-2）

假设N足够大，对于所有n?N的整数，幅值x(nT)很小，则式（4-2）变为

1?e?i?T X(?)?i??en?0N?1?i?nT x(nT) （4-3）

当??2?k/NT时，式（4-2）两边的值为

2?k1?e?i2?k/N X()?NTi2?k/NT?en?0N?1?i2?nk/N1?e?i2?k/Nx(nT)?X[k] （4-4）

i2?k/NT其中X[k]代表抽样信号x[n]?x(nT)的N点DFT。最后令??2?/NT，则上式变为

1?e?i2?k/NX[k]k?0,1,2, X(k?)?i2?k/NT,N?1 （4-5）

,?时的1首先用FFT算法求出X[k]，然后可用上式求出k?0,1,2,N

13

武汉工程大学毕业设计(论文)说明书

X(k?)。

应该强调的是，式（4-3）只是一个近似表示，计算得到的X(?)只是一个近似值。通过取更小的抽样间隔T，或者增加点数N，可以得到更精确的值。如果

??B时，幅度谱X(?)很小，对应于奈奎斯特抽样频率?s?2B，抽样间隔T选

择?/B比较合适。如果已知信号只在时间区间0?t?t1内存在，可以通过对

nT?t1时的抽样信号x[n]?x(nT)补零，使N足够大[12]。

将连续时间傅立叶变换进行数字近似，用函数fft（快速傅立叶算法）高效地计算这个近似值。很多信号都能用（4-1）式连续时间傅立叶变换（CTFT）来表示。利用MATLAB可以计算（CTFT）积分的数值近似。利用在密集的等间隔t的样本上的求和来近似这个积分，就可以用函数fft高效地计算这个近似值。所用的近似式是根据积分的定义得到的，即

????x(t)e?j?tdt?lim?x(n?)e?j?n?? （4-6）

??0n????对于一般信号，在足够小的τ下，上式右边的和式是对于CTFT积分的一个好的近似。若信号x(t)对于t?0和t?T为零，那么这个近似式就能写成

????x(t)e?j?tdt??x(t)e0T?j?tdt??x(n?)e?j?n?? （4-7）

n?0N?1式中T?n?，N为一整数。可以利用函数fft对一组离散的频率?k计算上式中的和式。如果N个样本x(n?)是存在向量x内的话，那么调用函数X=tau*fft(x)就可以计算出 X(k?1)???x(n?)en?0N?1?j?nk??X(j?k) （4-8）

N?2?k 0?k??N?2式中 ?k????2?k?2? N?1?k?N ??2?N?

以及N假设为偶数。为了计算高效，fft在负的频率样本之前先产生正频率样本。为了将频率样本置于上升的顺序，能用函数fftshift。为了将存入X中的X(j?k)的样本排列成使X(K?1)就是对于0?k?N?1，在????(2?kN?)上求得

14

武汉工程大学毕业设计(论文)说明书

的CTFT,可用X=fftshift(tau*fft(x))。

例4-1 利用FFT计算傅里叶变换 如图4-1所示的信号

?t?10?t?2 x(t)??0其它?其傅里叶变换为：

X(?)??cos(?)?sin(?)?i?2ei 2?利用下面的命令，可得到x(t)的近

似值和准确值。 图4-1 连续时间信号x(t)

N=input('Input N:'); T=input('Input T:'); %计算X(w)近似值 t=0:T:2;

x=[t-1 zeros(1,N-length(t))]; X=fft(x); gamma=2*pi/(N*T); k=0:10/gamma;

Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X; %计算真实值X(w) w=0.05:0.05:10;

Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w); plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact)); legend('近似值','真实值');

xlabel('频率(rad/s)');ylabel('|X|')

运行程序后输入N=128，T=0.1，此时??0.4909，得到实际的和近似的傅里叶变换的幅度谱如图4-2所示，此时近似值已经相当准确。通过增加NT可以增加更多的细节，减少T使得到的值更精确。再次运行程序后输入N=512，T=0.05，此时??0.2454，得到实际的和近似的傅里叶变换的幅度谱如图4-3所示。

15

武汉工程大学毕业设计(论文)说明书

图4-2 N=128，T=0.1时的幅度谱

图4-3 N=512，T=0.05时的幅度谱

16

武汉工程大学毕业设计(论文)说明书

4.2 FFT计算离散信号的线性卷积

已知两个离散时间信号x[n](n?0,1,2,M?1)与y[n](n?0,1,2,N?1)，取

L?M?N?1，对x[n]和y[n]右端补零，使得

x[n]?0,n?M?1,M?2, y[n]?0,n?N?1,N?2,,L?1

,L?1 （4-9）

利用FFT算法可以求得x[n]和y[n]的L点DFT，分别是X[k]和Y[k]，利用DTFT卷积性质，卷积x[n]*y[n]等于乘积X[k]Y[k]的L点DFT反变换，这也可以通过FFT算法得到。

例4-2 利用FFT计算线性卷积

已知x[n]?0.8nu[n]，其中u[n]为单位阶跃序列，信号y[n]如图4-4所示。由于当n?16时，x[n]很小，故M可以取为17；N取10，L?M?N?1?26。

利用下面的Matlab命令，可得到x[n]、y[n]的卷积图形如图4-4所示。 subplot(3,1,1); n=0:16; x=0.8.^n; stem(n,x);

xlabel('n');ylabel('x[n]'); subplot(3,1,2); n=0:15;

y=[ones(1,10) zeros(1,6)]; stem(n,y);

xlabel('n');ylabel('y[n]') subplot(3,1,3); L=26;n=0:L-1;

X=fft(x,L);Y=fft(y,L); Z=X.*Y;z=ifft(Z,L);

17

武汉工程大学毕业设计(论文)说明书

stem(n,z);

xlabel('n');ylabel('z[n]')

图4-4 信号x[n]、y[n]及其卷积z[n]=x[n]*y[n]

利用下面的Matlab命令，可得到信号x[n]、y[n]的幅度谱与相位谱如图4-5所示。

subplot(2,2,1); L=26;k=0:L-1;

n=0:16;x=0.8.^n;X=fft(x,L); stem(k,abs(X)); axis([0 25 0 5]);

xlabel('k');ylabel('|X[k]|') subplot(2,2,2); stem(k,angle(X)); axis([0 25 -1 1]);

xlabel('k');ylabel('Angle(X[k])(弧度)') subplot(2,2,3);

y=ones(1,10);Y=fft(y,L);

18

武汉工程大学毕业设计(论文)说明书

stem(k,abs(Y)); axis([0 25 0 10]);

xlabel('k');ylabel('|Y[k]|') subplot(2,2,4); stem(k,angle(Y)); axis([0 25 -3 3]);

xlabel('k');ylabel('Angle(Y[k])(弧度)')

图4-5 信号x[n]、y[n]的幅度谱与相位谱

4.3 FFT进行离散信号压缩

利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换[13]。

例4-3 对单位区间上的下列连续信号

1f(t)?t?cos(4?t)?sin(8?t)

2以fs?256Hz采样频率进行采样，将其离散化为28个采样值

f[n]?f(t)|t?nT?f(nT)?f(n/256),n?0,1,2,,255

用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图4-6a)所示，均方差为0.0429，相对误差

19

武汉工程大学毕业设计(论文)说明书

为0.0449；令绝对值最小的90%系数为0，得到重构信号图形如图4-6b)所示，均方差为0.0610，相对误差为0.0638。

a) 绝对值最小的80%系数为0的重构信号（FFT） b) 绝对值最小的90%系数为0的重构信号（FFT）

图4-6 用FFT压缩后的重构信号

相关Matlab程序如下: function wc=compress(w,r) %压缩函数compress.m %输入信号数据w,压缩率r %输出压缩后的信号数据 if(r<0)|(r>1)

error('r 应该介于0和1之间!'); end;

N=length(w); Nr=floor(N*r); ww=sort(abs(w)); tol=abs(ww(Nr+1)); wc=(abs(w)>=tol).*w;

function [unbiased_variance,error]=fftcomp(t,y,r)

%利用FFT做离散信号压缩

%输入时间t,原信号y,以及压缩率r if(r<0)|(r>1)

error('r 应该介于0和1之间!');

20

武汉工程大学毕业设计(论文)说明书

end; fy=fft(y);

fyc=compress(fy,r); %调用压缩函数compress.m yc=ifft(fyc);

plot(t,y,'r',t,yc,'b'); legend('原信号','重构信号');

unbiased_variance=norm(y-yc)/sqrt(length(t)); error=norm(y-yc)/norm(y); 输入以下Matlab命令： t=(0:255)/256;

f=t+cos(4*pi*t)+1/2*sin(8*pi*t);

[unbiased_variance,error]=fftcomp(t,f,0.8) unbiased_variance = 0.0429 error = 0.0449

如果用Harr尺度函数和Harr小波分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图4-7a)所示，均方差为0.0584，相对误差为0.0611；令绝对值最小的90%系数为0，得到重构信号图形如图4-7b)所示，均方差为0.1136，相对误差为0.1190。

a) 绝对值最小的80%系数为0的重构信号（Harr） b) 绝对值最小的90%系数为0的重构信号（Harr）

图4-7 用Harr小波压缩后的重构信号

21

武汉工程大学毕业设计(论文)说明书

相关Matlab程序如下

function [unbiased_variance,error]=daubcomp(t,y,n,r) %利用Daubechies系列小波做离散信号压缩 %输入时间t,原信号y,分解层数n,以及压缩率r

%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l^2误差 if(r<0)|(r>1)

error('r应该介于0和1之间!'); end;

[c,l]=wavedec(y,n,'db1');

cc=compress(c,r); %调用压缩函数compress.m yc=waverec(cc,l,'db1'); plot(t,y,'r',t,yc,'b'); legend('原信号','重构信号');

unbiased_variance=norm(y-yc)/sqrt(length(t)); error=norm(y-yc)/norm(y); 输入以下Matlab命令： =(0:255)/256;

f=t+cos(4*pi*t)+1/2*sin(8*pi*t);

[unbiased_variance,error]=daubcomp(t,f,8,0.8) unbiased_variance = 0.0584 error = 0.0611

结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好（甚至比小波还要好）。 4. 4 利用FFT对离散信号进行滤波

利用FFT算法对信号进行滤波的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值大于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行

22

武汉工程大学毕业设计(论文)说明书

逆傅里叶变换[14]。

例4-4 对被白噪声污染的信S?2?3cos?100?t??/6??1.5cos?150?t??/2?进行频谱分析，从中鉴别出有用的信号。要求：将信号的幅度换算成实际的幅度，信号的频率换算成实际的频率。 相关Matlab程序如下

fs=1000;%采样频率 N=512;%数据点数 n=0:N-1; randn(1,N);

t=0:1/fs:(N-1)/fs;%采样时间序列 f0=100;%信号频率

x=2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2); subplot(3,1,1); plot(t,x); xlabel('t');

ylabel('sin(2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2))'); title('时域信号');

Y=fft(x,N);%对信号进行FFT变换 magY=abs(Y);%求得FFT变换后的幅度

angY=angle(Y)*180/pi;%求得FFT变换后的相位 f=n*fs/N;%频率序列 subplot(3,1,2);

plot(f,magY);%画出幅频响应 xlabel('f'); ylabel('幅度'); title('N=512'); grid;

subplot(3,1,3);

plot(f,angY);%画出相频响应

23

武汉工程大学毕业设计(论文)说明书

xlabel('f'); ylabel('相位') title('N=512'); grid;

图4-8 滤波信号的相位，幅度

4.5 利用FFT提取离散信号中的最强正弦分量

这里最强是指在信号x[n],n?0,1,,N?1中某个正弦分量的振幅远大于其

它正弦分量的振幅。可以对x[n]求N点DFT来确定信号中是否有最强的正弦分量。如果信号x(t)是连续时间形式的，首先还要对其进行抽样，得到离散时间形式的信号x[n]?x(t)|t?nT?x(nT)，根据Nyquist定理，抽样间隔T应满T??/?max，其中?max是x(t)中的最大频率分量。要判断信号x[n]中是否包含最强正弦分量，采样数据至少要包含该分量一个完整周期的数据[15]。

例4-5 太阳耀斑数据的分析

太阳耀斑活动的周期是11年，这个事实可以通过提取耀斑数据的最强正弦 分量加以证实。耀斑数据可以从比利时皇家天文台(Royal Observatory of Belgium)太阳耀斑数据索引中心（Sunspot Index Data Center, SIDC）网站下载。网址是：http://sidc.oma.be/index.php3.

24

武汉工程大学毕业设计(论文)说明书

下载后的数据存放在文件“monthssn.dat”中，里面有四列数据，第一年是日期，第三列是太阳耀斑的平均数，第四列平滑后太阳耀斑的平均数，可以得到从1749年到当前年月（2006年4月）的耀斑数据。本次分析选择第三列1981年1月作为开始日期，2005年12月作为结束日期，共25年300个月份的数据。为此先把相关数据复制到Excel表格的第一列中，然后保存到Matlab所在目录下，并命名为“sunspotdata.csv”。然后输入以下命令，得到耀斑曲线如图9所示[16]。

spd=csvread('sunspotdata.csv',0,0,[0 0 299 0]); plot(spd); grid;

xlabel('月数');ylabel('耀斑平均数'); axis([0 300 0 200])

图4-9 1981年1月至2005年12月太阳耀斑的平均数

由上图可见，太阳耀斑的活动确实具有周期性，但周期的准确值不明显。可以通过数第一个峰值和第二个峰值之间的月份来估计周期的值。查验表中的数据得，第一个峰值为200.3，出现在第116个月（1990年8月），第二个峰值为170.1，出现在第235个月（2000年7月），所以周期是235-116=119月，和实际值132月比较接近。

25

武汉工程大学毕业设计(论文)说明书

下面利用DFT来分析spd[n]。首先，从图4-9中可以看出spd[n]从n?0到

n?300整个区间的平均值不为0，为了更方便地分析spd[n]，需要减去该均值，

得到结果如下：

1300 x[n]?spd[n]??spd[i] （4-12）

300i?1然后对x[n]进行傅里叶变换，得到X[k]。在Matlab中输入以下命令，得到

x[n]的幅度谱X[k]的图形如图10所示。

x=spd-sum(spd)/300;X=fft(x); k=0:299;

plot(k,abs(X));grid;

xlabel('k');ylabel('|X[k]|')

图4-10 x[n]的幅度谱

由图4-10可见，x[n]的幅度谱中有一个清晰的尖峰，这就证实x[n]中确实包含一个最强正弦分量。为了确定尖峰所对应的频率，用火柴棒图重画当k?0,

26

武汉工程大学毕业设计(论文)说明书

1,,10时X[k]的图形，在Matlab中输入以下命令，得到图形如图4-11所

示。

k=0:10;X=X(1:11); stem(k,abs(X));

xlabel('k');ylabel('|X[k]|')

图4-11 当k?0,1,,10时X[k]的火柴棒图

可以看出上图中有两个强谱分量k?2,3，频率分别为2?k/N?2??2/300弧度/月和2?k/N?2??3/300，周期分别为150月和100月。由于spd[n]的数据长度不是其中强正弦分量的整数倍，谱分量在k?3出现了泄露。要消除泄露，需要使数据的长度正好是其中强正弦分量的周期（太阳耀斑活动的周期，也即

11?12?132）的整数倍。为此，重新取spd[n]中从1至264之间的数据进行分析，

令

[n],?n1,2, y[n]?spd ,21264 z[n]?y[n]?y[i] （4-13）?264i?1然后对z[n]进行傅里叶变换，得到Z[k]。输入以下命令，得到k?0,1,,10时

Z[k]的火柴棒图如图4-12所示。

y=spd(1:264);z=y-sum(y)/264;Z=fft(z);

k=0:10;Z=Z(1:11);stem(k,abs(Z)); xlabel('k');ylabel('|Z[k]|')

27

武汉工程大学毕业设计(论文)说明书

图4-12 当k?0,1,,10时Z[k]的火柴棒图

注意到Z[k]的峰值仍然在k?2，但此时Z[2]比其附近的值大得多，这说明没有出现泄露，频率为??2?k/N?2??2/264弧度/月的正弦分量就是y[n]或者z[n]中唯一的最强正弦分量。该频率对应周期2?/??132月，正好等于11年。

k?2时Z[k]的值主要是由于太阳耀斑数据中的噪声引起的。然而和其它值相比，Z[1]、Z[3]和Z[4]的值相对较大，这说明除了周期为11年的最强分量，还包括其它正弦分量，此时太阳耀斑的活动并不是一个只包含单一频率4?/264的纯正弦运动。

28

武汉工程大学毕业设计(论文)说明书

5. 快速傅里叶变换在数字信号分析与处理的实际应用

5.1 快速傅里叶变换在喇曼光谱信号噪声平滑中的应用

电探测系统是光信号的转换、传输及处理的系统. 系统的各个部分在工作时总会受到一些无用信号的干扰,给光谱峰的检测判别及进一步的数据处理带来了不利因素.对光谱信号进行数字滤波,以获得更真实的光谱信息，显得格外重要. 目前最为通用和有效的信号滤波处理方法是快速傅立叶变换方法.纯水是一种较弱的喇曼散射介质,需要专用的喇曼散射光谱仪器才能获得高信噪比的喇曼光谱.我们以增强型的CCD 探头为探测器,结合普通的分光单色仪,在YAG 激光器532 nm 激光线的激励下获得低信噪比的纯水的喇曼光谱. 信噪比较差的喇曼光谱经过FFT 变换后,,用FFT 的逆变换将滤除噪声后的频谱信号转换成为光谱信号,最终获得信噪比较高的纯水的喇曼光谱[17]。

实验原理及结果如下：傅里叶变换的基本表达式为

nkx(k)??x(n)WN,(k?0,1...N?1 )n?0N?1 （5-1）

nkWN?exp[?2?kn]N （5-2）

式（5-1）中的x（n）(n=0,2,?N-1)是列长为N的输入序列，即实验采集到的时域上的切片数据;x(k)(=0,1,?N-1)是列长为N的输出序列，即经过傅里叶变换后的频域上的数据。

对数字化后的光谱信号而言, x(n)是一组离散的实数信号；而X(k)分为实部x(v)和虚部y(ν)2部分。x(ν)和y(ν)又可组成振幅谱A(ν)和相位谱P(ν)：

22A(V)?x(v)?y(v) （5-3）

p(v)?arctan

y(v)x(v) （5-4）

通过对式(5-3)和式(5-4)性能的考察,发现A(ν)和P(ν)中既含有目标信号的信息,也含有噪声的信息,如果二者所在的区域不同, 则可以通过傅里叶变换分析出噪声信息, 将之从捕获的信号中去除,从而达到噪声平滑的目的, 获得高信噪比的目标信号.

纯水普通喇曼散射的信号很弱,我们在532nm 脉冲激光泵浦液滴的条件下获

29

武汉工程大学毕业设计(论文)说明书

得其散射光谱.由于样品信号极其微弱,在将CCD 的增益调至最大时,获得如图1 所示的纯水的喇曼光谱. 光谱的信噪比值用如下方式估算:

设x(???n) 为含噪声图像y(???n)为消除噪声后的图像,图像的均方根误差为

???x(??n)-y(??n))n?1N2 （5-5） N信噪比定义为除噪声后的信号与均方根误差之比

N[????]/N SNR??Y(N?1? 计算出642. 86 ～ 643. 62 nm 光谱区的信噪比为

SN R ≈ 17.

图5-1.多通道光谱分析仪采集的含有噪声的纯水普通喇曼散射信号

30

（5-6）

武汉工程大学毕业设计(论文)说明书

图5-2 .傅里叶变换后的频谱图

对图5-2幅度谱纵轴取对数得图噪声幅度门限值低于2 ×105 ,经门限滤波处理,在频谱图中将幅度谱低于该门限值的频率成分去除,获得的频谱用FFT 的逆变换返回得到门限滤波曲线如图5-2所示.计算出642. 86～643. 62 nm 光谱区的信噪比为SN R≈484. 与图5-1相比,光谱的信噪比有了极大的改善.

从本实验可以得出：在光谱信号受到光子噪声调制的条件下,如果光谱信号的变化频率低于高频光子噪声的变化频率,则可以通过快速傅里叶变换,获得目标信号和噪声信号的频谱,进行低通滤波和门限滤波后,分别将具有高频和不同振幅的噪声信号去除,实现对弱光谱信号干扰噪声的抑制,得到高信噪比的光谱信号。快速傅里叶变换在效果上，减轻了噪声的干扰，同时计算也不会带来过于复杂的计算。

5.2.采用异步实现的快速傅里叶变换处理器

快速傅里叶变换(FFT)是数字信号处理领域一个重要的分析工具,广泛应用于雷达、通讯、图像处理、声纳和生物医学领域。已经开发出多种专用快速傅里叶变换处理器,大大提高了快速傅里叶变换的运算速度。异步集成电路具有功率效率高、电磁兼容性(EMC)好、功耗低和没有时钟歪斜(Skew)的特性,同时又具有潜在的高性能,以及便于系统模块化设计的优势[1]。异步集成电路运算的性能是平均性能,而不是最差性能。这样,当平均性能与最差性能差别较大时,异步集成电路有希望达到比同步电路更高的潜在性能。异步集成电路采用大量本地时序控

31

武汉工程大学毕业设计(论文)说明书

制信号来取代整体时钟,避免了当前在超大规模集成电路设计中遇到的时钟树设计和代价问题。

本实验原理及结果：异步实现的快速傅里叶变换处理器的结构如图5-3所示。处理器的控制由本地的握手信号控制,每个单元独立地工作,避免了同步电路中的时钟分配问题。处理器在输入数据准备好后开始工作,整个运算完成时产生一个完成信号。采用0.6μm 标准CMOS工艺,设计一个8点的异步快速傅里叶变换处理器。该处理器具有2×8 比特的输入,2×15 比特的输出,2×20 比特的内部运算精度。在电路设计完成之后,采用华晶2上华的0.6μm CMOS2P2M 混合电路工艺,建立了异步标准单元库,然后对异步快速傅里叶变换处理器进行了全定制设计。处理器的版图如图5-3所示。

图5-3 异步快速傅里叶变换处理器结构

功能仿真：用晶体管构成的电路网表描述每个单元(加法器、乘法器等) ,然后用Hspice 进行功能仿真。根据电路Hspice 仿真结果,通过抽象模型,建立每个单元的功能和延迟的逻辑模型。异步逻辑和运算模块的抽象过程比同步模块要复杂得多,因为同步模块只要用功能加上一个最差延迟就可以描述模块的功能性能模型。CMOS 的抽象过程就是用逻辑描述建立FFT的逻辑网表(带延迟) ,再用Verilog 进行逻辑仿真。

性能仿真：响应时间是异步集成电路性能分析时常用的度量标准[5]。响应时间是指请求信号到完成信号之间的延迟,它主要有两种类型:最差响应时间和平均响应时间。其中, 最差响应时间主要依赖于电路的结构和实现,而平均响应时间不仅与电路结构有关,还与输入的数据相关。文中采用Star2SimXT ,对整个

32

武汉工程大学毕业设计(论文)说明书

异步快速傅里叶变换处理器进行了电路仿真,得到芯片完成一次变换的最差响应时间为42.85 ns ,平均响应时间为31.15 ns ,功耗约为350.7mW。

从本实验可以得出：设计了一个异步的快速傅里叶变换处理器,该电路可以在异步逻辑控制下工作。性能分析表明,异步快速傅里叶变换处理器的平均性能较同步设计有优势。但是,异步集成电路完成信号的产生往往需要增加一部分电路。这不仅增加了芯片的面积,而且带来了一定的延迟,异步集成电路性能的优势能否实现,与这部分电路设计是否合理有很大的关系。另外,由于缺少成熟的EDA工具、算法和设计方法学的支撑,异步集成电路设计技术在超大规模集成方面还面临很多挑战，还需继续改进。

5.3 快速傅里叶算法在哈特曼夏克传感器波前重构算法中的应用

哈特曼夏克传感器因其波前测量实时性好等特点而广泛用于自适应光学系统中,随着应用研究的发展,哈特曼夏克波前测量传感器的空间分辨率也要相应提高。哈特曼夏克传感器测量的是波前相位斜率,需要经过波前复原求出相位值,复原的方法主要有区域法和模式法两类,为了满足实时性的要求,哈特曼夏克传感器的子孔径较少,测量的空间分辩率因此比干涉仪低。当增加哈特曼夏克传感器的子孔径数提高空间分辨率、提高测量精度时,区域法和模式法的运算量非常大,实时性降低,限制了高分辨率哈特曼夏克传感器在自适应光学系统等领域的进一步应用。针对实时性问题,提出了分块算法和迭代法进行波前重构。在区域法重构波前的基础上,应用快速傅里叶变换(FFT) 算法,提高波前复原算法的实时性,为高分辨率哈特曼夏克传感器在自适应光学技术及其它领域的应用作算法准备。

本实验原理及结果：快速傅里叶变换算法以其运算速度快、所需内存小而被广泛用于数字信号处理领域[9]。在求解由(1)式确定的线性方程组的过程中,需要实现方程系数矩阵的对角化,而这一过程可以通过快速傅里叶变换算法实现,从而实现(1) 式的快速求解。首先,不考虑区域中边界处的相位估计差分方程,在波前重构的区域内,即1≤i≤M -1,1≤j≤N -1,(1)式严格成立,并由它导出波前估计的矩阵方程组表示为

33

武汉工程大学毕业设计(论文)说明书

a0?1??2??1??q?1?A0?j??q?1??,(2?q?M?2)??M?2?A0?M?1??M?1,?4 -1 ??-1 4 -1 ?0??? ?-1 4??A0?? ?? 0 -1??? -1 4???? ??

(5-7)

(5-8)

对(5-8)式的矩阵AO作正交变换,得:

??````??q?D???????1jq?1?(2?q?M?2),?```???M?3?D?M?2??M?2?其中

``?q?Q`?q , ? Q`?qq?`D?1`??2??1`,

(5-9)

（5-10）

332(MN?NM) 应用快速傅里叶变换算法,乘法运算量可由直接作线性变换

次降为2(MNlog22N?MNlog22M)次,当哈特曼夏克传感器的子孔径数比较大时,运算速度可大幅度提高,从而提高哈特曼夏克传感器波前重构算法的实时性。在波前估计的计算式(5-8)中,只考虑了哈特曼夏克传感器区域内的估计点,需要知道区域边界处的相位值,才能准确求解（5-8）式,而哈特曼夏克传感器测量的是斜率值,给出的是诺依曼边界条件,需要作边界条件的近似求解,求得边界处的相位值。在边界上:

?0,j?1??0j??0j?gx,1j???M,j?1??MJ??MJ?gx,M,j,???i?1,0??i?1,0??j0?gy,i,1,??i?1,N??i,N??iN?gy,i,N??

(5-11)

34

武汉工程大学毕业设计(论文)说明书

由于实际被测的波前相位是连续光滑的曲面,则在边界上的相位点是封闭连续的曲线,设:

?0,0?0

则边界上的相位最小二乘解的矩阵表达式为 A??G (5-12)

其中,A为(5-8)式中AO的形式,对角线元素为2,维数[2(M +N)-1]×[2(M +N)-1],

从本实验可以得出：本文在应用区域法对波前进行最小二乘估计的过程中,应用快速傅里叶变换算法,在子孔径数较多的哈特曼夏克传感器的波前重构过程中,使算法的运算量大幅度降低,既节约处理系统的内存,又提高了波前重构的实时性,为解决高分辨率哈特曼夏克传感器实时性上的问题,在算法上提出了一种解决途径。从而可以在不降低哈特曼夏克传感器实时性、稳定性的前提下,进一步提高哈特曼夏克传感器的空间分辨率,提高测量精度。

35

武汉工程大学毕业设计(论文)说明书

致谢

此文得以完成，凝聚了许许多多老师、同学、朋友，亲人的心血和关爱！在我即将完成学业之际，谨向五年来给与我无私帮助、支持，关心和呵护过我的所有老师、同学、朋友、亲人致以最诚挚的谢意！

感谢华夏老师作为我的论文指导老师在本文的撰写过程中给予我大量的指导和帮助，花费了很多心血.尤其是在课题设计、研究方法、论文撰写等各个环节给予我的指导和帮助。

感谢我的同窗好友们，四年来朝夕共处的日子里，是他们给了我最大的温暖和感动，感谢他们在我论文写作过程中提出的宝贵建议和帮助。论文写作过程中借鉴和引用了许多学界前辈的观点和论据，向他们表示感谢！最后，特别感谢参加论文评审的各位老师！

36

武汉工程大学毕业设计(论文)说明书

参考文献

[1] 程佩青. 数字信号处理教程. 第2版[M]. 北京:清华大学出版社, 2001. [2] 张易知. 虚拟仪器的设计与实现[M]. 西安:西安电子科技大学出版社, 2002.

[3] 蒋正萍. 数字信号处理[M]. 北京:电子工业出版社, 2004.

[4] E. O. 布赖姆. 快速傅里叶变换[M ]. 柳 群译. 上海科学技术出版社, 1979.

[5] 付丽琴,桂志国,王黎明. 数字信号处理原理及实现[M] . 北京:国防工业出版社,2004.

[6] 叶卫平,方安平. 科技绘图及数据分析[M] . 北京:机械工业出版社,2004. [7] 蒋长锦 蒋勇. 快速傅里叶变换及 c 程序 [M].中国科技大学出版社.2004. [8] 王永仲.迂回相位编码的傅里叶变换计算全息图及其再现.红外技术. 2004. [9] 谭康. 泰勒公式家泰勒级数之妙用[J]. 高等数学研究. 2010.3. [10] 陆旦前.FFT算法的一种FPGA设计.现代电子技术. 2007. [11] 程乾生.《数字信号处理》.北京大学出版社.2003.

[12] Degrieck J, Verleysen P, Waele W. Optical measurement of target displacement and velocity in bird strikesimulation experiments[J]. Measurement Science and Technology, 2003, 14: 1-6.

[13] Mccarthy B D, Regan B J. Position measuring apparatus and method[P]. USA, US4885725,1989.12.

[14] Duhamel P, Hotmann H. Split-radix FFT algorithm. Electronic letters, 20(1): 14-16, 1984.

[15] Winograd S. On computing the discrete Fourier transform. Proc. Nat. Acad. Sci. USA, 73(Apr): 1005-1006, 1976.

[16] Duhamel P, Vetterli M. Fast Fourier transform: a tutorial review and a state of the art. Signal processing, 19: 259-299, 1990.

[17] Weihua Zheng, Kenli Li, and Keqin Li. A Fast Algorithm Based on SRFFT for Length N = q × 2m DFTs,NO. 2, FEBRUARY 2014.

37

本文来源：https://www.bwwdw.com/article/42mg.html