第六章 数字信号分析(I)-DFT与FFT

第六章数字信号分析（Ⅰ）——DFT与FFT●

计算机技术的发展是数字信号分析的基础。●

数字滤波优于模拟滤波，（1）速度快，例如，采用数字信号分析技术对于1024采样点进行A/D转换，仅需4～15μs，进行FFT运算须250ms，较快的只需数毫秒；一个蝶形FFT硬件运算只需2μs。（2）分辨力高，在高频段（50kHz)可达25Hz；在超低频段可达0.0025Hz。●数字信号分析一般包括：频谱分析与数字滤波等主要内容。前者又包含相关与统计分析等。第一节模拟信号离散化●

本章内容：主要介绍离散Fourier变换与快速Fourier变换的基本原理及应用。一、A/D与D/A转换1、A/D转换过程（1）采样（或称抽样）连续时间信号x(t)→离散采样信号x(n△t)。

量化误差呈等概率均匀分布，其概率密度函数p(R)=1/R，最大量化误差为±0.5R，其均方差为（2）关于量化误差问题或量化增量R取决于A/D转换器位数。例如，8位A/D，R=Vref/28。（3）编码——将离散幅值经过量化以后变为二进制数字D。信号幅值2、D/A转换过程译码——把数字信号RD恢复为有限幅值A的过程，即

D/A转换过程包括：译码与波形复原。波形复原——把离散幅值恢复为连续波形的过程，由保持电路实现。例如，零阶保持与一阶多角保持等。前者是在两个采样值之间，令输出保持上一个采样值的值；后者是在两个采样值之间，令输出为两个采样值的线性插值。经过保持变换构成的信号存在不连续点，用模拟低通滤波器消除输出波形的不连续点。

采样过程——是通过采样脉冲序列与连续时间信号x(t)相乘来完成的。根据采样脉冲序列的形状，分为理想脉冲采样与矩形脉冲采样。二、采样信号的Fourier变换1、时域采样（1）理想脉冲采样采样脉冲序列：采样信号：采样信号频谱沿频率轴每隔一个采样频率ωs，重复出现一次，即频谱产生周期延拓。幅值被Cn所加权，故频谱形状不变。（2）矩形脉冲采样周期矩形脉冲序列的傅立叶变换：故有：可见：Xs(ω)是X(ω)在以ωs为周期的重复过程中，其幅值按sinc(nωsτ/2)规律变化的函数。2、频域采样采样脉冲序列采样间隔则有可见：若频谱X(ω)被间隔为ω1的脉冲序列在频域中采样，则在时域中等效于x(t)以T1(=2π/ω1)为周期而重复。就是说，周期信号的频谱是离散的。根据因此，

上述分析证明：信号的时域与频域呈采样（离散）与重复（周期)关系。采样定理：ωs≥2ωm或fs≥2fm。因为时域采样间隔决定于fs，故又称为时域采样定理。三、采样定理1、频混现象F[xs(t)]为周期谱，其周期ωs=2π/Tsωs>2ωm,周期谱图相互分离ωs<2ωm,周期谱图相互重叠2、采样频率3、信号复原为了从Xs(ω)中无失真地选出X(ω)，用频域矩形窗函数H(ω)与Xs(ω)相乘，得

实现方法：将采样信号xs(t)通过传递函数为H(ω)的理想低通滤波器，则在滤波器输出端可得到频谱为X(ω)的连续信号x(t)。理想滤波器的H(ω)为：所以有：●连续信号可以展成正交采样函数（sinc(t)型函数）的无穷级数，级数的系数等于采样值x(nTs)。也可以说，若在采样信号xs(t)的每个采样值上画一个峰值为x(nTs)的sinc(t)型函数波形，则合成波形就是x(t)。所以，若xs(t)通过理想低通滤波器时，每个采样值产生一个脉冲响应，这些响应进行叠加就得到x(t)，从而达到由xs(t)恢复x(t)的目的。●

为什么称sinc(t)为内插函数？所谓内插是指从已知离散点的值，求在离散点之间另外一些点的值时，所采用的数学插值法。理论上讲，对不在取样时刻任意点的数值应该是无限加权样值的总和，但由于这里的内插函数是衰减的，因此，实际上可由该点附近的一组有限值之和而得到良好的逼近。●当ωs=2ωm，ωc=ωm，各个采样的冲激响应零点恰好落在采样时刻上。就采样点的数值而言，在这种情况下，各个冲激响应互相不产生“串扰”。●上述利用滤波器由采样信号恢复原信号的方法，又称为惠特克波形重构法或理想内插法。若x(t)是时域有限信号，并集中在-tm～tm的时间范围内。若在频域中以不大于1/2tm的频率间隔对X(ω)进行采样，则采样信号X1(ω)可唯一地表示原信号。●栅栏效应——频域采样后，只能获得采样点的频率成分，其余的频率成分一概被舍去。这犹如透过栅栏观赏光景，只能看到一部分，就可能使一部分有用的频率成分被漏掉，而丢掉了部分有用信息，此现象称为栅栏效应。4、频域采样定理第二节离散Fourier变换（DiscreteFourierTransform,DFT)●

