次端点镜像延拓问题
这里有个次端点镜像延拓的程序，不知道怎么用function imf = TryEmd(ip)%--------------------------------------------------------------------------
%
N=length(ip);
x=1:N;
=MaxMin(ip);
%--------------------------------------------------------------------------
%----函数说明(152)：
% 输入ip--待分解数据序列，一维列向量；
% 输出imf--分解后得到的IMF分量，二维数组表示；
%--------------------------------------------------------------------------
%----调用函数：
% Enovelope--次端点镜像法求解上、下包络线；(39)
% MaxMin--求局部极值点及其位置；(24)
% MaxMax--求数组极大值;(6)
% ZeroNum--求过零点数目；(17)
%--------------------------------------------------------------------------
%----参数设置：
PRE_yuzhi=0.005; %--偏差阈值；
NUM1=5; %--分解得到IMF分量的最多数目；
NUM2=800; %--提取IMF分量的最多筛分次数；
%--------------------------------------------------------------------------
flag=0;
I=0;
num1=0;
num2=0;
while(max_No(1)~=0 && min_No(1)~=0)
if(num1>=NUM1) break; end
while(1)
num2=num2+1;
=enovelope(x,max,max_No,min,min_No);
mean=(yy1+yy2)/2;
h1=ip-mean;
=MaxMin(h1);
if(max_No(1)==0 || min_No(1)==0)
flag=1;
break;
end
=enovelope(x,max,max_No,min,min_No);
mean=(yy1+yy2)/2;
hmax=MaxMax(mean);
N1=ZeroNum(h1);
tiaojian1=hmax-PRE_yuzhi;
tiaojian2=abs(length(max_No)+length(min_No)-N1);
if((tiaojian1<0 && tiaojian2<=1) || num2<=NUM2)
I=I+1;
imf(I,:)=h1;
break;
else
ip1=h1;
=MaxMin(ip1);
end
end
if(flag==1)
flag=0;
break;
else
ip=ip-imf(I,:);
num1=I;
=MaxMin(ip);
end
end
imf(I+1,:)=ip;
end
%%%%%%%%%%%%%%
function = Enovelope(x,max,max_No,min,min_No)
%--次端点镜像延拓法获得的上下包络线；
l=length(x);
l1=length(max);
N1=max_No(1)-1;
N2=length(x)-max_No(end);
l2=length(min);
N3=min_No(1)-1;
N4=length(x)-min_No(end);
max1_No(1)=1;
for i=2:(l1+1)
max1_No(i)=max_No(i-1)+N1;
end
max1_No(l1+2)=l+N1+N2;
max1(1)=max(1);
for i=2:(l1+1)
max1(i)=max(i-1);
end
max1(l1+2)=max(end);
min1_No(1)=1;
for i=2:(l2+1)
min1_No(i)=min_No(i-1)+N3;
end
min1_No(l2+2)=l+N3+N4;
min1(1)=min(1);
for i=2:(l2+1)
min1(i)=min(i-1);
end
min1(l2+2)=min(end);
x1=1:(l+N1+N2);
cs1=spline(max1_No,max1);%三次样条差值
y1=ppval(cs1,x1);%求出插值曲线
for i=1:l
y_up(i)=y1(i+N1);
end
x2=1:(l+N3+N4);
cs2=spline(min1_No,min1);%三次样条差值
y2=ppval(cs2,x2);%求出插值曲线
for
i=1:l
y_low(i)=y2(i+N3);
end
end
页:
[1]