在希尔伯特变换中我们已经知道,直接对解析信号做瞬属性的分析时,瞬时相位对时间的导数就是瞬时频率。但是瞬时频率直接根据解析信号这么按公式求,是没有物理意义的!因为"瞬时频率"会出现很多负频率,如下图所示为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),作者:即墨。 |
GMT+8, 2024-12-25 16:28 , Processed in 0.057563 second(s), 23 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.