经验模式分解（EMD）及希尔伯特-黄变换（HHT）简介及matlab实现

本文介绍过程涉及到两个链接工具包，可以自己网上搜索下载，以下提供了网盘下载的地址，因为作者主要做语音方面工作，所以后面的说明主要以说话人识别为例。（链接：https://pan.baidu.com/s/1LWzlEO6Vp7CqInjqPnV8_A
提取码：zga8 ）

一、经验模式分解（EMD）

关于经验模态分解的概念，网上有很多具体的讲解，这里就不进行细说，具体过程可以参考https://blog.csdn.net/lvsehaiyang1993/article/details/80875679。

目的：将任意信号（尤其是非平稳非线性时间序列信号）分解为线性稳态信号（IMF）
核心：：把任一个自由的信号分解成为若干个固有模态函数(IntrinsicModeFunction， IMF)和一个残余分量，每一个 IMF 代表了原信号不同频率段的振荡变化，反映信号的局部特征，而最后的残余分量则反映信号中的缓慢变化量。
过程：（该过程了解分解原理，知道思路即可，具体算法看看就好，毕竟有现成工具包）

判断h(t)是否为IMF的条件主要有两个：

在信号序列中，它的极值点与过零点的数目需要是一致的，或者相差必须少于一个
任意的时间点，信号序列的局部极大值与局部极小值分别构成的上包络线与下包络线的平均值须为 0，也就是说整个序列的波形必须是局部对称的。

IMF 特征：上包络和下包络是对称的，但大部分的信号并不符合这些条件。
IMF：是一个单分量信号，每时每刻有且仅有一个频率，即瞬时频率。

分解循环结束的条件有两个：残余分量r(t)是单调函数或者序列只有一个极值点

EMD分解举例

下图中：signal为原始信号，imf1~imf7为7个固有模态函数，即IMF分量，可以看到各个IMF分量分别对应不同的频率区间，res为剩余分量，可以看出其满足截至的条件。

二、希尔伯特变换（Hilbert Transform）

1、看了一些介绍，对Hibert变换的解释是说：相移90°，具体的意义是什么呢？
答：已知复信号的实部，求该复信号（虚部就是实部的希尔伯特变换）

实际接收信号是高频实信号（双边带），是基带信号经载波信号处理以后得到的，有用信息是基带信号是复信号（单边频谱）。复信号：求包络，求相位等操作方便

复信号：表示为a+ib，其中a为实信号输出，根据欧拉公式exp(x)=cosx+isinx，则b为a相移90°以后的信号，即a的希尔伯特变换，所以可以从实信号得到复信号

2、希尔伯特变换的条件是什么呢？
答：信号必须是窄带信号、不含噪声，不是复杂信号，不出现尖峰跳跃，数据尽量平稳，没有“骑波”，避免不可解释的负频率出现等等

带宽：信号频率的宽度，即信号所含谐波最高频率与最低频率之差
窄带信号：信号带宽远小于中心频率，也就是所含的频率分量区间尽量小
如果是实信号，其具有双边频谱，希尔伯特变换以后获得的瞬时频率、幅度和相位包括负值，无意义。

三、希尔伯特-黄变换（HHT）

希尔伯特-黄变换（HHT）就是先对信号进行EMD求固有模态函数IMFs，然后对各个IMF分量做希尔伯特变换，如下图所示。
因为希尔伯特变换的前提条件是：窄带信号、非复杂信号，平稳信号，但是生活中存在的大部分信号均不满足该条件，为了使用希尔伯特变换，必须将非线性平稳信号转换为平稳信号，因此，必须使用EMD进行模态的分解，再对分解得到的IMF进行希尔伯特变换，求得复信号以后再进行进一步的分析。

HHT后得到的是一种时频属性，既包括时域又包括频域的性质，是一种时频分析方法

HHT谱：信号的希尔伯特变换后做fft，表示信号幅值在整个频率段上随时间和频率的变化规律（有点类似语谱图的感觉）

HHT边谱：类似于傅里叶谱。对hilbert谱的做时间上的积分，表示信号幅值在整个频率段上随频率的变化情况，比傅里叶谱具有更高的频率分辨率。

HHT边际谱可以处理非平稳信号，如果信号中存在某一频率的能量出现，就表示一定有该频率的振动波出现，也就是说，边际谱能比较准确地反映信号的实际频率成分。

