一种为轴承故障诊断解调频段的新方法研究（Matlab代码实现）

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

轴承故障诊断对保证机械设备的安全十分重要。近年来，数据驱动的故障诊断方法得到了研究者的关注。

智能生产的发展加强了对早期故障诊断的需要，在机械设备中滚动轴承的运行状态往往直接影响机器能否正常工作。因此，精确的检测和诊断轴承的状态对于机械设备安全运行占有重要地位［

1］。 一直以来，研究者提出多种故障诊断算法对轴承故障进行分类。Cai 等［2］提出了能够进行诊断的面向对象贝叶斯网络复杂系统。Dong 等［3］使用 k － 最近邻方法来完成轴承的故障诊断。Konar 等［4］通过小波去噪对原始信号进行处理并提取特征，然后采用支持向量机进行轴承故障诊断。张钰等［5］从轴承振动信号中提取时域统计指标作为特征向量，然后利用随机森林进行诊断。这些方法不能直接从原始信号中提取特征，因此，研究直接从轴承振动信号中提取特征并自动识别轴承健康状态的诊断方法具有重要的研究价值。

📚2 运行结果

 

 

 部分代码：

%% Maximal Overlap Discrete Wavelet Packet Transform (modwpt)
% ==============================================================
tic
x = x - mean(x);

[wp,~,~] = modwpt(x(1:end),'db12','Fulltree',true,level);

for lev = 1 : level    
    jj = 0;
    for v = 0:lev-1
        jj = jj + 2^v;
    end
    jj = jj-1;   
    wpt(1:2^lev,1:length(x),lev) = wp(jj+1:jj+2^lev, 1:length(x));    
end

%% Plot the raw signal
% =======================
t = 0:1/Fs:length(x)/Fs-1/Fs;
figure('units','normalized','outerposition',[0.7 0.5 0.3 0.5])
plot(t,x)
set(gca,'Fontsize',15);set(gca,'FontWeight','Bold');set(gca,'Fontname','Times new roman')
xlabel('Time (s)'),xlim([0 max(t)])
ylabel('Acceleration (m/s^2)')
title('Raw time signal')

%% Compute Autogram   
% ====================
ss0 = floor(0.005*length(x));  % Removing the filter transient
if CUR == 0
    cor = 0.010;
elseif CUR == 1
    cor = 0.15;
end   
ss = floor(cor*length(x))+1;   % Removing a few first samples of AC
AAA = 1.55;

y = hilbert(x);
env = abs(y);
[acor,~] = xcorr((env).^2,'unbiased'); 
KUR(1,1) = Kurt((acor(length(x)+ss:floor(length(x)*AAA))),lmp);

