HHT，MEMD和信号处理

简介

本文系笔者在科研学习中学习到和总结到的有关HHT和MEMD在EEG信号处理上的可能用途和使用方法，有一些原理性的东西需要强大的数理知识基础，所以只谈谈对于这些知识的可行理解，欢迎大家批评指正，私聊或评论均可。谢谢！

HHT(Hilbert-Huang Transform) 希尔伯特-黄变换

希尔伯特-黄变换是由黄锷先生提出的一种对非线性、非稳态信号的解析方法，这极大地缓解了当时对非线性、非稳态信号束手无策的境况。从算法上来说，针对非稳态的振荡信号来说，Huang(EMD)可以分解出信号的本征模态(获得一些具有局部特征的信号波形)，这部分更具重要性，只是在使用谱分析手段分析数据时需要用到HT，所以，这可能是这个算法这样命名的由来吧。好了！闲话不多说，直接进入主题吧。

IMF(Intrinsic Mode Functions) 内涵(或者叫固有)模态函数

IMFs其实是在EMD分解过程中产生一些原始信号的特征模态的集合，你可以这样理解IMF，当所有的IMF叠起来(或者说拟合起来)，就变成了一个原始信号。IMF具有以下的特点：

1. IMF 在整个数据集上的极值点和过零点的数据只能相等或者相差为 1。
2. 在任意点，局部极大值的包络(envlope)和局部极小值的包络的均值为 0。
   它代表了与简单调和函数相对应的一种普遍简单的振荡模式。根据定义，IMF 是任何具有相同数目的极值点和零点交叉点的函数，其包络线对称于零。这一定义保证了 IMF 的希尔伯特转换表现良好。至于包络是什么,在下一小节会有具体提到。

EMD(Empirical Mode Decomposition) 经验模态分解

区别于其它几种常见的信号处理方法，如傅里叶变换(实现了时域信号到频域信号的互换)，小波变换(实现了从时域信号到时频域(又称作小波域)的变换)，关于小波变换，我这里自荐一下这篇blog浅谈一下小波变换和EEG信号处理 ，而EMD不仅实现了在时域信号上的信号特征提取，而且具备数据驱动性和自适应性，而不像小波变换一样，需要特定的小波函数或者小波基。

以下是一段关于EMD算法流程，大家先看一下，可以有一个大概的思路。
标准的EMD算法

Standard EMD algorithm:

1. 在原数据信号$x(t)$找出振荡波型的极值集合$x^{{}^{\prime }}(t)$。
2. 对极大值 (极小值也一样) 使用三次样条插值的方法获取振荡波形的上包络 $e_{max}(k)$(极小值是下包络 $e_{min}(k)$)；
3. 计算均值 $m(k)=\frac{1}{2}[e_{max}(k)+e_{min}(k)]$
4. 做差：$s(k)=x^{{}^{\prime }}(k)-m(k)$ 得一个振荡波形信号，若该信号满足 IMF 的停止原则，则将该信号 $s(k)$ 作为一个(或者说是第i个)IMF 分量$d_i(k)$，然后进入下一步；否则返回第二步继续进行计算；
5. 将得到的IMF分量$d_i(k)$作为新的信号$x(t)$继续进行分解运算，直至满足EMD的停止条件。

这里简单解释一下为什么是使用插值方法来获取包络线，因为一个自然信 号波形肯定是极其震荡的，这样就不能保证一个信号具有的相同数目的极大值或极小值了，那么就无法保证能够计算出m(k)了。所以需要使用插值方法对数据进行补足。

相信大家也都看到了我标粗的几个词，除了为什么要用三次样条插值外，另外两个也都还没有解释，先放在这，留个悬念，稍后在MEMD中仔细展开，但请一定注意，这两个停止原则是不同的，不能混为一谈。

这是一幅来自维基百科的EMD过程图，能够很好的帮助大家理解和记忆。当然，如果大家愿意去看维基百科理解HHT是更好的，也是更快的。

