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

第六章数字信号分析（）DFT与FFT,计算机技术的发展是数字信号分析的基础。,数字滤波优于模拟滤波，（1）速度快，例如，采用数字信号分析技术对于1024采样点进行A/D转换，仅需415s，进行FFT运算须250ms，较快的只需数毫秒；一个蝶形FFT硬件运算只需2s。（2）分辨力高，在高频段（50kHz)可达25Hz；在超低频段可达0.0025Hz。,数字信号分析一般包括：频谱分析与数字滤波等主要内容。前者又包含相关与统计分析等。,第一节模拟信号离散化,本章内容：主要介绍离散Fourier变换与快速Fourier变换的基本原理及应用。,一、A/D与D/A转换,1、A/D转换过程,（1）采样（或称抽样）,连续时间信号x(t)离散采样信号x(nt)。,量化误差呈等概率均匀分布，其概率密度函数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(ns/2)规律变化的函数。,2、频域采样,则有,可见：若频谱X()被间隔为1的脉冲序列在频域中采样，则在时域中等效于x(t)以T1(=2/1)为周期而重复。就是说，周期信号的频谱是离散的。,根据,因此，,上述分析证明：信号的时域与频域呈采样（离散）与重复（周期)关系。,采样定理：s2m或fs2fm。因为时域采样间隔决定于fs，故又称为时域采样定理。,三、采样定理,1、频混现象,Fxs(t)为周期谱，其周期s=2/Ts,s2m,周期谱图相互分离,sN/2时，代表了负频率点处理的结果。而虚部为奇函数，kN/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)的样值相乘，因此，共需要N1N2次乘法运算。,如果把线卷积改为求圆卷积，并借助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算法，依据的是维纳-辛钦关系，即自相关函数或互相关函数可以由功率谱密度或互谱密度函数来求得。,这种方法是一种迂回的方法，但它比直接时域计算方法快5100倍。当运用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)分布.,从概率统计学可知，卡埃平方分布,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个估计量的平均，得,此时变异系数为,分析：对于连续随机过程，样本总体长度为T(T=Nt)，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)相比，三角窗与汉宁窗分辨率低,不能分辨相邻谱线，但由于旁瓣衰减快,谱的分布区域窄而边沿清晰。,窗函数的选择应考虑被分析信号的性质与处理要求。若仅要求精确读出主瓣频率，而不考虑幅值精度，可选主瓣宽度较窄的矩形窗，例如测量物体的自振频率等；如果分析窄带信号，且有较强的干扰噪声，则选用旁瓣幅度小的窗函数，如汉宁窗、三角窗等；对于随时间按指数衰减的函数，可采用指数窗来提高信噪比。,一、最大熵谱估计的基本原理,本节将依据已经证明的最大熵定理（3-4）阐明这一方法的计算原理。,最大熵谱分析方法把信息熵的概念引入信号处理中，有时又称为现代时序谱分析方法。这是一种把自相关函数外推的方法。在分析过程中，没有固定的窗口函数。在每一步外推自相关函数中，使估计的相关函数包含过程的信息最多，即要求在过程的熵达到最大的条件下，确定未知的自相关函数值，借以达到谱估计的逼真和稳定度最好的目的。换句话说，就是采用谱熵为最大的准则来估计功率谱。,前面已经证明，N维高斯随机序列信源的相对熵为【参见式（3-15）】：,|R|(或记为detR)表示矩阵R的行列式。自相关矩阵是托布列兹矩阵，可写成,determinant行列式,第五节现代谱分析方法最大熵谱估计,根据维纳-辛钦(Wiener-Khintchine)关系式，功率谱密度与自相关函数之间互为Fourier变换关系。可以证明，自相关函数矩阵R的行列式detR的极限与功率谱密度函数S(f)的关系为：,熵率(Entropyrate)由上式可知，当序列N时，熵将发散，这时就不能用它来测量信息量。在此情况下可用熵率表示。熵率定义为,此式表明，信号序列的熵密度（熵率）是时间序列的功率谱密度的对数的Nyquist频率范围的积分。,当随机信号按Nyquist采样频率产生随机序列时，采样周期t=1/(2fN)。因此，可得到随机信号序列的熵率：,根据维纳-辛钦关系式，S(f)与自相关函数R(k)之间的关系为：,在维纳-辛钦关系约束下，当由观测到的随机序列数据估计自相关函数，并使每步估计的熵为最大时，即有,式中，功率谱密度S(f)应受到上面熵率导数等于零的条件约束，它的倒数可用截短的Fourier级数表示，即把1/S(f)展开为2m+1个有限级数之和，则得到,解Yule-Walker方程式，由式DFTx(n-m)=X(k)Wmk便得到MEM(最大熵法MaximumEntropyMethod)谱：,上式中，t为采样间隔；ak为预测滤波器的系数（或AR模型参数）；Pm为预测误差的方差（或白噪声序列方差）；a*k表示ak的复共轭。比较上式左右两边z的等幂项，得到如下方程式，即Yule-Walker(尤尔-沃克)方程式：,此式表明，MEM谱与参数ak和Pm有关。因此，如何根据观测数据的自相关函数，快速、精确地计算参数ak和Pm，是MEM分析方法的中心内容。不同的参数估计算法，就产生