滚动轴承故障的共振解调法

long9998

下面是我写的轴承故障信号处理程序，得到的故障特征不太明显 ，大家帮我看看，怎么改啊 ？附件是轴承数据
clc;
fs=8192;
N=8192;
n=0:N-1;
t=n/fs;
load data_03.txt;
s=data_03;
x=s';
%小波降噪
[c,l]=wavedec(x,3,'db3');
[thr,sorh,keepapp]=ddencmp('den','wv',x)
sd=wdencmp('gbl',x,'db3',3,thr,sorh,keepapp);
%带通滤波
[b,a]=butter(4,[0.125,0.5]);
[h,w]=freqz(b,a);
Sf=filter(b,a,sd);
%检波
sf=hilbert(Sf);
y=abs(sf);
y=y-mean(y);
%低通滤波
[c,d]=butter(4,0.5,'low');
y1=filter(c,d,y);
nfft=8192;
p=abs(fft(y1,nfft));%fft变换
plot((0:nfft/2-1)/nfft*fs,p(1:nfft/2));
%细化【50，200（左右）】
fa=50;
n=0:N-1;
%频移
b=2*n*pi*fa/fs;
q=y1.*exp(-i*b);
[b,a]=butter(4,0.125,'low');
y6=filter(b,a,q);
np=40;
y3=resample(y6,1,np);%重新采样，采样频率为fs/N
y2=abs(fft(y3));
w=fa:(nfft/np+fa);
figure
plot(w,y2)

beijing2009

为什么加一个带通滤波？

long9998

带通滤波是共振解调法的一个步骤，其作用是消除低频干扰。也可不通过任何滤波器直接通过解调、抗混滤波及谱分析

nkdtxf

顶一下，怎么没人回答楼主的问题呢，我也需要答案，还有啊，楼主怎么选择饿带通滤波器的中心频率呢？

zhlong

可以参考 文献  基于遗传算法和峰度最佳的滚动轴承故障诊断
               基于自适应复平移Morlet小波的轴承包络解调分析方法
               Bearing fault detection based on optimal wavelet filter and sparse code shrinkage

nkdtxf

fs=8192;
N=8192;楼主的，这句好像有问题啊。采样频率和采样点数不相等吧

xiangyu537

?????????????????????????

xiangyu537

怎末没人回答？？？？？？？？？？？

long9998

重新贴一个共振解调的程序！
%共振解调
clc;
clear;
close all;
%参数设置
%  sampleFreq   = 12000;
%  sampleLength = 8192;
sampleFreq   = 12000;%采样频率
sampleLength = 12000;%采样点数
%下载正常数据
% good      =  load('G:\CWRU_DATA\97.mat');
% data_good = good.X097_DE_time(1:sampleLength);
%下载内圈故障数据
% inner      = load('G:\CWRU_DATA\3001.mat');
% data_inner = inner.X056_DE_time(1:sampleLength);
%下载外圈故障数据
% outer      = load('G:\CWRU_DATA\198.mat');
% data_outer = outer.X198_DE_time(1:sampleLength);
%下载滚动体故障数据
% ball      = load('G:\CWRU_DATA\3005.mat');
% data_ball = ball.X048_DE_time(1:sampleLength);
%产生仿真信号
t=0:1/sampleFreq:(sampleLength-1)/sampleFreq;
MySignal=sin(2*pi*500*t).* (sin(2*pi*5*t)+1);
%下载齿轮数据
% t=0:1/sampleFreq:(sampleLength-1)/sampleFreq;
% gear=load('G:\齿轮试验台2130\TXT_data\Bearing9\Equipment1-911.txt');
%故障诊断课轴承数据
% t=0:1/sampleFreq:(sampleLength-1)/sampleFreq;
% bear=load('C:\Documents and Settings\jmnmu\My Documents\移动文件夹\Data5\0\09-06-16-23-06-18 .txt');
%信号赋值
Data1 = MySignal;
%==========================================================================
%直接进行FFT变换
%==========================================================================
for n=1:sampleLength
   if abs(Data1(n))>10;
    Data(n)=0;
else Data(n)=Data1(n);
   end
end
fft_result = abs(fft(Data)) * 2 / sampleLength;
%画图的坐标变换
time_plot_s = 0:1/sampleFreq:(sampleLength-1) / sampleFreq;
fft_plot_Hz = sampleFreq*(1:sampleLength/2)/sampleLength;
figure;
subplot(221)
plot(time_plot_s,Data);
title('时域波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;
subplot(222)
plot(fft_plot_Hz,fft_result(1:sampleLength/2));
title('频域波形');
ylabel('振幅/m/s^2');
xlabel('频率/Hz');
grid on;
N=length(Data1);
Xrms1=sqrt(sum(Data1.^2)/N)
beta1=sum(Data1.^4)/N
kv1=beta1/(Xrms1^4)
%==========================================================================
%预处理：去除趋势项
%==========================================================================
m=1;
data=Data';
a=polyfit(t,Data,m);
y=Data-polyval(a,t);
%==========================================================================
%预处理：小波降噪
%==========================================================================
[c,l]=wavedec(y,3,'sym8');
[thr,sorh,keepapp]=ddencmp('den','wv',y);
Data=wdencmp('gbl',y,'sym8',3,thr,sorh,keepapp);
%==========================================================================
%预处理：自相关
%==========================================================================
%  Lag=12000;
% [c,lags]=xcorr(y,Lag,'unbiased');
%  Data=c(1:Lag);

%==========================================================================
%共振解调法
%==========================================================================
%带通滤波
subplot(223);
plot(time_plot_s,Data);

[b,a]        = butter(4,[0.2,0.8]) ;
filter_data  = filter(b,a,Data);
%包络
envelop_hil = hilbert(filter_data);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
subplot(223)
plot(time_plot_s,envelop_abs);
title('包络时域波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;
subplot(224)
plot(fft_plot_Hz,10*envelop_fft(2:sampleLength/2+1));
title('包络频域波形');
ylabel('振幅/m/s^2');
xlabel('频率/Hz');
grid on;
N=length(envelop_fft);
Xrms=sqrt(sum(envelop_fft.^2)/N)
beta=sum(envelop_fft.^4)/N
kv=beta/(Xrms^4)

200810100659

很不错。谢谢分享

xjd

谢谢9楼分享程序，我想问一下带通滤波的中心频率如何选取？现在我在用LabVIEW做共振解调，共振解调频率有问题，得到的包络频谱图不理想。

mazhiping12311

谢谢分享，但是我还是不懂啊！

王晓华

回复 9 # long9998 的帖子

你做的包络是Hilbert有没有能量算子的程序呢？急求！

lragy

顶一下！！！！！！！！虽然现在看不太懂 ，小弟刚选这个方向，以前是学电路的！

BuDingYu

正在做这方面的研究 哈哈 顶一下