伤城风絮 发表于 2015-4-2 19:30

HHT信号分析程序是在整不明白,求大神指导

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 dedian
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);
= 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();
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
=hhspectrum(IMF);% IMF(1,:)
=toimage(Aa,fa,tt,length(tt));
disp_hhs(im,tt1/sf);
plot_hht_3d(IMF,1024,4000);%三维hilbert图
%%%%%%%%%bianji 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);
怎么解决

伤城风絮 发表于 2015-4-2 19:32

毕业设计急!!!各位大神求帮忙

qqeell 发表于 2015-4-3 14:19

电脑不行,算不动

wenchao89 发表于 2015-4-4 19:01

论坛不错。。。
页: [1]
查看完整版本: HHT信号分析程序是在整不明白,求大神指导