轴承matlab处理程序

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

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

3、样点数fs=10000;%采样频率f=(0:N-1)*fs/N;%进行对应的频率转换G201p=abs(fft(G201);%进行fft变换，G201P为G201进行fft变换后结果G2011P=abs(fft(G201l);%进行fft变换，G2011P为G2011进行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.5x104),处理后避免了其对周围小峰值产生的负面影响，便于频域分析。消除趋势项(原理公式见报告P

4、10)使用最小二乘法,命令窗口输入：t=(0:1/fs:(N-1)/fs);%离散时间列向量G201x=Po1yfit(t,G201,6);%计算多项式待定系数向量G201x=G201-po1yva1(G201x,t);%用G201减去多项式系数生成的趋势项,G201x即为消除趋势项后的数据subp1ot(2,1,1),p1ot(G201);subp1ot(2,1,2),p1ot(G201x);%显示G201与G201x得到以下图形：与前面零均值化处理中做频域图的方法一样，做出G201与G201x的频谱图G201P与G201xp,得到图形如下：从时域图形和频域图形上看，消除趋势项与零均值化处理

5、的功能相似。不过，需要注意的是，它更重要的消除趋势项，因为本数据中的多项式趋势项很小，所以没有明显的变化。1.3平滑处理(原理公式见报告P11)使用五点三次平滑，命令窗口输入：a=G201;fork=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;forj=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

6、*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得到以下图形：与前面零均值化处理中做频域图的方法一样，做出G201与G201Ph的频谱图G201p与G201Php，得到图形如下：从时域图形上看，平滑处理使图形变得平滑，去除毛刺，从频域图形上看，高频部分明显变少变小，而低频部分基本无变化。因为故障的频率主要集中在低中频

7、部分，这样处理后不仅对故障的分析无影响，而且去除部分噪音，减少干扰。1.4滤波处理(原理公式见报告P13)%使用巴特沃斯滤波器进行滤波，命令窗口输入：wp=2400;%通带截至频率2400hzws=2800;%阻带截至频率2800hzrp=2;%通带波动系数rs=60;%阻带波动系数N,wn=buttord(wp/(fs/2),ws/(fs/2)jpjs,z);%建立巴特沃斯滤波器num,den=butter(N,wn);%建立数字滤波器H,W=freqz(num,den);%分析滤波器的幅频特性plot(W*fs/(2*pi),abs(H);grid;%巴特沃斯滤波器频率响应图得到巴特沃斯滤

8、波器频率响应图：继续输入：G2011b=filtfilt(num,den,G201);%G2011b为G201滤波后的数据subp1ot(2,1,1),p1ot(G201);subp1ot(2,1,2),p1ot(G2011b);%显示G201与G2011b得到以下图形：与前面零均值化处理中做频域图的方法一样，做出G201与G2011b的频谱图G201P与G2011bp，得至U图形如下：在时域内不能明显的看出处理前后的区别。但从频域图可以看出，2500Hz后的频率几乎不存在。因为低通滤波器的通带截至频率为2400hz,阻带截至频率为2800hz。可见滤波效果是很好的。以上介绍了一些数据预处理的

9、方法，鉴于本文采集的原始信号数据较好，故只做零均值化这一项处理。.时域特征值提取（原理公式见P15）命令窗口输入：G201m=sum(G201l)/20000;%G201m为均值，G2011为零均值化处理后结果，下同G201f=sum(G201l-G201m).A2);%G201f为方差G201rms=sqrt(sum(G201l.A2)/20000);%G201rms均方根值G201peak=(max(G201l)-min(G201l)/2;%G201peak为峰值G201c=G201peak/G201rms;%G201c为峰值因子G201k=sum(G201l.A4)/(G201rms.A

10、4)*20000);%G201k为峭度系数G201s=(G201rms*20000)/sum(abs(G201l);%G201s为波形因子G201cl=G201peak/(sum(sqrt(abs(G201l)/20000).A2;%G201cl裕度因子G201i=(G201peak*20000)/sum(abs(G201l);%G201i脉冲因子由此得到G201的时域特征值根据前述方法一次得到G202G2010，Z201Z2010的时域特征值，建立表格状态样本时域特征值均值(10)方差均方根值RMS峰值peak峭度系数K峰值因子C裕度因子CL脉冲因子I波形因子SG2016.85222340.

11、800.34212.270113.3236.635718.22512.4641.8785499G20222.34522605.740.36102.354914.2176.524218.82112.6271.9355027G203-32.38542902.360.38092.488613.5326.532618.93212.6361.9343030G20415.32902630.680.36272.480313.8756.838819.02812.9641.8957684故G20515.89442510.690.35432.337913.2636.598518.574012.5271.8985障

12、21轴G206-3.73632647.010.36382.393613.5386.579318.24012.4291.8892承285G207-2.06752379.660.34492.287112.3886.630517.56212.1191.8278230G208-4.51022548.620.35702.451214.0456.866619.81913.2831.9346759G2098.08802496.800.35332.337112.6306.614617.375712.0821.826643G20104.40802871.860.37892.316711.7416.113816.

13、91311.3861.8625746Z2015.24191940.060.31151.58504.32035.08898.18286.73971.3244Z20227.71791805.180.30041.50304.46845.00278.06126.63511.3263Z203-23.98241698.730.29141.37644.62554.72287.70336.31551.3372Z2042.58211677.680.28961.73994.88596.00739.68087.97701.3279正Z2053.42741890.520.30751.52314.50354.95397

14、.94036.55221.3226常Z20628.42331688.290.29051.32473.92824.55947.21485.96471.3082轴Z20716.67021629.540.28541.46184.58805.12138.23386.79211.3262承Z208-17.69651605.220.28331.34904.43664.76187.67706.32011.3272Z20920.48481714.370.29281.45734.66104.97748.09836.64661.3354Z2010-4.03201790.060.29921.68775.19085.

15、64129.37957.64671.3555列出时域参数的数字表后可以简单分析，故障轴承和正常轴承在方差，峰值，峭度系数，裕度因子，脉冲因子，波形因子差别较为明显，而在均值，均方根值，峰值因子差别不明显。.频域特征值提取(原理公式见P18)频域参数命令窗口输入：fori=2:20000G201g(i)=(G201l(i)-G201l(i-1)/(1/10000);endfori=2:20000G201gg(i)=G201g(i)*G201l(i);endG201msf=(sum(G201g).A2)/(4*(piA2)*sum(G201l.A2);%G201msf为均方频率G201fc=(su

16、m(G201gg)/(2*pi*sum(G201lC2);%G201fc重心频率G201vf=G201msf-G201fc.A2;%G201vf为频率方差由此得到G2011的频域参数。根据前述方法一次得到G202G2010,Z201Z2010的时域特征值，建立表格状态样本频域参数重心频率频率方差均方频率故G201727.4390768683.09001297850.5751障G202732.7359755254.24021292156.1458轴G203783.7554830793.71491445066.2883承G204772.6559858228.74341455225.8325G205

17、708.5887743009.38451245107.3676G206807.02478766520240G207776.64488252423696G208752.5714842258.30101408621.9914G209839.0520952037.65141656045.8589G2010774.6435842014.28501442086.8989Z2011947.72046177234.38169970849.1539Z2022063.44766752047.589611009863.5570Z2032124.788969696

18、66.443511484394.4028Z2042129.98197038297.661611575120.4785正Z2052049.52726681185.600010881747.3107常Z2062047.97916659325.170810853543.6293轴Z2072151.69967045509.332011675320.4950承Z2082181.39907213367.433511971868.8764Z2092171.50927046245.891511761698.1452Z20102088.47246637350.777710999067.8211从上表可以看出，频

19、域参数的特征值重复性和差异性都是比较良好的。4.2傅里叶变换（原理公式见P20）将G2011和Z2011（Z2011为Z201零均值化后数据）的fft变换后的G2011P与Z2011P做出，程序如下：G201lp=abs(fft(G201l,16384);G201lp=G201lp(1:8192,1);Z201lp=abs(fft(Z201l,16384);Z201lp=Z201lp(1:8192,1);subplot(2,1,1),plot(G201lp);subplot(2,1,2),plot(Z201lp);如下图所示：能够区分两个状态且能代表自己频谱的区域有：点（326,1）、区域（2

20、5603000）、点（3278,）用区域（63KT6646）、区域（68507300）用变1标记。从故障轴承抽样数据对比图形可以看出,各个特征值的性重复较好。G201ffttz=G201lp(326,1),sum(G201lp(2560:3000,1),G201lp(3278,1),sum(G201lp(6310:6646,1),sum(G2011P(6850:7300,1);%G201ffttz为G2011进行fft变换后提取的特征值建立表格：状态样本FFT频域特征值(326,1)(2560：3000,1)(3278,1)(6310：6646,1)(6850：7300,1)故G20124.0

21、37016457.4591223.41115511.60235070.6419障G20236.131815337.1272216.35055283.45645046.5340轴G20328.267118178.6398213.62015642.92055046.5340承G20411.124618178.6398163.98695662.49845687.7720G205114.039916154.9008168.20395562.42985151.0799G20636.935119195.7369179.82815528.26136344.3349G20711.462918551.45391

22、62.56425491.62415443.5653G208141.138515024.3193148.25905700.04945710.0924G20937.458620519.3069148.70935955.78156390.3742G201092.740220333.6468138.53706092.24495809.4190Z201247.171610356.38888.211015539.769111057.7255Z202179.159810083.433627.494314772.372010632.8833Z203187.585210018.71816.226915396.4

23、91310692.6573正Z204212.92359502.587717.208415035.716910475.6605常Z205132.327710125.071418.669915691.313911033.8520轴Z206210.835610029.603417.639215133.589610612.5788承Z207205.152910220.055918.711914764.755710577.8874Z208227.172010414.902020.944014186.354211097.2232Z209173.456810161.368810.417015232.0201

24、10877.8867Z2010198.617310366.74376.862415144.132911205.66324.3功率谱处理(原理公式见P22)采用Welch平均周期法，采样频率为10000Hz,长度为16384点，分段时每段长度为4096,相邻两段重叠的点数为2048，因此分成了7段，窗函数为缺省。命令窗口输入以下程序：p,f=spectrum(G201l,4096,2048,fs);G201gl=p(:,1);%采用Welch平均周期法，G201功率谱处理结果将G201gl和Z201gl显示出来，得到：能够区分两个状态且能代表自己频谱的区域有：点（82,1）、区域（660739）

25、、点（820,1）、点（1473,1）、点（1616,1）、点（1639，1）点（1778，1）。从故障轴承抽样数据对比图形可以看出,各个特征值的性重复较好。对正常轴承数据随机抽取Z203gl、Z204gl、Z209gl、Z2010gl数据对比图形如下:KWP忙uisn_.T一-值北m-*皿T：iSK.I9N-XW9：_V:1iniKiri?向VUXiU-一LJCEJ|YO.5ffi6_-1一.”EE：E-.KlAtfiV:1*1?rL-Y:1BTI_i_1-L-JKH：UN出；,产Sbs_L-.-i.1411yliKimtV:1皿.1,KKBy：l.05，JC:1E5we图3-15正常轴承重

26、复性功率谱从正常轴承抽样数据对比图形可以看出，各个特征值的性重复较好。将G201G2010,Z201Z2010的傅里叶变换特征值与特征区域提取出来，建立表格：状Welch平均周期法功率谱特征值样本态（82，1）（820，1）（1473,1）（1616,1）（1639，1）（1778，1）（660：739，1）G2010.36430.56650.21450.01790.04310.00878.8665G2020.23510.50680.18650.02460.04410.012810.8937G2030.26840.52910.12290.01150.03530.011012.0233故G204

27、0.27000.34990.07680.01560.06600.008010.2235障G2050.42580.36920.13480.01790.04520.00919.7017轴G2060.22950.41230.11890.01260.08690.014714.1895承G2070.15160.32130.02270.01630.04390.006511.9705G2080.29760.28780.04740.01300.05170.00747.6765G2090.21950.32900.10010.02170.03190.009415.106G20100.36300.28870.157

28、10.02510.04680.011413.565Z2010.65380.00910.02171.22961.31050.24273.8739Z2020.36680.01870.01171.43032.31200.31313.5024Z2030.47740.00680.00971.27641.44310.26823.4486正Z2040.56580.00970.01601.44681.97070.20482.9251常Z2050.46120.01110.00950.42011.42860.22223.3279轴Z2060.59870.00920.02771.23900.94610.27473.

29、2920承Z2070.43350.01320.01340.89882.08890.19683.3916Z2080.55980.01360.01630.54911.89510.16723.5790Z2090.41640.00950.02321.02071.47640.19303.5333Z20100.39260.01640.01631.28561.42980.18733.8456从上表可以看出，故障轴承和正常轴承功率谱的特征值重复性和差异性都是比较良好的。.时频分析法一一小波包络(原理公式见P24)500)H75)H:(卜500)H75)H:(卜5加口加在这里选择(30)频段的信号为例进行重构，

30、再包络解调，程序如下：wpt=wpdec(G201l3db4);%小波包进行3层分解，分解用的小波基为db4cz=wprcoef(wpt,30);%对(30)节点即(0875Hz)进行重构G201xb=abs(hilbert(cz);%G201xb为G2011的包络分析G201xbp=abs(fft(G201xb,4096);%G201xbp为G201xb进行FFT变换后数据故障轴承G202与正常轴承Z202的(30)节点的包络谱：图3-16节点（30）G202与Z202小波包包络谱对比发现两者并没有明显的区别，很难从此节点提取出特征值，说明故障频率不在0875Hz频率段。再对（31）节点的样

31、本1数据进行分析：图3-17节点（31）G202与Z202小波包包络谱从图中可以得出在点（2，1）、点（513,1）和点（1025，1）处的值个状态的差异性比较好，因此初步确定这三点为其特征值。接下来对（32）、（33）几个节点进行了包络谱分析：图3-18节点(32)G202与Z202小波包包络谱图3-18节点(33)G202与Z202小波包包络谱在众多的包络谱分析中发现其规律是：他们都是在包络谱的(2,1)、(513,1)、(1025,1)点出现峰值并且各个状态有区别，即差别性较好，并且通过重复性检验，其可以作为特征值。因此将这三点作为特征值，并取(31)、(32)、(33)三个节点作为此种

32、方法分析的特征值提取节点。特征值提取的程序如下：wpt=wpdec(G201l,3,db4);fork=1:3cz=wprcoef(wpt,3k);hom=abs(hilbert(cz);czf=abs(fft(hom,4096);G201xbp(1,3*(k)-2)=czf(2,1);G201xbp(1,3*(k)-1)=czf(513,1);G201xbp(1,3*(k)=czf(1025,1);end得到的G201xbp为4096*9的矩阵，第一行就是所求的各节点3个点的特征值。12345678911140.82375.985610,479257,2915I9.4964口.913662,

33、201524,80815.8094依此建立如下表格：状态样本小波包包络谱特征值（31）（32）（33）（2，1）(513,1)(1025,1)（2，1）(513,1)(1025,1)（2，1）(513,1)(1025,1)G201140.82375.985510.479257.291519.45640.915652.201524.80815.8094G202146.3216.16529.797766.867420.45920.796290.958228.26856.26617G203168.2919.481512.902764.289221.36120.6669103.759931.21457

34、.3365故G204187.25917.203512.036855.834720.15920.754997.767229.56787.1333障G205120.15774.97749.117938.364721.64820.663064.718724.17255.8821轴G206126.1719.998613.977838.176219.99110.629763.416833.85788.0195承7G207172.10008.891213.615048.372116.89800.834875.363429.52816.8948G208129.8306.317313.272052.23101

35、9.16160.765483.543128.15307.92583G209116.45794.113611.304852.985220.01131.237453.141328.90417.5804G2010135.15729.374115.721574.642928.54920.9465107.163344.01199.8800Z20146.381212.393721.06194.710110.46040.295412.422318.86958.1549Z20223.954111.761418.22852.190710.23340.379916.762117.05978.3838Z20313.

36、48868.445916.57215.43278.06850.345113.014517.45928.2029正Z20493.863514.109119.01434.954610.50310.452517.713519.01018.4352常Z20539.649412.643519.611511.75429.97510.350420.583019.69528.4686轴Z20638.831911.601818.92456.05379.57180.346915.048718.10158.2054承Z20727.221112.947618.80511.98989.88270.462919.4168

37、16.95098.4433Z20813.718510.653916.82794.08689.21350.36089.594918.36949.2244Z20963.891213.132218.80909.11319.85580.39784.604518.66667.9271Z201019.222713.034419.74573.828710.06120.353624.758616.12317.9186.特征值归一化(原理公式见报告P29)程序如下：fori=1:33forj=1:20gy(i,j)=(tz(i,j)-min(tz(i,：)/(max(tz(i,：)-min(tz(i,：);%t

38、z为原特征值矩阵，gy为归一化后的特征值矩阵。33为原始特征值个数，20为数据样本个数。endend得到下列表格：故障轴承G201G20100.64530.900000.78470.79400.47110.49860.45840.66560.60510.56710.771310.79060.69810.80310.59700.72730.68730.9765时0.60250.796110.81350.72750.82480.63110.75510.71720.9795域0.81230.885110.99290.87050.91840.82690.96790.86980.8523特0.91311

39、0.93340.96680.90730.93400.82230.98340.84580.7594征0.89990.85160.85520.98800.88380.87550.897710.89080.6737值0.87360.92080.92960.93730.90120.87480.820910.80610.76940.88810.91030.91150.95630.89660.88330.840810.83580.74080.909110.99810.93660.94100.92620.82830.99860.82640.8836频域0.01280.01640.05100.043500.0

40、6680.04620.02990.08860.0448参0.00400.00190.01360.017800.02070.01270.01530.03230.0153数0.00490.00440.01860.019600.02640.01710.01520.03830.0184FFT谱0.05470.10590.072600.43600.10930.00140.55080.11160.3458特0.63130.52960.78750.78750.60380.87980.82140.501210.9831征10.96750.95490.72640.74580.79930.71980.65400.

41、65600.6092值0.021900.03450.03640.02680.02350.02000.04000.06460.07770.0039000.10410.01700.21070.06450.10770.21820.12390.42350.16630.23260.23580.54600.155100.29070.13520.4209功10.89330.93320.61300.64750.72450.56190.50210.57570.5037率10.86340.55320.32830.61120.53370.06440.18490.44200.72谱0.00450.009100.002

42、90.00450.00080.00330.00100.00710.009特0.00490.00540.00150.01500.00580.02410.00530.008700.0065征0.00720.02050.01470.00490.00850.026700.00290.00950.0160值0.48780.65420.74690.59920.55630.92480.74260.390110.87350.73280.76440.890810.61390.64850.91280.66950.59260.7002小0.18730.20530.53700.30910.08640.58880.47

43、800.220500.5262波0.11400.05690.31690.244400.40690.37650.34780.18310.5529包0.76120.89300.85750.74110.50070.49810.63840.69150.70191络0.55600.60500.64900.59030.66300.58210.43110.54160.58311谱0.65840.53160.39440.48780.39020.35490.57270.499010.6912特0.46410.84200.96680.90840.58610.57340.68990.76970.47331征0.31

44、140.43550.54110.48210.28860.63590.48070.43140.45831值00.11220.37500.32520.01790.54290.26660.51990.43511正常轴承Z201Z20100.61880.98840.13820.57500.588910.80670.24160.86950.46630.25810.15420.07210.05590.21990.06400.018700.08410.1425时0.28890.17520.08300.06450.24800.07380.021500.09730.1629域0.22360.15320.0444

45、0.35670.170500.11780.02090.11390.3119特0.03810.05250.06780.09310.055900.0640.04940.07120.1227征0.22950.19210.07080.62760.171000.24350.08770.18120.4689值0.07680.06710.03880.19560.057600.08080.03670.07010.17170.10590.09160.04790.27490.080300.11300.04860.09320.22980.02580.02890.04620.03140.023000.02870.03

46、030.04340.0754频域0.84130.91990.96160.96510.91050.90940.979810.99330.9369参0.83990.92870.96230.97290.91780.91440.974110.97420.9110数0.81350.91030.95460.96300.89840.89570.972410.98040.9093FFT谱10.71190.74760.85490.51350.84610.82200.91530.68770.7943特0.07750.05270.046800.05650.04780.06510.08280.05980.0784征0

47、.00910.097900.05060.05730.05250.05750.06780.01930.0029值0.98540.91170.97170.937010.94640.91100.85540.95570.94740.97600.90700.91670.88150.97210.90370.89810.98240.94681功10.42850.64870.82480.61650.89030.56130.81280.52730.4799率0.00410.021300.00520.00770.00430.01140.01210.00480.0172谱0.05950.01070.00100.03

48、1700.08880.01900.03320.06680.0332特0.84870.9890.881310.28470.85520.61820.37460.70310.8877征0.560810.61890.85030.61260.40090.90220.81720.63350.6131值0.770410.85360.64680.70350.87480.62070.52410.60830.58970.07790.0470.043000.03310.03010.03830.05370.04990.07560.18930.060200.46250.15050.14580.07900.00130.2

49、9010.0330小0.82840.76510.433410.85340.74910.88380.65430.90230.8925波10.76280.62410.82860.87860.82100.81110.64550.81140.8898包0.03740.00280.04740.04080.13440.055900.02890.09800.0253络0.11680.105700.11890.09310.07340.08860.05590.08730.0973谱00.08970.05280.16680.05840.05470.17780.06940.10870.0618特0.07620.11850.08200.12780.15580.10180.14440.048700.1965征0.09850.03360.04790.10350.12810.07090.02970.08050.09120值0.57620.63240.58800.64510.65330.58860.64720.83890.52020.51827.建立BP神经网