球磨机信号分析和关键参数预报系统—磨音信号分析方法毕业论文.doc

号：例如TK172 编号：例如BY 15 13/14/2 -沈阳化工大学本科毕业论文题 目：球磨机信号分析和关键参数预报系统磨音信号分析方法 院 系：信息工程学院 专 业：通信工程 班 级：通信1201 学生姓名：王一帆 指导教师：赵立杰 论文提交日期： 年 月 日 论文答辩日期： 年 月 日论文独创性声明本人郑重声明：此处所提交的学位论文是本人在导师指导下，在沈阳化工大学本科毕业环节独立进行研究工作所取得的成果。据本人所知，论文中除已注明部分外不包含他人已发表或撰写过的研究成果。对本文的研究工作做出重要贡献的个人和集体，均已在文中以明确方式注明。本声明的法律结果将完全由本人承担。作者签名： 日期： 年 月 日关于学位论文使用授权的说明本论文的研究成果归沈阳化工大学所有，本论文的研究内容不得以其它单位的名义发表。本学位论文作者和指导教师完全了解沈阳化工大学有关保留、使用学位论文的规定，即：学校有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅；本人授权沈阳化工大学可以将论文的全部或部分内容编入有关数据库进行检索、交流，可以采用影印、缩印或其他复制手段保存论文和汇编本学位论文。（保密的论文在解密后应遵循此规定）作者签名： 导师签名： 日期： 年 月 日摘要球磨机的作用是粉碎并细碎已碎物料，广泛应用于水泥，建筑和耐火材料，化肥，金属选矿以及玻璃陶瓷等各种生产线，对各种尺寸矿石和可磨性物料进行干式或者湿式粉磨加工。由于选矿生产过程具有非线性，多变性，时变性，大滞后，强耦合的特点，各个粉磨要素生产指标不能在线测量，致使选矿工程师不能实时调整整个生产过程，这引起研究人员对选矿生产指标预报算法的研究。算法应用到工业现场需要软件平台的支撑，但是国内外针对选矿生产指标预测软件的研发相对较少，虽然算法研究人员也会根据自己的需求开发预测软件，但是这些软件的人机交互性能友好性，易用性较差，通常难以满足其他算法研究人员。为了实现磨矿过程球磨机负荷参数的在线检测，采用MATLAB和C#.混合编程方式，开发实现了球磨机振动，振声信号分析和关键参数预报的软测量系统。该系统基于集合经验模态分解（Ensemble Empirical Mode Decomposition,EEMD)技术和区间偏最小二乘（interval partial least-squares,IPLS）技术提取与磨机负荷参数密切相关的本征模态函数（intrinsic mode functions,IMF）频域特征，构建基于本征模态函数特征空间的选择性集成模型，实现磨机负荷参数的测量。MATLAB的Deploytool工具将信号EEMD分解，IMF频谱变换，IPLS特征选择，关键参数模型训练和模型预测一系列m函数编译生成DLL程序集合，在C#.net编程环境中，通过调用上述程序实现球磨机信号分析和关键参数预报软件系统的快速开发。系统测试结果表明该系统能够有效地选择简体振动和振声信号IMF频谱特征，系统准确性和可靠性较高，对改进磨矿过程控制和优化具有重要意义。关键词：磨机负荷，EEMD,特征选择，软测量，混合编程 AbstractBall mill is the key equipment used to crush has broken material. It is widely used in cement, building materials, refractory material, fertilizer, metal processing and glass ceramics and other production industries, for different size of ore and materials may be grinding of dry or wet powder grinding. Ball mill is widely used in grinding ore and other materials, is now widely used in mineral dressing, building materials and chemical industry, etc. Through vibration acceleration signal of industrial practice to collect data, parameters of the mill load forecast development. Accurate judgment through the data analysis the working state of the ball mill.Because the mineral processing production process with nonlinear, variability, time-varying, big lag, the characteristics of strong coupling, concentrate grade production indexes such as cannot online measurement, the mineral processing engineers can real-time adjust the whole production process, this caused the researchers study of ore dressing production index prediction algorithm. Algorithm is applied to the industrial field needs the support of software platform, but for beneficiation production index prediction software research and development at home and abroad is relatively small, although the algorithm, the researchers will also according to your requirements development prediction software, but these software friendly human-computer interaction performance, ease of use is poorer, the researchers usually difficult to meet other algorithm. In order to realize the on-line measurement of grinding process parameters of ball mill load, using MATLAB and c # hybrid programming way, development has realized the ball mill vibration, acoustic signal analysis and key parameters of vibration prediction of soft measurement system. The system is based on the collection of Empirical Mode Decomposition (Ensemble Empirical Mode Decomposition, the EEMD) technology and interval partial least squares (interval partial further - squares, IPLS) technology parameter extraction and mill load is closely related to the intrinsic Mode function (the intrinsic Mode functions provides, the IMF) frequency domain feature, build based on intrinsic Mode function feature space selective integration model, realize the mill load parameter measurement. MATLAB EEMD Deploytool tool to signal decomposition, the IMF spectrum transform, IPLS feature selection, the key parameters of model training and model prediction series of m function compiled DLL procedures set in the C#.net programming environment, by calling the above signal analysis program to realize ball mill and the rapid development of key parameters prediction software system. System test results show that the system can effectively select simplified vibration and acoustic signal of vibration spectrum characteristics, the IMF system accuracy and high reliability, to improve the grinding process control and optimization is of great significance.Key words: mill load, EEMD, feature selection, the soft measurement, mixed programming目录第一章 绪论11.1引言11.2球磨机噪声机理2第二章 信号分解EEMD4 2.1 EMD42.2 EEMD分解52.3数值信号分解实验6第三章 球磨机振声信号EEMD分解 93.1 球磨机振声信号采集9 3.2. 信号分解模块63.3 振声信号图形9 3.4振声信号EEMD分解结果 10第四章 振声信号IMF特征提取114.1基于 IMF 频谱特征的集成建模方法124.2 IMF分量选择134.3分解信号的PCA主元分析144.4PCA特征提取方法15第五章 结论和展望16第6章 全文总结13参考文献19致谢语20沈阳化工大学学士学位论文参考文献第一章 绪论1.1引言磨矿分级作业是选矿过程的一个重要环节，而其中对球磨机运行状态的检测与控制，又是磨矿分级作业的核心部分。球磨机运行状态的一个重要参考就是球磨机负荷，球磨机负荷是指球磨机内瞬时的全部装载量，包括新给矿量、循环负荷、加水量和加球个数等。随着自动化技术的发展，很多先进的控制方法被应用于球磨机负荷的检测与控制中，但未取得较好的控制效果。由于线性系统理论已趋于成熟，所以目前最好的球磨机负荷检测与控制方法依然是基于数学模型的控制。最早的球磨机预测模型是 Bond模型;1972年，Wickham在忽略了球磨机负载率和出料的影响前提下提出了球磨机的混合理想模型；1977年,Austin提出了基于动力学的球磨机模型；1987年，Leung提出了撞击和磨损方程用于描述磨矿过程；Man通过对Wickham模型的改进，提出了基于混合理想模型的比例改变建模方法。上述模型的提出极大地推进了对球磨机负荷的研究，但是这些模型都是针对不同情况而建立的，如有的模型用于球磨机负荷的预测，有的模型给出用磨损方程描述球磨机负荷。球磨机负荷可靠检测是选矿过程节能降耗和高效稳定运行的必要条件。磨机负荷的检测通常采用磨机运行时外部响应信号振动及振声间接获得。目前磨机负荷参数检测系统实现方式包括嵌入式仪表或者软件工具包。目前这些开发方式在矿石来料不稳定、多变情况下由于模式固定、适应性较差。为了有效解决球磨机负荷监测问题，目前国际研究热点是通过具有较高灵敏度和较强抗干扰性筒壁振动信号，来监测球磨机运行。充分提取和利用这种非线性非平稳信号蕴含的大量丰富有价值信息，通过特征选择和优化技术提取振动和振声信号集合经验模态分解(EEMD)本征模态函数 IMFs 特征，构建基于IMFs 特征空间的集成模型，并在过程漂移时自适应更新集成模型结构和参数，适应时变非线性磨矿过程复杂工况变化，是保证球磨机负荷集成模型泛化性、可信度和适应性的有效手段。因此，需要开发具有模块 化结构和实验功能的球磨机信号分析和磨机负荷监测软件系统。由于球磨机外部响应信号振动、振声变量时域特征难以提取，频域信息具有较大的冗余成分，国内外提出了基于频谱特征选择的磨机负荷检测方法，并对此进行了大量的研究。目前针对振动、振声信号的处理效果比较好的是多传感器信息融合的磨机负荷集成建模。近年来，混合编程技术运用越来越成熟，可以集合两个或者多个软件的优点进行混合编程，大大提高软件的开发效率。MATLAB 作为一种高级语言，具有极高的编程效率，但它也存在界面开发能力弱、执行速度慢、算法和数据不利于保密以及依赖 MATLAB 编译环境不能独立运行等一些局限性。C#.NET开发的系统界面友好、执行速度快、易维护、可升级等特点，能有脱离编译环境而运行，而且算法和数据保密性强，但是在工程计算方面与 MATLAB 相比就相形见绌了。因此实现 MATLAB与C#的混合编程，结合两者的优点，必将提高软件的开发效率以及软件的性能，也可以工作人员带来很大的方便。针对球磨机负荷难以直接测量问题，利用 MATLAB 与 C#混合编程开发了一套球磨机信号分析和关 键参数预报系统。利用Matlab 的 Deploytool 工具将信号 EEMD 分解、IMF 频谱变换、iPLS 特征选择、 关键参数模型训练和预测一系列 m 函数编译为 C#.NET dll 程序集，供 C#.net 调用，并由 C#.NET 完成前端界面开发和端口操作任务。该系统充分利用混合编程优势互补，提高了编程效率。1.2球磨机噪声机理球磨机噪声主要分以下几种：球磨机运行时筒体与钢球和物料之间碰撞产生的机械噪声；电动机轴承运转时产生的噪声，电动机转子不平衡引起机壳振动产生的声辐射；排粉机的空气动力噪声和通风管道振动噪声等。上述噪声中，应用磨音电耳检测的磨音指的是筒体与钢球和物料之间碰撞产生的机械噪声，其余噪声均属于干扰噪声。磨音本质上属于机械噪声中的撞击噪声，主要包括钢球与钢球、钢球与物料及筒体与钢球之间的撞击产生的撞击噪声。其发声机制包括以下几种：撞击瞬间，物体间的高速流动空气所引起的喷射噪声；撞击瞬间，圆筒、钢球、物料产生的突然变形，导致在该附近激发产生的压力冲击噪声；撞击瞬间，物料破碎形成的外向辐射的压力脉冲噪声；撞击后引起的受撞部件结构共振所激发的结构共振噪声。以上几种发声机制中，结构共振噪声的影响最强，维持时间最长。球磨机噪声的共振结构主体为筒体，当球磨机负荷较小时，磨音主要来源于钢球之间及钢球和衬板的摩擦碰撞，声音较大且频率较高，听起来清脆；当球磨机负荷较大时，磨音主要来源于物料和钢球及物料和衬板之间的摩擦碰撞，声音较小且频率较低，听起来沉闷。第二章 信号分解EEMD2.1EMDEMD，全称为经验模态分解(EMD Empirical Mode Decomposition )算法是Hilbert-Huang变换(HHT)的核心算法。经验模态分解(EMD)算法是通过算法过程定义的，而并非由确定的理论公式定义的，所以对其进行准确的理论分析非常困难，我们目前只能借助大量的数字仿真试验不断对其性能进行深入的研究。 EMD算法的目的在于将性能不好的信号分解为一组性能较好的本征模函数（IMFIntrinsic Mode Function ），且IMF须满足以下两个性质：(1)信号的极值点(极大值或极小值)数目和过零点数目相等或最多相差一个；(2)由局部极大值构成的上包络线和由局部极小值构成的下包络线的平均值为零。EMD算法的计算步骤如下：（1）找出原数据序列X(t)的所有极大值点和极小值点，将其用三次样条函数分别拟合为原序列的上和下包络线；上下包络线的均值为m1；将原数据序列减去m1可得到一个减去低频的新序列h，即h1=X(t)m1；一般h1不一定是平稳数据序列，为此需对它重复上述过程。如h1的包络均值为m11，则去除该包络平均所代表的低频成分后的数据序列为h11，即h11=h1m11重复上述过程，这样就得到第一个本征模函数分量c1，它表示信号数据序列最高频率的成分。（2） 用X(t)减去c1，得到一个去掉高频成分的新数据序列r1；对r1再进行上述分解，得到第二个本征模函数分量c2；如此重复直到最后一个数据序列rn不可被分解，此时，rn代表数据序列X(t)的趋势或均值。在算法中的极值点是指一阶导数为零的点。2.2EEMD分解EEMD是针对EMD方法的不足,提出了一种噪声辅助数据分析方法。EEMD分解原理为:当附加的白噪声均匀分布在整个时频空间时,该时频空间就由滤波器组分割成的不同尺度成分组成。当信号加上均匀分布的白噪声背景时,不同尺度的信号区域将自动映射到与背景白噪声相关的适当尺度上去。当然,每个独立的测试都可能会产生非常嘈杂的结果,这是因为每个附加噪声的成分都包括了信号和附加的白噪声。既然在每个独立的测试中噪声是不同的,当使用足够测试的全体均值时,噪声将会被消除。全体的均值最后将会被认为是真正的结果,随着越来越多的测试，附加的噪声被消除了，唯一持久稳固的部分是信号本身。全称为Ensemble Empirical Mode Decomposition (集合经验模态分解)(Wu and Huang, 2009)，是EMD（经验模分解）(Huang et al. 1998; Huang and Wu, 2008)的改进算法，有效的解决了EMD的混频现象。原始信号通过EEMD算法流程磨音振动信号分解系统2.3数值信号分解实验x=0:0.01:10;Y=sin(x);Nstd=0.2;NE=50;allmode=eemd(Y,Nstd,NE);plot(allmode)EEMD振声数据分解M文件 % EEMD振声数据分解clearclcfilepath=D:ZS_dataD=dir(filepath);kk=length(D);for i=3:kk ifilename=D(i).namefile=strcat(filepath,filename);load(file);dataX=y(:,1);fs=8000; %采样频率for i=10:13nn1=i*fs+1; %第10秒开始取值nn2=nn1+4*fs*60/53-1; X=dataX(nn1:nn2); %第10秒开始 取4个运转周期end%figure(1)%plot(X);title(4个周期振动信号);ylabel(振动幅值);xlabel(采样点数)Xmean=mean(X);Xstd=std(X);outliers=abs(X-Xmean)8*Xstd; %剔除离群点数据nout=sum(outliers)local=find(outliers=1);X(local)=(X(local-1)+X(local+1)/2;%figure(2)%plot(X);title(振动信号预处理);ylabel(振动幅值);xlabel(采样点数)% EEMD分解Nstd=0.2;NE=50;allmodel=eemd(X,Nstd,NE); eval(save D:ZS_IMFZS_IMF_ filename(1:8) allmodel) % figure(3)% subplot(5,1,1);plot(allmodel(:,1);title(IMF1);% subplot(5,1,2);plot(allmodel(:,2);title(IMF2);% subplot(5,1,3);plot(allmodel(:,3);title(IMF3);% subplot(5,1,4);plot(allmodel(:,4);title(IMF4);% subplot(5,1,5);plot(allmodel(:,5);title(IMF5);% figure(4)% subplot(5,1,1);plot(allmodel(:,6);title(IMF6);% subplot(5,1,2);plot(allmodel(:,7);title(IMF7);% subplot(5,1,3);plot(allmodel(:,8);title(IMF8);% subplot(5,1,4);plot(allmodel(:,9);title(IMF9);% subplot(5,1,5);plot(allmodel(:,10);title(IMF10);% figure(5)% subplot(5,1,1);plot(allmodel(:,11);title(IMF11);% subplot(5,1,2);plot(allmodel(:,12);title(IMF12);% subplot(5,1,3);plot(allmodel(:,13);title(IMF13);% subplot(5,1,4);plot(allmodel(:,14);title(IMF14);% subplot(5,1,5);plot(allmodel(:,15);title(IMF15);% figure(6)% subplot(5,1,1);plot(allmodel(:,16);title(IMF16);% subplot(5,1,2);plot(allmodel(:,17);title(IMF17);N0=8000 ; %; 数据长度fs=8000; nfft=N0;window=hann(1024);noverlap=512;range=oneside; % twosideFs=fs;allmodel=allmodel;n_imp,p_imf=size(allmodel); for jj=1:1:n_imp p_data_imf=allmodel(jj,:); for ii=1:1:3 % (Tn-5) %从第8秒起始,到倒数第五秒的数据 nn1=(ii-1)*N0+1; %计算取数起始 1 nn2=nn1+N0-1; %计算取数终止 32768 xn=p_data_imf(nn1:nn2); % 取数个数57962 Pxx_welch,f=pwelch(xn,window,noverlap,nfft,Fs,range); Pall_welch(:,ii)=Pxx_welch ; % 数据存放在A中，按列。不同列表示不同时刻数据。 end mPall_welch,nPall_welch=size(Pall_welch); % 求全部周期的功率值 for j=1:mPall_welch t_welch=mean(Pall_welch(j,:); %求A中每行的均值。每行代表不同的频率点 B_welch(j,1)=t_welch; %均值存放在 B，按行存放。 end ZS_IMF_FFT(:,jj)=B_welch; % 含有全部 16个IMF分量的PSD end %B_welch_all= B_welch_all; eval(save D:ZS_IMFZS_IMF_FFT_ filename(1:8) ZS_IMF_FFT); %imf_NVOICE01频谱数据 end38第3章 . 球磨机振声信号EEMD分解3.1球磨机振声信号采集本次试验所用的球磨机筒体尺寸为460mm460mm（如图2.2所示）。实验数据来源于XMQL 420450格子型球磨机数据采集系统。该球磨机可以随时停机、改变磨机内矿物量和钢球量，实验方便。试验球磨机转速为57rpm. 实验用球磨机内壁上安装有锰钢衬板。实验过程中使用30mm、20mm和15mm三种锰钢钢球。磨机内所装钢球的最大重量为80kg，矿物处理量为10kg/h。实验中处理的矿石为经过破碎的铜矿石，颗粒尺寸小于6mm。球、料和水混合均匀，研磨1分钟。其中振动加速度信号采样频率为51,200 Hz，声音信号采样频率表8KHz。 针对磨机筒体振动和振声信号组成复杂、信息存在冗余性和互补性、现有球磨机负荷测量模型大多 假设振动信号具有线性平稳性，同时传统单模型结构的稳定性和泛化性有限，本文提出基于集合经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)和区间偏最小二乘(Interval Partial Least Squares, iPLS)本征模态频谱特征提取的球磨机负荷参数集成建模方法。筒体振动和振声信号采用非线性、非平稳 EEMD 自适应分解技术获取反映球磨机负荷参数变化的本征模态函数(IMFs)，利用区间偏最小二乘 iPLS 方法去掉本征模态函数频谱中不相关和冗余成分，提取本征模态函数频谱特征，融合时域电流信号，建立磨机负荷参数随机向量泛函连接网络 RVFL(Random Vector Functional Link)集成模型，通过对多个 RVFL 成员模型输出加权平均获得磨矿浓度、填充率和料球比三个球磨机负荷参数估计值。3.2. 信号分解模块 集合经验模态分解是在经验模态分解(EMD)的基础上发展起来的一种信号处理方式，它是EMD方法的改进。它能够自适应的将信号分解为不同时间尺度的IMF集合，每个IMF中只包含一种振荡形式，代表了信号的本征波动模式，这些波动模式的产生与球磨机的负荷密切相关。 基于EEMD的球磨机振动信号分解步骤如下：首先在振动信号 Xv ( t) 中加入振幅为 的白噪声得到： (1)寻找筒体振动信号Yv(t)的极大值点和极小值点；采用样条插值法拟合Yv(t)的上下包络线，计算上下包络线的均值曲线mv1(t) ，得到振动信号Yv(t)与均值曲线的差值hv1(t) (2)hv1(t)为原始信号，重复步骤(2)和(3)得到 (3)反复进行 k 次，直到得到第一个 IMF 分量 hVk (t) ，定义第一个振动信号 IMF 分量 (4)原始振动信号Yv(t)减去第一阶 IMF 分量得到剩余部分 rV1 (t) ，令Yv(t)=rV1 (t)，重复步骤(2)，(3) 和(4)，得到第二阶 IMF 分量 。 重复步骤(2)(4)n 次可得到多个 IMF 分量，当 小于给定的某个值或成为一个单调函数的时候则 整个过程停止。得到第一次筛分结果L (5) 随机产生新的白噪声 ，返回步骤(1)； 将上述步骤(1)(7)重复 N 次，得各阶 IMF 集合 (6) 相应阶的 IMF 平均值作为振动信号该阶的 IMF 结果： (7) 则原振动信号Xv ( t)可表示为 IMF 与振动信号分解残差的形式： (8) 按照上述分解步骤，最终球磨机的振动信号Xv ( t)会被分解为一组时间尺度不同的本征模态函数。 振声信号 XA ( t ) 的 EEMD 分解方法与振动信号基本相同，最终可以得到振声信号的本征模态函数，振声信号可以表示为： (9) 其中 rA ( t ) 为振声信号分解的残差。3.3振 声 数 据 图 像振声信号频谱实验程序：function p_data,Pxx_all_welch_sound=BallMill_WaveFFT(filename)%p_data 振声时域数据数据%Pxx_all_welch_sound振声频谱数据%p_data,Pxx_all_welch_sound=BallMill_WaveFFT(NVOICE01.mat)FileLength=length(filename);FileType=filename(FileLength-2:FileLength);if FileType=wav data,Fs,bits=wavread(filename); % Fs=8000elseif FileType=mat eval(load filename );else disp(data file type not defined)endy=data;y=y(:,1);%单声道 %8000*60ymean=mean(y);yvar=var(y);siglength=length(y);data=y;m,n=size(data);p=128*19; %6784;q=5*125; %1875;dataN= resample(data(:,1),p,q);m,n=size(dataN);mm=floor(32768); % 数据处理点数k=ceil(m/mm); % ceil(x) : 大于x 的最小整数 ,数据的分组数kmm2=m-mm*(k-1); % 数据点的尾数 小于 m1的点数data_ycl=; % 变量名与文件名要一致，便于读取，计算尺寸j=1;for j=1:k if j4*e*delta ; %剔除异常数据 for kkk=1:m1 %补入值 if erastdataN(kkk)=1 if kkk=1 dataNi(kkk)= dataNi(kkk+1); elseif kkk=m1 dataNi(kkk)= dataNi(kkk-1); else dataNi(kkk)= (dataNi(kkk-1)+dataNi(kkk+1)/2; end end end dataNmeanNi=mean(dataNi); %中心化处理 dataNi=dataNi-e*dataNmean; else %e=ones(m2,1); dataNi=dataN(k-1)*mm+1):m),:); % 121-139 m2,n2=size( dataNi); % 39, 6 e=ones(m2,1); dataNmean=mean(dataNi); %求平均 delta=std(dataNi); %求标准方差 erastdataN=abs(dataNi-e*dataNmean)4*e*delta ; %剔除异常数据 for kkkk=1:m2 %补入值 if erastdataN(kkkk)=1 if kkkk=1 dataNi(kkkk)= dataNi(kkkk+1); elseif kkkk=m1 dataNi(kkkk)= dataNi(kkkk-1); else dataNi(kkkk)= (dataNi(kkkk-1)+dataNi(kkkk+1)/2; end end end dataNmeanNi=mean(dataNi); %中心化处理 dataNi=dataNi-e*dataNmean; end data_ycl=data_ycl;dataNi; %endm_ycl,n=size(data_ycl);% *% 以上为对数据进行重采样和剔出，处理周期为32768点% *PP=;FF=;SS=;p_data=data_ycl(1:1:end,1); % 从第1个整周期开始取数据T0=57/60; % 时域信号的采样周期Delta_f=1/T0;% 频率分辩率N0=32768 ; %; 数据长度fs=floor(N0*57/60); %采样频率 28945T=1/fs; %时域信号的采样间隔t=0:T:T0;freq=0:Delta_f:fs;% 设置参数 main for welchnfft=N0;window=hann(1024);noverlap=512;range=oneside; % twosideFs=fs;m10,n10=size(p_data);Tn=floor(m10/N0); %Tn一个周期for ii=8:(Tn-3) %从第8秒起始,到倒数第五秒的数据 nn1=(ii-1)*N0+1; %计算取数起始 1 nn2=nn1+N0-1; %计算取数终止 32768 xn=p_data(nn1:nn2); % 取数个数57962 Pxx_welch,f=pwelch(xn,window,noverlap,nfft,Fs,range); Pall_welch(:,ii-7)=Pxx_welch ; % 数据存放在A中，按列。不同列表示不同时刻数据。 mPall_welch,nPall_welch=size(Pall_welch); % 求全部周期的功率值 for j=1:mPall_welch t_welch=mean(Pall_welch(j,:); %求A中每行的均值。每行代表不同的频率点 B_welch(j,1)=t_welch; %均值存放在 B，按行存放。 end Pxx_all_welch_sound=Pall_welch; Pxx_all_welch_sound(:,nPall_welch+1)=B_welch;endPxx_all_welch_sound=Pxx_all_welch_sound(1:6000,:);figure(1)plot(p_data);xlabel(Time Samples Number);ylabel(Acoustic Signal)title(filename Time Domain Curve);LenData=length(p_data);Xaxis=axis;Xaxis(2)=LenData;axis(Xaxis);saveas(gca,Acoustic_Time_Signal,jpeg)figure(2)Ymax=max(max(Pxx_all_welch_sound(100:6000,:);plot(Pxx_all_welch_sound);xlabel(Frequncy domain (Hz);ylabel(Acoustic Signal Power Spectrum)LenWelch=length(Pxx_all_welch_sound);Xaxis=axis;Xaxis(2)=LenWelch;Xaxis(4)=Ymax;axis(Xaxis);title(filename Frequency Domain Curve);saveas(gca,Acoustic_Spectrum,jpeg)close all force实验结果：3.4 振声信号EEMD分解结果经验模态分解（Empirical Mode Decomposition,EMD）法是美国国家宇航局美籍华人黄锷（N. E. Huang）等人于1998年创造性地提出的一种新型自适应信号时频处理方法，特别适用于非线性非平稳信号的分析处理。经验模态分解7方法根据信号的时间尺度特征进行对信号进行分解并且不需要事先设置基函数。和傅里叶分解与小波分解方法相比是有很大优势的，傅里叶分解与小波分解方法需要建立先验性的谐波基函数和小波基函数。因此，EMD 方法在理论上可以应用于任何类型的信号的分解， 因而在处理非平稳及非线性数据上，具有非常明显的优势,适合于分析非线性、非平稳信号序列，具有很高的信噪比。经验模态分解可以使复杂的信号分解为有限个本征模函数（Intrinsic Mode Function，简称IMF），其中分解出的各个本征模态函数（IMF）中分别包含了原始信号对于不同时间尺度的局部特征信号。在进行EMD分解时，所求解的固有模态函数(IMF)要求满足下述条件：l 要求在整个信号长度上，过零点和极值点的数目必须相等或至多只相差为一点；l 在任意时刻，上包络线和下包络线