for lev = 1:level    
    for i = 1:2^(lev)
        y = hilbert(wpt(i,ss0:length(x)-ss0,lev));
        env = abs(y);
        [acor,~] = xcorr((env).^2,'unbiased');  
        KUR(lev+1,i) = Kurt((acor(length(x)+ss:floor(length(x)*AAA)))',lmp);
    end     
end

T = 2^level;
for lev = 0:level
    TE = T/2^lev; 
    for i = 1:2^lev    
        color((i-1)*TE+1:(i)*TE,lev+1) = KUR(lev+1,i);         
    end       
end

[M,~] = max(KUR(:));
disp('----')
disp(['Autogram -- Max Kurtosis = ',num2str(M)])
[i, j] = find(KUR==M);
LM1 = i-1;
JM1 = j;
disp(['level = ',num2str(LM1)])
disp(['node = ',num2str(JM1)])

%% Plot Autogram
% ================
figure('units','normalized','outerposition',[0.0 0.5 0.3 0.5])
imagesc(0:0.5*Fs/2^level:0.5*Fs,0:1:level,color')
set(gca,'Fontsize',16);set(gca,'FontWeight','Bold')
set(gca,'Fontname','Times new roman')
xlabel('Frequency (Hz)','FontSize', 16);ylabel('Level','FontSize', 18)
colorbar
title(['Autogram - K_{max}= ',sprintf('%0.1f',M),' @ level ',num2str(LM1),', Bw= ',num2str(Fs/2^(LM1+1)),'Hz, f_c= ',num2str(j/2^LM1*Fs*0.5-Fs/2^(LM1+2)),'Hz'],'FontSize', 13)
toc                                     

%% Compute defect frequencies for the CWRU bearings
% ========================================================
fd = 0; fs = 0; fr = 0; asdf = 0;

if Spe ~= 0
speedr = [1797 1772 1750 1730]; 
ratio = [5.415 3.585 0.3983 2.357; 4.947 3.053 0.3816 1.994] ;
fr = speedr(Spe)/60;
if Drive_end == 1
    j = 1;
else
    j = 2;
end
fd = fr * ratio(j, Dn);
if Dn == 1
    asdf = 1;
    fs = fr;
elseif Dn == 4
    asdf = 1;
    fs = ratio(j,3)*fr;
else
    asdf = 0;
    fs = 0;
end
end

%% Squared Envelope Spectrum (SES)
% ==================================
NN = 4000; % length of the SES

% Manually set the node and level for demodulation
% JM1 = 15; % node
% LM1 = 5;  % level

if LM1 == 0
    c0 = x;
else
    c0 = squeeze(wpt(JM1,:,LM1)');
end

c = c0(ss0:end-ss0);
y = hilbert(c);
c = abs(y);

[acor,~] = xcorr(((c).^2),'unbiased');
XA111 = (acor(length(c)+ss:floor(length(c)*AAA-1)));
Han = hanning(length(c));
nfft = 2*ceil(length(c)/2);
c = c';
FT = abs(fft(((c).^2).*Han',nfft));
FT11 = FT;
Le = length(FT)-1;
XA111 = XA111 - min(XA111);
MEAN = movmean(XA111,10000);

% Plot SES
% ===========

if UL == 1
    figure('units','normalized','outerposition',[0.3 0 0.4 1])
    subplot(3,1,1)
    PLA(fd,fs,fr,Fs,Le,NN,FT11,asdf,1,Dn)
    title('SES (node with the highest kurtosis)','FontSize', 12)
    set(gca, 'XLabel', [])
else
figure
    PLA(fd,fs,fr,Fs,Le,NN,FT11,asdf,1,Dn)
    title('SES (node with the highest kurtosis)','FontSize', 12)
end

% Plot the spectra after Upper and Lower threshold 
% ==================================================

if UL == 1
        poi = XA111 - UP*MEAN;
        poiU = poi.*(poi>0);
        poi = XA111 - 1*MEAN;
        poiL = poi.*(poi<0);
        
    poiU = poiU';
    poiL = poiL';
    Han = hanning(length(XA111));
    FTpoiL = abs(fft(poiL(1:length(XA111)).*Han',nfft));
    FTpoiU = abs(fft(poiU(1:length(XA111)).*Han',nfft));
    
    subplot(3,1,2)
    PLA(fd,fs,fr,Fs,Le,NN,FTpoiL,asdf,0,Dn)
    title('SES (after lower threshold)','FontSize', 12)
    set(gca, 'XLabel', [])
    subplot(3,1,3)
    PLA(fd,fs,fr,Fs,Le,NN,FTpoiU,asdf,0,Dn)
    title('SES (after upper threshold)','FontSize', 12)
end

%% Compute CSES and Average CSES
% ===================================
clear Z X Y aaf 

N0 = 10;
for L = 1: level+1

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1][1]赵志宏,李春秀,窦广鉴等.基于MTF-CNN的轴承故障诊断研究[J].振动与冲击,2023,42(02):126-131.DOI:10.13465/j.cnki.jvs.2023.02.015.

[2]Moshrefzadeh A, Fasana A an effective approach for selecting the optimal demodulation band in rolling element bearings diagnosis

🌈4 Matlab代码实现