轴承matlab处理程序文件.docx
《轴承matlab处理程序文件.docx》由会员分享，可在线阅读，更多相关《轴承matlab处理程序文件.docx（33页珍藏版）》请在冰点文库上搜索。
轴承matlab处理程序文件
1．数据导入matlab
1.1启动Matlab软件
1.2点击
载入故障数据中的G2015，Workspace窗口出现：
1.3取第一组数据G201，命令窗口输入：
G201=G2015（1:
1:
20000）;
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轴的平移。
再看看频域变化，命令窗口输入：
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的频谱图
得到下面图形：
从频域图可以明显看出，零均值后消除
处出现一个由直流分量产生的大谱峰（将近达到
），处理后避免了其对周围小峰值产生的负面影响，便于频域分析。
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与G201x
得到以下图形：
与前面零均值化处理中做频域图的方法一样，做出G201与G201x的频谱图G201p与G201xp，得到图形如下：
从时域图形和频域图形上看，消除趋势项与零均值化处理的功能相似。
不过，需要注意的是，它更重要的消除趋势项，因为本数据中的多项式趋势项很小，所以没有明显的变化。
1.3平滑处理（原理公式见报告P11）
使用五点三次平滑，命令窗口输入：
a=G201';
fork=1:
2
b
（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-2
b（j）=（-3*（a（j-2）+a（j+2））+12*（a（j-1）+a（j+1））+17*a（j））/35;
end
b（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;
end
G201ph=a';%G201ph为五点三次平滑法处理的数据
subplot（2,1,1）,plot（G201）;subplot（2,1,2）,plot（G201ph）;%显示G201与G201ph
得到以下图形：
与前面零均值化处理中做频域图的方法一样，做出G201与G201ph的频谱图G201p与G201php，得到图形如下：
从时域图形上看，平滑处理使图形变得平滑，去除毛刺，从频域图形上看，高频部分明显变少变小，而低频部分基本无变化。
因为故障的频率主要集中在低中频部分，这样处理后不仅对故障的分析无影响，而且去除部分噪音，减少干扰。
1.4滤波处理（原理公式见报告P13）
%使用巴特沃斯滤波器进行滤波，命令窗口输入：
wp=2400;%通带截至频率2400hz
ws=2800;%阻带截至频率2800hz
rp=2;%通带波动系数
rs=60;%阻带波动系数
[N,wn]=buttord（wp/（fs/2）,ws/（fs/2）,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与G201lb的频谱图G201p与G201lbp，得到图形如下：
在时域内不能明显的看出处理前后的区别。
但从频域图可以看出，2500Hz后的频率几乎不存在。
因为低通滤波器的通带截至频率为2400hz，阻带截至频率为2800hz。
可见滤波效果是很好的。
以上介绍了一些数据预处理的方法，鉴于本文采集的原始信号数据较好，故只做零均值化这一项处理。
3．时域特征值提取（原理公式见P15）
命令窗口输入：
G201m=sum（G201l）/20000;%G201m为均值，G201l为零均值化处理后结果，下同
G201f=sum（（G201l-G201m）.^2）;%G201f为方差
G201rms=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*20000）/sum（abs（G201l））;%G201i脉冲因子
由此得到G201的时域特征值
根据前述方法一次得到G202~G2010，Z201~Z2010的时域特征值，建立表格
状态
样本
时域特征值
均值（
）
方差
均方根值RMS
峰值peak
峭度系数K
峰值因子C
裕度因子CL
脉冲因子I
波形因子S
故障轴承
G201
6.8522
2340.80
0.3421
2.2701
13.3234
6.6357
18.2259
12.4649
1.8785
G202
22.3452
2605.74
0.3610
2.3549
14.2170
6.5242
18.8212
12.6277
1.9355
G203
-32.3854
2902.36
0.3809
2.4886
13.5320
6.5326
18.9323
12.6360
1.9343
G204
15.3290
2630.68
0.3627
2.4803
13.8756
6.8388
19.0288
12.9644
1.8957
G205
15.8944
2510.69
0.3543
2.3379
13.2632
6.5985
18.5740
12.5271
1.8985
G206
-3.7363
2647.01
0.3638
2.3936
13.5382
6.5793
18.2408
12.4295
1.8892
G207
-2.0675
2379.66
0.3449
2.2871
12.3882
6.6305
17.5623
12.1190
1.8278
G208
-4.5102
2548.62
0.3570
2.4512
14.0457
6.8666
19.8195
13.2839
1.9346
G209
8.0880
2496.80
0.3533
2.3371
12.6304
6.6146
17.3757
12.0823
1.8266
G2010
4.4080
2871.86
0.3789
2.3167
11.7417
6.1138
16.9134
11.3866
1.8625
正常轴承
Z201
5.2419
1940.06
0.3115
1.5850
4.3203
5.0889
8.1828
6.7397
1.3244
Z202
27.7179
1805.18
0.3004
1.5030
4.4684
5.0027
8.0612
6.6351
1.3263
Z203
-23.9824
1698.73
0.2914
1.3764
4.6255
4.7228
7.7033
6.3155
1.3372
Z204
2.5821
1677.68
0.2896
1.7399
4.8859
6.0073
9.6808
7.9770
1.3279
Z205
3.4274
1890.52
0.3075
1.5231
4.5035
4.9539
7.9403
6.5522
1.3226
Z206
28.4233
1688.29
0.2905
1.3247
3.9282
4.5594
7.2148
5.9647
1.3082
Z207
16.6702
1629.54
0.2854
1.4618
4.5880
5.1213
8.2338
6.7921
1.3262
Z208
-17.6965
1605.22
0.2833
1.3490
4.4366
4.7618
7.6770
6.3201
1.3272
Z209
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.6877
5.1908
5.6412
9.3795
7.6467
1.3555
列出时域参数的数字表后可以简单分析，故障轴承和正常轴承在方差，峰值，峭度系数，裕度因子，脉冲因子，波形因子差别较为明显，而在均值，均方根值，峰值因子差别不明显。
4．频域特征值提取（原理公式见P18）
4.1频域参数
命令窗口输入：
fori=2:
20000
G201g（i）=（G201l（i）-G201l（i-1））/（1/10000）;
end
fori=2:
20000
G201gg（i）=G201g（i）*G201l（i）;
end
G201msf=（sum（（G201g）.^2））/（4*（pi^2）*sum（G201l.^2））;%G201msf为均方频率
G201fc=（sum（G201gg））/（2*pi*sum（G201l.^2））;%G201fc重心频率
G201vf=G201msf-G201fc.^2;%G201vf为频率方差
由此得到G201l的频域参数。
根据前述方法一次得到G202~G2010，Z201~Z2010的时域特征值，建立表格
状态
样本
频域参数
重心频率
频率方差
均方频率
故障轴承
G201
727.4390
768683.0900
1297850.5751
G202
732.7359
755254.2402
1292156.1458
G203
783.7554
830793.7149
1445066.2883
G204
772.6559
858228.7434
1455225.8325
G205
708.5887
743009.3845
1245107.3676
G206
807.0247
876652.1313
1527941.0240
G207
776.6448
825242.1517
1428419.3696
G208
752.5714
842258.3010
1408621.9914
G209
839.0520
952037.6514
1656045.8589
G2010
774.6435
842014.2850
1442086.8989
正常轴承
Z201
1947.7204
6177234.3816
9970849.1539
Z202
2063.4476
6752047.5896
11009863.5570
Z203
2124.7889
6969666.4435
11484394.4028
Z204
2129.9819
7038297.6616
11575120.4785
Z205
2049.5272
6681185.6000
10881747.3107
Z206
2047.9791
6659325.1708
10853543.6293
Z207
2151.6996
7045509.3320
11675320.4950
Z208
2181.3990
7213367.4335
11971868.8764
Z209
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（1:
8192,1）;
subplot（2,1,1）,plot（G201lp）;subplot（2,1,2）,plot（Z201lp）;
如下图所示：
能够区分两个状态且能代表自己频谱的区域有：
点（326，1）、区域（2560～3000）、点（3278，1）、区域（6310~6646）、区域（6850~7300）用
标记。
对故障轴承数据随机抽取G202fft、G206fft、G207fft、G209fft数据对比图形如下：
从故障轴承抽样数据对比图形可以看出，各个特征值的性重复较好。
对正常轴承数据随机抽取Z203fft、Z204fft、Z206fft、Z208fft数据对比图形如下：
图3-12正常轴承重复性FFT谱
从正常轴承抽样数据对比图形可以看出，各个特征值的性重复较好。
利用以下程序将G201~G2010，Z201~Z2010的傅里叶变换特征值与特征区域提取出来：
G201ffttz=[G201lp（326,1）,sum（G201lp（2560:
3000,1））,G201lp（3278,1）,sum（G201lp（6310:
6646,1））,sum（G201lp（6850:
7300,1））];%G201ffttz为G201l进行fft变换后提取的特征值
建立表格：
状态
样本
FFT频域特征值
（326，1）
（2560：
3000，1）
（3278，1）
（6310：
6646，1）
（6850：
7300，1）
故障轴承
G201
24.0370
16457.4591
223.4111
5511.6023
5070.6419
G202
36.1318
15337.1272
216.3505
5283.4564
5046.5340
G203
28.2671
18178.6398
213.6201
5642.9205
5046.5340
G204
11.1246
18178.6398
163.9869
5662.4984
5687.7720
G205
114.0399
16154.9008
168.2039
5562.4298
5151.0799
G206
36.9351
19195.7369
179.8281
5528.2613
6344.3349
G207
11.4629
18551.4539
162.5642
5491.6241
5443.5653
G208
141.1385
15024.3193
148.2590
5700.0494
5710.0924
G209
37.4586
20519.3069
148.7093
5955.7815
6390.3742
G2010
92.7402
20333.6468
138.5370
6092.2449
5809.4190
正常轴承
Z201
247.1716
10356.3888
8.2110
15539.7691
11057.7255
Z202
179.1598
10083.4336
27.4943
14772.3720
10632.8833
Z203
187.5852
10018.7181
6.2269
15396.4913
10692.6573
Z204
212.9235
9502.5877
17.2084
15035.7169
10475.6605
Z205
132.3277
10125.0714
18.6699
15691.3139
11033.8520
Z206
210.8356
10029.6034
17.6392
15133.5896
10612.5788
Z207
205.1529
10220.0559
18.7119
14764.7557
10577.8874
Z208
227.1720
10414.9020
20.9440
14186.3542
11097.2232
Z209
173.4568
10161.3688
10.4170
15232.0201
10877.8867
Z2010
198.6173
10366.7437
6.8624
15144.1329
11205.6632
4.3功率谱处理（（原理公式见P22）
采用Welch平均周期法，采样频率为10000Hz，长度为16384点，分段时每段长度为4096，相邻两段重叠的点数为2048，因此分成了7段，窗函数为缺省。
命令窗口输入以下程序：
[p,f]=spectrum（G201l,4096,2048,[],fs）;
G201gl=p（:
1）;%采用Welch平均周期法，G201功率谱处理结果
将G201gl和Z201gl显示出来，得到：
能够区分两个状态且能代表自己频谱的区域有：
点（82，1）、区域（660～739）、点（820，1）、点（1473,1）、点（1616,1）、点（1639，1）点（1778，1）。
对故障轴承数据随机抽取G202gl、G204gl、G206gl、G208gl数据对比图形如下：
图3-14故障轴承重复性功率谱
从故障轴承抽样数据对比图形可以看出，各个特征值的性重复较好。
对正常轴承数据随机抽取Z203gl、Z204gl、Z209gl、Z2010gl数据对比图形如下：
图3-15正常轴承重复性功率谱
从正常轴承抽样数据对比图形可以看出，各个特征值的性重复较好。
将G201~G2010，Z201~Z2010的傅里叶变换特征值与特征区域提取出来，建立表格：
状态
样本
Welch平均周期法功率谱特征值
（82，1）
（820，1）
（1473,1）
（1616,1）
（1639，1）
（1778，1）
（660：
739，1）
故障轴承
G201
0.3643
0.5665
0.2145
0.0179
0.0431
0.0087
8.8665
G202
0.2351
0.5068
0.1865
0.0246
0.0441
0.0128
10.8937
G203
0.2684
0.5291
0.1229
0.0115
0.0353
0.0110
12.0233
G204
0.2700
0.3499
0.0768
0.0156
0.0660
0.0080
10.2235
G205
0.4258
0.3692
0.1348
0.0179
0.0452
0.0091
9.7017
G206
0.2295
0.4123
0.1189
0.0126
0.0869
0.0147
14.1895
G207
0.1516
0.3213
0.0227
0.0163
0.0439
0.0065
11.9705
G208
0.2976
0.2878
0.0474
0.0130
0.0517
0.0074
7.6765
G209
0.2195
0.3290
0.1001
0.0217
0.0319
0.0094
15.106
G2010
0.3630
0.2887
0.1571
0.0251
0.0468
0.0114
13.565
正常轴承
Z201
0.6538
0.0091
0.0217
1.2296
1.3105
0.2427
3.8739
Z202
0.3668
0.0187
0.0117
1.4303
2.3120
0.3131
3.5024
Z203
0.4774
0.0068
0.0097
1.2764
1.4431
0.2682
3.4486
Z204
0.5658
0.0097
0.0160
1.4
升级会员 