【语音信号处理五】希尔伯特黄变换

【参考资料】
【1】https://www.cnblogs.com/qiweiwang/archive/2010/12/20/1911736.html
【2】https://zhuanlan.zhihu.com/p/74413218
【3】https://www.pianshen.com/article/9161269068/
【4】《希尔伯特变换与信号的包络_瞬时相位和瞬时频率》 https://wenku.baidu.com/view/9aa65710580216fc700afd57.html

1. 概述

所谓HHT 希尔伯特黄变换是将信号进行EMD分解，将分解完的每个IMF分量用Hibert变换得到其瞬时频谱，最终进行信号的时频分析。因此在本笔记中重点记录hilbert变换和EMD分解。

2. Hilbert变换与瞬时频率

2.1 复信号的定义

已知有欧拉公式：

e j θ = c o s θ + j s i n θ e^{j \theta} = cos \theta + j sin \theta ejθ=cosθ+jsinθ

由公式 θ = ω t = 2 π T t \theta = \omega t = \dfrac{2 \pi}{T}t θ=ωt=T2πt 得到下述复变函数公式：

e j ω = c o s ω t + j s i n ω t e^{j\omega} = cos \omega t + j sin \omega t ejω=cosωt+jsinωt 从而反推得到：

c o s ω t = e j ω + e − j ω cos \omega t = e^{j\omega} + e^{-j\omega} cosωt=ejω+e−jω

由于上面和式中如果根据欧兰公式展开，其虚部是互相抵消的，因此我们可以将实信号表示成：

c o s ω t = R e s { e j ω } cos \omega t = Res\{ e^{j\omega} \} cosωt=Res{ejω} ,即 e j ω e^{j\omega} ejω 这个信号的实部，此时我们称为 e j ω e^{j\omega} ejω为其复信号

2.2 实信号和复信号的频谱关系

推导详见【参考文献4】，结论为：

称 q ( t ) q(t) q(t)为 x ( t ) x(t) x(t)的复信号，其中 Q ( f ) Q(f) Q(f)和 X ( f ) X(f) X(f)分别为其频谱，这有：

