|
基于改进的EMD的数据包络分析:
代码下载链接:https://pan.baidu.com/s/1WqcjVf2c50CYxy7-htP8ag
提取码:5ebi
- function imf = emd2(x)
- % Empiricial Mode Decomposition (Hilbert-Huang Transform)
- % imf = emd(x)
- % Func : findpeaks
- 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 = getspline2(x1); % 上包络线
- s2 = -getspline2(-x1); % 下包络线
- x2 = x1-(s1+s2)/2; % 此处的x2为文章中的h
-
- sd = sum((x1-x2).^2)/sum(x1.^2);
- x1 = x2;
- end
-
- imf = [imf, x1'];
- x = x-x1;
- end
复制代码 样条拟合改进如下:
- function s = getspline2(x)
- % 三次样条函数拟合成元数据包络线
- N = length(x);
- p = findpeaks(x);
- % s = spline([0 p N+1],[0 x(p) 0],1:N);
- if(length(p)>2)
- x1 = p(1); y1 = x(x1);
- x2 = p(2); y2 = x(x2);
- x3 = 0; y3 = (y2-y1)./(x2-x1+eps).*(x3-x1)+y1;
- x4 = p(end); y4 = x(x4);
- x5 = p(end-1); y5 = x(x5);
- x6 = N+1; y6 = (y5-y4)./(x5-x4+eps).*(x6-x4)+y4;
- s = spline([0 p N+1],[y3 x(p) y6],1:N);
- else
- s = spline([0 p N+1],[0 x(p) 0],1:N);
- end
复制代码 代码运行效果:
https://www.bilibili.com/video/B ... id_from=333.999.0.0
|
|