轴承matlab处理程序文档格式.docx
- 文档编号:7256712
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:29
- 大小:622.06KB
轴承matlab处理程序文档格式.docx
《轴承matlab处理程序文档格式.docx》由会员分享,可在线阅读,更多相关《轴承matlab处理程序文档格式.docx(29页珍藏版)》请在冰点文库上搜索。
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:
%显示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,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,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,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
4
7
5
G202
70
G203
6
20
60
3
G204
G205
9
40
1
G206
8
G207
23
G208
G209
G2010
正常轴承
Z201
50
Z202
30
Z203
Z204
Z205
Z206
Z207
Z208
Z209
Z2010
列出时域参数的数字表后可以简单分析,故障轴承和正常轴承在方差,峰值,峭度系数,裕度因子,脉冲因子,波形因子差异较为明显,而在均值,均方根值,峰值因子差异不明显。
4.频域特征值提取〔原理公式见P18〕
4.1频域参数
fori=2:
20000
G201g(i)=(G201l(i)-G201l(i-1))/(1/10000);
G201gg(i)=G201g(i)*G201l(i);
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的频域参数。
频域参数
重心频率
频率方差
均方频率
从上表可以看出,频域参数的特征值重复性和差异性都是比较良好的。
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:
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〕
90
10
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〕
从上表可以看出,故障轴承和正常轴承功率谱的特征值重复性和差异性都是比较良好的。
5.时频分析法——小波包络〔原理公式见P24〕
在这里选择(30)频段的信号为例进行重构,再包络解调,程序如下:
wpt=wpdec(G201l,3,'
db4'
%小波包进行3层分解,分解用的小波基为db4
cz=wprcoef(wpt,[30]);
%对〔30〕节点即〔0~875Hz〕进行重构
G201xb=abs(hilbert(cz));
%G201xb为G201l的包络分析
G201xbp=abs(fft(G201xb,4096));
%G201xbp为G201xb进行FFT变换后数据
故障轴承G202与正常轴承Z202的〔30〕节点的包络谱:
图3-16节点〔30〕G202与Z202小波包包络谱
比照发现两者并没有明显的区别,很难从此节点提取出特征值,说明故障频率不在0~875Hz频率段。
再对(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)三个节点作为此种方法分析的特征值提取节点。
特征值提取的程序如下:
cz=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);
得到的G201xbp为4096*9的矩阵,第一行就是所求的各节点3个点的特征值。
依此建立如下表格:
小波包包络谱特征值
〔31〕
〔32〕
〔33〕
〔2,1〕
〔513,1〕
〔1025,1〕
6.特征值归一化〔原理公式见报告P29〕
程序如下:
fori=1:
33
forj=1:
gy(i,j)=(tz(i,j)-min(tz(i,:
)))/(max(tz(i,:
))-min(tz(i,:
)));
%tz为原特征值矩阵,gy为归一化后的特征值矩阵。
33为原始特征值个数,20为数据样本个数。
得到以下表格:
故障轴承G201~G2010
时
域
特
征
值
80
频域
参
数
FFT谱
功
率
谱
小
波
包
络
正常轴承Z201~Z2010
值
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 轴承 matlab 处理 程序