Q ( f ) = { 2 Q ( f ) f > 0 0 f < 0 Q(f)=\left\{ \begin{aligned} 2Q(f)& & f > 0 \\ 0 & & f < 0 \end{aligned} \right. Q(f)={2Q(f)0f>0f<0

针对这种关系，我们可以定义一个滤波器,其频谱为 H 1 ( f ) = { 2 f > 0 0 f < 0 H_1(f)=\left\{ \begin{aligned} 2 & & f > 0 \\ 0 & & f < 0 \end{aligned} \right. H1(f)={20f>0f<0

同时得到其对应的时间函数： h 1 ( t ) = δ ( t ) + i 1 π t h_1(t) = \delta(t) + i \dfrac{1}{\pi t} h1(t)=δ(t)+iπt1
备注：这步不知道怎么推出来的?

2.3 Hilbert变换的定义

由2.2存在公式（频域的乘法为时域的卷积），我们有：

q ( t ) = h 1 ( t ) ∗ x ( t ) q(t) = h_1(t) * x(t) q(t)=h1(t)∗x(t)

q ( t ) = （ δ ( t ) + i 1 π t ） ∗ x ( t ) q(t) = （\delta(t) + i \dfrac{1}{\pi t}） * x(t) q(t)=（δ(t)+iπt1）∗x(t) 得到

q ( t ) = x ( t ) + i 1 π t ∗ x ( t ) q(t) = x(t) + i \dfrac{1}{\pi t} * x(t) q(t)=x(t)+iπt1∗x(t)

令 x ~ ( t ) = 1 π t ∗ x ( t ) = 1 π ∫ − ∞ ∞ x ( τ ) t − τ d τ \widetilde{x}(t) = \dfrac{1}{\pi t} * x(t) = \dfrac{1}{\pi} \int_{-\infty}^{\infty} \dfrac{x(\tau)}{t - \tau} d \tau x (t)=πt1∗x(t)=π1∫−∞∞t−τx(τ)dτ

我们称 x ~ ( t ) \widetilde{x}(t) x (t) 为Hilbert变换

这里的Hilbert变换相当于做了一个滤波，频谱滤波为： H 1 ( f ) = { i f > 0 − i f < 0 H_1(f)=\left\{ \begin{aligned} i & & f > 0 \\ -i & & f < 0 \end{aligned} \right. H1(f)={i−if>0f<0

举例：

2.4 解析信号、包络、瞬时相位和瞬时频谱

利用实信号 x ( t ) x(t) x(t)及其hilbert变换 x ~ ( t ) \widetilde{x}(t) x (t)构成一个解析信号：

q ( t ) = x ( t ) + x ~ ( t ) q(t) = x(t) + \widetilde{x}(t) q(t)=x(t)+x (t)

假设有一工程中常用的窄带信号：

x ( t ) = a ( t ) c o s ( ω t + φ ( t ) ) x(t) = a(t) cos(\omega t + \varphi(t)) x(t)=a(t)cos(ωt+φ(t)) 及其hilbert变换 x ~ ( t ) = a ( t ) s i n ( ω t + φ ( t ) ) \widetilde{x}(t) = a(t) sin(\omega t + \varphi(t)) x (t)=a(t)sin(ωt+φ(t))推导详见参考文献4 略

我们令 θ ( t ) = ω t + φ ( t ) \theta(t) = \omega t + \varphi(t) θ(t)=ωt+φ(t) 则有如下一些核心定义：

包络： e ( t ) = x ( t ) 2 + x ~ ( t ) 2 e(t) = \sqrt{x(t)^2 + \widetilde{x}(t)^2} e(t)=x(t)2+x (t)2

瞬时相位： θ ( t ) = a r c t g x ~ ( t ) x ( t ) \theta(t) = arctg \dfrac{\widetilde{x}(t)}{x(t)} θ(t)=arctgx(t)x (t)

瞬时频率： u ( t ) = d θ ( t ) d t = d d t a r c t g x ~ ( t ) x ( t ) u(t) = \dfrac{d \theta (t)}{dt} = \dfrac{d}{dt} arctg \dfrac{\widetilde{x}(t)}{x(t)} u(t)=dtdθ(t)=dtdarctgx(t)x (t) 是瞬时相位的导数

备注：F(t) = sin(2πft + φ)：f就是频率；2πft + φ 就是相位；对时间的导数就是2πf，频率是相位的变换速度

2.5 python代码

import numpy as np
import pylab as pl
import scipy.signal as signal
import matplotlib.pyplot as plt
from scipy import fftpack
import math
from tftb.processing import inst_freqt     = np.arange(0, 0.1, 1/20000.0)
x     = 4*t*np.sin(2*np.pi*200*t + 3.0) + 0.1
hx    = fftpack.hilbert(x)
e     = np.sqrt(x**2 + hx**2)
u1, _ = inst_freq(hx)t = []
for i in range(len(x)):t.append(math.atan(hx[i]/x[i]))u = []
for i in range(len(t)-1):u.append(t[i+1] - t[i])plt.figure(figsize=(10,7))
plt.subplot(3, 1, 1)
plt.plot(x,'r') #hilbert变换
plt.subplot(3, 1, 1)
plt.plot(hx,'g') #hilbert变换
plt.subplot(3, 1, 1)
plt.plot(e,'b') #包络曲线
plt.subplot(3, 1, 2)
plt.plot(t,'r') #瞬时相位
plt.subplot(3, 1, 3)
plt.plot(u,'r') #瞬时频率
plt.show()

备注：这里是按照自己的理解去实现离散Hilbert变换的瞬时相位和频率，但和 inst_freq 似乎不太对？？？

3. 经验模态分解EMD分解

3.1 EMD分解的基本步骤：

步骤1：

找出原始信号 x ( t ) x(t) x(t)局部极大值，用三次B样条构筑上包络曲线
找出原始信号 x ( t ) x(t) x(t)局部极小值，用三次B样条构筑下包络曲线

步骤2：

计算上包络曲线和下包络曲线的均值 m 1 m_1 m1
计算原始曲线与包络均值的差： h 1 = x ( t ) − m 1 h_1 = x(t) - m_1 h1=x(t)−m1

步骤3：

判断 h 1 h_1 h1是否满足IMF条件

备注：IMF条件需要满足两点：
2.1 经过分解得到的时间序列，其极值点数量与过0点数量相差不超过1个
2.2 经过插值方式拟合的极值包络线在分析时间序列内的均值为0

如果不满足IMF条件，这重复步骤1、步骤2，使得k步迭代后的 h 1 , k ( t ) h_{1,k}(t) h1,k(t)作为第一个IMF分量 c 1 ( t ) = h 1 , k ( t ) c_1(t)=h_{1,k}(t) c1(t)=h1,k(t)

步骤4：

计算剩余信号 r 1 ( t ) = x ( t ) − c 1 ( t ) r_1(t) = x(t) - c_1(t) r1(t)=x(t)−c1(t)

步骤5：

针对剩余信号 r ( t ) r_(t) r(t)重复步骤1~3,求得其余的IMF分量 r ( t ) … r N ( t ) r_(t) \dots r_{N}(t) r(t)…rN(t)
最终原始信号被分解为IMF分量和剩余信号的和，即 x ( t ) = ∑ n = 1 N c n ( t ) + r N ( t ) x(t) = \sum_{n=1}^N c_n(t) + r_N(t) x(t)=∑n=1Ncn(t)+rN(t)

备注：在其他资料中，步骤3的k步迭代是否满足IMF条件，采用公式：

要求SD的值在0.2~0.3之间

3.2 python代码


import pyhht
from pyhht.visualization import plot_imfst = np.arange(0, 0.1, 1/20000.0)
x = 4*t*np.sin(2*np.pi*200*t + 3.0) + 0.1
decomposer = pyhht.emd.EMD(x)
imfs = decomposer.decompose()
plot_imfs(x, imfs)