Hilbert-Huang希尔伯特黄变换（HHT）的Matlab实现

代码中的plot_hht函数实现了Hilbert-Huang变换(HHT)，HHT就是Hilbert-Huang Transform（希尔伯特黄变换），美国航天航空局黄鳄教授发明的，大概的过程是：先对信号进行经验模态分解（Empirical mode decomposition －EMD），得出本征模态函数（IMF intrinsic mode function），再对本征模态函数进行希尔伯特变换，从而过进一步得该信号的希尔伯特谱、时频能量谱等，以便对信号进行分析，据介绍对非线性及非平稳信号有较好的分析和处理效果。

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）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件      958087  2007-01-25 14:59  HHT.pdf
     文件         965  2008-04-24 13:36  plot_hht.m
     文件         893  2008-03-13 11:49  emd.m
     文件         146  2007-08-26 00:42  findpeaks.m
     文件      102892  2007-04-28 15:57  Hum.wav