HHT原理

Hiibert Huang变换是由Huang等人于1998年提出来的一种信号分析方法，它主要由两个部分组成:经验模型分解(Empirical Mode Decomposition, EMD)和希尔伯特变换（Hilbert Transform，HT），其中EMD是核心。

经验模式分解方法是一种自适应的、高效的数据分解方法。由于这种分解是以局部时间尺度为基础的，因此它适应于非线性、非平稳过程。通过分解，原始的数据序列可用IMF分量 $c_i(t)$以及剩余项 $r_n(t)$表示：
$$x(t)=\sum _{i=1}^n{c}_i(t)+r_n(t)\text{\hspace{1em}(1)}$$
EMD 将信号$x(t)$分解为 n 个 IMF，对每个 IMF 分 量 即$c_i(t)$作HT，继而可求取每个IMF的瞬时频率和瞬时幅值信息。这类本征模态函数的瞬时频率（Instantaneous Frequency， IF）有着明确的物理意义。
$$y_i(t)=\frac{1}{\pi }{\int }_{-\infty }^{\infty }\frac{c_i(t^{\prime })}{t-t^{\prime }}d{t}^{\prime }\text{\hspace{1em}(2)}$$
$c_i(t)$和$y_i(t)$构成解析信号$z(t)$：
$$z_i(t)=c_i(t)+j{y}_i(t)=a_i{e}^{j{\theta }_i(t)}\text{\hspace{1em}(3)}$$
其中$a_i(t)=\sqrt{c_i(t{)}^2+y_i(t{)}^2}$代表瞬时幅值，${\theta }_i=arctan(y_i(t)/c_i(t))$代表瞬时相位，${\omega }_i(t)=d{\theta }_i(t)/dt$代表瞬时频率。

由瞬时幅值$a_i(t)$和瞬时频率${\omega }_i(t)$可将信号表示为：
$$x(t)=\sum _{i=1}^n{a}_i(t)e^{j\int {\omega }_i(t)dt}\text{\hspace{1em}(4)}$$
上式省略了剩余项$r_n(t)$，因为$r_n(t)$幅值小，不是一个单调函数就是一个常数，对信号分析和信号提取没有实质影响。在时间-频率面上画出每个IMF 以其幅值加权的瞬时频率曲线，这个时间-频率分布谱 图就是Hilbert谱，记为 $H(\omega ,t)$。由式(4)可以看出，Hilbert谱其实就是傅里叶变换 的一种扩展。与傅里叶变换中的常数幅值和固定频率相比较，式(4)具有时变的幅值和频率，它更能反应出信号的非线性和非平稳等特征信息。

基于HHT的基音检测方法

Hilbert-Huang变换适用于非线性非平稳信号处理， 基于HHT的基音检测方法流程图如下：

（1）带通滤波。由于人类的基音频率范围为60-500Hz，所以先将原始信号经过一个60-900Hz的带通滤波器做初步处理。

（2）分帧。虽然HHT方法不需要对语音信号做短时平稳的假设，因而不需要对语音信号做分帧加窗的处理，但语音数据的长度太长会影响EMD分解的效率，所以一般还是必须对语音信号分帧， 只是分帧的目不再是为了保证帧内数据的短时平稳。

（3）清浊音判定，得出需要进行基音检测的浊音部分。

（4）EMD分解。EMD分解是一种自适应的、具有高时频分辨能 力的数据分解方法. 通过 EMD分解,信号在每一个时间局部的频率结构被分离出来，即IMF分量。

（5）Hilbert变换。通过HT能够得出各个IMF分量的瞬时频率和瞬时振幅。

（6）计算信号的瞬时能量。可以由式(5)计算:
$$E(t)=\int H^2(\omega ,t)d\omega \text{\hspace{1em}(5)}$$
（7）对瞬时能量求一阶导数。在声门脉冲发生时刻,瞬时能量一定是增加的,，瞬时能量的一阶导数一定大于某个正数，因此，理论上我们就可以通过确定瞬时能量一阶导数的局部极大值点来确定声门脉冲发生的时刻,从而确定基音周期。

MATLAB编程实现

程序中用到了emd和hht函数，这在MATLAB2018a之后才被官方引进，本实验用的版本是2020a

单帧处理

首先我们选出一帧浊音帧进行处理，该帧长为100ms，用y0表示，它的时域波形如下：

