希尔伯特-黄变换说明及程序.docx

Hilbert-Huang 变换希尔伯特-黄转换希尔伯特-黄转换(Hilbert-Huang Transform),由台湾中央研究院院士黄锷 (Norden E. Huang)等人提出，将欲分析资料分解为本质模态函数(intrinsicmode functions, IMF)，这样的分解流程称为经验模态分解(Empirical Mode Decomposition, EMD)的方法然后将 IMF 作希尔伯特转换(Hilbert Transform)， 正确地获得资料的瞬时频率此方法处理对象乃针对非稳态与非线性讯号与其 他数学转换运算(如傅立叶变换)不同，希尔伯特-黄转换算是一种应用在数据 资料上的算法，而非理论工具质模态函数(IMF)任何一个资料，满足下列两个条件即可称作本质模态函数1. 局部极大值(local maxima)以及局部极小值(local minima)的数目之和必须 与零交越点(zero crossing)的数目相等或是最多只能差1，也就是说一个极值 后面必需马上接一个零交越点2. 在任何时间点，局部最大值所定义的上包络线(upper envelope)与局部极小 值所定义的下包络线，取平均要接近为零。

因此，一个函数若属于IMF，代表其波形局部对称于零平均值此类函数类似于 弦波(sinusoid-like)，但是这些类似于弦波的部分其周期与振幅可以不是固 定因为，可以直接使用希尔伯特转换，求得有意义的瞬时频率经验模态分解(EMD)EMD算法流程图建立IMF是为了满足希尔伯特转换对于瞬时频率的限制条件之前置处理，也是一 种转换的过程我们将IMF来做希尔伯特转换可以得到较良好的特性，不幸的是 大部分的资料并不是IMF，而是由许多弦波所合成的一个组合如此一来，希尔 伯特转换并不能得到正确的瞬时频率，我们便无法准确的分析资料为了解决非线性(non-linear)与非稳态(non-stationary)资料在分解成IMF时所遇到的 困难，便发展出EMD经验模态分解是将讯号分解成IMF的组合经验模态分解是借着不断重复的筛选 程序来逐步找出IMF以讯号•、'U ：为例，筛选程序的流程概述如下：步骤1 :找出,「厂中的所有局部极大值以及局部极小值，接着利用三次样条(cubic spline)，分别将局部极大值串连成上包络线与局部极小值串连成下包络 线步骤2 :求出上下包络线之平均，得到均值包络线；山'「。

步骤3 :原始信号「"与均值包络线相减，得到第一个分量，\ '"hi (t) = S (t) — 7H1 (f)步骤4 :检查％^是否符合IMF的条件如果不符合，则回到步骤1并且将当作原始讯号，进行第二次的筛选亦即姻(f) = hi (t) 一 (t)重复筛选上次hk(t) = h— (i) - mk (i)直到0、仃符合IMF的条件，即得到第一个IMF分量「I仃:，亦即Q (£) = hk (f)步骤5 :原始讯号减去、仃:可得到剩余量「\'"，表示如下式n(i) = s@)—勺⑴步骤6 :将「I 口:当作新的资料，重新执行步骤1至步骤5,得到新的剩余量 ,七门：如此重复、次侦£)=『1⑴一血(C々(£) = R (土)— 4 (t)以⑴=以一1 (£) 一 % (t)当第”个剩余量，，已成为单调函数(monotonic function)，将无法再分解IMF 时，整个EMD的分解过程完成原始讯号:可以表示成"个IMF分量与一个平 均趋势(mean trend)分量的组合，亦即1（"=立"）+『"） fc = l如此一来，原始资料便分解成n个IMF和一个趋势函数，我们便可将IMF做希尔 伯特转换来进行瞬时频率的分析。

结论傅立叶变换是将一个讯号分解成无限多个弦波来分析资料，但是希尔伯特-黄转 换则是将一个讯号分解成数个近似于弦波的讯号（周期、振幅不固定）和一个趋 势函数来做分析两者各有其优缺点，整理如下优点：1. 避免复杂的数学运算2. 可分析频率会随时间变化的讯号3. 较适于分析气候、经济等具有趋势的资料4. 可以找出一个函数的趋势缺点：1. 缺乏严谨的物理意义2. 需要复杂的递回，运算时间反而比短时距傅立叶变换要长3. 希尔伯特转换未必能正确计算出本质模态函数之瞬时频率4. 无法使用快速傅立叶变换5. 只有在特例（组合较简单的资料）时使用希尔伯特-黄转换较快2.Matlab 代码 EMD 分解代码(emd.m)function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform) % imf = emd(x)% Func : findpeaksx = 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;endimf(end+1} = x1;x = x-x1;endimf(end+1} = x;% FUNCTIONSfunction u = ismonotonic(x)u1 = length(findpeaks(x))*length(findpeaks(-x));if u1 > 0, u = 0;else, u = 1;endfunction u = isimf(x)N = length(x);u1 = sum(x(1:N-1).*x(2:N) < 0);u2 = length(findpeaks(x))+length(findpeaks(-x));if abs(u1-u2) > 1, u = 0;else, u = 1endfunction s = getspline(x)N = length(x);p = findpeaks(x);s = spline([0 p N+1],[0 x(p) 0],1:N);2.2 找极值代码(findpeaks.m)function n = findpeaks(x)% Find peaks.% n = findpeaks(x)n = find(diff(diff(x) > 0) < 0);u = find(x(n+1) > x(n));n(u) = n(u)+1;2.3绘制时-频曲线以及尺度分解代码(plot_hht.m)function plot_hht(x,Ts)% Plot the HHT.% plot_hht(x,Ts)% Ts : time interval (sec)% :: Syntax% The array x is the input signal and Ts is the sampling period.% Func : emd% Get HHT.imf = emd(x);for k = 1:length(imf)b(k) = sum(imf{k}.*imf{k});th = angle(hilbert(imf{k}));d{k} = diff(th)/Ts/(2火pi);end[u,v] = sort(-b);b = 1-b/max(b);% Set time-frequency plots.N = length(x);c = linspace(0,(N-2)*Ts,N-1);for k = v(1:2)figure, plot(c,d{k},'k.','Color',b([k k k]),'MarkerSize',3);set(gca,'FontSize',8,'XLim',[0 c(end)],'YLim',[0 1/2/Ts]);xlabel('Time'), ylabel('Frequency');end% Set IMF plots.M = length(imf);N = length(x);c = linspace(0,(N-1)*Ts,N);for k1 = 0:4:M-1figurefor k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]); endxlabel ('Time');end。