现在，我先简单的填一下前面的一个坑：关于EMD和IMF的停止原则

1. IMF的停止条件： 通常来说，IMF的停止条件有很多种，最为常见的有四种，而其中最常见的就是由EMD算法的发明人黄锷先生提出来的H98和H99;
   H98: H98是由黄锷先生于1998年提出的基于柯西收敛准则的一种方法，这种方法可能因为计算的复杂性，较之于H99没那么常用；
   H99: 同样是由黄锷先生于1999年提出的，这种停止原则在预设最大计算次数，同时在最大计算次数以内，当满足极值点数目和零点数目相等或相差为1时结束，这一句话是不是很熟悉，因为这就是IMF的性质之一！
   2.EMD的停止条件： EMD最终停止运算是当最后的IMF有单调性的时候，才会停止运算。
   下面的实验过程中，为了方便对照，我并没有严格按照EMD的停止条件进行实验，而是预设了一个最大的IMF数量，所以各位需要引起注意。

EMD 在EEG信号的应用——寻找各个模态的特征

以下是我使用PyEMD工具包对波恩大学癫痫脑电数据集进行EMD分解的结果：分两种情况，但两种情况都是在预设IMF数量的前提下进行的，而非严格的符合EMD的停止条件。

这里的IMFs的数量为7，最后一个IMF并不满足单调性条件，只是因为在使用EMD分解时，使用的是预设最大IMFs数量，而不是EMD停止满足条件。

而这一幅图(下图)是经过数据切片后得到的数据进行EMD操作，刚好在IMFs数量为7(也可以说是6，最后一个一般称为残差项(residual))。

具体结果，可通过下面的python代码进行实现。

代码如下

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD
import  matplotlib as mpl
import os


mpl.use('TKAgg')
#不使用这个的话，PyCharm可能会由bug
# Reference:https://blog.csdn.net/qq_51248682/article/details/129396134

def GenerateData(path):
    #生成数据
    eeg_signal = np.loadtxt(path)
    return eeg_signal
def EMDFunction(signal):
    emd = EMD()
    #对EEG信号进行emd分解
    imfs = emd(signal,max_imf=6)

    #绘图
    plt.figure(figsize=(7,30))
    plt.subplot(len(imfs) + 1, 1, 1)
    # subplot(a,b,c) 参数a,b,c分别表示 行 列 图像的位置(1,2,3,...)
    plt.plot(signal)
    plt.title('Raw Data')

    for i ,imf in enumerate(imfs):
        plt.subplot(len(imfs)+1,1,i+2)
        plt.plot(imf)
        plt.title('IMF{}'.format(i+1))

    plt.tight_layout()
    plt.show()
def main():
    #获取path
    path = "D:\\PythonCode\\FirstPaper_SDP\\UBonn_eegdatasets_sampled\\A_Z\\Z001.txt"
    EEG_data = GenerateData(path)

    # 切片
    # EEG_data = EEG_data[:250]
    # print(np.array(EEG_data).shape)

    EMDFunction(EEG_data)

if __name__ == "__main__":
    main()

另外，需要注意的是，笔者在download python工具包PyEMD的时候，会出现有两个分析包重叠的情况，即用于地球测绘的pyemd和PyEMD冲突的情况，这里可能需要进入conda的envs目录下修改自己下好的用于EMD分析的数据包的目录名字，改为PyEMD即可。

HT(Hilbert Transform)变换和 Gabor变换