DFT并非泛指对任意离散信号取傅里叶积分或傅里叶级数，而是专指适用于计算机计算的FT。这是因为：●

DFT表达式的导出方法：

（1）从离散时间序列的Z变换基础上导出，即有限长序列的离散傅立叶变换可解释为它的Z变换在单位圆上的采样；原因：对x(t)进行FT或IFT运算时，无论在时域或在频域，都需进行包括(-∞,∞)区间的积分运算。而要实现DFT，则必须做到：1）把连续信号（时域或频域）改造为离散数据；2）把计算范围收缩到一个有限区间；3）实现正、逆傅立叶变换运算。

（2）把DFT作为连续信号傅立叶变换的一种特殊情况来导出。其物理概念较前者更清楚。一、离散Fourier变换关系式1、时域采样2、时域截断对连续时间信号进行DFT,一般包括时域采样,时域截断,频域采样三个步骤.

皱纹问题：因矩形函数有突变阶跃点，在时域截断后，反映在频域上产生皱纹。此即由于Gibbs现象产生的能量泄漏效应。

★采样信号经过截断处理后，虽然在时域为有限长的离散样本，但频域内仍为连续函数，若实现逆变换，还必须改造频域函数为有限离散值。3、频域采样采样脉冲序列：令频域采样脉冲序列为δ1(f)，1）根据频域采样定理（f0≤1/2tm)，选取采样间隔f0=1/T0(T0为时域截断信号分布区间，它相当于2tm)。2）又根据FT的对称性，对应的时域函数为：采样信号：被δ1(f)采样后的频域采样信号为：逆傅立叶变换表明：是周期为T0的离散函数，每个周期内有N个离散点。δ函数的卷积)特性积分是在一个周期T0内进行，即r=0δ函数的筛选(积分)特性由于是周期函数，所以其傅立叶变换也是等间隔脉冲序列，即T0=NTs这表明：是

的傅立叶变换，在频域内（k区间(-∞,∞))是一个被Ck所加权的周期性离散脉冲序列，每一个周期内有k=N个样本点。

Ck表示频谱中的一个周期内N个采样点的复数值，以下用X(kf0)来表示同理可证明X(kf0)的逆傅立叶变换：通常将上两式写成：

这就是DFT，它通过连续FT，将N个时域采样点与N个频率采样点联系起来，是连续FT的一种特殊情况.

●此式物理意义是十分明确的，在变换过程中，都仅仅涉及到了处理区间(0≤n≤N-1;0≤k≤N-1)的N个x(n)和N个X(k)值。实质上是表明了两个N维矢量的相互线性变换（映射）。其中，x(n)可以分解为N个谐波复指数序列，每个谐波分量的频率为f0或kf0，复振幅为(1/N)X(k)。同样，X(k)被分解为N个复指数序列之和，每个分量的复振幅是x(n)，频率是nf0。●实际信号x(n)在非处理区间，可能是零(时域有限信号)，可能是周期的(周期信号)，可能是非零、非周期的(随机信号)，甚至是不定的，但只要x(n)在所关心的处理区间是确定的，就可沿用上述DFT与IDFT关系。二、DFT与FT的关系DFT与FT之间是一个近似，因为DFT需要采样与截断。而近似的程度是被分析波形的函数。1、频域有限的周期信号，时域截断长度等于周期时●时域截断波形的FT为X(f)*δ0(f)*U(f)，与X(f)相比，产生严重畸变。●当这一畸变图形被δ1(f)采样后[图(g)]，畸变却避免了。这是因为频域采样脉冲的间隔为1/T0，在这些采样频率坐标点1/T0,2/T0，…处，图(e)中的实线除了在±1/T0点有数值外，其余点都是零。

●●x(t)经采样、截断后的幅值关系设x(t)的时域幅值为A，★直接作FT变换以后频域幅值为A/2。x(t)经时域采样、截断再卷积以后([x(t)δ0(t)u(t)]*δ1(t))其幅值为AT0，不是原来的A。x(t)在频域内的X(f)经过卷积及频域采样后([X(f)*δ0(f)*U(f)]δ1(f))其幅值为AT0/2Ts，也不是原来X(f)的A/2。因此，如果希望用DFT来计算FT，就必须对DFT的变换偶对间的常数因子作一些调整，即将离散时间函数乘以因子Ts/T0，这样就得到频率函数所要求的A/2的常数因子，所以:●以上分析表明，对频域有限的周期信号，当截断长度等于其周期（或周期的整数倍）时，DFT与FT之间的差别仅仅是一个比例因子T0/Ts。因此，欲使DFT与FT等价，就必须：（1）时间函数x(t)是周期性的；（2）x(t)是频域有限信号；（3）采样频率至少是x(t)的上限频率的2倍；（4）截断函数u(t)必须正好在x(t)的一个周期内(或整数倍周期)是非零的.★更正2、频域有限的周期信号，时域截断长度不等于周期时

