<>HHT原理

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

经验模式分解方法是一种自适应的、高效的数据分解方法。由于这种分解是以局部时间尺度为基础的,因此它适应于非线性、非平稳过程。通过分解,原始的数据序列可用IMF分量
ci(t)c_i(t)ci​(t)以及剩余项 rn(t)r_n(t)rn​(t)表示:
x(t)=∑i=1nci(t)+rn(t)(1) x(t)=\sum_{i=1}^n c_i(t)+r_n(t)\tag{1} x(t)=i=1∑n​ci​(
t)+rn​(t)(1)
EMD 将信号x(t)x(t)x(t)分解为 n 个 IMF,对每个 IMF 分 量 即ci(t)c_i(t)ci​(t)
作HT,继而可求取每个IMF的瞬时频率和瞬时幅值信息。这类本征模态函数的瞬时频率(Instantaneous Frequency, IF)有着明确的物理意义。
yi(t)=1π∫−∞∞ci(t′)t−t′dt′(2) y_i(t)=\frac{1}{π}\int_{-\infty}^\infty
\frac{c_i(t')}{t-t'}dt'\tag{2}yi​(t)=π1​∫−∞∞​t−t′ci​(t′)​dt′(2)
ci(t)c_i(t)ci​(t)和yi(t)y_i(t)yi​(t)构成解析信号z(t)z(t)z(t):
zi(t)=ci(t)+jyi(t)=aiejθi(t)(3) z_i(t)=c_i(t)+jy_i(t)=a_i e^{j\theta
_i(t)}\tag{3}zi​(t)=ci​(t)+jyi​(t)=ai​ejθi​(t)(3)
其中ai(t)=ci(t)2+yi(t)2a_i(t)=\sqrt{c_i(t)^2+y_i(t)^2}ai​(t)=ci​(t)2+yi​(t)2​
代表瞬时幅值,θi=arctan(yi(t)/ci(t))\theta _i=arctan(y_i(t)/c_i(t))θi​=arctan(yi​(t)/ci
​(t))代表瞬时相位,ωi(t)=dθi(t)/dt\omega _i(t)=d\theta _i(t)/dtωi​(t)=dθi​(t)/dt代表瞬时频率。

由瞬时幅值ai(t)a_i(t)ai​(t)和瞬时频率ωi(t)\omega _i(t)ωi​(t)可将信号表示为:
x(t)=∑i=1nai(t)ej∫ωi(t)dt(4) x(t)=\sum_{i=1}^n
a_i(t)e^{j\int\omega_i(t)dt}\tag{4}x(t)=i=1∑n​ai​(t)ej∫ωi​(t)dt(4)
上式省略了剩余项rn(t)r_n(t)rn​(t),因为rn(t)r_n(t)rn​(t)
幅值小,不是一个单调函数就是一个常数,对信号分析和信号提取没有实质影响。在时间-频率面上画出每个IMF 以其幅值加权的瞬时频率曲线,这个时间-频率分布谱
图就是Hilbert谱,记为H(ω,t)H(ω,t)H(ω,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)=∫H2(ω,t)dω(5) E(t)=\int H^2(\omega,t)d\omega\tag{5} E(t)=∫H2(ω,t)dω(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("中值滤波后基音跟踪曲线");
输入的语音信号为“我正在学习语音信号处理”,执行结果如下:

技术
©2019-2020 Toolsou All rights reserved,
SpringBoot JpaRepository 数据库增删改查关于Bellman-Ford算法的个人理解(精华)2020年6月26日 C#类库 异常处理帮助类C#/.NET 系统优化专题(redis第六篇 数据结构【List】)vue使用THREE.js创建一个可以控制的立方体SpringBoot 与JPA结合中 JpaRepository 里自定义查询(精华)2020年6月26日 C#类库model PageInput(精华)2020年7月13日 微信小程序 页面间通信【JAVA】【华为校园招聘笔试-软件】2020-09-09 mysql无备份恢复