关于希尔伯特变换，其实可以理解为一种谱分析的手段，这可能也是黄锷先生为什么会在自己的算法中加上它的原因。
希尔伯特变换会使实信号转变成虚信号，而且两个信号之间的相位差为$\frac{\pi }{2}$，两个信号会具有相同的瞬时频率，那么，Gabor就通过三角变换，将相位差变成了$\sin$和$\cos$变换，这就有了如下的关系。
假设$$y(t)=H(x(t)),H表示希尔伯特变换,即y(t)表示x(t)的希尔伯特变换$$
而，对应的解析信号(经过Gabor变换)为$$a(t)=x(t)+j\cdot y(t){=}^{欧拉公式}A(t)e^{j\theta (t)}$$
则，对应的赋值为
$$A(t)=\sqrt{x^2(t)+y^2(t)}$$
对应的相位为
$$\theta (t)=\arctan (\frac{y(t)}{x(t)})$$
所以，对应的瞬时频率为
$$\omega (t)=\frac{d\theta (t)}{dt}$$
关于HT的一般定义如下：
下面是关于HT变换的公式定义及相关性质。

函数$u(t)$的希尔伯特变换定义为
$$H(u)(t)=\frac{1}{\pi }p\cdot v\cdot {\int }_{-\infty }^{+\infty }\frac{u(t)}{t-\tau }d\tau$$
其中，$p\cdot v\cdot$表示柯西主值(the Cauchy principal value)，也就是瑕积分的积分值。也可以用$P$替换。

上式也可以转换成
$$H(u)(t)=\frac{2}{\pi }\underset{ϵ\to 0}{\lim }{\int }_{ϵ}^{\infty }\frac{u(t+\tau )-u(t-\tau )}{2\tau }d\tau$$

从定义上来看，某函数的希尔伯特变换就是该函数和函数$h=\frac{1}{\pi t}$的卷积

希尔伯特变换可以看作是在单位圆上将信号旋转$\frac{\pi }{2}$所致，所以，希尔伯特变换存在逆变换即为$H^3$。可以这样理解，希尔伯特变换将信号转换成了解析信号的形式输出，而解析信号的额模即为包络线。

希尔伯特变换与傅里叶变换(Fourier Transform)

在频率域上，对函数$f$进行希尔伯特变换，实质为在$f$的傅里叶变换结果$F(\omega )$前乘一个系数$H(\omega )=-jsgn(f)$

即相当于把$F(\omega )$的所有正频率向后移动$\frac{\pi }{2}$相位,所有负频率向前移动$\frac{\pi }{2}$相位，即如下：
KaTeX parse error: Undefined control sequence: \notag at position 23: …) = \hat{f}(t),\̲n̲o̲t̲a̲g̲ ̲\\ = f(t) \a…
根据卷积定理，时域上的卷积等于频域上的相乘。
F { f ( t ) ∗ g ( ( t ) ) } = F { f ( t ) } ⋅ F { g ( t ) } F\{f(t) \ast g((t))\} = F\{f(t)\} \cdot F\{g(t)\} F{f(t)∗g((t))}=F{f(t)}⋅F{g(t)}
KaTeX parse error: Undefined control sequence: \notag at position 125: …ay} \right. \̲n̲o̲t̲a̲g̲ ̲
其中H表示希尔伯特变换，F表示傅里叶变换。

HSA(Hilbert Specturm Analysis) 希尔伯特谱分析

希尔伯特谱分析其实是基于两个式子
由
$$\omega =\frac{d\theta (t)}{dt},这个式子已经在前面完成了推导。$$

和希尔伯特变换可以推导出
$$X(t)=\sum _{j=1}^n{a}_j(t)exp(i\int {\omega }_j(t)dt)$$
这样，$X(t)$其实就是一个关于时间和频率的函数，那么就可以在时频图上去分析信号了，这就是希尔伯特谱分析了，这也是希尔伯特·黄变换的一个重要部分。

震荡信号的简单谱分析(HHT)

HHT具有以下的特征：

1. HHT反映的信号的频谱特征随时间变化的规律。
2. HHT可以对局部特征进行反映，这点主要得益于EMD的作用。EMD可以自适应地进行时频局部化分析，有效提取原信号的特征信息。
   以下是关于希尔伯特-黄变换分析低频震荡信号的一个例子，主要参考于HHT算法实现（基于Python）。代码我就不贴了，需要的话直接点上面的链接获取就行。
   以下是一个结果。
   

