weixin 发表于 2022-10-12 15:52

HHT--希尔伯特黄变换

在希尔伯特变换中我们已经知道,直接对解析信号做瞬属性的分析时,瞬时相位对时间的导数就是瞬时频率。但是瞬时频率直接根据解析信号这么按公式求,是没有物理意义的!因为"瞬时频率"会出现很多负频率,如下图所示为HT求解的瞬时频率!希尔伯特-黄变换 (HHT) 可以解决这一问题,对非稳态信号做时频分析!
简单来说,理解HHT变换,主要是经验模式分解 (EMD)+希尔伯特谱分析两部分,目的是获得信号中具有实际物理意义的瞬时频率分量,进而实现高分辨率的时频分析。

1998年,Norden E. Huang(黄锷:中国台湾海洋学家)等人提出了EMD方法,并引入了Hilbert谱的概念和Hilbert谱分析的方法,美国国家航空和宇航局 (NASA) 将这一方法命名为Hilbert-Huang Transform,即希尔伯特-黄变换。

HHT主要内容包含两部分,第一部分为EMD;第二部分为Hilbert谱分析 (Hilbert Spectrum Analysis,简称HSA)。简单说来,HHT处理非平稳信号的基本过程是:

首先,利用EMD方法将给定的信号分解为若干本征模态函数 (IMF),利用这类函数的局部特性,可以使得函数在任意一点的瞬时频率都有意义。

然后,对每一个IMF进行HT,得到相应的Hilbert谱,即将每个IMF表示在联合的时频域中。

每一个固有模态函数可以表示为:

对其进行希尔伯特变换,把采样点时间、瞬时频率、瞬时振幅放在3维空间里做时频谱分析了。最后,汇总所有IMF的Hilbert变换就会得到原始信号的希尔伯特谱H(w,t),下图为HHT图。

HHT的优势
它的关键就在于"经验模态分解"这一步。因为这种信号的分解,完全是"方法适应于数据"的,即不同的数据,方法会根据数据的情况做变换和调整,而不是始终用固定的那几种基函数(小波变换、短时傅里叶变换)。这种自适应思想下的工具,一定是非常优越的!

经验模态分解何时结束?
经验模态分解就是为了分解出一系列的固有模态函数,所以当剩余数据已经无法再分出来固有模态函数时,经验模态分解这一步就彻底结束了。

如何判断剩余数据无法再分解了?
剩余数据极大值或极小值点数目有一个为0时(均值条件没法满足了),或剩余数据已经是单调时(没法做包络线了),就可以结束了。

以下是 EMD算法的Matlab源程序,仅供参考!
%主函数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%非主函数,被调用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%非主函数,被调用<br>%判断x是否单调,返回0代表不是单调,返回1代表是单调function u = ismonotonic(x)    u1 = length(findpeaks(x))*length(findpeaks(-x));%如果最大/最小值有一个为0即可判断程序满足退出条件了    if u1 > 0      u = 0;    else      u = 1;          endend %非主函数,被调用。判断当前x是不是真IMFfunction 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));%求极值点个数    ifabs(u1-u2) > 1      u = 0;    else      u = 1;           endend%非主函数,被调用,作用是获得x的包络线function s = getspline(x)    N = length(x);    p = findpeaks(x);    s = spline(,,1:N);end

来源:CAE 进行时微信公众号(ID:gh_bbbde3e9698b),作者:即墨。

页: [1]
查看完整版本: HHT--希尔伯特黄变换