基于Matlab的1_3倍频程分析

1、用Matlab语言实现噪声的1/3倍频程分析摘 要: 在声学测量研究中，1/3倍频程谱反映了声源的能量分布情况.本文基于Matlab软件开发平台，实现了对高斯白噪声的1/3倍频程分析，验证了该算法的正确性，具有精度高,性能稳定的特点。关键词: 故障诊断；传感器优化布置；高斯白噪声；功率谱；1/3倍频程0 引言“设备故障诊断（Condition Monitoring and Faults Diagnosis）”是近十几年发展起来的一门新兴技术，包含两方面的内容：一是对设备的现场运行状态进行监测；二是在出现故障情况时对故障进行分析与诊断。这二者是密不可分和相互关联的。掌握设备现在的状况及信息，预知

2、和预测有关故障或异常的程度，分析故障产生的原因，判断故障发展趋势及其对将来的影响，从而找出必要的对策或解决方法，是设备故障诊断的功能。 运用设备诊断技术所取得的经济效益是明显的。据日本资料报道，采用诊断技术后，事故率减少 75%，维修费用可降低 25-50%。英国对两千个工厂进行的调查表明，采用诊断技术后，维修费用每年可节约 3 亿英磅。目前我国的机器设备总值约为 8000 多亿元，每年用于设备大修、小修和处理故障的费用一般占固定资产原值的 3%5%，采用诊断技术改进维修方式和方法后，一年取得的经济效益可达数百亿元。因而减少停机时间而创造的社会效益将非常巨大。显然，设备故障诊断与监测技术对企业

3、的正常生产经营是必不可少的，必须把它作为企业管理与发展的一个重要内容。设备故障诊断一般分两个阶段四个步骤实施。两个阶段为状态监测和故障诊断。故障诊断的四个步骤为：信号检测、特征提取（信号处理）、状态识别和诊断决策。其具体内容为：（l）信号检测：按不同诊断目的选择最能表征工作状态的信号。这种工作状态信号称为初始模式。（2）特征提取（信号处理）：将初始模式向量进行信号处理、变换，去掉冗余信息，提取故障特征、形成待检模式；（3）状态分析：将待检模式与样式模式（故障档案）对比和状态分类，判断出故障类型。（4）诊断决策：根据判别结果采取相应对策，对机械设备及工作进行必要的预测及修正。 1 传感器优化布置

4、结构健康监测检测是近年来发展起来的结构无损检测技术。它利用结构的某些信息,运用一定的数学方法,来判定结构是否损伤以及损伤的位置和程度。近10年来,国内外学者普遍认同的损伤评估方法是试验模态分析法。进行模态实验的第1步就是获得被测结构激励和响应的时域信号,而传感器的配置方案是首先要确定的。不适当的传感器配置将影响识别参数的精度,而且传感器本身需要一定的成本,与其配套使用的数据采集和处理设备的代价也都较高。从经济方面考虑,希望采用尽可能少的传感器。因此,确定传感器的最佳数目,并将它们配置在最优位置,具有重要的实用价值。 要进行传感器的优化配置,首先要确定合理的、能反映设计要求的优化配置准则。目前发

5、展起来的优化准则很多,其中,基于识别误差最小准则的方法是使用较多的1种方法;如:Kammer提出的有效独立法;模型缩减准则也是1种常用的准则,但这种方法只能保证低阶模态的精度;另外还有可控可观度准则,模态应变能准则等。Came等认为模态置信度MAC(Modal Assurance Criterion)矩阵是评价模态向量交角的1个很好的工具。其次,传感器的优化配置还必须选用适当的优化方法。近几年发展起来的随机类方法主要有模拟退火算法和遗传算法,这种方法不易陷入局部最优解,但现有的随机搜索技术可靠性并不高。目前使用最多的是序列法中的逐步削减法,它每次从剩余传感器的可选位置中去掉1个或多个对目标函数

6、贡献最小或较小的可选位置,直到剩余最优可选位置为止。逐步累积法与削减法相反,它是不断的从剩余可选位置中选取1个最优的加入到优化配置中,直至达到最优的数目为止。本文在此基础上提出了修正的逐步累计法,算例表明,修正的逐步累计法具有较高的计算精度和优化结果。1.1 基本理论 1.1.1基于QR分解的传感器配置原理 根据模态叠加原理,系统的响应可以表示为: (1-1)其中其中,u为物理坐标,u,为第i阶模态向量,为模态向量矩阵,为第i阶模态坐标,q。s代表传感器数量,m为所需识别的模态个数。式(1)的最小二乘解为 (1-2)如果考虑测量噪声,式(1)应改写为 (1-3)其中,v代表方差为的高斯分布白噪

