轴承matlab处理程序.doc

1、1 数据导入 matlab1.1 启动 Matlab 软件1.2 点击 载入故障数据中的 G2015，Workspace 窗口出现：1.3 取第一组数据 G201，命令窗口输入：G201=G2015(1:1:20000);2. 数据预处理在测试中由数据采集所得的原始信号，在分析前需要进行预处理，以提高数据的可靠性和真实性，并检查信号的随机性，以便正确地选择分析处理方法。预处理工作主要包括三个方面：一是除去信号中的外界干扰信号和剔除异常数据，如趋势项和异点；二是对原始数据进行适当的平滑或拟合；三是对原始信号的特性进行检验。当然这些处理工作不是全部必需的，可以选项或两项内容，当认为原始信号获取工作

2、十分可靠或原始数据简单可以直接判断的情况下，也可以不进行这些预处理工作。以下所做数据预处理，故障轴承以 G201 为例，正常轴承以 Z201 为例，观察原始数据经过不同方法做处理前后的变化。1.1 零均值化处理（原理公式见报告 P8）命令窗口输入：G201l=G201-sum(G201)/20000;%G201l 为零均值处理后的数据。 “20000”为采样点数。sum 为求和语句subplot(2,1,1),plot(G201);subplot(2,1,2),plot(G201l);%显示 G201 与 G201l得到下面图形：从时域图形上看，是波形整体在 Y 轴的平移。再看看频域变化，命令

