声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1622|回复: 3

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

[复制链接]
发表于 2015-4-2 19:30 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

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);
怎么解决

回复
分享到:

使用道具 举报

 楼主| 发表于 2015-4-2 19:32 | 显示全部楼层
毕业设计急!!!各位大神求帮忙
发表于 2015-4-3 14:19 | 显示全部楼层
电脑不行,算不动
发表于 2015-4-4 19:01 | 显示全部楼层
论坛不错。。。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-19 17:29 , Processed in 0.054151 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表