7、声,这里假设测量噪声相互独立并且对各个传感器测量信号的统计特性相同,则与q的协方差2为: （1-4）其中,Q称为Fisher信息矩阵。当Q取极大值时,协方差P最小,就能够得到较好的估计。所以必须使Q的某1种范数最大,这里选取常用的2-范数Q2。由于 （1-5）所以 (1-6)因此,以上对Q的要求可通过的选择来实现。根据矩阵理论,列主元QR分解是选取矩阵列向量组具有较大范数子集的1种简捷有效的方法。设有限元模型所得的模态向量矩阵对应于可测自由度的子集为5,5。一般有m<n,并且r()=m,即矩阵列满秩。由于列主元QR分解选择的是列向量组的子集,所以,进行的列主元QR分解: （1-7）其中,

8、E为置换矩阵,则中对应于1 ,2, n的行(即自由度)就是的行向量组中具有较大范数的子集(其中代表R矩阵的第i列)。1.1.2模态置信度MAC矩阵 由结构动力学原理可知,结构各阶模态向量在节点上的值形成了1组正交向量,但由于量测自由度远小于结构模型的自由度并且受到测试精度和测量噪声的影响,测得的模态向量已不可能保证其正交性,在极端的情况下甚至会由于向量间的空间交角过小而丢失重要的模态。因此,在选择测点时有必要使量测的各模态向量保持较大的空间交角,从而尽可能地把原来模型的特性保留下来。Kammer的EI法事实上从另1个角度起到了这个作用。Came等认为模态置信度MAC矩阵是评价模态向量空间交角的

9、1个很好的工具,其公式表达如下: （1-8）其中和分别为第i阶和第j阶模态向量。由式(1-8)可以看出,MAC矩阵考虑的是模态向量矩阵列空间的度量特性。而根据1.1.1中的讨论,利用QR分解的传感器配置实际上讨论的是模态向量矩阵行空间的特性,所以虽然这种传感器配置可以保证q估计值的质量,但并不一定能够得到良好的MAC,必须采取新的措施来提高MAC以满足振型匹配的要求。由式(1-8)可以看出,MAC矩阵的非对角元(ij)代表了相应2模态向量的交角状况。换句话说,当MAC阵的某一非对角元 (ij)=1时,表明了第i向量与第j向量交角为0,2向量不可分辨;而当 (ij)=0时,则表明第i向量与第j向

10、量相互正交,2向量较易识别,故测点的布置应力求使MAC的非对角元向最小化发展。2 噪声的物理量度 2.1.1噪声的基本概念各种频率和声强杂乱无序组合的声音叫做噪声。声音是由物体振动引起的。物体振动通过媒质中传播引起人耳或其它接收器的反应，就是声音。振动的物体是声音的声源，产生噪声的物体或机械设备等成为噪声源。振动在弹性介质中以波的形式传播，这种弹性波叫做声波。当振源频率在2020000Hz之间时，人的耳朵可以感受到它。当振源频率低于20Hz或高于20000Hz时，人耳无法听到。低于20Hz的波动叫次声波，高于20000Hz的波动叫超声波。人耳日常听到的声音，通常来自空气所传播频率在202000

11、0Hz之间的声波。多年来噪声信号都是作为有害信号在系统设计或改进时加以抑制。在多数以控制为目的的测量中，传感器的设计都是尽量排除噪声信号。其实，噪声信号中带有大量设备运行状态信息，通过对噪声信号进行分离和分析，可以了解设备运行状态，对设备进行状态检测和故障诊断。近场噪声分析方法是一种非接触测量方法，它是利用声音传感器接近被测部件来测取噪声信号，通过分析噪声信号进行频谱分析，比较被测信号频谱特性与正常运行状态下的频谱特性来诊断机械故障。2.1.2 声压和声压级当没有声波时，空气处于静止状态时，其压强为大气压强P0，当有声压存在时，局部空气产生压缩或膨胀，在压缩的地方压强增加，在膨胀的地方压强减少

12、，这样就在原来大气压上又迭加了一个压强的变化量。这个迭加上去的压强变化与静压强的差值称为声压，用Pe表示。声压可正可负是围绕钧值摆动的量，其均值恒为零，通常仪器检测的声压为声压的均方根值Pe，称为有效声压。 （Pa） （2-1）其中T为周期，根据简谐波，Pm为声压幅值。声压的大小表示了声波的强弱。由于正常人耳能听到的最弱声音的声压和能使人耳感到疼痛的声音的声压大小之间相差一百万倍，表达和应用起来很不方便。同时，实际上人耳对声音的感受也不是线性的，它不是正比与声压绝对值的大小，而是同它的对数近似成正比。因此如果将两个声音的声压之比用对数的标度来表示，那么不仅应用简单，而且也接近于人耳的听觉特性。

