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

Hilbert-Huang 变换 希尔伯特-黄转换希尔伯特-黄转换（Hilbert-Huang Transform），由台湾中央研究院院士黄锷（Norden E. Huang）等人提出，将欲分析资料分解为本质模态函数（intrinsic mode functions, IMF），这样的分解流程称为经验模态分解(Empirical Mode Decomposition, EMD)的方法。然后将IMF作希尔伯特转换(Hilbert Transform)，正确地获得资料的瞬时频率。此方法处理对象乃针对非稳态与非线性讯号。与其他数学转换运算（如傅立叶变换）不同，希尔伯特-黄转换算是一种应用在数据资料上的算法，而非理论工具。质模态函数(IMF)任何一个资料，满足下列两个条件即可称作本质模态函数。 局部极大值(local maxima)以及局部极小值(local minima)的数目之和必须与零交越点(zero crossing)的数目相等或是最多只能差1，也就是说一个极值后面必需马上接一个零交越点。 在任何时间点，局部最大值所定义的上包络线(upper envelope)与局部极小值所定义的下包络线，取平均要接近为零。因此，一个函数若属于IMF，代表其波形局部对称于零平均值。此类函数类似于弦波（sinusoid-like），但是这些类似于弦波的部分其周期与振幅可以不是固定。因为，可以直接使用希尔伯特转换，求得有意义的瞬时频率。经验模态分解(EMD)EMD算法流程图建立IMF是为了满足希尔伯特转换对于瞬时频率的限制条件之前置处理，也是一种转换的过程。我们将IMF来做希尔伯特转换可以得到较良好的特性，不幸的是大部分的资料并不是IMF，而是由许多弦波所合成的一个组合。如此一来，希尔伯特转换并不能得到正确的瞬时频率，我们便无法准确的分析资料。为了解决非线性（non-linear）与非稳态（non-stationary）资料在分解成IMF时所遇到的困难，便发展出EMD。经验模态分解是将讯号分解成IMF的组合。经验模态分解是借着不断重复的筛选程序来逐步找出IMF。以讯号为例，筛选程序的流程概述如下:步骤 1: 找出中的所有局部极大值以及局部极小值，接着利用三次样条(cubic spline)，分别将局部极大值串连成上包络线与局部极小值串连成下包络线。步骤 2: 求出上下包络线之平均，得到均值包络线。步骤 3: 原始信号与均值包络线相减，得到第一个分量。步骤 4: 检查是否符合IMF的条件。如果不符合，则回到步骤1并且将当作原始讯号，进行第二次的筛选。亦即重复筛选次直到符合IMF的条件，即得到第一个IMF分量，亦即步骤 5: 原始讯号减去可得到剩余量，表示如下式步骤 6: 将当作新的资料，重新执行步骤1至步骤5，得到新的剩余量。如此重复次.当第个剩余量已成为单调函数(monotonic function)，将无法再分解IMF时，整个EMD的分解过程完成。原始讯号可以表示成个IMF分量与一个平均趋势(mean trend)分量的组合，亦即如此一来，原始资料便分解成n个IMF和一个趋势函数，我们便可将IMF做希尔伯特转换来进行瞬时频率的分析。结论傅立叶变换是将一个讯号分解成无限多个弦波来分析资料，但是希尔伯特-黄转换则是将一个讯号分解成数个近似于弦波的讯号（周期、振幅不固定）和一个趋势函数来做分析。两者各有其优缺点，整理如下优点：1.避免复杂的数学运算2.可分析频率会随时间变化的讯号3.较适于分析气候、经济等具有趋势的资料4.可以找出一个函数的趋势缺点：1.缺乏严谨的物理意义2.需要复杂的递回，运算时间反而比短时距傅立叶变换要长3.希尔伯特转换未必能正确计算出本质模态函数之瞬时频率4.无法使用快速傅立叶变换5.只有在特例（组合较简单的资料）时使用希尔伯特-黄转换较快2.Matlab代码/matlabcentral/fileexchange/196812.1 EMD分解代码（emd.m）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 imfend+1 = x1; x = x-x1;endimfend+1 = x; % FUNCTIONS function u = ismonotonic(x) u1 = length(findpeaks(x)*length(findpeaks(-x);if u1 0, u = 0;else, u = 1; end function u = isimf(x) N = length(x);u1 = sum(x(1:N-1).*x(2:N) 1, u = 0;else, u = 1; end function 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) 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(imfk.*imfk); th = angle(hilbert(imfk); dk = diff(th)/Ts/(2*pi);endu,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,dk,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 =