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

质模态函数(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/19681

2.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

imf{end+1} = x1;

x = x-x1;

end

imf{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);