最近一直在做HHT，大家看看，我把我最近做的贴出来

ihod

%设置参数
clc;clear;
load('IR007_0.mat');
fs=12000;
N=1024;
t=1/12000:1/12000:N/12000;
x=X105_DE_time(1:1024,1)';
%EMD分解，展示各imf分量
imf=emd(x);
x1=imf(1,:);
x2=imf(2,:);
x3=imf(3,:);
x4=imf(4,:);
x5=imf(5,:);
x6=imf(6,:);
x7=imf(7,:);
x8=imf(8,:);
x9=imf(9,:);
figure(1);
subplot(4,1,1);plot(t,x1),title('imf1');
subplot(4,1,2);plot(t,x2),title('imf2');
subplot(4,1,3);plot(t,x3),title('imf3');
subplot(4,1,4);plot(t,x4),title('imf4');
xlabel('时间(time)t/s'),ylabel('幅值/mm');
figure(2);
subplot(5,1,1);plot(t,x5),title('imf5');
subplot(5,1,2);plot(t,x6),title('imf6');
subplot(5,1,3);plot(t,x7),title('imf7');
subplot(5,1,4);plot(t,x8),title('imf8');
subplot(5,1,5);plot(t,x9),title('res');
xlabel('时间(time)t/s'),ylabel('幅值/mm');

%计算imf的行列维数
[m,n]=size(imf);

for i=1:m-1
%--------------------------------------------------------------------
%计算各IMF的方差贡献率
%定义：方差为平方的均值减去均值的平方
%均值的平方
%imfp2=mean(c(i,:),2).^2
%平方的均值
%imf2p=mean(c(i,:).^2,2)

%各个IMF的方差
mse(i)=mean(imf(i,:).^2,2)-mean(imf(i,:),2).^2;
end;
mmse=sum(mse);%总的方差
for i=1:m-1
%方差百分比，也就是方差贡献率
mseb(i)=mse(i)/mmse*100;
end;
figure(3);
bar([1:m-1],mseb);%可以改进，最好能做出柱状图
title('各imf分量的方差贡献率');
xlabel('方差'),ylabel('百分比');



%原始信号的Hilbert变换
hx=hilbert(x);
xr=real(hx);xi=imag(hx);
%计算瞬时振幅
A=sqrt(xr.^2+xi.^2);
figure(4);plot(t,A);title('瞬时振幅')
%计算瞬时相位
P=angle(hx);
figure(5);plot(t,P);title('瞬时相位')
%计算瞬时频率
dt=diff(t);
dx=diff(P);
sp=dx./dt;
figure(6);plot(t(1:N-1),sp);title('瞬时频率')

%imf1的Hilbert变换
xn1=hilbert(imf(1,:));
xr1=real(xn1);
xi1=imag(xn1);
%imf1的瞬时振幅
A1=sqrt(xr1.^2+xi1.^2);
figure(7);
subplot(2,1,1);plot(t,A1);
xlabel('时间（t）');ylabel('瞬时振幅');title('imf1')
%imf2的Hilbert变换
xn2=hilbert(imf(2,:));
xr2=real(xn2);
xi2=imag(xn2);
%imf2的瞬时振幅
A2=sqrt(xr2.^2+xi2.^2);
subplot(2,1,2);plot(t,A2);
xlabel('时间（t）');ylabel('瞬时振幅');title('imf2')
%imf1的瞬时相位
P1=atan2(xi1,xr1);
figure(8);
subplot(2,1,1);plot(t,P1);
xlabel('时间（t）'); ylabel('瞬时相位');title('imf1')
%imf2的瞬时相位
P2=atan2(xi2,xr2);
subplot(2,1,2);plot(t,P2);
xlabel('时间（t）'); ylabel('瞬时相位');title('imf2')
%imf1瞬时频率
xh1=unwrap(P1);
fs=1000;
xhd1=fs*diff(xh1)/(2*pi);
figure(9);
subplot(2,1,1);plot(t(1:1023),xhd1);
xlabel('时间（t）');ylabel('瞬时频率');title('imf1')
%imf2瞬时频率
xh2=unwrap(P2);
fs=1000;
xhd2=fs*diff(xh2)/(2*pi);
subplot(2,1,2);plot(t(1:1023),xhd2);  
xlabel('时间（t）');ylabel('瞬时频率');title('imf2')