●时域内：信号周期为T’0，截断长度为T0，且T’0≠T0。频域内：X(f)*δ0(f)*U(f)，注意这个卷积图形是频率间隔为1/T’0的脉冲的卷积，而sinc(t)型函数的零点是1/T0,2/T0，…。卷积以后的叠加图形在1/T0,2/T0，…等点处并不是零。这表明：由于截断后的波形不是周期的整数倍，故其积分平均值不为零。当这一波形被频率间隔为1/T0的脉冲作频率采样以后，零点不再恰好与每个样本点（脉冲采样点）相重合，故引起DFT与FT之间的差异。产生这一现象的原因解释：从时域看，不按周期的整数倍对x(t)进行采样和截断，会产生具有间断点的周期函数(图中(g))，这些剧烈变化，将在频域中产生附加的频率成分。从频域看，时域截断等效于sinc(t)型函数与X(f)的单个脉冲的卷积，结果频率函数不再是一个脉冲，而是频率的连续函数，在原来脉冲的位置上，这个连续函数具有局部最大值，还有一系列称为旁瓣的峰值。这些旁瓣在频域抽样后就造成附加的频率成分，这就是所谓泄漏效应，因为时域中截断是必须的，所以泄漏效应是DFT所固有的。3、时域有限而频域无限信号

●这类信号在时域采样后必然会产生频混现象，因此采样时间间隔的选择必须使频混现象减小到允许的限度以内。

●如果对有限长的波形进行采样，选择的采样点数N正好等于它的样本点总数，则不必在时域进行截断。截断被省略，经时域采样后的函数作傅立叶变换以后[图(c)]，与频域采样脉冲δ1(f)相乘，这个乘积在时域等效于图(c)和(d)所示的时间函数的卷积。最后所产生的波形是周期的，周期由原函数的N个样本确定，所以它是原信号x(t)的复制品。这个周期函数的傅立叶变换即为图(e)所示的频域采样后的波形。

●对于这类函数，如果选择N等于时域有限信号的样本点数，则误差仅由频混造成。如果选择采样间隔足够小，就可以减少由频混引起的误差。因此，在这种情况下DFT与FT可很好的一致起来。(a)(b)(c)(d)(e)4、一般周期信号

●对于时域无限频域无限的一般周期信号，采样后存在着频混效应，如果选取截断函数精确地等于周期的整数倍，那么将不产生泄漏效应。这种情况下变换的误差源主要是频混效应，若时域截断不等于周期，则产生泄漏效应.

●实际情况中所遇到的信号，往往既不是有限时间，也不是有限带宽或周期性的信号。这类任意信号在作离散傅立叶变换中，存在着频混和泄漏效应。减小采样间隔可以减小频混效应；改善截断函数，即选择合适的窗函数，可以减小或抑制泄漏现象。5、任意信号★可见，如果处理恰当，在许多应用中都可以用DFT来得到本质上和FT等价的结果。其中，值得记住的一个重要概念是，DFT意味着在时域、频域两方面都周期化，而时域函数的N个样本点则表示周期化以后所形成的新周期函数的一个周期。★最后应指出，DFT实际上是建立了函数x(t)的N个时间样本点与N个频率样本点之间的互换关系，利用这一关系可从x(n)计算X(k)，也可以从X(k)计算x(n)。这种相互关系，在数据处理的数字量分析法中经常遇到，例如，从数据的相关函数计算功率谱密度函数，也可以从功率谱密度函数计算相关函数。三、DFT的性质（与连续傅立叶变换的性质类同）1、线性2、时移特性有限长序列x(n)的移位x(n)x(n)有限长序列的圆移位(或循环移位)

●时移特性——若将x(t)沿时间轴位移t0，则其FT要乘以因子e-j2πft0，即与此类同，如果DFT[x(n)]=X(k)，则3、频移特性DFT：若DFT[x(n)]=X(k)

则DFT[x(n)ej2πln/N]=X(k-l)

或

DFT[x(n)W-ln]=X(k-l)

FT：若F[x(t)]=X(f)

则F[x(t)ej2πf0t]=X[(f-f0)]

时间函数x(n)乘以指数项ej2πln/N，则DFT向右圆移l单位。4、离散卷积