13、这种用对数标度来表示的声压称为声音级，它用dB表示。即 （2-2）其中，Lp为声压级，单位为分贝（dB），P0是参考声压级，国际上规定P020x10Pa，这就是人耳刚能听到的1KHz的纯音的声压值。当采用声压级的概念后，听阈与痛阈的声压比从100万倍的变化范围变成0120dB的变化，所以“级”的大小能衡量声音的相对强弱，声压级被定义成声压平方比值的对数值。由于声压平方的比值也与声功率成比例，因此声压级与声功率级就联系起来了。2.1.3 声强与声强级在声波传播方向上单位时间内垂直通过单位面积的声能量，成为声音的强度或简称声强，用I表示，单位W/。声强的大小可用来衡量声音的强弱，声强越大，我们听到

14、的声音越强；声强越小，我们听到的声音越轻。声强与离开声源的距离有关，距离越远，声强越小。根据定义，瞬时声强I可表示为I（t）P（t）*u（t） （2-3）式中P（t）为瞬时声压，u（t）为质点的振速，由于u（t）是矢量，所以I（t） 也是矢量，即代表声能的传播方向。与声压一样，声强用“级”来表示，即声强级,它的单位也是分贝（dB）,即 （2-4）其中为基准声强，W/，它相当于人耳能听到的最弱的声音的强度。声强级与声压级的关系是：式中为空气密度，c为空气中声速。媒质的c随媒质的温度和气压而改变。如果在测量条件时恰好c400，则。对一般情况，声压级与声强级相差一个修正项，此值是比较小的。3 幅度谱

15、分析方法噪声控制中的数字化技术是指对噪声信号的检测,变换,记录,传输和分析处理采用数字编码信号,通过数字计算方法对信号进行滤波,检测,快速傅立叶变换等,用以提高噪声信号的处理速度和精度。3.1傅立叶变换随机振动,可以认为是一个具有零均值的各态历经的平稳高斯过程1,"因此现场测试所得的振动时程曲线可以视为许多不同频率的正弦波叠加的结果,即可以用傅立叶级数的形式表示"有限时间内非周期连续函数的傅立叶变换为正变换:逆变换:实际的采样信号总是离散的并且信号长度是有限的, 对等时距的离散点进行傅立叶分析称为离散点的傅立叶变换: (n=0,1,2,N-1)式中,N为采样数据量。在噪声测

16、量中，采用1/3倍频程谱分析能详细的反映出噪声源的频谱特性，便于较全面的了解声源产生激励和提出最佳的降噪对策。本文就是基于上述背景和目的，应用Matlab语言实现了1/3倍频程分析功能，并将结果以柱状图的形式可视化，给人以直观的视觉感受。3.2 1/3倍频程的实现一般噪声的频率分布宽阔，在实际的频谱分析中，不需要也不可能对每个频率成分进行具体分析，为了便于观察噪声信号宏观上的能量分布，忽略信号频率或相位信息微小变化对观察结果的影响，把2020000Hz的声频范围分为几个段落，每个段落称为频带或频程，频带的中心频率为：，和为该频带的上限和下限频率。频带的划分也有规定，一般规定集中m倍频程，m由下

17、式确定：当m1时，称为倍频程，类推，没1/3时，称为1/3倍频程。在此频带（）的频谱，称为1/3倍频程谱。按照1/3倍频程的方法，可将声频范围分为更多的频带，便于较仔细地研究。对于噪声来讲，由于1/3倍频程能够很好地体现噪声带宽地能量分布情况，为噪声控制提供参数，采取合理地措施，从而达到降噪地目的，所以在噪声分析过程1/3倍频程地分析显得尤为重要。我国对1/3倍频程地频率范围及其中心频率都做了规定，程序中按照上、下限的规定将功谱中的频率划分成多个频带,在每个频带中将所有的功率谱幅值相加、平均,这样就借助功率谱实现了 1/3 倍频程功能,最终结果以柱状图的形式显示。这里需要注意的是，在求和的过程

18、中要用原始的测量值（噪声为Pa），如果要使用声压级dB（A）,需要转化为Pa，噪声dB（A）与Pa的关系：dB=20xlog（Pa/20-5）（在空气中），所以在叠加时，一定要采取对数加法，否则，就会导致算法错误。1/3倍频程的实现流程如图1所示（请见下页）：三分之一倍频程的中心频率和频率范围如表1所示：读入数据进行功率谱分析将功率谱幅值dB转换为帕（Pa） Pa20x按1/3倍频程的频率范围对功率谱域进行划分，在22.4Hz到18000Hz内形成29个频段。在每个频段内将所有的功率谱密度相加，平均，并以柱状图形式显示结果。生成信号图1 1/3倍频程的实现流程表1 1/3倍频程的中心频率和频率

19、范围 (Hz)中心频率频率范围中心频率频率范围2522.42880071090031.52835.5100090011204035.545125011201400504556160014001800635671200018002240807190250022402800100901123150280035501251121404000355045001601401805000450056002001802246300560071002502242808000710090003102803551000090001120040035545012500112001400050045056016000

