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')