基于改进的EMD的数据包络分析

Halcom

基于改进的EMD的数据包络分析：
代码下载链接：https://pan.baidu.com/s/1WqcjVf2c50CYxy7-htP8ag
提取码：5ebi

1. function imf = emd2(x)
   
2. % Empiricial Mode Decomposition (Hilbert-Huang Transform)
   
3. % imf = emd(x)
   
4. % Func : findpeaks
   
5. 
   
6. x   = transpose(x(:));  % 转置
   
7. imf = [];
   
8. while ~ismonotonic(x) % 当x不是单调函数，分解终止条件
   
9.    x1 = x;
   
10.    sd = Inf; % 均值
    
11.    % 直到x1满足IMF条件，得c1
    
12.    while (sd > 0.1) || ~isimf(x1) % 当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件
    
13.       s1 = getspline2(x1);   % 上包络线
    
14.       s2 = -getspline2(-x1); % 下包络线
    
15.       x2 = x1-(s1+s2)/2;    % 此处的x2为文章中的h
    
16.       
    
17.       sd = sum((x1-x2).^2)/sum(x1.^2);
    
18.       x1 = x2;
    
19.    end
    
20.    
    
21.     imf  = [imf, x1'];
    
22.     x = x-x1;
    
23. end

复制代码

样条拟合改进如下：

1. function s = getspline2(x)
   
2. % 三次样条函数拟合成元数据包络线
   
3. N = length(x);
   
4. p = findpeaks(x);
   
5. % s = spline([0 p N+1],[0 x(p) 0],1:N);
   
6. if(length(p)>2)
   
7.     x1 = p(1);  y1 = x(x1);
   
8.     x2 = p(2);  y2 = x(x2);
   
9.     x3 = 0;  y3 = (y2-y1)./(x2-x1+eps).*(x3-x1)+y1;
   
10.     x4 = p(end);    y4 = x(x4);
    
11.     x5 = p(end-1);  y5 = x(x5);
    
12.     x6 = N+1;  y6 = (y5-y4)./(x5-x4+eps).*(x6-x4)+y4;
    
13.     s = spline([0 p N+1],[y3 x(p) y6],1:N);
    
14. else
    
15.     s = spline([0 p N+1],[0 x(p) 0],1:N);
    
16. end

复制代码

代码运行效果：
https://www.bilibili.com/video/B ... id_from=333.999.0.0