EMD拓展——MEMD(Multivariate Empirical Mode Decomposition) 多元经验模态分解

这个算法来自于一篇10年的Paper(Multivariate empirical mode decomposition BY N. REHMAN* AND D. P. MANDI)，虽然很老了，但在很多地方都被用到，很值得一学。这个算法比对每个维度的信号使用EMD算法具有以下的优点：

1. 各个维度获取到的IMFs数量相同，不存在一个维度多或者一个维度少的情况；
2. 在计算均值线$m(t)$的时候，使用的是全部信号的上下包络线，考虑到各个通道之间可能会存在的影响，很具有借鉴性。
   首先让大家看个结果，然后再具体的讲述一下算法过程及需要理解的前置知识。
   以下是使用MEMD对一个随机生成的六元信号进行MEMD分解的结果：
   原始六元信号
   
   Channel1结果
   
   Channel2结果
   
   Channel3结果
   
   Channel4结果
   

Channel5结果

Channel6结果

从运行结果上来看，各个方向上(channels)，在经过MEMD分解后得到的最终的IMFs的数量最终是一致的。那么为什么会有这样的结果呢？
下面将从算法的角度分析一下MEMD算法，需要了解的前置知识主要有MC方法和多元样条插值的方法。

MEMD算法

在已知采样获得的方向向量 x k θ = { x 1 k , x 2 k , … , x n k } x^θ_k =\{x_1^k,x_2^k,…,x_n^k \} xkθ​={x1k​,x2k​,…,xnk​}, 和代表n元信号的n个成分(没有明确说n个成分是什么(从仿真来看，可能就是n个通道)，但正好我们需要处理的EEG,MEG信号都是多通道的)的n维输入向量 v ( t ) t = 1 T = { v 1 ( t ) , v 2 ( t ) , … , v n ( t ) } v(t)_{t=1}^T=\{v_1 (t),v_2 (t),…,v_n (t)\} v(t)t=1T​={v1​(t),v2​(t),…,vn​(t)}.
算法如下：

1. 将输入信号$v(t{)}_{t=1}^T$投影到方向向量$x_k^{\theta }$上,得k个投影$p_k^{\theta }(t){\}}_{k=1}^K$；此时值(并不是复信号等类似信号了)是可比的。
2. 通过多元样条插值方法获得K条包络线$e_k^{\theta }(t){\}}_{k=1}^K$
3. 通过K条包络线计算均值$m(t)=\frac{1}{K}{\sum }_{k=1}^K{e}_k^{\theta }(t)$ 。(用在脑电、脑磁信号上似乎很合理，信号之间具有联系)
4. 获取‘detail’参数$d(t)=x(t)-m(t)$,如果$d(t)$满足停止原则,就对$x(t)-d(t)$继续上一步步骤，否则对$d(t)$。

那么问题来了？采样的方向向量如何获得呢？多通道信号如何获得其包络线呢？下面我就论文中提到的相关方法做一个简单的介绍和分析。

Quasi-MC(Monte Carlo) 拟蒙特卡罗方法

相信大家都多少学习过，或者听过蒙特卡罗方法(Monte Carlo)，蒙特卡罗方法其实就是使用均匀分布的思想，求解一些直接求解很难得问题，最著名的就是美国曼哈顿原子弹计划中求解高阶玻尔兹曼方程的问题。那么拟蒙特卡罗和蒙特卡洛有什么区别呢？为什么要学习蒙特卡洛算法呢？这篇Paper为什么要使用QMC呢？我们带着问题来继续下面的部分。

MC 和 QMC

