马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
clc
clear
sf=22800;
Z=360;
dearta=0.1;
A=xlsread('15kw转速.xls');
B=xlsread('15kw信号.xls');
% ttt=1/sf:1/sf:(length(B)/sf);%消除趋势项
% M=4; %趋势项阶数
% G=polyfit(ttt',B,M); %计算多项式待定系数向量G
% B=B-polyval(G,ttt'); %用x减去多项式系数a生成的趋势项
pb=2800;% ziji de dian
pa=floor(pb/2280*720);
B1=B(pb:pb+2285);
A1=A(pa-35:pa+720-1-35);
A11=A(pa-35:pa+720-1);
%%%%%%%%%%5 %去除奇点
for i=3:length(A1)
if abs(abs(A1(i)-A1(i-1))-abs(A1(i-1)-A1(i-2)))>=10
A1(i)=A1(i-1)+(A1(i-1)-A1(i-2));
end
end
%%%%%%%% jie bi gen zong fen xi
T=60./(A1.*Z); %由瞬时转速 求出每个脉冲间隔
R=(T.*Z)/360; %每个齿段内柴油机转过单位度数所需时间
figure(2)
bar(R,0.2) %柱形图表示
R2=0:360/360:720;%每个齿的角度
b=R2(1:end-1);
c=R2(2:end);
RR=(c+b)/2; %确定每个齿的中间曲轴转角度数
figure(4)
plot(RR,R);%插值前的角度时间图
Rx=linspace(R2(1),R2(720),720/dearta);%对角度的线性插值新序列
Tx=interp1(RR,R,Rx,'spline');%对时间的三次样条插值
figure(5)
plot(Rx,Tx);%插值后的角度时间图
t=1/sf;
for nj=1:length(Tx)
T0=sum(Tx(1:nj)*dearta);
% T0=sum(Tx(1:nj));
k1=T0/t;
k2=fix(k1);
if k2==0
a(nj)=0;
else
a(nj)=B1(k2)+(B1(k2+1)-B1(k2))/t*(T0-k2*t);
end
% a(nj)=B(k2)+(B(k2+1)-B(k2))/t*(T0-k2*t);
end
% subplot(2,1,1)
figure(6)
plot(Rx,a);%与角度对应的振动信号
fid = a-mean(a);
N=length(fid);
t=1/sf:1/sf:N/sf;
%%%%%%EMD fen jie
IMF = emd(fid);
[m,n] = size(IMF);
figure(7)
for i = 1:10
subplot(10,1,i);
plot(Rx,IMF(i,1:end));
end
%计算IMF相关性系数
for i=1:m;
a=corrcoef(IMF(i,:),fid);
xg(i)=a(1,2);
end
for i=1:m-1
mse(i)=mean(IMF(i,:).^2,2)-mean(IMF(i,:),2).^2;
end;
mmse=sum(mse);
for i=1:m-1
mse(i)=mean(IMF(i,:).^2,2)-mean(IMF(i,:),2).^2;
%方差百分比,也就是方差贡献率
mseb(i)=mse(i)/mmse*100;
end;
%显示各个IMF的方差和贡献率
for i=1:m-1
disp(['imf',int2str(i)]) ;disp([mse(i) mseb(i)]);
end;
len = size(IMF,1);
for k = 1:len
len1 = length(IMF(k,:));
e(k) = sum(IMF(k,:).*IMF(k,:))/len1;% 时域均方值能量
end
E=sum(e);
for k=1:len;
p(k)=e(k)/E;
H(k)=-(p(k)*log(p(k)));%能量熵
end
Hen=sum(H);
%%%%%%%%%san wei pu
[Aa,fa,tt]=hhspectrum(IMF);% IMF(1,:)
[im,tt1]=toimage(Aa,fa,tt,length(tt));
disp_hhs(im,tt1/sf);
plot_hht_3d(IMF,1024,4000);%三维hilbert图
%%%%%%%%%bian ji pu
for k=2:size(im,1)
bjp(k)=sum(im(k,:))*(1/sf)*(sf/N);
end
f2=(0:N-3)/N*(sf/2);
figure(10)
plot(f2,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
xlabel('频率 f/Hz');
ylabel('幅值');
title('边际谱');
程序如上,但运行时总出现
??? Out of memory. Type HELP MEMORY for your options.
Error in ==> disp_hhs at 67
im =10*log10(im/M);
Error in ==> SHICE at 104
disp_hhs(im,tt1/sf);
怎么解决
|