【故障诊断】频域多点峰度重复瞬变提取轴承故障诊断研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客 

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

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

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

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

💥1 概述

本文从故障滚动体轴承或齿轮振动信号中提取重复瞬变 （RT）。在我们的测试中，它可用于干扰抑制，特别是脉冲噪声，复合故障诊断和某些情况下的早期故障诊断。

轴承故障诊断对保证机械设备的安全十分重要 。 近年来，数据驱动的故障诊断方法得到了研究者的关 注。 与传统的依赖于专家经验的故障特征提取方法不同，深度学习方法可以实现端到端自动故障特征提取与分类 。

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

📚2 运行结果

 

 

  部分代码：

% CHECK INPUT VALUES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = length(x);
N2 = log2(N) - 7;
if nlevel > N2
   error('Please enter a smaller number of decomposition levels');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SET THE UPPER AND LOWER OF FREQUENCY HUNTING ZONE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fcf_lower=FCF-5;
fcf_upper=FCF+5;
lower=floor(fcf_lower*N/Fs);
upper=ceil(fcf_upper*N/Fs);
fprintf('The lower frequency of frequency hunting zone: %4.4f\nThe upper frequency of frequency hunting zone: %4.4f\n',fcf_lower,fcf_upper);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FAST COMPUTATION OF THE MKURTGRAM 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
% Analytic generating filters
N = 16;            fc = .4;
h = fir1(N,fc).*exp(2*1i*pi*(0:N)*.125);
g = fir1(N,fc).*exp(2*1i*pi*(0:N)*.375);

N = fix(3/2*N);
h1 = fir1(N,2/3*fc).*exp(2*1i*pi*(0:N)*.25/3);
h2 = h1.*exp(2*1i*pi*(0:N)/6);
h3 = h1.*exp(2*1i*pi*(0:N)/3);

