Hilbert-Huang变换可以高内聚信号特征自适应的将信号分解成若干固有模态函数。更适用于非线性非平稳信号的处理。
缺点：
1、存在端点延拓；
2、分解判据确定；
3、Hilbert解调固有局限性。
一、介绍
Hilbert-Huang变换是经验模态分解（EMD）和Hilbert时频谱的统称。
步骤：
1、将信号用EMD分解为若干固有模态函数（IMF）；
2、对每个IMF分量进行Hilbert变换得到瞬时频率和瞬时幅值；
3、得到信号的时频分布。
1.1 固有模态函数（IMF）
IMF必须满足条件：
（1）曲线的极值点和零点的数目相等或至多相差1；
（2）在曲线的任意一点，包络的最大极值点和最小极值点的均值等于0。
1.2 EMD原理
瞬时频率
瞬时频率对多数信号不适用。
由于大多数非平稳信号不直接满足IMF条件，Huang提出了假设：任何复杂信号都是由一些相互独立的IMF分量组成；每个IMF分量可以是线性的，也可以是非线性的。
EMD过程：
（1）确定信号所有的局部极大值和极小值点；
（2）用三次样条函数对所有极大值点和极小值点分别进行插值运算，拟合出上下包络并求出平均值m1，计算
（3）如果h1是IMF，则h1是函数x(t)的第1个IMF分量；若h1不是IMF，则重复前两个步骤k此得到
h1k满足IMF条件，就是第1个IMF，记c1。
（4）从x(t)中减去c1，得到残差
重复（1）~（3），得到c2，以此类推，当rn为单调函数不能再提取IMF分量时，循环结束。得到n个IMF分量
终止判据
SD可以控制迭代的次数。一般取0.2~0.3比较合适。
1.3 EMD的局限性
1、原信号与其Hilbert变换构成解析信号，然后用解析信号的模和相位导数作为原信号的估计值，产生估计误差；
2、在调制信号两端出现较大误差。
1.4 matlab实现
imf=end(x)
二、Hilbert-Huang变换仿真实例
例1 对仿真信号进行EMD分解。
clc;
clear;
fs=1024;
n=1024;
t=0:1/fs:(n-1)/fs;
x=sin(2*pi*250*t)+sin(2*pi*100*t)+sin(2*pi*50*t)+0.1*randn(1,length(t));
imf=emd(x);
[m,w]=size(imf);
figure(1)
subplot(411)
plot(t,x);
title('原始波形');
subplot(412)
plot(t,imf(:,1));
ylabel('IMF1')
subplot(413)
plot(t,imf(:,2));
ylabel('IMF2');
subplot(414)
plot(t,imf(:,3));
ylabel('IMF3');
xlabel('t/s')计算前3个IMF分量
对前3个IMF分量进行fft变换
nfft=512;
fz=(1:nfft/2)*fs/nfft;
figure(2)
for(i=1:3)
subplot(3,1,i)
y(:,i)=fft(imf(:,i),nfft);
imffft(:,i)=abs(y(:,i));
plot(fz,imffft(1:nfft/2,i)/nfft);
xlim([0 nfft-1]);
end
从fft变换图可以看出，前三个IMF分量刚好对应频率250，100和50.
将信号进行小波包分解
T=wpdec(x,3,'db1');
Y(:,1)=wprcoef(T,[3,0]);
Y(:,2)=wprcoef(T,[3,1]);
Y(:,3)=wprcoef(T,[3,2]);
Y(:,4)=wprcoef(T,[3,3]);
Y(:,5)=wprcoef(T,[3,4]);
Y(:,6)=wprcoef(T,[3,5]);
Y(:,7)=wprcoef(T,[3,6]);
Y(:,8)=wprcoef(T,[3,7]);
figure(3)
for(i=1:8)
subplot(8,1,i)
plot(t,Y(:,i));
xlim([0 (n-1)/fs]);
end
利用前三个IMF分量对信号进行重构
figure(4)
xx=imf(:,1)+imf(:,2)+imf(:,3);
plot(t,x,t,xx,'-*')
title('原始信号&重构信号')
legend('原信号','重构信号')
从上图可以看出，重构信号与原信号基本重合。
三、Hilbert-Huang在机械故障诊断中的应用
3.1 EMD-AR谱提取柴油机活塞、活塞销故障特征
AR模型是时间序列分析中最基本的数学模型。不仅可以用来故障诊断，还可以对故障进行早期预测。还可以分析短样本信号，克服了Hilbert分离算法存在加窗效应的问题。
EMD先将复杂的非稳态振动信号分解为单分量信号，然后用AR模型进行分析。
AR模型主要用于平稳过程。
例2 利用EMD-AR谱方法分析柴油机活塞销、活塞磨损故障
采样频率12.8kHz，转速为1800r/min，采样点数为16384.
clc;
clear;
close all
fs=25600;
n=16384;
t=0:1/fs:(n-1)/fs;
s1=load('Sig1.txt');
s2=load('Sig2.txt');
s3=load('Sig3.txt');
figure(1)
subplot(311)
plot(t,s1);
title('正常信号');
axis tight
subplot(312)
plot(t,s2);
title('活塞销磨损信号')
axis tight
subplot(313)
plot(t,s3);
title('活塞磨损波形');
axis tight
xlabel('t/s')显示时域信号
EMD分解
imf1=emd(s1);
imf2=emd(s2);
imf3=emd(s3);
figure(2)
for(i=1:6)
subplot(6,1,i)
plot(t,imf1(:,i));
xlim([0 (n-1)/fs]);
end显示s1前6个分量
对前6个imf分量进行AR谱分析
N=1024;
for(i=1:6)
xpsd=pburg(imf1(:,i),10,N);
xpsd1(i,:)=10*log10(xpsd(1:400)+0.000001);
end
for(i=1:6)
xpsd=pburg(imf2(:,i),10,N);
xpsd2(i,:)=10*log10(xpsd(1:400)+0.000001);
end
for(i=1:6)
xpsd=pburg(imf3(:,i),10,N);
xpsd3(i,:)=10*log10(xpsd(1:400)+0.000001);
end
ss=N/2;
dd=(0:1/ss:1)*fs/2;
d=dd(1:400);
figure(3)
for(i=1:6)
subplot(3,2,i)
plot(d,xpsd1(i,:),'b-',d,xpsd2(i,:),'b--',d,xpsd3(i,:),'b-.');
legend('正常','活塞销磨损','活塞磨损');
ylabel('EMD-AR谱/db');
end
将前6个imf分量的AR谱累加
q1=sum(xpsd1);
q2=sum(xpsd2);
q3=sum(xpsd3);
figure(4)
plot(d,q1,'b-',d,q2,'b--',d,q3,'b-.');
legend('正常','活塞销磨损','活塞磨损');
xlabel('频率/Hz');
ylabel('EMD-AR谱/dD')
例3 将例2方法应用到轴承故障
clc;
clear;
close all
fs=25600;
n=16384;
t=0:1/fs:(n-1)/fs;
s1=load('Sig1.txt');
s2=load('Sig2.txt');
s3=load('Sig3.txt');
figure(1)
subplot(311)
plot(t,s1);
title('正常信号');
axis tight
subplot(312)
plot(t,s2);
title('活塞销磨损信号')
axis tight
subplot(313)
plot(t,s3);
title('活塞磨损波形');
axis tight
xlabel('t/s')时域图
EMD分解
imf1=emd(s1);
imf2=emd(s2);
imf3=emd(s3);
figure(2)
for(i=1:6)
subplot(6,1,i)
plot(t,imf1(:,i));
xlim([0 (n-1)/fs]);
end
N=1024;
for(i=1:6)
xpsd=pburg(imf1(:,i),10,N);
xpsd1(i,:)=10*log10(xpsd(1:400)+0.000001);
end
for(i=1:6)
xpsd=pburg(imf2(:,i),10,N);
xpsd2(i,:)=10*log10(xpsd(1:400)+0.000001);
end
for(i=1:6)
xpsd=pburg(imf3(:,i),10,N);
xpsd3(i,:)=10*log10(xpsd(1:400)+0.000001);
end
ss=N/2;
dd=(0:1/ss:1)*fs/2;
d=dd(1:400);
figure(3)
for(i=1:6)
subplot(3,2,i)
plot(d,xpsd1(i,:),'b-',d,xpsd2(i,:),'b--',d,xpsd3(i,:),'b-.');
legend('正常','活塞销磨损','活塞磨损');
ylabel('EMD-AR谱/db');
end
q1=sum(xpsd1);
q2=sum(xpsd2);
q3=sum(xpsd3);
figure(4)
plot(d,q1,'b-',d,q2,'b--',d,q3,'b-.');
legend('DE','BA','FE');
xlabel('频率/Hz');
ylabel('EMD-AR谱/dD')