20、14000180006305607104 程序结果此次程序验证信号为高斯白噪声,结果如下所示:1 功率谱函数2 柱状图显示5 实验结论由白噪声的频率特性可知,在相等带宽内能量是相同的。用倍频程分析，则中心频率增加一倍，能量就增加一倍，即3dB；而用1/3倍频程进行分析时，中心频率增加1/3倍，能量就增加1/3倍，即1dB；由上图可知，该算法的正确性（由于低频部分的fft的线速很少，与理论有一些偏差，而中高频部分很好的反映了这一特性）。与传统的方法相比，不但计算速度显著提高了，而且计算精度也提高了。6实验体会此次科研训练为期32个学时，由于对专业知识的匮乏，几经波折，终于完成。通过这次科研训练，

21、学到了不少东西：首先，最大的体会是，自己的基础很薄弱，对很多东西都不了解，碰到了很多问题。因此，在以后的学习中，需要加强对专业知识的学习。其次，我对科研的过程有了一个大体的了解，并且对科学精神有了一个明确的认识，做学问要一丝不苟，不可敷衍了事，不可弄虚作假。第三，学会了查找资料的方法，我现在明白了如果需要学习知识，除了图书馆里藏书和期刊之外，网络数据库也是一个不错的选择，比如中国期刊网，Web of Science，等；并且网络的作用也很大，比如一些专业的论坛，可以去论坛浏览，寻找自己需要的东西，还可以发帖问问题，等等。第四，学会了一些软件的使用。学会了一些Mat lab编制一些简单的程序，使

22、用Ansys对一个模型进行应力，应变分析和模态分析。第五，对实验内容有一个深入的了解。了解了声压级，功率谱分析和1/3倍频程第六，学会了独立解决问题，在程序编制过程中遇到不少问题，出现了几次的反复和错误，通过查阅图书，上网搜集资料，逐一解决。附：Matlab程序带注释%生成高斯白噪声n=8192;y=randn(1,n); y=y-mean(y); y=y/std(y); a=0; b=sqrt(1); y=a+b*y; % 采样频率fs=51200;%频率范围f=22.4,28,35.5,45,56,71,90,112,140,180,224,280,355,450,560,710,900,

23、1120,1400,1800,2240,2800,3550,4500,5600,7100,9000,11200,14000,18000;%中心频率fm=25,31.5,40,50,63,80,100,125,160,200,250,315,400,500,630,800,1000,1250,1600,2000,2500,3150,4000,5000,6300,8000,10000,12500,16000;n3=length(fm);deltf=fs/n;ff=0:fs/n:fs*(n/2-1)/n;yf=do_fft(y);figure,plot(ff,yf);title('高斯白噪声

24、的功率谱');ylabel('功率/Pa');xlabel('频率/Hz');yn=yf.*yf;for i=1:29 f1=f(i); f2=f(i+1); n1=floor(f1/deltf); n2=floor(f2/deltf); temp=0; for j=n1:1:n2 temp1=yn(j)/20-5; yn(j)=20*10(temp1); temp=temp+yn(j); end energy(i)=temp; temp2=log10(energy(i)/(20*10(-5); ydb(i)=temp2*20/deltf; endfig

25、urebar(ydb);title('1/3倍频程分析结果');ylabel('分贝/dB');xlabel('频率/Hz');参考文献【1】 马大猷，噪声与振动手册【M】，机械工业出版社，2002。9【2】 张沛上，噪声防治【M】，电子工业出版社，1988.10【3】 王济，胡晓，Matlab在振动信号处理中的应用，中国水利出版社，2006.1【4】 陈花玲，机械工程测试技术，机械工业出版社，2002.1【5】 王建春等，噪声测量中1/3倍频程与倍频程频谱的关系与计算【J】，噪声与振动控制，1996（5）：39-41【6】 成峰等，基于Matl

26、ab的1/3倍频程FIR数字滤波器设计【J】，福州大学学报（自然科学版），2003，31（2）：160-165【7】 杨永红，一种基于Measurement Studio 计算动态信号1/3倍频程谱的方法【J】，北京工商大学学报（自然科学版），2003，21（3）：37-39【8】 黄晶晶，雷勇，基于VC+的1/3倍频程设计与实现，电子测量技术，2006，29（6）135-136【9】 刘福强，张令弥，作动器传感器优化配置的研究进展【J】，力学进展，2001，3（2）：380【10】 Kammer DC.Sensor placement foron-orbit modalidentification and CorrelationoflargeStructuresJ,Control and Dynamics,1991,14(2):251-259【11】 Amir Fijany and Farrokh Vatan，A New Efficient Algorithm for Analyzing and Optimizingthe System of Sensors，IEEE