首先，我先向大家讲一个关于MC方法的故事，以方便大家理解：
蒙特卡洛算法的原理只是通过随即均匀采样估算某类数据，而最早的蒙特卡洛思想出现在布
丰投针计算圆周率 π 的实验中，直到二战时期，由于美国曼哈顿原子弹计划的需要，在计算中子在原子弹中的增殖和
扩散过程中，需要求解高维玻尔兹曼方程 (高阶偏微分方程 (或者说是在高维空间中的积分))。

假设要计算积分$\int f(x)g(x)dx$,可以这样理解，这个积分其实是满足g(x)特征的f(x)分布的期望值$E_{g(x)}[x]$，而根据大数定律，当满足f(x)的样本数量足够多的时候，样本$g(x_1)\dots g(x_n)$的期望值$E[g(x)]$与$E_{g(x)}[x]$近似相等。即需要获取足够多满足f(x)的样本，通过f(x)可以得到其cdf $F(x)\in [0,1]$,那么可以在[0,1]上做均匀采样获取$i,i=F(x)$,那么$x=F^{-1}(i)$。

但如果存在$f(x)$也难积分的情况呢？这就需要找到一个分布$q(x),s.t.mq(x)\ge f(x)$,并且$q(x)$可积分。 但此时又存在一个情况，有些地方的采样可能满足$f(x)$但不一定满足$q(x)$,就需要用到接受拒绝采样，以$p=\frac{f(x)}{mq(x)}$接受采样结果。否则则以$1-p$的概率拒绝采样。

上面就是一个标准的蒙特卡洛的接受拒绝采样的过程。

下面是一个使用MC方法求解问题的一个例子。

求解积分${\int }_0^1\sqrt{(1-x^2)}$
这个问题显然可以通过三角代换求解出积分值为$\pi$，但如果使用MC方法呢？(当然这里仅是为了验证MC方法的正确性，在实际问题中不会存在这么简单的积分)
$${\int }_0^1\sqrt{1-x^2}={\int }_0^{1}dx{\int }_0^{f(x)}dy=P(Y\le f(X))$$
即对 X 在 (0,1) 上采样，对 Y 在 (0,1) 上采样，满足式子
$$x^2+y^2\le 1$$
最后，用满足条件的样本点除以总的样本点个数，即可获得相应的积分近似值。

代码实现如下

import random
import math


def montecarlo(N):
    i = 0
    count =0
    while i<=N:
        x = random.random()
        y = random.random()
        if pow(x,2) + pow(y,2) <1:
            count+=1
        i+=1
    pi = 4*count /N
    print(pi)
    print(math.pi)

if __name__ == "__main__":
    montecarlo(80000)
# »3.14975  程序输出结果
# »3.141592653589793  math.pi 值

MC扩展——MH算法(Metropolis-Hastings)

MH算法其实之上是MCMC(Markov Chain Monte Carlo)中的一种，而MH算法具备以下优点：
马尔可夫-蒙特卡罗方法 (MCMC) 较一般的 MC(Accept-Reject) 方法具有以下的优点：

1. 基于马尔可夫过程，采样速度更快；
2. 采样的接受率相对来说更高，MH 的采样接受率还未达到 1，Gibbs 采样接受率为 1.
   从马尔可夫过程出发，一个马尔可夫过程必定满足以下关系：
   $$P(x_t|x_{t-1},x_{t-2},\dots ,x_1)=P(x_t|x_{t-1})$$
   假设最终采样的结果序列服从一个分布$\pi (x)$ , 而序列产生的过程是依照 Markov Chain 的，所以，序列中前后两个采
   样结果必定满足一个概率分布 $K(x^{t-1}\to x^t)$，为以示区别，令前后两个采样结果分别表示为 $x$ 和$x^{\ast }$，两者服从同一个分布 $\pi (x)$,所以，一定存在以下关系
   $${\pi }_t(x)={\int }_{\pi }{\pi }_{t-1}(x)K(x\to x^{\ast })dx$$
   而使上式成立的一个充分条件是