依据运算方式不同，可分为线卷积与圆卷积。（1）时域线卷积或已知x(n)和h(n)两有限序列，卷积如图示x(n)1234h(n)12344812163691224681234302011420114y(n)nnn（2）时域圆卷积用圆卷积方法计算上例，公式如下：y(n)=x(n)

*

h(n)圆卷积符号

可见，圆卷积与线卷积所得结果不同。这是因为线卷积过程中，经反折再向右平移，在左端将依次留出空位。而圆卷积过程中，经反折的圆移序列向右移去的样值又从左端循环出现，这使得两种情况下相乘叠加而得之数值不同。01234

为解决圆卷积与线卷积结果不同的问题，将x(n)和h(n)都适当地补一些零值，以扩展其长度。那么在作圆卷积时，向右移去的零值，从左端出现仍取零值，这样就与线卷积的情况相同.补零扩展的条件为：

x(n)的样本点数h(n)的样本点数432121000438300000432110004324000000432132100041262000043210432100012620004321004321000830004321432100016941000432100043210004000x(m)x(m)x(m)x(m)x(m)x(m)x(m)h(0-m)h(1-m)h(2-m)h(3-m)h(4-m)h(5-m)h(6-m)411411202030●若选取的L不够长，圆卷积将首尾交叠混淆，其结果与线卷积不一致（这可看做是一种混叠现象）。取L≥N1+N2-1，可避免这一现象。●

圆卷积可利用快速傅立叶变换技术，实现快速卷积。因此，对于有限长序列求线卷积的问题，可转化为圆卷积来求解，以便利用快速傅立叶变换技术。（3）离散时域卷积定理运用这一定理，可对两个时域周期序列x(n)与h(n)分别计算离散傅立叶变换，再将结果相乘，然后计算乘积的离散傅立叶逆变换，即可得两个时域周期序列的卷积。这一定理为用快速傅立叶变换计算时域卷积提供了依据。（4）离散频域卷积定理即两个周期为N的时域周期采样函数，它们的乘积的离散傅立叶变换等于它们的离散傅立叶变换的卷积。5、离散相关定理变换对称为离散相关定理，即两个周期为N的时域周期序列，它们的时域离散相关的离散傅立叶变换等于它们的离散傅立叶变换的乘积。运用这个定理，可以等效地在频域中确定相关性。连续函数时域卷积定理：离散时域卷积定理：6、巴什瓦定理对于离散信号，时域功率和频域功率之间的关系，由下式给出：若x(n)为实序列第三节快速Fourier变换（FastFourierTransform,FFT)

●

FFT是一种减少DFT计算时间的算法。在FFT出现之前，虽然DFT为离散信号的分析从理论上提供了变换工具，但是很难实现，因为计算时间很长。例如，对采样点N=1000，DFT算法运算量约需200万次，而FFT仅约需1.5万次。

●

FFT方法于1965年由美国库利-图基(J.W.Cooley-J.W.Tukey)首先提出，曾被认为是信号分析技术的划时代的进步。一、FFT算法的基本原理1、DFT的计算量N×N对称方阵，Wnk=[Wnk]TN×N对称方阵，W-nk=[W-nk]T计算量：每计算一个X(k)值，需进行N次复数相乘和(N-1)次复数相加。当计算X(0)，X(1)，…共N个X(k)值时，需要N2次复数相乘，N(N-1)次复数相加。矩阵[W]与[x(n)]相乘过程中存在不必要的重复运算，这是简化运算的关键.以N=4为例讨论如下：分析：（1）不必要的计算①W0=1，②

WN/2=[e-j2π/N]N/2=-1①

Wnk的周期性，即Wnk=Wn(k+N)=Wk(n+N)，当N=4，有W2=W6，W1=W9等（2）可利用的特性②

Wnk的对称性，即W(nk+N/2)=-Wnk，当N=4，有W3=-W1，W2=-W0等原计算式周期性简化对称性简化这就是库利-图基FFT算法的基本思想。

2、减小运算工作量的途径

FFT算法有多种变型，其算法很多，但每种变型的建立，多是考虑了被分析数据的特性，或者利用计算机特性，或者利用专用计算机FFT硬件特性等。本节以基2FFT算法作为讨论的起点，因为它包含了FFT算法的基本要素，运算过程比较单纯，适于人们学习。3、FFT计算方法

●基2算法设一个点序列x(n)，要求采样点数N=2M，M为正整数。

●基2算法的出发点把N点DFT运算分解为两组N/2点的DFT运算，即把x(n)按n为偶数和n为奇数分解为两部分。将x(n)的DFT计算分为奇、偶两部分以2r表示偶数n，2r+1表示奇数n，r=0,1,2,…,(N/2-1)一个N点的DFT被分解为两个N/2点的DFT

