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

在我们正式开始讲解Hilbert-Huang Transform之前,不放先来了解一下这一伟大算法的两位发明人和这一算法的应用领域

Section I 人物简介

  希尔伯特:公认的数学界“无冕之王”,1943年去世于瑞士苏黎世。除此之外,自不必过多介绍。

  黄锷:1937年出生于湖北省;1975年进入NASA(美国国家宇航局);美国国家工程院院士。

Section II Hilbert-Huang的应用领域

  医学领域:探测心率不齐、登革热的扩散、血压的变化

  交通领域:探测公路桥梁安全

  安全领域:辨识发言者的身份

  地理领域:地震工程

  航天领域:卫星资料分析

  在了解了这一伟大发明的背景后,下面我们要正式的开始入手希尔伯特-黄变换了,我将尝试以尽可能简要的语言向大家介绍这一发明,并尽可能的避免不必要的数学推导。

Section III Hilbert-Huang的算法详细介绍

    如下图所示,在希尔伯特-黄的运算步骤中,原始脑电信号/其他时序信号被作为Huang的算法的输入,在经过huang的算法处理过后被当做Hilbert的输入进行处理。这便是Hilbert-Huang最简单明了的运算步骤。在这里为了继续往下的讲解更加方便,我们先来介绍两个概念。上文中,我们提到了“huang的算法”,在正式的书面语言中,我们并不这么称呼它,而是将“huang的算法”称为EMD(Empirical mode decomposition,经验模式分解)。而另外一个概念IMF在这里直接讲解或许会使大家晕头转向(或许有人注意到图中的IMF后面有一个‘s’,而这里却没有加‘s’,对英语只有基础了解的人也应猜到IMF不止一个)。

当将到这里的时候,大部分人或许会萌生出一个念头——“难道Huang仅仅是对Hilbert的锦上添花吗”。好吧,至少本人当初就是这样想的,毕竟Hilbert比huang更早出名,而且Hilbert是数学史上公认的“大牛”,哦,不对,是“大王”。用当前时兴的话来说就是“huang有可能抱了Hilbert的大腿”。但当我真正了解了这一伟大的发明之后,我才彻彻底底打消了这个十分愚蠢的念头。

我个人并不喜欢吊人胃口,这里把结论说在前面“Huang的算法几乎是Hilbert使用的前提条件,Hilbert Transform则是Hilbert-Huang算法的精要所在”(注意句中出现了“几乎”一词)。下面我就给大家讲一下这句话的由来。比如我们造了一款叫做“榔头”的手机,“榔头”手机对用户的使用提出了下列要求:1.晚上不能使用。2.下雨天不能打。3.室内不能打。4.室外的偏远郊区也不能打。实际上,Hilbert正是这样一款“榔头”手机,它对用户的使用提出了近乎苛刻的要求。Hilbert变换算法要求输入信号只能是线性稳态的。请注意这里是两个词“线性”“稳态”。无论是在自然界还是在人类社会中,绝大部分的信号要么是“线性非稳态”,要么是“非线性稳态”,要么干脆是“非线性非稳态”。我们关心的重点——EEG信号正是这样一类“非线性非稳态”的信号。这也就导致了绝大部分信号不能够愉快的进入Hilbert的“碗里”来。此时,Huang的EMD算法起到了这样的作用,它能够将所有的时域信号转化为“线性稳态”,解了Hilbert算法的软肋。

首先,我们先说一说Huang的EMD算法。为了讲解清晰起见,我将对照下图予以讲解:

上图中,深蓝色的线条是EEG信号(截取自瑞士联邦理工学院DEAP数据库 s01 trail1 channel1的前200个数据点)。图中,红线上的红点是该EEG signal的极大值点,绿线上的绿点是该signal的极小值点。我们分别为极大值点和极小值点做三次包络线做好的包络线分别是红色包络线和绿色包络线两条线。为这两条线做出均值线即为图中围绕y=0轴(注意y=0轴的位置,并非是图中的坐标系的x轴,x轴所代表的线是y=-10)震荡的浅蓝色均值线。之后,我们将原始EEG信号减去均值线,得到疑似IMF线(图中未标出)(这里终于出现IMF了,可是我还是无法让你直观理解,大家先暂且忍忍,强行记忆一下)。之后,我们对疑似IMF进行判断(需要同时满足两个条件,下面讲)。如果满足条件,则疑似IMF升级为正式IMF。然后将原始信号减去正式IMF的结果赋值给原始信号,说白了就是让这一IMF从原始信号里“滚蛋”。就好比蒙面歌王的某一期的获胜者一样,都赢了,不滚干嘛,难道还要和加时赛选手(减完后剩下的原始信号,也即新原始信号)在一起吗?另一个方面,如果疑似IMF未能通过检验,则将当前IMF作为原始信号,并回到做极值点的包络线那一步重新开始。现在讲一下两个重要的条件:

                     条件1:均值线(总得有很多数构成吧)的平均值趋近于0(一般和0做差<0.1)

                     条件2:原始信号的极值点个数(包括极大值点个数+极小值点个数)和原始信号同y=0的交点个数之差不能大于1(小于等于1)