$$DetailBalanceEquation：\pi (x^{\ast })K(x^{\ast }\to x)=\pi (x)K(x\to x^{\ast })$$
只需将细致平衡方程(Detail Balance Equation)带入条件即可验证。
而MH(Metropolis-Hastings)算法即如下：
输入：抽样的目标分布的密度函数p(x),函数f(x) \
输出：p(x)的随机样本$x_{m+1},x_{m+2},\dots ,x_n$，函数样本均值$f_{mn}$；

1. 初始化 $x^{(0)}$
2. for i=0 to N-1
   u $\sim$ U(0,1)
   $x^{\ast }\sim q(x^{\ast }|x^{(i)})$
   if $u<\min (1,\frac{\pi (x^{\ast })q(x|x^{\ast })}{\pi (x)q(x^{\ast }|x)})$
   $x^{(i+1)}=x^{\ast }$
   else
   $x^{(i+1)}=x^{(i)}$

其中$q(x,x^{\ast })$表示建议分布，通常有以下两种形式：

1. Metropolis选择，假设建议分布是对称的，即对任意的$x$和$x^{{}^{\prime }}$有
   $$q(x,x^{{}^{\prime }})=q(x^{{}^{\prime }},x)$$
   那么，接受分布$\alpha (x,x^{{}^{\prime }})$可以简化成 min ⁡ { 1 , p ( x ′ ) p ( x ) } \min \{ 1, \frac{p(x^{'})}{p(x)}\} min{1,p(x)p(x′)​}
2. 独立抽样。假设$q(x,x^{{}^{\prime }})$与当前状态x无关，即$q(x,x^{{}^{\prime }})=q(x^{{}^{\prime }})$。
   那么，接受分布为 α ( x , x ′ ) = min ⁡ { 1 , w ( x ′ ) w ( x ) } \alpha(x,x^{'}) = \min \{1,\frac{w(x^{'})}{w(x)} \} α(x,x′)=min{1,w(x)w(x′)​}，
   其中$w(x^{{}^{\prime }})=p(x^{{}^{\prime }})/q(x^{{}^{\prime }}),w(x)=p(x)/q(x)$。

Halton序列与Hammersley序列

Halton序列是一种低差异度的无偏采样点生成算法。它是基于一个简单的观察：为了获得一个均匀分布的样本点集合，我们可以使用一组不同的质数作为基础，将每个维度映射到一个序列中。

具体来说，在Halton序列中，我们选择一组互不相同的质数$x_1,\dots ,x_n$作为基础，其中n是数据点的维度。则在1维Halton序列中第i个样本为
$$r_i^x=\sum _{k=0}^s\frac{a_k}{x^{k+1}}$$
其中m是预先定义的整数，$a_k$是一个介于0和$x_i$之间的整数，用于确定该维度上的n个样本点中的第k个点。
$$i=\sum _{k=0}^s{a}_k{x}^k$$

其获取到的采样序列为$(r_i^{x_1},r_i^{x_2},\dots ,r_i^{x_n})$

Hammersley序列
这种采样的个数n是事先知道的，其获取的第i次样本为
$(i/n,r_i^{x_1},r_i^{x_2},\dots ,r_i^{x_{n-1}})$

两者采样结果对比如下

(a)表示序列算法采样图，(b)表示Hammersley序列算法结果图。

关于Halton和Hammersley可以参考如下链接加深理解低差异序列 (low-discrepancy sequences)之Halton序列均匀产生多维随机数的介绍与实现。

Multivariate Spine Interpolation 多元样条插值

样条插值在每个间隔中使用低阶多项式，并选择多项式以使它们平滑地吻合在一起。结果函数被称为样条曲线。
多元插值是函数变量多余一个时所使用的插值方法，因为变量在空间坐标系中，所以，这种插值方法又叫做空间插值 (spatial interpolation)。

从某种意义上来说，插值方法其实是一种变相的拟合，不过这种拟合是一种及其过度的拟合，需要拟合到所有数据的特征，然后从拟合函数上找到除样本点外的一些数据点(或者称为未知点)。

一般的插值方法包括但不仅限于以下几种：
$\circ$ 线性插值:

线性插值是用一系列首尾相连的线段依次连接相邻各点，每条线段内的点的高度作为插值获得的高度值。利用插值点与其他数据点(区间内)的斜率值与区间端点的梯度值进行插值；

一般有如下形式，设$x_1,x_2$为插值区间，$x_k$为插值点，
有
$$\frac{x_1-x_2}{y_1-y_2}=\frac{x_1-x_k}{y_1-y_k}$$

$\circ$ 多项式插值(拉格朗日插值)

对n个样本点过拟合出一个最高次可能为n-1的多项式，一次性拟合，就比如下图：

使用插值后得到的一个5阶多项式为
$$f(x)=-0.0001521x^6-0.003130x^5+0.07321x^4-0.3577x^3+0.2255x^2+0.9038x$$
而这种方法在加入新的数据点时，会使得函数更新一次，并且可能会导致龙格现象。

所以可以所用样条插值的方法，例如使用三次样条插值：
$$f(x)=\{\begin{array}{l}-0.1522x^3+0.9937x,ifx\in [0,1],\\ -0.01258x^3-0.4189x^2+1.4126x-0.1296,ifx\in [1,2],\\ 0.1403x^3-1.3359x^2+3.2467x-1.3623,ifx\in [2,3],\\ 0.1579x^3-1.4945x^2+3.7225x-1.8381,ifx\in [3,4],\\ 0.05375x^3-0.2450x^2-1.2756x+4.8259,ifx\in [4,5],\\ -0.1871x^3+3.3673x^2-19.3370x+34.9282,ifx\in [5,6].\end{array}$$

MEMD分解多通道信号

最后，基于MEMD算法进行相关实验，即可获得本部分最初的实验结果。
MEMDPython算法代码来自于GithubMultivariate Empirical Mode Decomposition algorithm
以下是我写的例子代码：

# -*- coding: utf-8 -*-
"""
@author: zl
"""

import numpy as np
import matplotlib.pyplot as plt
from MEMD_all import memd


from scipy.io import wavfile

# Define six audio signals
sig1 = np.random.uniform(-1, 1, size=200)
sig2 = np.random.uniform(-1, 1, size=200)
sig3 = np.random.uniform(-1, 1, size=200)
sig4 = np.random.uniform(-1, 1, size=200)
sig5 = np.random.uniform(-1, 1, size=200)
sig6 = np.random.uniform(-1, 1, size=200)

# Combine the signals into a single six-channel signal
six_channel_signal = np.array([sig1, sig2, sig3, sig4, sig5, sig6])
print(six_channel_signal.shape)
#
# fig,axis = plt.subplots(nrows=6,figsize=(10,15))
# # print(six_channel_signal[0].shape)
# for i in range(6):
#     axis[i].plot(six_channel_signal[i])
#     axis[i].set_title('Channel{}'.format(i+1))
#
# plt.tight_layout()
# plt.show()

# inp = np.loadtxt("D:/PythonCode/FirstPaper_SDP/UBonn_eegdatasets_sampled/E_S/S001.txt")
imf = memd(six_channel_signal)
imf_x = imf[0,0,:]
imf_y = imf[:,1,:]
imf_z = imf[:,2,:]
print(np.array(imf_x).shape)
#
for j in range(6):
	fig,axi = plt.subplots(nrows=9,figsize=(10,15))
		for i in range(9):
			imf_i = imf[i,j,:]
			axi[i].plot(imf_i)
			axi[i].set_title('IMF{}'.format(i+1))
			plt.tight_layout()
			plt.show()

最后，就能获得本部分开头的实验结果了。

最后，关于HHT即相关应用的简单讲述到此结束了，如果其中有不足的地方或者有什么疑问的话，欢迎大家私聊或留言。谢谢！