有关emd分解的matlab源程序

把 G.Rilling 写的 EMD.m注释汉化了，同时根据个人理解添加了部分注释，最后附上一个EMD分解＋HHT谱＋边际谱的例子。例子程序中所用到的函数是G.Rilling提供的

% EMD 计算经验模式分解
%
%
%   语法
%
%
% IMF = EMD（X）
% IMF = EMD（X...‘Option_name‘Option_value...）
% IMF = EMD（XOPTS）
% [IMFORTNB_ITERATIONS] = EMD（...）
%
%
%   描述
%
%
% IMF = EMD（X） X是一个实矢量，计算方法参考[1]，计算结果包含在IMF矩阵中，每一行包含一个IMF分量，
% 最后一行是残余分量，默认的停止条件如下[2]：
%
%   在每一个点 mean_amplitude < THRESHOLD2*envelope_amplitude （注：平均幅度与包络幅度的比值小于门限2）
%   &
%   mean of boolean array {（mean_amplitude）/（envelope_amplitude） > THRESHOLD} < TOLERANCE 
%  （注：平均幅度与包络幅度比值大于门限的点数占信号总点数中的比例小于容限）
%   &
%   |#zeros-#extrema|<=1 （注：过零点和极值点个数相等或者相差1）
%
% 这里 mean_amplitude = abs（envelope_max+envelope_min）/2 （注：平均幅度等于上下包络相互抵消后残差的一半的绝对值，理想情况等于0）
% 且 envelope_amplitude = abs（envelope_max-envelope_min）/2 （注：包络幅度等于上下包络相对距离的一半，理想情况等于上下包络本身的绝对值）
% 
% IMF = EMD（X） X是一个实矢量，计算方法参考[3]，计算结果包含在IMF矩阵中，每一行包含一个IMF分量，
% 最后一行是残余分量，默认的停止条件如下[2]：
%
%   在每一个点 mean_amplitude < THRESHOLD2*envelope_amplitude（注：平均幅度与包络幅度的比值小于门限2）
%   &
%   mean of boolean array {（mean_amplitude）/（envelope_amplitude） > THRESHOLD} < TOLERANCE
%  （注：平均幅度与包络幅度比值大于门限的点数占信号总点数中的比例小于容限）
%
% 这里平均幅度和包络幅度的定义与前面实数情况下类似
%
% IMF = EMD（X...‘Option_name‘Option_value...） 设置特定参数（见选项）
%
% IMF = EMD（XOPTS） 与前面等价，只是这里OPTS是一个结构体，其中每一个域名与相应的选项名称一致。
%
% [IMFORTNB_ITERATIONS] = EMD（...） 返回正交指数
%                       ________
%         _  |IMF（i:）.*IMF（j:）|
%   ORT = \ _____________________
%         /
%         -       || X ||^2        i~=j
%
% 和提取每一个IMF时进行的迭代次数。
%
%
%   选择
%
%
%  停止条件选项:
%
% STOP: 停止参数 [THRESHOLDTHRESHOLD2TOLERANCE]
% 如果输入矢量长度小于 3 只有第一个参数有效，其他参数采用默认值
% 默认值: [0.050.50.05]
%
% FIX （int）: 取消默认的停止条件，进行  指定次数的迭代
%
% FIX_H （int）: 取消默认的停止条件，进行  指定次数的迭代，仅仅保留 |#zeros-#extrema|<=1 的停止条件，参考 [4]
%
%  复 EMD 选项:
%
% COMPLEX_VERSION: 选择复 EMD 算法（参考[3]）
% COMPLEX_VERSION = 1: “algorithm 1“
% COMPLEX_VERSION = 2: “algorithm 2“ （default）
% 
% NDIRS: 包络计算的方向个数 （默认 4）
% rem: 实际方向个数 （根据 [3]） 是 2*NDIRS
% 
%  其他选项:
%
% T: 采样时刻 （线性矢量） （默认: 1:length（x））
%
% MAXITERATIONS: 提取每个IMF中，采用的最大迭代次数（默认：2000）
%
% MAXMODES: 提取IMFs的最大个数 （默认: Inf）
%
% DISPLAY: 如果等于1，每迭代一次自动暂停（pause）
% 如果等于2，迭代过程不暂停 （动画模式）
% rem: 当输入是复数的时候，演示过程自动取消
%
% INTERP: 插值方法 ‘linear‘ ‘cubic‘ ‘pchip‘ or ‘spline‘ （默认）
% 详情见 interp1 文档
%
% MASK: 采用 masking 信号，参考 [5]
%
%
%   例子
%
%
% X = rand（1512）;
%
% IMF = emd（X）;
%
% IMF = emd（X‘STOP‘[0.10.50.05]‘MAXITERATIONS‘100）;
%
% T = linspace（0201e3）;
% X = 2*exp（i*T）+exp（3*i*T）+.5*T;
% IMF = emd（X‘T‘T）;
%
% OPTIONS.DISLPAY = 1;
% OPTIONS.FIX = 10;
% OPTIONS.MAXMODES = 3;
% [IMFORTNBITS] = emd（XOPTIONS）;
%
%
%   参考文献
%
%
% [1] N. E. Huang et al. “The empirical mode decomposition and the
% Hilbert spectrum for non-linear and non stationary time series analysis“
% Proc. Royal Soc. London A Vol. 454 pp. 903-995 1998
%
% [2] G. Rilling P. Flandrin and P. Goncalves
% “On Empirical Mode Decomposition and its algorithms“
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03 Grado （I） June 2003
%
% [3] G. Rilling P. Flandrin P. Goncalves and J. M. Lilly.
% “Bivariate Empirical Mode D