频谱分析论文.docx

基于加窗FFT的频谱分析摘要随着计算机和微电子技术的飞速发展，基于数字信号处理的频谱分析已经应用到各个领域并且发挥着重要作用。信号处理方法是当前机械设备故障诊断中重要的技术基础之一，分析结果的精确程度是诊断成功与否的关键因素。研究频谱分析是当前主要的发展方向之一。数字信号处理基本上从两个方面来解决信号的处理问题：一个是时域方法，即数字滤波；另一个是频域方法，即频谱分析。FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域，这样有助于对信号进行分析。Matlab是具有很强的科学计算和图形显示界面的软件系统。本文主要通过用MATLAB编程来对信号进行频谱分析，并分析加矩形窗、汉宁窗、哈明窗等窗函数在不同截断时间长度及不同采样频率下对信号频谱分析的影响。关键词：频谱分析、FFT、窗函数、泄漏误差1 绪论1.1 引言在数字信号处理中，离散傅里叶变换(DiscreteTime Fourier Transform，DFT)是常用的变换方法，它在数字信号处理系统中扮演着重要角色。由离散傅里叶变换发现频率离散化，可以直接用来分析信号的频谱、计数滤波器的频率响应，以及实现信号通过线系统的卷积运算等，因而在信号的频谱分析方面有很大的作用。由于DFT的运算量太大，即使是采用计算机也很难对问题进行实时处理，所以经过很多学者的不懈努力，便出现了通用的快速傅里叶变换（FFT）。快速傅里叶变换（Fast Fourier Transform，FFT）并不是与离散傅里叶变换不同的另一种变换，而是为了减少DFT计算次数的一种快速有效的算法。对FFT算法及其实现方式的研究是很有意义的。目前，FFT己广泛应用在频谱分析、匹配滤波、数字通信、图像处理、语音识别、雷达处理、遥感遥测、地质勘探和无线保密通讯等众多领域。在不同应用场合，需要不同性能要求的FFT处理器。在很多应用领域都要求FFT处理器具有高速度、高精度、大容量和实时处理的性能。因此，如何更快速、更灵活地实现FFT变得越来越重要。1.2 频谱分析的发展频谱分析在生产实践和科学研究中获得日益广泛的应用。例如，在声纳系统中，为了寻找海洋水面船只或潜艇，需要对噪声信号进行谱分析，以提供有用信息，判断舰艇运动速度、方向、位置、大小等。因此对谱分析方法的研究，受到普遍注意和重视，是当前信号处理技术中一个十分活跃的课题。1965年库利首次提出了快速傅里叶变换(FFT)算法，FFT和频谱分析很快发展成为机械设备故障诊断、振动分析、无线电通信、信息图象处理和自动控制等多种学科重要的理论基础。然而长期的应用和近年来的理论分析表明：经快速傅立叶变换得到的离散频谱，频率、幅值和相位均可能产生较大误差，即使加其他窗时，也不能完全消除此影响。从70年代中期，有关学者开始致力于频谱校正理论的研究以期解决离散频谱误差较大的问题。1975年John CBurges等从事电学领域研究工作的学者采用插值法对加矩形窗的离散化频谱进行校正，解决了电学中的离散高次谐波参数的精确测量问题。1996年，余佳兵，史铁林等提出了采用复调制细化谱分析将已产生频谱干涉的密集频率成分分离开，消除干涉，再用比例法进行校正以解决密集频率成分的离散频谱的校正问题。2001年，徐培民、杨积东、闻邦椿提出了自动识别和修正离散频谱中两临近谱峰参数的方法，不仅能识别间距不到一个频率分辨率的两个密集频率成分，而且能识别峰间距为1-6个频率的临近谱峰参数。1.3 本文主要研究内容本实验的目的主要是分析矩形窗、汉宁窗、哈明窗、布莱克曼窗函数对信号矩形频谱分析，分析各窗函数对频谱分析的影响。用MATLAB编出程序，把各种函数的图像进行对比，然后在频率不变的情况下改变截断时间，以及在截断时间不变的情况下，改变频率，观察信号的变化。2离散傅里叶变换2.1离散傅里叶变换的定义设x(n)是一个长度为M的有限长序列，则x(n)的N点离散傅立叶变换为：X(k)=DFTx(n)=,k=0,1,.,N-1。 (1)X(k)的离散傅里叶逆变换为：x(n) =IDFTX(k)= ,k=0,1,.,N-1 (2)式中，N称为DFT变换区间长度，NM。2.2快速傅里叶变换（FFT）快速傅里叶变换（FFT）是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅里叶变换的算法进行改进获得的，它对离散傅里叶变换并没有新的发现。有限长序列x(n)及其频域表示X(k)可由以下离散傅立叶变换得出 （8） （9）其中。式（8）称为离散傅立叶正变换，式（9）称为离散傅立叶逆变换，x(n)与X(k)构成了离散傅立叶变换对。根据上述公式，计算一个X(k)，需要N次复数乘法和N-1次复数加法，而计算全部X(k)( )，共需要次复数乘法和N(N-1)次复数加法。实现一次复数乘法需要四次实数乘法和两次实数加法，一次复数加法需要两次实数加法，因此直接计算全部X(k)共需要4次实数乘法和2N(2N-1)次实数加法。当N较大时，对实时信号处理来说，对处理器计算速度有十分苛刻的要求，于是如何减少计算离散傅里叶变换运算量的问题变得至关重要。为减少运算量，提高运算速度，就必须改进算法。计算DFT过程中需要完成的运算的系数里，存在相当多的对称性。通过研究这种对称性，可以简化计算过程中的运算，从而减少计算DFT所需的时间。如前所述，N点的DFT的复乘次数等于。显然，把N点的DFT分解为几个较短的DFT，可是乘法的次数大大减少。另外，旋转因子具有明显的周期性和对称性，其周期为：其对称性表现为： 或FFT算法就是不断地把长序列的DFT分解成几个短序列的DFT，并利用的周期性和对称性来减少DFT的运算次数。具有以下固有特性：(1)的周期性：(2)的对称性： (3)的可约性：另外，。利用的上述特性，将x(n)或X(k)序列按一定规律分解成短序列进行运算，这样可以避免大量的重复运算，提高计算DFT的运算速度。算法形式有很多种，但基本上可以分为两大类，即按时间抽取(Decimation In Time，DIT)FFT算法和按频率抽取(Decimation In Frequency，DIF)FFT算法。3 用FFT加窗进行频谱分析3.1相关的MATLAB功能函数简介(1)图像显示函数：plot，x轴和y轴均为线性刻度；(2)图形生成函数：figure，生成图形窗口；(3)FFT函数：fft（S），其中S为加窗后的一个帧信号。因为实时信号FFT的频域样值关于中点（即采样频率的1/2）对称，所以fft(S)给出的矩阵数据只有前半部分有用。3.2 相关窗函数简介(1)矩形窗：wRn=RN(n)，能产生矩形波形，MATLAB中用boxcar(N)表示；(2)汉宁窗：wHn=0.51-cos2nN-1RNn，汉宁窗的幅度函数由三部分相加，旁瓣互相对消，能使能量更集中在主瓣中。在MATLAB中用函数hanning(N)表示；(3)哈明窗：wHmn=0.54-0.46cos2nN-1RN(n),这种改进的升余弦窗能量更集中在主瓣中，主瓣的能量约占99.96%，MATLAB中用hamming(N)表示；(4)布莱克曼窗：wBln=0.42-0.5cos2nN-1RN(n)，这样其幅度函数由五部分组成，他们都是移位不同，且幅度也不同的函数，使旁瓣再进一步抵消，旁瓣峰值幅度进一步增加。MATLAB中用Blackman(N)表示。3.3 分析与结论(1)矩形窗Tp=0.1s,Fs=1000hzTp=2s,Fs=1000hzTp=0.1s,Fs=1500矩形窗的主瓣比较集中，缺点是旁瓣较高，并有负旁瓣，导致变换中带进了高频干扰和泄漏，甚至出现负频现象。由上面图像可以看出增大截断时间Tp，即矩形窗口加宽，则窗谱将被压缩变窄（pi/T），旁瓣的影响减小。(2)汉宁窗Tp=0.1s,Fs=1000hzTp=1s,Fs=1000hzTp=0.1,Fs=1500hz汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄漏观点出发，汉宁窗优于矩形窗，但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降。当截断时间Tp增加时旁瓣减小，频率分辨率下降。(3)哈明窗Tp=0.1s,Fs=1000hzTp=1s,Fs=1000hzTp=0.1s,Fs=1500hz哈明窗与汉宁窗都是余弦窗，不过哈明窗能使旁瓣达到更小。分析表明，哈明窗的衰减速度比汉宁窗的慢。(4)布莱克曼窗Tp=0.1s,Fs=1000hzTp=1s,Fs=1000hzTp=0.1s,Fs=1500hz布莱克曼窗的幅度谱主瓣宽度大约是矩形窗的三倍。旁瓣峰值增加。当截断时间一定时，增加采样频率频谱变得光滑。参考文献1高西全，丁玉美.数字信号处理.西安：西安电子科技大学出版社，2008.2程佩青，数字信号处理教程.北京：清华大学出版社，2003.3MATLAB 7.0从入门到精通，求是科技.人民邮电出版社，2006.4张威.MATLAB基础与编程入门.西安：西安电子科技大学出版社，2006.5郑阿奇.MATLAB实用教程M.北京：电子工业出版社.2009.6郭仕剑，王宝顺，贺志国，杨可心.MATLAB7.X数字信号处理M。北京：中国邮电出版社.2006 7赵桂芳等基于DSP的快速傅立叶变换的实现J黄石理工学院学报20078乔瑞萍，崔涛，张芳娟.TMS32054xDSP原理及应用M.西安：西安电子科技大学出版社.2005：9杨宇，叶宇风，王洪.基于DSP的实时信号频谱分析模块设计J.电子测量技术.2006年4月第二期 10 汪安民.TMS320C54x DSP实用技术M.清华大学出版社，2002附录：MATLAB程序：矩形窗clear;Fs=1000;T=1/Fs;Tp=0.1;N=Tp*Fs;w=100*pi;n=1:N;Xn=24*sin(w*n*T)+12*sin(3*w*n*T+40*pi/180)+6*sin(5*w*n*T+60*pi/180)+3*sin(7*w*n*T+(80*pi)/180)+1.5*sin(9*w*n*T+120*pi/180);Xn=Xn/max(abs(Xn);Xk=fft(Xn,4096);fk=Fs*(0:4095);plot(fk,abs(Xk)/max(abs(Xk);xlabel(Hz);ylabel(幅值);title(矩形窗);汉宁窗clear;Fs=1000;T=1/Fs;Tp=0.1;N=Tp*Fs;w=100*pi;n=1:N;Xn=24*sin(w*n*T)+12*sin(3*w*n*T+40*pi/180)+6*sin(5*w*n*T+60*pi/180)+3*sin(7*w*n*T+(80*pi)/180)+1.5*sin(9*w*n*T+120*pi/180);Xn=Xn/max(abs(Xn),wn=hanning(N);Xn1=Xn.*wn;Xk=fft(Xn1,4096);fk=Fs*(0:4095);plot(fk,abs(Xk)/max(abs(Xk);xlabel(Hz);ylabel(幅值);title(汉宁窗);哈明窗clear;Fs=1000;T=1/Fs;Tp=0.1;N=Tp*Fs;w=100*pi;n=1:N;Xn=24*sin(w*n*T)+12*sin(3*w*n*T+40*pi/180)+6*sin(5*w*n*T+60*pi/180)+3*sin(7*w*n*T+(80*pi)/180)+1.5*sin(9*w*n*T+120*pi/180);Xn=Xn/max(abs(Xn),wn=hamming(N);Xn1=Xn.*wn;Xk=fft(Xn1,4096);fk=Fs*(0:4095);plot(fk,abs(Xk)/max(abs(Xk);xlabel(Hz);ylabel(幅值);title(哈明窗);布莱克曼窗clear;Fs=1000;T=1/Fs;Tp=0.1;N=Tp*Fs;w=100*pi;n=1:N;Xn