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代码 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 imf{end+1} = x1; x = x-x1;endimf{end+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) < 0);u2 = length(findpeaks(x))+length(findpeaks(-x));if abs(u1-u2) > 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) < 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-1 figure for k2 = 1:min(4,M-k1), subplot(4,1,k2), plot(c,imf{k1+k2}); set(gca,'FontSize',8,'XLim',[0 c(end)]); end xlabel('Time');end。