●必须注意

G(k)和H(k)只有N/2个点，k=0,1,2,…,N/2-1。而X(k)却需要N个点，k=0,1,2,…,N-1。如果以G(k)和H(k)表达全部X(k)，应利用G(k)和H(k)的两个重复周期，由周期性可知●加权系数WkN为：●计算X(k)的全部关系式：G(k)和H(k)可分别看成是序列x(2r)与x(2r+1)的N/2点DFT

此式表明，一个N点的DFT可分解成两个N/2点的DFT

●采用蝶形流程图，以N=4为例，说明上式的计算过程。由G(0),G(1),H(0),H(1)计算X(0),X(1),X(2),X(3)蝴蝶结

一个蝴蝶结包括两次复数乘法，两次复数加法，但其中有重复。H(0)与W04相乘以及与-W04相乘可以改成只与W04相乘，再分别加减，这样就使运算量减少至只有一次复数乘法和两次复数加(减)法。同理，用第二个蝴蝶结计算X(1)和X(3)，也只有一次乘法，两次加(减)法。这样，由G(k)和H(k)计算X(k)的过程中，包含N/2个蝴蝶结运算，共N/2次复数乘法，和N次复数加(减)法。●由x(0),x(1),x(2),x(3)计算G(0),G(1),H(0),H(1)时的计算量矩阵表示W12=e-j2π/2=-1=-W02实现N=4的FFT计算全过程，所需乘法运算次数为2×N/2(4次)，加法运算次数为2×N(8次)。显然，这比用DFT直接计算所需(乘法N2=16，加法N(N-1)=12）运算工作量大为减小。2次乘法，4次加减法●N=23=8的FFT计算量分三级蝶形运算，每级需乘法N/2次，加法N次。全部运算3×N/2=12次乘法，3×N=24次加法。而直接DFT运算量N2=64次乘，N(N-1)=56次加。●N=2M的FFT计算量运算分解为M级蝶形图，每级包含N/2次乘，N次加，故算法的工作量为复数乘法复数加法●直接DFT与FFT算法所需乘法次数的比较。当N=256，二者比值为64；当N=1024，其比值为204.8二、FFT算法的应用

FFT是实施DFT的一种快速算法，FFT的应用实质上是DFT的应用。FFT可直接用来处理离散信号数据,也可用于对连续时间信号分析的逼近。【实例1】单边指数函数的DFT函数傅立叶变换1、对连续时间变换的逼近FT图DFT图令β=1，a=1令β=1，a=1，利用DFT方法，按下列步骤得到离散变换图：①时域采样，设N=32，Ts=0.25，得到每点样值x(nTs)；②计算

可见，DFT是FT的逼近，其实部是偶函数，在频域k=N/2点对称，在k>N/2时，代表了负频率点处理的结果。而虚部为奇函数，k>N/2处是负频率处理结果。【实例2】方波的谐波分析将DFT用于方波的谐波分析，需要计算以得到各次谐波系数值。对方波在时域的采样点数N=32。在k=N/2点处对称。可以看出，低次谐波比较逼近，而高次谐波有误差，这是由于频率混叠效应所致。虽然可以通过提高采样频率来减少这一现象，但不可能完全避免，因为周期方波为时域无限、频域无限信号。2、卷积运算（1）快速卷积运算过程长度为N1的序列x(n)和长度为N2的序列h(n)卷积，其结果y(n)长度为N1+N2-1●卷积运算中，每个x(n)的样值必须与每个h(n)的样值相乘，因此，共需要N1×N2次乘法运算。

●如果把线卷积改为求圆卷积，并借助FFT技术，可减少运算量。●快速卷积运算过程实现快速卷积算法中，由于利用了DFT分析，即时域或频域都是周期性的离散数据，当对他们作卷积运算时，将出现一种周期数据之间的叠带求和现象，给计算结果带来一种所谓的环绕误差。下面分析其产生原因和避免方法。第一步：利用FFT算法计算两信号的DFT第二步：在各频率点处两信号的变换相乘第三步：运用IFFT算法，计算变换式乘积的反变换

实现这一过程共需两次FFT和一次IFFT运算（相当于三次FFT运算），此外，完成X(k)与H(k)两序列相乘，需作N次乘法。在一般的有限冲激响应（FIR）数字滤波器中，由h(n)求H(k)这一步是预先设计好的，数据已置于存贮器中，故实际只需两次FFT的运算量。如果假定N1=N2=N，则全部复数乘法运算次数为。可见，随N值增大，计算量显著减少，故圆卷积的方案可以快速完成卷积运算。（2）圆卷积的环绕误差●环绕误差产生原因

图①为两个非周期离散序列的卷积。采用直接线卷积或补零圆卷积方法很容易求得其计算结果[图中(g)]。两个非周期离散序列的卷积①两个周期离散序列的圆卷积②

原因：当采用DFT分析方法，上述两个非周期离散信号被改造为时域、频域相对应的周期离散信号，导致卷积结果与图①不同（见图2）。主要区别：当h(n-m)向右移动时，h(n-m)另一周期的一部分进入到求和区域，导致错误的计算结果，被称为“环绕误差”或“叠带效应”。●避免环绕误差的方法

方法：对x(n)与h(n)分别在尾部填补N1(x(n)的样点数)与N2(h(n)的样点数）零值点，即使其周期加倍。如果x(n)与h(n)的长度相等，则都加长N点。采用补点方法后所得计算结果如图。

补点后的副作用：对x(n)与h(n)进行补点后，避免了环绕误差的同时，使的DFT(或FFT)算法所需容量加倍，在各DFT表示式中，必须用2N来代替N，在各函数的尾部补填足够的零值使有效周期为2N，这就能够对两个含有N点的函数进行正常的卷积运算。两个含有N点的非周期性函数的离散卷积给出一个具有2N-1点新函数，当在原函数上填加N个零点后，所得圆卷积的周期为2N，比原来的非周期信号的卷积多一个点，每周内多一个附加零点。如果这两函数相当靠近，但长度不等，则首先将短函数的尾部补零使与长函数的长度相等，然后再补零到2N-1，这也即对较短信号补充的零点数总共超过了N个。总结以上各点，快速卷积过程可按如下步骤进行：（a）用补零法修正x(n)和h(n)，以避免环绕误差的出现。（b）用FFT算法计算两个修正后的函数的DFT，得到X(k)与H(k)（c）将X(k)与H(k)相乘，得到（d）利用FFT算法，计算出Y(k)的IDFT，即3、相关运算相关函数的数字计算方法有时域直接计算与FFT快速算法。直接计算方法是依据下述定义进行的，即互相关函数（在数字信号分析中，一般用符号r）自相关函数：这种计算方法与卷积运算相类同（卷积多一个时间反折），也是一个乘、加序列，所需计算量很大。★更正★更正

●相关函数的FFT算法，依据的是维纳-辛钦关系，即自相关函数或互相关函数可以由功率谱密度或互谱密度函数来求得。

①这种方法是一种迂回的方法，但它比直接时域计算方法快5～100倍。

②当运用FFT方法计算相关函数时，也必须注意到环绕误差的影响，它类似于圆卷积中的误差。解决的方法也是对时间序列x(n)与y(n)补零扩展。第一步：对x(n)和y(n)作FFT分析，得到复频谱X(k)与Y(k)第二步：对X(k)与Y(k)作共轭乘积，得到或当x(n)与y(n)相同时，得到自功率谱密度:

第三步：作IFFT分析，从功率谱密度获得相关函数，即

★本节内容：随机信号的谱分析与谱估计技术，即对功率谱密度的传统估计方法。

★DFT与FFT是信号处理的重要工具，尤其是DFT的基本概念与算法，为信号频谱分析的各种应用铺平了道路。

★传统谱分析方法，是基于Fourier变换的谱分析方法，包括相关函数法与周期图法。第四节谱分析与谱估计

★

谱分析与谱估计在生产实践与科学研究中获得了日益广泛的应用。【例1】在声纳系统中，通过对噪声信号进行谱分析，以判断水面舰艇或潜艇的运动速度、方向、位置、大小等；【例2】对飞机、轮船、汽车、汽轮机、电机、机床等主体或部件进行实际运动的谱分析，可以提供设计数据和检验设计效果，或者寻找振源和诊断故障。

1）相关函数法(又称BT法)1958年由布莱克曼-图基(Blackman-Tukey)提出。它是通过统计分析，从时域上先求信号的自相关函数，再作Fourier变换，求得功率谱估计值。2）周期图法它是直接将数据进行Fourier变换，再取其幅度平方，得到信号的功率谱密度。此法是在1965年由Cooley-Tukey提出的FFT方法问世以后，被用于谱估计。★特别注意一、周期图法作功率谱估计

★自相关函数作为时移的函数是最能较完整地表征随机信号的特定统计平均量值的。而随机信号的功率谱密度，正是自相关函数的Fourier变换。

★对于随机信号而言，其自身的Fourier变换是不存在的，只能用功率谱密度来表征它的统计平均谱特性。根据相关定理与维纳-辛钦关系式（参见图2-5）易于证明随机信号序列x(n)的功率谱密度：1)“泄漏”问题以上两种传统方法本质上是一样的，都认为“有限长的数据段，可以看作是无限长的取样序列给予开窗截断后的结果”。不论是数据开窗，还是自相关函数开窗，在频率域内都会发生“泄漏”现象，即功率谱主瓣内的能量泄漏到旁瓣内。这样，弱信号的主瓣很容易被强信号的旁瓣所淹没或歪曲，造成谱的模糊与失真。

2）所有旁瓣抑制技术，都是以损失谱分辨率为代价。3）在谱分析应用中，频率分辨率与低旁瓣一样是个重要指标，有时甚至更重要。因此，解决高分辨率与低旁瓣的矛盾是谱分析中的一个重点问题。

●

为什么叫周期图法？

由于序列x(n)的离散傅里叶变换X(k)具有周期函数的性质，故把它称为长度为N的实平稳随机信号序列x(n)的周期图。

●

周期图法作谱估计时，存在的两个主要问题：功率谱密度的统计变异性和能量泄漏。前者是统计误差，是由于在功率谱测量中收集到的数据数量有限，而引起的不确定度。后者是谱分析中所固有的，它将造成估计的偏度误差。因此，实际应用中对周期图法进行修改，以尽量减小估计误差。其方法是：（1）采取平均化处理，以减少统计变异性；（2）采用窗处理，以减少泄漏。其估计值周期图法计算功率谱密度流程

统计变异性产生原因——对于有限长数据的处理，由于数据的概率性质，而产生了统计性误差。二、平均化处理（以减少统计变异性）一般用变异系数（Coefficientofvariation，或称为标准化标准偏差）来表征谱估计的变异性质，其定义为谱估计的均方差均值当用周期图法作谱估计时，周期图是一个复数，故而有可以证明，实部XR(k)与虚部XI(k)是等方差和零均值的两个不相关的随机变量。由于Fourier变换是线性运算，如果被分析数据x(n)是高斯分布，则XR(k)与XI(k)也是高斯随机变量。|X(k)|2=X2R(k)+X2I(k)表明，谱估计值是两个独立高斯变量的平方和，它相当于具有两个自由度的卡埃平方(Chi-Square)分布.从概率统计学可知，卡埃平方分布自由度为n的卡埃平方变量n为独立变量数，即自由度

zi服从正态分布的独立随机变量1、谱估计的变异性变量χ2的均值和方差：变异系数

显然，用式作谱估计时，相当于具有自由度n=2（实部和虚部），其估计的变异系数εr=1。这表明估计的相对误差达到100%，即估计的变异性和被估计量一样大，这样的估计是不合用的。为此可采用平均化处理方法来提高估计精度.平均周期图的方法——是将序列x(n)分段，求各段周期图，再进行平均。2、平均化处理方法设序列x(n)[或x(t)]总体长度为N（或T），将其分为q段，每段长度为Nq（或Te），对每段数据Nq作谱估计，得到

，此时有如果各段频率分量的实部XRj(k)与虚部XIj(k)是互为独立的随机变量，那么有此时χ2分布的自由度n=2q。将各段估计谱在对应的频率上作q个估计量的平均，得

N此时变异系数为

●分析：对于连续随机过程，样本总体长度为T(T=N△t)，△t为采样间隔，分段长度Te(T/q)，那么分析带宽Be=1/Te，故有。这一关系式与模拟分析方法中所得到的结论是一致的，欲提高谱估计精度，必须同时考虑到样本总体长度T与频率分析带宽Be。三、窗口函数因不可能对无限长信号进行分析，故必须截断。假定截断区间为(-T,T)，因为对|τ|>T时的Rx(τ)值假定为零，所得到的估计谱为近似谱，即正弦信号的真谱与估计谱之间的关系1、真谱与估计谱根据维纳-辛钦定理，理论谱密度的定义为：两者间总是有误差例如，正弦信号的真谱与估计谱之间的关系。估计谱是真谱与窗谱的卷积，即理想窗谱W(ω)应为δ函数。此时估计谱与真谱完全相同。但要得到δ函数形窗谱，其时域窗口必然是无限宽，这不可能。实际窗口为有限宽，窗谱W(ω)为sinc(t)型函数,故有主瓣能量泄漏到旁瓣.

2、窗函数★窗函数设计的目的——改善窗谱形状。

★

窗函数的基本要求——窗谱的主瓣要窄且高，以提高分辨率；旁瓣应小，正负交替接近相等，以减小泄漏或负谱现象。

★加窗的其他作用——可抑制噪声，提高频率辨识能力。例如，1）冲击测量或脉冲激振时，实际有用信号延续时间很短。如果采样时间较长，可利用窗口控制，避免在余留样本时间内，无实际信号输入时而有噪声混入。2）在系统辨识中，小阻尼系统的响应衰减缓慢。如在较短时间内截断，会丢掉有用信息；若加指数窗口，使其增加衰减速度，既保留了有用信息，又可防止噪声信号混入。★改善窗谱形状的基本思想——改善截断处的不连续状态。Gibbs现象的研究已经表明：时域内的间断，反映到频域，必然发生振荡现象；反之，频域内的间断，反映到时域，也同样发生振荡现象。

1）幂窗：采用时间变量某种幂次的函数，如矩形、三角形、梯形或其他时间t的高次幂2）三角函数窗：应用正弦或余弦函数等组合成复合函数，如汉宁窗、海明窗等3）指数窗：采用指数时间函数，如e-at，如高斯窗等★窗函数主要类型（1）矩形窗（时间变量的零次幂窗）

优点是主瓣比较集中；缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负谱现象。（2）三角窗（也称费杰（Fejer）窗，是幂窗的一次方形式）特点：主瓣宽约为矩形窗的两倍，旁瓣小，且无负旁瓣（3）汉宁(Hanning)窗（又称升余弦窗）汉宁窗w(t)●汉宁窗的谱窗为3个矩形窗的频谱之和，即3个sinc(t)型函数之和。括号中的两项相对于第一个谱窗向左、右各移动了π/T，从而使旁瓣互相抵消，消去高频干扰和漏能.

●汉宁窗与矩形窗的谱图对比图(a)为W(ω)－ω关系；图(b)为相对幅度(相对于主瓣衰减)－logω关系。可以看出，汉宁窗主瓣加宽（第一个零点在2π/T处）并降低，旁瓣则显著减小。第一个旁瓣衰减-32dB，而矩形窗第一个旁瓣衰减-13dB。此外，汉宁窗的旁瓣衰减速度也较快，约为60dB/10oct，而矩形窗为20dB/10oct。由以上比较可知，从减小泄漏观点出发，汉宁窗优于矩形窗。但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨力下降。octave——倍频程（4）海明(Hamming)窗（也是余弦窗的一种，为改进的升余弦窗）

海明窗与汉宁窗都是余弦窗，只是加权系数不同。海明窗的加权系数使旁瓣更小。分析表明，海明窗的第一旁瓣衰减为-42dB。海明窗的频谱也是由3个矩形时窗的频谱合成，但其旁瓣衰减速度为20dB/10oct，比汉宁窗衰减速度(约为60dB/10oct)慢。（5）高斯窗（是一种指数窗）a为常数，决定了函数曲线衰减的快慢。适当选取a值，可使截断点（T为有限值）处的函数值比较小，则截断造成的影响就比较小。高斯窗谱无负的旁瓣，主瓣较宽，故频率分辨力低。第一旁瓣衰减达-55dB。高斯窗函数常用来截断一些非周期信号，如指数衰减信号等。

●除以上常用窗函数外，还有多种窗函数，如帕仁(Parzen)窗、布莱克曼(Blackman)窗、凯塞(Kaiser)窗等。●五种典型窗函数的性能特点

【例】矩形窗、三角窗、汉宁窗分析同一信号数据。(a)为被分析信号的真谱；(b)和(c)是用两种时宽(T/2与T)的矩形窗分析的结果。可见，时域窗口窄，分辨率低，相邻的两谱线不能分辨。(d)和(e)是分别用三角窗与汉宁窗分析的结果，与图(c)相比，三角窗与汉宁窗分辨率低,不能分辨相邻谱线，但由于旁瓣衰减快,谱的分布区域窄而边沿清晰。窗函数类型-3dB带宽等效噪声带宽旁瓣幅度(dB)旁瓣衰减速度(dB/10oct)矩形三角形汉宁海明高斯0.89B1.28B1.20B1.30B1.55BB1.33B1.23B1.36B1.64B-13-27-32-42-55-20-60-60-20-20(a)(b)(c)(d)(e)

●窗函数的选择应考虑被分析信号的性质与处理要求。▲若仅要求精确读出主瓣频率，而不考虑幅值精度，可选主瓣宽度较窄的矩形窗，例如测量物体的自振频率等；▲如果分析窄带信号，且有较强的干扰噪声，则选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；▲对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。一、最大熵谱估计的基本原理

●

本节将依据已经证明的最大熵定理（3-4）阐明这一方法的计算原理。

●最大熵谱分析方法把信息熵的概念引入信号处理中，有时又称为现代时序谱分析方法。这是一种把自相关函数外推的方法。在分析过程中，没有固定的窗口函数。在每一步外推自相关函数中，使估计的相关函数包含过程的信息最多，即要求在过程的熵达到最大的条件下，确定未知的自相关函数值，借以达到谱估计的逼真和稳定度最好的目的。换句话说，就是采用谱熵为最大的准则来估计功率谱。前面已经证明，N维高斯随机序列信源的相对熵为【参见式（3-15）】

：|R|(或记为det[R])表示矩阵[R]的行列式。自相关矩阵是