3、窗口输入：N=20000; %采样点数fs=10000; %采样频率f=(0:N-1)*fs/N; %进行对应的频率转换G201p=abs(fft(G201); %进行 fft 变换，G201p 为 G201 进行 fft 变换后结果G201lp=abs(fft(G201l); %进行 fft 变换，G201lp 为 G201l 进行 fft 变换后结果subplot(2,1,1),plot(f(1:N/2),G201p(1:N/2);subplot(2,1,2),plot(f(1:N/2),G201lp(1:N/2); %显示G201 与 G201p 的频谱图得到下面图形：从频域图可以明显看

4、出，零均值后消除 处出现一个由直流分量产生的大谱峰0（将近达到 ） ，处理后避免了其对周围小峰值产生的负面影响，便于频域分析。4105.1.2 消除趋势项（原理公式见报告 P10）使用最小二乘法，命令窗口输入：t=(0:1/fs:(N-1)/fs); %离散时间列向量G201x=polyfit(t,G201,6); %计算多项式待定系数向量G201x=G201-polyval(G201x,t); %用 G201 减去多项式系数生成的趋势项,G201x 即为消除趋势项后的数据subplot(2,1,1),plot(G201);subplot(2,1,2),plot(G201x);%显示 G201

5、 与 G201x得到以下图形：与前面零均值化处理中做频域图的方法一样，做出 G201 与 G201x 的频谱图 G201p 与G201xp，得到图形如下：从时域图形和频域图形上看，消除趋势项与零均值化处理的功能相似。不过，需要注意的是，它更重要的消除趋势项，因为本数据中的多项式趋势项很小，所以没有明显的变化。1.3 平滑处理（原理公式见报告 P11）使用五点三次平滑，命令窗口输入：a=G201;for k=1:2b(1)=(69*a(1)+4*(a(2)+a(4)-6*a(3)-a(5)/70;b(2)=(2*(a(1)+a(5)+27*a(2)+12*a(3)-8*a(4)/35;for j

6、=3:N-2b(j)=(-3*(a(j-2)+a(j+2)+12*(a(j-1)+a(j+1)+17*a(j)/35;endb(N-1)=(2*(a(N)+a(N-4)+27*a(N-1)+12*a(N-2)-8*a(N-3)/35;b(N)=(69*a(N)+4*(a(N-1)+a(N-3)-6*a(N-2)-a(N-4)/70;a=b;endG201ph=a; %G201ph 为五点三次平滑法处理的数据subplot(2,1,1),plot(G201);subplot(2,1,2),plot(G201ph);%显示 G201 与 G201ph得到以下图形：与前面零均值化处理中做频域图的方法

7、一样，做出 G201 与 G201ph 的频谱图 G201p与 G201php，得到图形如下：从时域图形上看，平滑处理使图形变得平滑，去除毛刺，从频域图形上看，高频部分明显变少变小，而低频部分基本无变化。因为故障的频率主要集中在低中频部分，这样处理后不仅对故障的分析无影响，而且去除部分噪音，减少干扰。1.4 滤波处理（原理公式见报告 P13）%使用巴特沃斯滤波器进行滤波，命令窗口输入：wp=2400; %通带截至频率 2400hzws=2800; %阻带截至频率 2800hzrp=2; %通带波动系数rs=60; %阻带波动系数N,wn =buttord(wp/(fs/2),ws/(fs/2)

8、,rp,rs,z);%建立巴特沃斯滤波器num,den=butter(N,wn);%建立数字滤波器H,W=freqz(num,den);%分析滤波器的幅频特性plot(W*fs/(2*pi),abs(H);grid;%巴特沃斯滤波器频率响应图得到巴特沃斯滤波器频率响应图：继续输入：G201lb=filtfilt(num,den,G201);% G201lb 为 G201 滤波后的数据subplot(2,1,1),plot(G201);subplot(2,1,2),plot(G201lb);%显示 G201 与 G201lb得到以下图形：与前面零均值化处理中做频域图的方法一样，做出 G201 与

9、 G201lb 的频谱图 G201p 与G201lbp，得到图形如下：在时域内不能明显的看出处理前后的区别。但从频域图可以看出，2500Hz 后的频率几乎不存在。因为低通滤波器的通带截至频率为 2400hz，阻带截至频率为 2800hz。可见滤波效果是很好的。以上介绍了一些数据预处理的方法，鉴于本文采集的原始信号数据较好，故只做零均值化这一项处理。3 时域特征值提取（原理公式见 P15）命令窗口输入：G201m=sum(G201l)/20000; %G201m 为均值，G201l 为零均值化处理后结果，下同G201f=sum(G201l-G201m).2); %G201f 为方差G201rms

10、=sqrt(sum(G201l.2)/20000); %G201rms 均方根值G201peak=(max(G201l)-min(G201l)/2; %G201peak 为峰值G201c= G201peak/G201rms; %G201c 为峰值因子G201k=sum(G201l.4)/(G201rms.4)*20000); %G201k 为峭度系数G201s=(G201rms*20000)/sum(abs(G201l); %G201s 为波形因子G201cl=G201peak/(sum(sqrt(abs(G201l)/20000).2; %G201cl 裕度因子G201i=(G201peak

11、*20000)/sum(abs(G201l); %G201i 脉冲因子由此得到 G201 的时域特征值根据前述方法一次得到 G202G2010，Z201Z2010 的时域特征值，建立表格时域特征值状态样本 均值（）710方差均方根值RMS峰值peak峭度系数K峰值因子C裕度因子CL脉冲因子I波形因子SG201 6.8522 2340.80 0.3421 2.2701 13.32346.6357 18.225912.46491.8785G202 22.3452 2605.74 0.3610 2.3549 14.21706.5242 18.821212.62771.9355G203 -32.385

12、4 2902.36 0.3809 2.4886 13.53206.5326 18.932312.63601.9343G204 15.3290 2630.68 0.3627 2.4803 13.87566.8388 19.028812.96441.8957G205 15.8944 2510.69 0.3543 2.3379 13.26326.5985 18.574012.52711.8985G206 -3.7363 2647.01 0.3638 2.3936 13.53826.5793 18.240812.42951.8892G207 -2.0675 2379.66 0.3449 2.2871

13、12.38826.6305 17.562312.11901.8278G208 -4.5102 2548.62 0.3570 2.4512 14.04576.8666 19.819513.28391.9346故障轴承G209 8.0880 2496.80 0.3533 2.3371 12.63046.6146 17.375712.08231.8266G2010 4.4080 2871.86 0.3789 2.3167 11.74176.1138 16.913411.38661.8625Z201 5.2419 1940.06 0.3115 1.5850 4.3203 5.0889 8.1828 6

14、.7397 1.3244Z202 27.7179 1805.18 0.3004 1.5030 4.4684 5.0027 8.0612 6.6351 1.3263Z203 -23.9824 1698.73 0.2914 1.3764 4.6255 4.7228 7.7033 6.3155 1.3372Z204 2.5821 1677.68 0.2896 1.7399 4.8859 6.0073 9.6808 7.9770 1.3279Z205 3.4274 1890.52 0.3075 1.5231 4.5035 4.9539 7.9403 6.5522 1.3226Z206 28.4233

15、1688.29 0.2905 1.3247 3.9282 4.5594 7.2148 5.9647 1.3082Z207 16.6702 1629.54 0.2854 1.4618 4.5880 5.1213 8.2338 6.7921 1.3262Z208 -17.6965 1605.22 0.2833 1.3490 4.4366 4.7618 7.6770 6.3201 1.3272Z209 20.4848 1714.37 0.2928 1.4573 4.6610 4.9774 8.0983 6.6466 1.3354正常轴承Z2010 -4.0320 1790.06 0.2992 1.6

16、877 5.1908 5.6412 9.3795 7.6467 1.3555列出时域参数的数字表后可以简单分析，故障轴承和正常轴承在方差，峰值，峭度系数，裕度因子，脉冲因子，波形因子差别较为明显，而在均值，均方根值，峰值因子差别不明显。4 频域特征值提取（原理公式见 P18）4.1 频域参数命令窗口输入：for i=2:20000G201g(i)=(G201l(i)-G201l(i-1)/(1/10000);endfor i=2:20000G201gg(i)=G201g(i)*G201l(i);endG201msf=(sum(G201g).2)/(4*(pi2)*sum(G201l.2); %

17、G201msf 为均方频率G201fc=(sum(G201gg)/(2*pi*sum(G201l.2); %G201fc 重心频率G201vf=G201msf-G201fc.2; %G201vf 为频率方差由此得到 G201l 的频域参数。根据前述方法一次得到 G202G2010，Z201Z2010 的时域特征值，建立表格频域参数状态样本重心频率 频率方差 均方频率G201 727.4390 768683.0900 1297850.5751G202 732.7359 755254.2402 1292156.1458故障轴 G203 783.7554 830793.7149 1445066.28

18、83G204 772.6559 858228.7434 1455225.8325G205 708.5887 743009.3845 1245107.3676G206 807.0247 876652.1313 1527941.0240G207 776.6448 825242.1517 1428419.3696G208 752.5714 842258.3010 1408621.9914G209 839.0520 952037.6514 1656045.8589承G2010 774.6435 842014.2850 1442086.8989Z201 1947.7204 6177234.3816 99

19、70849.1539Z202 2063.4476 6752047.5896 11009863.5570Z203 2124.7889 6969666.4435 11484394.4028Z204 2129.9819 7038297.6616 11575120.4785Z205 2049.5272 6681185.6000 10881747.3107Z206 2047.9791 6659325.1708 10853543.6293Z207 2151.6996 7045509.3320 11675320.4950Z208 2181.3990 7213367.4335 11971868.8764Z20

20、9 2171.5092 7046245.8915 11761698.1452正常轴承Z2010 2088.4724 6637350.7777 10999067.8211从上表可以看出，频域参数的特征值重复性和差异性都是比较良好的。4.2 傅里叶变换（原理公式见 P20）将 G201l 和 Z201l（Z201l 为 Z201 零均值化后数据）的 fft 变换后的 G201lp 与 Z201lp 做出，程序如下：G201lp=abs(fft(G201l,16384);G201lp=G201lp(1:8192,1);Z201lp=abs(fft(Z201l,16384);Z201lp=Z201lp

21、(1:8192,1);subplot(2,1,1),plot(G201lp); subplot(2,1,2),plot(Z201lp);如下图所示：能够区分两个状态且能代表自己频谱的区域有：点（326，1） 、区域（2560 3000） 、点（3278 ，1 ） 、区域（63106646） 、区域（68507300）用 标记。对故障轴承数据随机抽取 G202fft、G206fft、G207fft、G209fft 数据对比图形如下：从故障轴承抽样数据对比图形可以看出，各个特征值的性重复较好。对正常轴承数据随机抽取 Z203fft、Z204fft、Z206fft、Z208fft 数据对比图形如下：