基础版EMD分解函数Matlab程序

**
不调用matlab自带emd(x)函数，将其内容简化为如下部分
EMD分解基础步骤可以参见：意念回复：经验模态分解（EMD） 原始程序为百度搜索，结合ChatGPT后给出相应注释。


% 日期：2023.06.07
% 注释：调用子程序EMD分解基本结构（原始文件来源未知，根据内容添加注释）


clc
clear all

% 定义输入信号
Ts = (1/512);
Fs = 1/Ts;
t=0:Ts:1;
v1 = 0.5*cos(2*20*pi*t);
v2 = cos(2*30*pi*t);
input = v1+v2;
% x = awgn(input,1);
x=input;

%绘制出入信号时域图
figure(1)
plot(t,x);hold on;
grid on
xlabel('t');ylabel('x');title('初始信号');
legend('0.5cos(40x*pi)+cos(60x*pi)');

%求取本证模态函数imf并绘制其时图像及频谱
imf = emd(x);

[M, N]=size(imf);
figure(2);
subplot(M+1,1,1);
plot(x);hold on;
ylabel('x')
set(gca,'xtick',[])
axis tight;
title('合成信号x做EMD分解');

for i=2:M
    subplot(M+1,1,i);
    plot(imf(i-1,:),'r');hold on;xlim([100 400]);
    ylabel (['IMF ' num2str(i-1)]);
    set(gca,'xtick',[])
    xlim([1 length(x)])
end

subplot(M+1,1,M+1)
plot(imf(M,:));hold on;
ylabel(['RS ' num2str(M)])
xlim([1 length(x)])

%计算相关系数并绘制
 for i = 1:1:M
     cc(i)=min(min(corrcoef(imf(i,:),x)));
 end
figure(3)
plot(cc,'-g<','LineWidth',1.5,'MarkerEdgeColor','b','MarkerFaceColor','b','MarkerSize',5);
set(gca,'XGrid', 'on', 'YGrid', 'on');
legend('CC');
xlabel('IMF');
ylabel('相关系数');

% 重构信号
cg_ev=imf(1,:)+imf(2,:);
figure(4);
plot(t,x);
hold on
plot(t,cg_ev,'r');
xlabel('t/s');ylabel('幅值');legend('原信号','重构信号');
%


function imf = emd(x)

x= transpose(x(:));%转置为行矩阵
imf = [];

while ~ismonotonic(x) %当x不是单调函数，分解终止条件
    x1 = x;
    sd = Inf;%均值

    %直到x1满足IMF条件，得c1
    % 两个条件：）
    %（1）在任意时刻，由局部极大值点形成的上包络线和由局部极小值点形成的下包络线的平均值为零，
    %     即上、下包络线相对于时间轴局部对称
    %（2）在整个时程内极值点个数与过零点个数相等或最多相差1

    while (sd > 0.1) || ~isimf(x1) %当标准偏差系数sd大于0.1或x1不是固有模态函数时，分量终止条件

        s1 = getspline(x1);%上包络线
        s2 = -getspline(-x1);%下包络线

        x2 = x1-(s1+s2)/2;%此处的x2为中间信号，即为链接中的C1,1，(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;% 余量
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FUNCTIONS
function u = ismonotonic(x)
%u=0表示x不是单调函数，u=1表示x为单调的
u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0
    u = 0;
else
    u = 1;
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = isimf(x)
%u=0表示x不是固有模式函数，u=1表示x是固有模式函数
N  = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);% 确定过零点个数：受限判断相邻两个点的成绩，小于0则表示一正一副，置1，反之置0.
u2 = length(findpeaks(x))+length(findpeaks(-x));% 确定极值点个数
if abs(u1-u2) > 1 % 在整个时程内极值点个数与过零点个数相等或最多相差1
    u = 0;
else
    u = 1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function s = getspline(x)
%三次样条函数拟合成元数据包络线
N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);%spline(x,y,xq)，x 坐标； y - x 坐标处的函数值； xq - 查询点
% y处前后取0的作用是使包络线在起始点与末尾处不发散。
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function n = findpeaks(x)
% 寻找极值点
n    = find(diff(diff(x) > 0) < 0);
u    = find(x(n+1) > x(n));
n(u) = n(u)+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%