那么这样一个程序什么时候可以循环结束呢,答案是,当某一次IMF被发现是单调函数或者是缺少极大/小值点即可让程序结束。下图是程序流程图。

空口无凭,我们处理一段真正的脑电试试看(程序会在之后给出)

图中共有4*2个图,位于(1,1)这个位置的是脑电的原始信号。之后从(1,2)->(4,2)均为IMF。其中,除了(1,1)自身,每一副图都是(1,1)的一个IMF(现在知道什么是IMF了吧)。通过观察不难发现。一个典型的IMF分量的上下包络线肯定是对称的。最后一个IMF(4,2)被称为余项用r表示。观察即可知该IMF分量没有极小值点(端点除外),所以程序才会结束。通常来讲,别的书上会这样用数学公式告诉你:

其中ξ(t)就是原始信号,IMFi就是K个固有模态函数。rK就是原始信号减完IMF后剩下的余项。

下面是求解EMD算法的Matlab源程序,特此声明,本程序为我本人在网上找到的,除了注释外,其他版权皆归属原作者,由于不清楚原作者是谁,未能标出,如果侵犯权利,请联系我删除源码。

%非主函数,被调用
function n = findpeaks(x)%用于寻找极值点,该函数只会求极大值
%	Find peaks.
%   n = findpeaks(x)
    n = find(diff(diff(x)>0)<0);%一阶导数大于0二阶导数小于0的点
    u = find(x(n+1)>x(n));
    n(u) = n(u) + 1;
end
%非主函数,被调用%判断x是否单调,返回0代表不是单调,返回1代表是单调
function u = ismonotonic(x)
    u1 = length(findpeaks(x))*length(findpeaks(-x));%如果最大/最小值有一个为0即可判断程序满足退出条件了
    if u1 > 0
        u = 0;
    else
        u = 1;
    end
end 
%非主函数,被调用。判断当前x是不是真IMF
function u = isimf(x)
    N  = length(x);
    u1 = sum(x(1:N-1).*x(2:N) < 0);%求x与y=0轴交点的个数
    u2 = length(findpeaks(x))+length(findpeaks(-x));%求极值点个数
    if abs(u1-u2) > 1
        u = 0;
    else
        u = 1;
    end
end
%非主函数,被调用,作用是获得x的包络线
function s = getspline(x)
    N = length(x);
    p = findpeaks(x);
    s = spline([0 p N+1],[0 x(p) 0],1:N);
end
%主函数
function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% imf = emd(x)
% Func : findpeaks
    x = transpose(x(:));
    imf = [];
    while ~ismonotonic(x)
        x1 = x;
        sd = Inf;
        while (sd > 0.1) || ~isimf(x1)
            s1 = getspline(x1);
            s2 = -getspline(-x1);
            x2 = x1-(s1+s2)/2;

            sd = sum((x1-x2).^2)/sum(x1.^2);
            x1 = x2;
        end

        imf(end+1,:) = x1;
        x = x-x1;
    end
    imf(end+1,:) = x;
end

Section IV Hilbert算法的介绍

  在上一章中,我们介绍了EMD算法,在这一部分中,我会介绍Hilbert算法,这一节有些许数学趣味,对数学趣味不感兴趣的直接跳到应用部分。

 由最后一步可以知道,当频率大于0时,相位向左移90度;反之,向右移90度。这便是希尔伯特变换。

一般来讲,对于原始信号x(t)的希尔伯特变换H[x(t)],通常被写为

z(t)=x(t)+j H[x(t)]

其中,x(t)被称为复信号z(t)的实部,H[x(t)]被称为复信号z(t)的虚部, z(t)被称为x(t)的解析信号

一般情况下,matlab会将z(t)给出,而不直接给出原始信号的希尔伯特变换,所以需要使用imag函数求解z(t)的虚部,这才是真正的希尔伯特变换。

  

时间: 2024-10-12 07:26:36

Hilbert-Huang Transform(希尔伯特-黄变换)的相关文章

如何使用 css3 transform 属性来变换背景图?

