emd分解matlab代码2016a

由于EMD分解算法比较复杂，需要多个函数配合使用，下面是EMD分解的MATLAB代码示例：

function [IMF,residue]=emd(x)

%emd函数的输入为待分解的信号x，输出为IMF和残差residue

%其中IMF为分解出的各个固有模态函数，residue为最终的残差信号

%下面是emd函数的具体实现：

N=length(x);

imf_num=0;%初始化IMF的数量

residue=x;%将初始信号赋值给残差变量

while true

imf_num=imf_num+1;%IMF数量加1

%对残差信号进行迭代分解

c=residue;%将残差信号赋值给中间变量c

while true

%计算信号的上下包络线

u=extrema(c,'max');

d=extrema(c,'min');

%计算平均值

m=(u+d)/2;

%计算均值包络线与原信号的差

h=c-m;

%判断是否为IMF

if isimf(h)

IMF(imf_num,:)=h;%将分解出的IMF存入IMF矩阵中

break;

end

%如果不是IMF，则将h作为下一次迭代的输入信号

c=h;

end

%将分解出的IMF从残差信号中减去

residue=residue-IMF(imf_num,:);

%判断是否满足停止条件

if isstop(residue)

break;

end

end

%最后剩下的残差信号即为residue

end

%以下是一些辅助函数

function [yext]=extrema(y,mode)

%计算信号的上下包络线

len=length(y);

yext=zeros(1,len);

for i=2:len-1

if strcmp(mode,'max')%计算上包络线

if y(i)>y(i-1)&&y(i)>y(i+1)

yext(i)=y(i);

end

elseif strcmp(mode,'min')%计算下包络线

if y(i)<y(i-1)&&y(i)<y(i+1)

yext(i)=y(i);

end

end

end

%处理端点

if strcmp(mode,'max')

yext(1)=y(1);

yext(len)=y(len);

elseif strcmp(mode,'min')

yext(1)=y(1);

yext(len)=y(len);

end

%插值处理

idx=find(yext~=0);

yext=yext(idx);

x=idx;

yext=interp1(x,yext,1:len,'pchip');

end

function [flag]=isimf(x)

%判断信号是否为IMF

%需要满足以下两个条件：

%1.信号的局部极值点的个数必须等于或者比信号的零交叉点的个数少1

%2.信号的极值点的幅值的平均值为0

n=0;

while true

n=n+1;

%计算信号的零交叉点

zc=find(x(1:end-1).*x(2:end)<0);

%计算信号的局部极值点

ext=find((diff(x(1:end-1))>0)&(diff(x(2:end))<0));

%判断是否满足条件1

if abs(length(ext)-length(zc))<=1

break;

end

%将信号减去平均值

x=x-mean(x);

%判断是否满足条件2

if n>=10

break;

end

end

%计算极值点幅值的平均值

ext_val=x(ext);

ave_ext_val=mean(abs(ext_val));

%判断是否满足条件2

if ave_ext_val>0.2*std(abs(x))

flag=false;

else

flag=true;

end

end

function [flag]=isstop(x)

%判断是否满足停止条件

%需要满足以下两个条件：

%1.分解出的IMF数量不能超过原信号的一半

%2.最终的残差信号不能为0

if length(x)<=1

flag=true;

return;

end

if length(x)==2

flag=all(x==0);

return;

end

if length(x)==3

flag=all(diff(x)==0);

return;

end

if length(x)>3

%计算信号的零交叉点

zc=find(x(1:end-1).*x(2:end)<0);

%计算信号的局部极值点

ext=find((diff(x(1:end-1))>0)&(diff(x(2:end))<0));

%计算分解出的IMF数量

imf_num=length(ext)-length(zc)+1;

%判断是否满足条件1

if imf_num>=ceil(length(x)/2)

flag=true;

else

flag=false;

end

%计算最终的残差信号

residue=x-sum(IMF(1:imf_num,:));

%判断是否满足条件2

if sum(abs(residue))<=0.2*sum(abs(x))

flag=true;

else

flag=false;

end

end

end

在使用时，只需将待分解的信号作为输入参数传入emd函数中即可，例如：

x=1:100;

IMF=residue=emd(x);

其中，IMF为分解出的各个IMF，residue为最终的残差信号。需要注意的是，EMD分解算法比较耗时，对于较长的信号需要分段分解才能得到正确的结果。