求助：轴承外圈故障包络谱分析（带数据）

dongyun.wang

刚刚编的轴承外圈故障包络分析，转速为1730，轴承外圈的故障的特征频率为103Hz,进行了取出趋势项，小波去噪，想得到
带有故障频率特征的细节包络谱，得到的分析结果却不理想。请高手帮忙分析一下。圣诞快乐 ！！数据在附件里
clc;
clear;
close all;
sampleFreq   = 12000;%采样频率
sampleLength = 2048;%采样点数
t=0:1/sampleFreq:(sampleLength-1)/sampleFreq;
Data1 = load('D:\MATLAB7\work\outer1.txt');
%==========================================================================
%直接进行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;
%==========================================================================
%预处理：去除趋势项
%==========================================================================
m=1;
data=Data';
a=polyfit(t,Data,m);
y=Data-polyval(a,t);
%==========================================================================
%预处理：小波降噪
%==========================================================================
[c,l]=wavedec(y,3,'sym3');
[thr,sorh,keepapp]=ddencmp('den','wv',y);
Data=wdencmp('gbl',y,'sym3',3,thr,sorh,keepapp);
%==========================================================================
%预处理：自相关
%==========================================================================
%  Lag=12000;
% [c,lags]=xcorr(y,Lag,'unbiased');
%  Data=c(1:Lag);

%==========================================================================
%共振解调法
%==========================================================================
%db10小波进行4层分解
%一维小波分解
[c,l] = wavedec(Data,4,'db4');
%重构第1-4层细节系数
d4 = wrcoef('d',c,l,'db4',4);
d3 = wrcoef('d',c,l,'db4',3);
d2 = wrcoef('d',c,l,'db4',2);
d1 = wrcoef('d',c,l,'db4',1);
envelop_hil = hilbert(Data);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
%==========================================================================
%d1的包络
%==========================================================================
figure(2)
envelop_hil = hilbert(d1);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
subplot(221)
plot(fft_plot_Hz,envelop_fft(2:1025));
title('d1包络波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;
%==========================================================================
%d2的包络
%==========================================================================
envelop_hil = hilbert(d2);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
subplot(222)
plot(fft_plot_Hz,envelop_fft(2:1025));
title('d2包络波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;
%==========================================================================
%d3的包络
%==========================================================================
envelop_hil = hilbert(d3);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
subplot(223)
plot(fft_plot_Hz,envelop_fft(2:1025));
title('d3包络波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;
%==========================================================================
%d4的包络
%==========================================================================
envelop_hil = hilbert(d4);
envelop_abs = abs(envelop_hil);  
%fft变换
envelop_fft = abs(fft(envelop_abs))*2 /sampleLength;
subplot(224)
plot(fft_plot_Hz,envelop_fft(2:1025));
title('d4包络波形');
ylabel('振幅/m/s^2');
xlabel('时间/s');
grid on;

liu222bing

回复 1 # dongyun.wang 的帖子

你这个用的也是国外网站的数据吧？

dongyun.wang

是的 数据是没问题的 不知道编程错在哪里
我作出的图大，没办法传上来。

hustyoung

数据有用，能否给个出处！！！！！！！

dongyun.wang

是美国Case Western Reserve Univisity的轴承数据，很早以前下的数据，具体网址不记得了了，得上网搜一下

duguzi

应该设法把图贴出来，这样大家看着也直观一些，你觉得哪里不对也应该说明，别人看程序也能有针对性。

dongyun.wang

回复 6 # duguzi 的帖子

F:\QQ.BMP
我的理想的d1的包络谱中有103赫兹及其倍频应该很明显地显示出来。
但是我的图中没有显示

duguzi

数据是有问题的，采样频率12000，采样点数2048，频率分辨率为12000/2048=5.8Hz，过于粗略，你要寻找103Hz，d1包络谱中有99.6Hz，近似相等，差别不到一个频率分辨率，可以认为是同一频率。
另我用自己的试验台振动信号处理了一下，该程序没有大问题，按你的处理过程，最终结果和直接Hilbert包络解调几乎一致。不要怀疑自己的算法了。

留个邮箱，可以把我的数据给你

dongyun.wang

回复 8 # duguzi 的帖子

dongyun.wang@163.com 谢谢

kiefer0107

你好，楼主。谢谢你的帖子，我从你的程序中学到了很多东西。问一个和你题目不相关的问题。
为什么你作频谱图时要从第2个点开始做图？
plot(fft_plot_Hz,envelop_fft(2:1025));

dongyun.wang

因为第一个点的值太高了，以至于后面的值都看不清了，你可以试一试，所以把第一个点去掉了

duguzi

回复 10 # kiefer0107 的帖子

主要是因为原始信号经过Hilbert变换再取绝对值后，即原始程序中的
envelop_hil = hilbert(d4);
envelop_abs = abs(envelop_hil);  
此时envelop_abs变量中全是大于等于0的数字，即该变量中包含有一较大的直流分量，
再进行fft变换，包络谱中在0Hz便有一较大的值，以至于其他有效非0Hz的频率成分被0Hz的
直流分量压制，在包络谱中显得比较小了，难以分辨，故作者从第二个值开始绘制。
其实在fft变换之前将envelop_abs先减掉均值，再进行fft即可。
即原始程序中可稍微修改一句：
%fft变换
envelop_fft = abs(fft(envelop_abs-mean(envelop_abs)))*2 /sampleLength;

yazhoulong

,支持数据，如果有声发射数据就好了

miwueshine

这几天看到大家的分析 学习了很多东西 谢谢大家的分析

yeguchan