%计算HHT的时频谱
[A, fa, tt] = hhspectrum(imf);
if size(imf,1) > 1
    [A,fa,tt] = hhspectrum(imf(1:end-1, :));
else
    [A,fa,tt] = hhspectrum(imf);
end
[E,tt1]=toimage(A,fa,tt,length(tt));   
disp_hhs(E,[],fs) %二维图形显示HHT视频谱，E是求得的HHT谱

for i=1:m-1
    faa=fa(i,:);
[FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图
figure(11);
surf(FA,TT1,E)
title('HHT时频谱三维显示')
end
   
%画边际谱
E=flipud(E);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
f = (0:N-3)/N*(fs/2);

figure(12)
plot(f,bjp);



有几个问题：
1）我用的数据是西储大学轴承数据，采样频率是12k，内圈故障，但是边际谱出来的故障频率不对啊，我计算故障频率应该是162hz左右，图像怎么不对，问题出在哪？
2）我自己觉得我对边际谱的理解不深，到底边际谱是做什么的？应该单独对一个imf分量求边际谱还是对所有的imf分量都求？
3）谁帮我看看我的imf分解怎么样？我自己感觉还可以啊？
4）程序有问题吗，哪些地方需要改进？

希望大家能帮我，共同进步，谢谢，祝大家每天开心！

ihod

回复 2 # summerxt404 的帖子

谢啦，不容易啊，为什么没人帮我看啊

zhangnan3509

西储的轴承数据用HHT处理，很多效果都不好，从上面的IMF看就知道了，调制很明显。

ihod

回复 3 # zhangnan3509 的帖子

谢谢，那我现在应该怎么办，换数据吗?还是用别的方法？换数据的话，我也没有别的采集数据啊，HHT最近也是看了一段，稍有了解，但是很浅。您有什么办法吗？

ihod

怎么没人啊

ihod

我有几个问题，关于在画边际谱这一块的
%画边际谱
E=flipud(E);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;
end
f = (0:N-3)/N*(fs/2);

figure(12)
plot(f,bjp);

横坐标是f，但是f = (0:N-3)/N*(fs/2);
也就是说f与采样频率fs直接相关，而我的采样频率是12000HZ，但是我的故障频率162HZ根本根本不能分辨出来啊，怎么回事

tamujin

回复 1 # ihod 的帖子

楼主你好，我现在也做了一下，发现边际谱也没有看到故障频率。而且你的几个问题我也有疑惑，请问你现在有进展了吗？

ihod

回复 7 # tamujin 的帖子

你是交大的吗？我也在交大

tamujin

回复 8 # ihod 的帖子

是啊，校友你好，多多交流啊:handshake

dingdingysu

也遇到了这样的问题，顶一下！！！

daxue123

怎么都赶在一块了，头疼

znas0707

那个三维的频谱我怎么就出不来呢？

午后虫鸣

很巧，距离LZ发帖，到今天正好两年了，目前我也在研究HHT。前辈啊

cwb

你好，可以解释一下toiamge得到的结果中E到底是什么意思吗？它的行和列是受toiamge函数中的参数个数控制的，我一直都没搞明白为什么默认的行是400行。。。如果加上tt和length(tt)的话又出错，，，
还有为什么对E的每一行进行累加得到的就是边际谱，难道E的每一行就代表某个频率在整个时间跨度上的幅度值吗？
如能得到解答，感激不尽，谢谢！！！

cwb

为何要计算方差贡献率呀？