HHT--希尔伯特黄变换

weixin

在希尔伯特变换中我们已经知道，直接对解析信号做瞬属性的分析时，瞬时相位对时间的导数就是瞬时频率。但是瞬时频率直接根据解析信号这么按公式求，是没有物理意义的！因为"瞬时频率"会出现很多负频率，如下图所示为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;   

        end

end　

%非主函数，被调用。判断当前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;   

       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

来源：CAE 进行时微信公众号（ID：gh_bbbde3e9698b），作者：即墨。