执行以下语句：

emd(y0);

执行后就会显示y0的IMF分量如下：

通过以下语句就可以完成y0的Hilbert-Huang变换

imf=emd(y0);
[hs,f,t,imfinsf,imfinse]=hht(imf,fs); % imfinsf 瞬时频率 imfinse 瞬时能量

单独输入hht(imf,fs,'FrequencyLimits',[0 1600])，就可以显示y0的HHT谱

接下来对瞬时能量imfinse求一阶导数

deg=[0;diff(imfinse(:,1))]; % 瞬时能量差分
plot(deg);

结果如下：

继续进行阈值化处理，并找到局部最大值，这些局部最大值发生的时刻就认为是声门脉冲发生时刻

for i=1:length(deg)
    if(deg(i)<0.1e-3)
        deg(i)=0;
    end
end
deg=medfilt1(deg,3);
TF=islocalmax(deg,'MinSeparation',fs*0.005);
t=linspace(0,dT_ms/1000,length(deg));
plot(t,deg,t(TF),deg(TF),'r*')

结果如下：

最后求出这些局部最大值的时间差即可，这个时间差就可以认为是基音的周期，可计算上述帧信号的基音频率为104.9180Hz。可见使用HHT方法可以较好地提取浊音信号中的声门脉冲，接下来应用HHT对一段完整语音进行基音检测

基音检测

程序如下：

clear all;
clc;
[y,fs]=audioread('filename.wav');
n=length(y);
t=0:1/fs:(n-1)/fs;
hf=HF;
lf=LF;
y=filter(hf.Numerator,1,y); % 高通 60Hz
y=filter(lf.Numerator,1,y); % 低通 900Hz
figure()
subplot(3,1,1)
plot(t,y);
xlabel("时间/s");
ylabel("幅值");
title("原始信号");
% 分帧
dT_ms=50; % 一帧信号的时长
w=floor(dT_ms/1000*fs); % 一帧信号的样点数
overlap=w/2; % 帧移

frame=floor(n/(overlap)); % 帧数
Y=zeros(frame,w); % 存储分帧加窗后的信号
y=[y;zeros((frame+1)*overlap-n,1)]; % 补零
for i=1:frame
   interval=(i-1)*(overlap)+(1:w);
   Y(i,:)=y(interval).*hamming(w); % 分帧加窗
end

% 清浊音判定
voice=zeros(frame,1);
for i=1:frame
   y0=Y(i,:);
   imf=emd(y0);
   [hs,f,t]=hht(imf,fs);
   hs=full(hs);
   F=zeros(w,1);
   for m=1:size(hs,2)
       for n=1:size(hs,1)
            F(m)=F(m)+hs(n,m)*f(n);
       end
   end
   voice(i)=mean(F);
   if(voice(i)<10e-3)
       voice(i)=0;
   else
       voice(i)=1;
   end
end

% 基音检测
pitch=zeros(frame,1);
for i=1:frame
   y0=Y(i,:);
   imf=emd(y0); % EMD分解
   [hs,f,t,imfinsf,imfinse]=hht(imf,fs); % HT变换得到瞬时能量imfinse
   if(voice(i)==1) % 如果是浊音
        deg=[0;diff(imfinse(:,1))]; % 瞬时能量差分
        for k=1:length(deg)
            if(deg(k)<1e-6)
                deg(k)=0;
            end
        end
        deg=medfilt1(deg,3);
        TF=islocalmax(deg,'MinSeparation',fs*0.005); % 寻找局部最大值所在位置
        T=find(TF==1);
        T=diff(T);
        T=mean(T);
        T=T/fs;
        ff=1/T;
        pitch(i)=ff;
   else
       pitch(i)=0;
   end
end

t=linspace(0,(length(y)-1)/fs,length(pitch));
subplot(3,1,2)
plot(t,pitch);
xlim([min(t),max(t)]);
ylim([0,200]);
xlabel("时间/s");
ylabel("频率/Hz");
title("基音跟踪曲线");

pitch=medfilt1(pitch,3); % 中值滤波
subplot(3,1,3)
plot(t,pitch);
xlim([min(t),max(t)]);
ylim([0,200]);
xlabel("时间/s");
ylabel("频率/Hz");
title("中值滤波后基音跟踪曲线");

输入的语音信号为“我正在学习语音信号处理”，执行结果如下：