MKwav = MK_wpQ(x,h,g,h1,h2,h3,nlevel,upper,lower);        % multipoint kurtosis of the complex envelope
MKwav = MKwav.*(MKwav>0);                            % keep positive values only!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GRAPHICAL DISPLAY OF RESULTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
Level_w = 1:nlevel;    Level_w = [Level_w;Level_w+log2(3)-1];    
Level_w = Level_w(:); Level_w = [0 Level_w(1:2*nlevel-1)'];
freq_w = Fs*((0:3*2^nlevel-1)/(3*2^(nlevel+1)) + 1/(3*2^(2+nlevel)));
imagesc(freq_w,1:2*nlevel,MKwav),colorbar,[I,J,M] = max_IJ(MKwav);
xlabel('frequency [Hz]'),set(gca,'ytick',1:2*nlevel,'yticklabel',round(Level_w*10)/10),ylabel('level k')
fi = (J-1)/3/2^(nlevel+1);   fi = fi + 2^(-2-Level_w(I));
title(['Mkurt_{max}=',num2str(round(100*M)/100),' @ level ',num2str(fix(10*Level_w(I))/10),', Bw= ',num2str(Fs*2^-(Level_w(I)+1)),'Hz, f_c=',num2str(Fs*fi),'Hz'])
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SIGNAL FILTERING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r = [];
test = input('Do you want to filter out repetitive transients signal from the Mkurtgram? (yes = 1 ; no = 0): ');
while test == 1
   fi = input(['    Enter the optimal carrier frequency (btw 0 and ',num2str(Fs/2),') where to filter the signal: ']);
   fi = fi/Fs;
   lev = input(['    Enter the optimal level (btw 0 and ',num2str(nlevel),') where to filter the signal: ']);
   r = Find_wav_kurt(x,h,g,h1,h2,h3,lev,fi,Fs);
   test = input('Do you want to keep on filtering out repetitive transients? (yes = 1 ; no = 0): ');
end 
end
% END OF THE MAIN ROUTINE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LIST OF SUBROUTINES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function K = Mkurt(x,upper,lower)
% Compute the multipoint kurtosis of signal x
   E=abs(x).^2;
   E=E-mean(E);
   N=length(E);
   y=fft(E,N);
   y=abs(y)*2/N;
   L=round(N/2);
   y=y(1:L);
   range=lower:1:upper;
   M=length(range);
   k = 1;
   T = zeros(L,M);
   for i=range
       for j=i:i:L
           T(j,k)=1;
       end
       k=k+1;
   end
   MK =zeros(M,1);
   for i=1:M
%        MK(i)=sum(((T(:,i).*y).^4))*(sum(T(:,i).^2))^2/((sum(y.^2))^2 * sum(T(:,i).^8));
       MK(i)= nanmean((T(:,i).*y).^4)/nanmean(y.^2)^2;
   end
   K=max(MK);
end 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [I,J,M] = max_IJ(X)
% Returns the row and column indices of the maximum in matrix X.
[temp,tempI] = max(X);
[M,J] = max(temp);
I = tempI(J);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = binary(i,k)
% Returns the coefficients of the binary expansion of i: 
% i = a(1)*2^(k-1) + a(2)*2^(k-2) + ... + a(k)
if i>=2^k
   error('i must be such that i < 2^k !!')
end
a = zeros(1,k);
temp = i;
for l = k-1:-1:0
   a(k-l) = fix(temp/2^l);
   temp = temp - a(k-l)*2^l;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LIST OF SUBROUTINES FOR THE WAVELET PACKET KURTOGRAM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function K = MK_wpQ(x,h,g,h1,h2,h3,nlevel,upper,lower)
% Generate Mkurtgram through 1/3-binary tree of filter-banks and multipoint
% kurtosis
level = nlevel;
x = x(:);                 % shapes the signal as column vector if necessary

   Kd3 = Kd3*ones(1,Long);
   KQ = [Ka1 Ka2 Ka3 Kd1 Kd2 Kd3; KaQ KdQ];
end

if level == nlevel    
   K1 = Mkurt(x,upper,lower);        % multipoint kurtosis of the raw signal is computed here
   K = [K1*ones(1,length(K));K];
   
   [a1,a2,a3] = TBFB(x,h1,h2,h3);
   Ka1 = Mkurt(a1(Lh:end),upper,lower);

a = a(2:2:N);
a = a(:);
% highpass filter
d = filter(g,1,x);
d = d(2:2:N);
d = d(:);
end
% ------------------------------------------------------------------------
function [a1,a2,a3] = TBFB(x,h1,h2,h3)
% Trible-band filter-bank (ternary decomposition).
% (Subroutine of MK_wpQ_local')
N = length(x);
% lowpass filter
a1 = filter(h1,1,x);
a1 = a1(3:3:N);
a1 = a1(:);
% passband filter
a2 = filter(h2,1,x);
a2 = a2(3:3:N);
a2 = a2(:);
% highpass filter
a3 = filter(h3,1,x);
a3 = a3(3:3:N);
a3 = a3(:);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = Find_wav_kurt(x,h,g,h1,h2,h3,Sc,Fr,Fs)
% Filter the optimal band of the original signal
level = fix(Sc) + (rem(Sc,1)>=0.5)*(log2(3)-1);
Bw = 2^(-level-1);
freq_w = (0:2^level-1)/(2^(level+1))+Bw/2;
[~,J] = min(abs(freq_w-Fr));
fc = freq_w(J);
i = round((fc/Bw-1/2));
if rem(level,1) == 0
   acoeff = binary(i,level);
   bcoeff = [];
   temp_level = level;
else
   i2 = fix(i/3);
   temp_level = fix(level)-1;
   acoeff = binary(i2,temp_level);
   bcoeff = i-i2*3;
end
acoeff = acoeff(end:-1:1);
c = MK_wpQ_filt(x,h,g,h1,h2,h3,acoeff,bcoeff,temp_level);
spec = input('    Do you want to see the envelope spectrum? (yes = 1 ; no = 0): ');
figure
t = (0:length(x)-1)/Fs;
tc = linspace(t(1),t(end),length(c));
subplot(1+spec,1,1)
plot(tc,abs(c))
xlabel('Time [s]'),xlim([tc(1) tc(end)])
ylabel('Amplitude [g]')
title('Envelope of filtered signal')
if spec == 1
   env = abs(c).^2;
   env=env-mean(env);
   nfft=length(env);
   S=fft(env,nfft);


end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = MK_wpQ_filt_local(x,h,g,h1,h2,h3,acoeff,bcoeff,level)
% (Subroutine of function MK_wpQ_filt)
[a,d] = DBFB(x,h,g);    % perform one analysis level into the analysis tree
N = length(a);                       
d = d.*(-1).^(1:N)';
if level == 1
   if isempty(bcoeff)
      if acoeff(level) == 0
         c = a(length(h):end);
      else
         c = d(length(g):end);
      end
   else
      if acoeff(level) == 0
         [c1,c2,c3] = TBFB(a,h1,h2,h3);
      else
         [c1,c2,c3] = TBFB(d,h1,h2,h3);
      end
      if bcoeff == 0

🎉3 参考文献

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

[1]Y. Liao, P. Sun, B. Wang, L. Qu, Extraction of repetitive transients with frequency domain multipoint kurtosis for bearing fault diagnosis, Measurement Science and Technology, 29 (2018) 055012, 10.1088/1361-6501/aaae99.

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

🌈4 Matlab代码实现