本文和大家分享的主要是使用 css3 transform 属性来变换背景图相关内容,一起来看看吧,希望对大家学习css3有所帮助. 使用 css3 transform 属性可以轻易的旋转,倾斜,缩放任何元素.目前即使没有任何前缀也可以在绝大部分浏览器上很好的使用 . 如果你要在黑莓浏览器或者 UC 浏览器使用这个属性, 你需要加 -webkit- 前缀. #myelement { -webkit-transform: rotate(30deg); transform: rotate(30deg)

(转)View Transform(视图变换)详解

原文作者讲得太好了,唯有这篇让我对视图矩阵了解的清晰了很多. ---------------------------------------------------------------------------- 什么是View Transform 我们可以用照相机的原理来阐释3D图形的绘制过程,想象一下,我们在摄影的时候都需要做哪些工作,大致可分为如下几个步骤 摆放好待拍摄的物品,或者人物. 调整好拍摄角度. 调整焦距. 拍摄. 好了,来分析一下,上面的第一步就相当于世界变换了,将一个模型置

talib 中文文档(十四):Math Transform Functions 数学变换

Math Transform Functions ACOS - Vector Trigonometric ACos 函数名:ACOS 名称:acos函数是反余弦函数,三角函数 real = ACOS(close) ASIN - Vector Trigonometric ASin 函数名:ASIN 名称:反正弦函数,三角函数 real = ASIN(close) ATAN - Vector Trigonometric ATan 函数名:ASIN 名称:数字的反正切值,三角函数 real = ATA

形象易懂讲解算法I——小波变换

https://zhuanlan.zhihu.com/p/22450818?refer=dong5 最早发于回答:能不能通俗的讲解下傅立叶分析和小波分析之间的关系? - 咚懂咚懂咚的回答现收入专栏. 从傅里叶变换到小波变换,并不是一个完全抽象的东西,可以讲得很形象.小波变换有着明确的物理意义,如果我们从它的提出时所面对的问题看起,可以整理出非常清晰的思路. 下面我就按照傅里叶-->短时傅里叶变换-->小波变换的顺序,讲一下为什么会出现小波这个东西.小波究竟是怎样的思路.(反正题主要求的是通俗形

拉普拉斯变换-工程数学笔记

Laplace transform 通过拉普拉斯变换,可将以时间t为自变量的函数f(t)转化为以复数s为自变量的函数F(s),其逆变换称为拉普拉斯逆变换,即将F(s)变换为f(t),具体变换为: 常用的拉普拉斯变换如下: 当多个函数相乘时: 示例如下:

transform的2D部分,嗯…就这个标题了。

上一次写了transition的内容,这次就写拼写很类似的另外一个属性transform好了……我英语差这件事就不要吐槽了,下面是正文,真的: transition是过渡,transform是变换. transform分为2D变换和3D变换,简直碉堡了,其实3D变换就是比2D变换多了1D,可以简单这么理解,具体是不是等下次说3D的时候再说,这次只说2D. 在2D转换里我们可以实现斜切(skew),缩放(scale),旋转(rotate)以及位移(translate)元素的效果(还有一个矩阵-ma

WP8.1 UI 编程 六、变换特效和三维特效

1. 变换特效 变换原理:是二维变换矩阵 M11 M12 0 M21 M22 0 OffsetX OffsetY  1 WP只支持仿射变换,因此矩阵右边是0.0.1. (x,y,1)乘矩阵得到(x1,y1,1),新坐标为(x1,y1). 即:坐标(x,y)经矩阵变换后,新坐标为(x*M11 + y*M21 + OffsetX,x*M12 + y*M22 + OffsetY). WP提供了很多Transform类以变换对象,只需应用到UIElement的RenderTransform属性即可. 列

unity3D游戏开发之Transform的坐标变换注意事项

Transform是unity的核心类之一.表示的是物体的平移,旋转和缩放. 而position和localPosition, 分别表示的是,transform的位置是世界空间,和父空间的描述. 注意是 父空间,并不是自身空间. 注意到这点后,在空间的变换时就会省心很多了. 如果想搞清楚transform.position的变换过程,可以这样来测试: Java代码 //父空间转世界. Debug.Log(transform.position); Debug.Log(transform.paren

H.264 Transform

变换是视频.图像编码的核心部分.目前所采用的变换算法都是从傅里叶变换演变而来.单纯的变换并不会导致视频(图像)的码率变小,反而会增大.但是非常巧妙的一点是:变换把图像从空域转换成的时域,把由色块组成的图像变为由基准色调与图像细节组成:低频代表图片的基准色调,高频代表图像细节,类比电路中的基频与谐波.变换会使得图像的低频系数集中于某一点(左上角),频率向右下角递增.一般来说,4x4大小的图像大多只是颜色平缓的色块,不会有太多的细节,因此低频系数会较大,而高频系数较小.另外,人的眼睛对于高频系数,即