👨🎓个人主页:研学社的博客
💥💥💞💞欢迎来到本博客❤️❤️💥💥
🏆博主优势:🌞🌞🌞博客内容尽量做到思维缜密,逻辑清晰,为了方便读者。
⛳️座右铭:行百里者,半于九十。
📋📋📋本文目录如下:🎁🎁🎁
目录
💥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.