www.pudn.com > HHT_power-system_power-quality_disturbances-detect > emdyy.m, change:2010-04-26,size:1133b


function imf = emd(x)    %EMD 
 
 
x   = transpose(x(:));    %转置 
imf = []; 
while ~ismonotonic(x)     %当x不是单调函数,分解终止条件 
   x1 = x; 
   sd = Inf;               %直到x1满足IMF条件,得c1 
   while (sd > 0.1) || ~isimf(x1)   %当标准偏差系数sd大于0.1或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; 
 
 
 
function u = ismonotonic(x) 
 
u1 = length(findpeaks(x))*length(findpeaks(-x)); 
if u1 > 0, u = 0; 
else      u = 1; end            %u=0表示x不是单调函数,u=1表示x为单调的 
 
 
 
 
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    %u=0表示x不是固有模式函数,u=1表示x是固有模式函数 
 
function s = getspline(x)       %三次样条函数拟合成元数据包络线 
 
N = length(x); 
p = findpeaks(x); 
s = spline([0 p N+1],[0 x(p) 0],1:N);