傅里叶变换只能处理平稳信号。傅里叶变换不能得到瞬时频率，即不能得到某个时刻的频率值。

HHT谱区别于傅里叶谱的地方，其得到的谱信息既包括时域，又包括频域，傅里叶变换只包括频域，没有对应的时域信息。

四、基于EMD的语音信号处理（以说话人识别为例）

语音信号是非稳态信号，只有在短时才具有稳态性，所以在信号处理的时候，假设信号具有稳态性，所以通过分帧（帧长10—30ms）将信号转换为稳态信号，但是分帧的过程会因为加窗类型和窗长的不同，对原始信号产生影响。EMD的作用与加窗分帧类似，其本质就是将非稳态信号转换为稳态信号，然后进行后续的处理，其中的IMF就是不同的稳态信号，对应的是不同的频率区间，一般由高到低排列，直到出现残余量r，过程终止。
语音信号经EMD过程，分解为各个频率范围的IMF分量。在后续的信号处理过程中，可以根据需求对各个IMF分量进行加权后求和，得到新的语音信号，然后可以完成进一步的操作。下图为基于EMD的说话人特征参数提取的过程，可做参考。

五、希尔伯特-黄变换（HHT）的matlab实现

1、EMD分解

[x,fs]=audioread('C:\Users\如初\Desktop\matlab程序\s01_1_1.wav')
y=x(:,1)%双通道，只读取单声道的数据
y=y(501:1000)%选取501-1000之间的500个点做分析
t= 1:length(y)
imf=emd(y)%进行EMD分解各个固有模态函数IMF（i）
emd_visu(y,t,imf)%模态分解图（图1、2、3分别对应模态分解图，信号合成图1，信号合成图2）

图1为EMD后各个IMF分量以及剩余信号r的波形，图2和图3均为各个IMF分量加权求和以后得到的新的信号，其中图2 是从imf1到imf7的顺序，图3是从imf7到imf1的顺序。

2、HHT谱图及HHT边谱图

%HHT谱图
[A,f,tt] = hhspectrum(imf)%imf各个分量分别进行希尔伯特变换，然后分别做fft，此处的f是归一化的频率
[E,t,Cenf]=toimage(A,f)%合成HHT谱图（类似语谱图）
disp_hhs(E)%绘制HHT谱图%HHT边谱图
for k=1:size(E,1)bjp(k)=sum(E(k,:))/fsendfigure(5)plot(Cenf(1,:)*fs,bjp)xlabel('频率 / Hz')ylabel('幅值')title('HHT边谱图')

3、瞬时频率和瞬时幅度（包络）

%各个imf对应的波形
for i = 1:size(imf,1)
figure(6)
subplot(3,1,1)
plot((1:500)/fs,imf(i,:))%横坐标/fs，是由采样点转换为时间的标度
hold on
end
legend('1','2','3','4','5','6','7','8');
xlabel('时间 t/s');ylabel('波形');title('imf各个分量波形')%瞬时频率：各个imf分量的瞬间频率
for i = 1:size(imf,1)
figure(6)
subplot(3,1,2)
plot(t/fs,f(i,:)*fs)%f(i,:)*fs是因为之前f输出是归一化的频率，需要*fs才是原始的频率
legend('imf(i)')
hold on
end
legend('1','2','3','4','5','6','7','8');
xlabel('时间 t/s');ylabel('幅值');title('imf各个分量瞬时包络')%瞬时幅值（包络）：各个imf分量的瞬间频率
for i = 1:size(imf,1)
figure(6)
subplot(3,1,3)
plot(t/fs,A(i,:))
hold on
end
legend('1','2','3','4','5','6','7','8');
xlabel('时间 t/s');ylabel('幅值');title('imf各个分量瞬时包络')

六、参考文献

1、陈雨莺.基于 EMD 的说话人特征参数提取方法研究[D]. 湘潭大学,2015

2、信号处理：希尔伯特-黄变换

3、Hilbert-Huang Transform: matlab 希尔伯特-黄变换matlab实现 :

4、希尔伯特变换和瞬时频率问题–连载（二）

5、希尔伯特变换的本质是什么？

以上均为个人学习的总结，如有不妥，还望指教！