基于EMD的包络谱进行故障诊断matlab程序实例

源地址为http://blog.sina.com.cn/s/blog_574d08530100rbe1.html

刚开始接触滚动轴承故障诊断通常都是一头雾水。

其实只要按部就班就可以了。

滚动轴承故障诊断分为数据采集、数据处理和故障辨识（或故障预测）。

一接到故障诊断这个课题，你首先要明白，这三个步骤中，你想搞哪块。一般在其中一块有突破，基本上硕士就可以毕业了。

下面介绍的是用EMD和包络解调进行数据处理，然后人工进行故障辨识。

故障数据：http://download.csdn.net/source/3107403

该故障数据是美国凯斯西储大学轴承数据中心的，数据使用方法见：http://blog.sina.com.cn/s/blog_574d08530100p3pr.html

matlab程序：（内圈故障、外圈故障）

%转速： 1750转/分
%转频：29.16（29）
%采样率：12K
%轴承型号：6205
%inner ring:5.4152
%outer ring:3.5848
�ge train:0.39828
%rolling element:4.7135
%外圈故障：104.57
%内圈故障：157.94（158）
%滚动体双故障：137.49
%滚动体单故障：68.75
%保持架外圈故障：11.62（12）【外圈静止，内圈转动】
%-----------------
%转速1750的6025轴承的深度7的内圈点蚀故障
%驱动计数端的内圈故障，故障明显，基于EMD的包络解调有效
%风扇计数端的内圈故障，故障效果不好，基于EMD的包络解调不是很有效
%基础计数端的内圈故障，故障效果不好，基于EMD的包络解调无效，只能看到转频，故障频率不明显
load 1750I7.mat;%内圈故障数据
x=X107_DE_time;%驱动计数段的内圈故障
fs=12000;%采样率
N=10240;%采样点数（100倍）
M=0;%采样数据段的起始位置
n=M:N-1;
t=n/fs;%信号时间序列
X=X107_DE_time(1:N);%装载 驱动计数端的内圈故障数据
%X=X107_FE_time(1:N);%装载 风扇计数端的内圈故障数据
%X=X107_BA_time(1:N);%装载 基础计数端的内圈故障数据
%X=X107_DE_time(1:N)-X107_BA_time(1:N);
y=X';%信号幅值序列
k_in=kurtosis(y);%峭度系数，正常轴承为3左右
figure;%画原始信号时域和频域图
subplot(211);plot(t,y);title('原始信号时域波形');
subplot(212);hua_fft1(y,fs);title('原始信号频谱');
figure;%原始信号的包络谱
subplot(211);hua_baoluo1(y,fs,1);title('原始信号包络谱');
subplot(212);hua_baoluo1(y,fs,1,500);title('原始信号部分频段包络谱');
imf=emd1(y);%经验模态分解
figure;%前三个IMF分量
subplot(311);plot(t,imf(1,:));title('IMF1时域波形图');
subplot(312);plot(t,imf(2,:));title('IMF2时域波形图');
subplot(313);plot(t,imf(3,:));title('IMF3时域波形图');
figure;%前三个IMF分量频谱
subplot(311);hua_fft(imf(1,:),fs,1);title('IMF1频谱');
subplot(312);hua_fft(imf(2,:),fs,1);title('IMF2频谱');
subplot(313);hua_fft(imf(3,:),fs,1);title('IMF3频谱');
figure;%前三个IMF分量选择频段内的包络谱
xf1=0;%需要查看的包络谱频率段起点频率
xf2=1000;%需要查看的包络谱频率段终止频率
subplot(311);hua_baol(imf(1,:),fs,1,xf1,xf2);title('IMF1包络谱');
subplot(312);hua_baol(imf(2,:),fs,1,xf1,xf2);title('IMF2包络谱');
subplot(313);hua_baol(imf(3,:),fs,1,xf1,xf2);title('IMF3包络谱');

从前三个IMF分量都能找到内圈故障频率谱线，并且在其倍频处也能找到

%----------

%外圈故障数据
load 1750o7-3.mat;
x=X146_DE_time;%驱动计数段的内圈故障
fs=12000;%采样率
N=10240;%采样点数（100倍）
M=0;%采样数据段的起始位置
n=M:N-1;
t=n/fs;%信号时间序列
X=X146_DE_time(1:N);%装载 驱动计数端的内圈故障数据
%X=X107_FE_time(1:N);%装载 风扇计数端的内圈故障数据
%X=X107_BA_time(1:N);%装载 基础计数端的内圈故障数据
%X=X107_DE_time(1:N)-X107_BA_time(1:N);
y=X';%信号幅值序列
k_in=kurtosis(y);%峭度系数，正常轴承为3左右
figure;%画原始信号时域和频域图
subplot(211);plot(t,y);title('原始信号时域波形');
subplot(212);hua_fft1(y,fs);title('原始信号频谱');
figure;%原始信号的包络谱
subplot(211);hua_baoluo1(y,fs,1);title('原始信号包络谱');
subplot(212);hua_baoluo1(y,fs,1,500);title('原始信号部分频段包络谱');
imf=emd1(y);%经验模态分解
figure;%前三个IMF分量
subplot(311);plot(t,imf(1,:));title('IMF1时域波形图');
subplot(312);plot(t,imf(2,:));title('IMF2时域波形图');
subplot(313);plot(t,imf(3,:));title('IMF3时域波形图');
figure;%前三个IMF分量频谱
subplot(311);hua_fft(imf(1,:),fs,1);title('IMF1频谱');
subplot(312);hua_fft(imf(2,:),fs,1);title('IMF2频谱');
subplot(313);hua_fft(imf(3,:),fs,1);title('IMF3频谱');
figure;%前三个IMF分量选择频段内的包络谱
xf1=0;%需要查看的包络谱频率段起点频率
xf2=1000;%需要查看的包络谱频率段终止频率
subplot(311);hua_baol(imf(1,:),fs,1,xf1,xf2);title('IMF1包络谱');
subplot(312);hua_baol(imf(2,:),fs,1,xf1,xf2);title('IMF2包络谱');
subplot(313);hua_baol(imf(3,:),fs,1,xf1,xf2);title('IMF3包络谱');

从包络谱中可以看到外圈故障频率的谱线