www.gusucode.com > 局部均值分解源代码 难得的matlab程序代码源码 > lmd/example_lmd.m

    %测试zhaochun.m的性能
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
h=y;
a=1;
wucha1=0.001;
t=1:length(h);
%retry
h_pos=position(h);%极值点位置序列
h_pos=position(h);%极值点位置序列
tpoint=t(h_pos);%极值点时间值
hpoint=h(h_pos);%极值点幅度值
subplot(3,1,1);plot(t,h);
hold on;plot(tpoint,hpoint,'r+');hold off;
hpoint=bianjie(h,hpoint,1);%端点处理后的极值点，多出了2个极值点
[mi,ai]=find_miai(hpoint);%找出极值点之间的均值函数和包络函数
mi_point=junzhi(mi);%所有极值点处的均值序列（幅值）-纵坐标（点序列）
ai_point=junzhi(ai);
lmi_point=mi(1);%左端点的均值
rmi_point=mi(length(mi));%右端点的均值
lai_point=ai(1);%左端点的包络
rai_point=ai(length(ai));%右端点的包络
%
mi_point_d=link(lmi_point,rmi_point,mi_point);%连接端点均值及所有极值点处的 均值 （带端点的均值序列）（点序列）
ai_point_d=link(lai_point,rai_point,ai_point);
tpoint_d=link(t(1),t(length(t)),tpoint);
mi_fun=chadian1(tpoint_d,mi,mi_point_d);%包含端点和极值点和普通点的 均值序列-平缓前的均值序列
ai_fun=chadian1(tpoint_d,ai,ai_point_d);%
subplot(3,1,2);plot(t,mi_fun);
hold on;plot(t,ai_fun,'r');
kmax=max(diff(tpoint_d));%找出时间跨度最大的  相邻几点 间的 距离
kmax1=uint16(kmax/3);
kmax1=double(kmax1);
jiou=uint8(rem(kmax1,2));
if kmax1<3
    move_k=3; % b<3,滑动跨度=3，
elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1
            move_k=(kmax1-1);
   else move_k=kmax1; %b>=3,b是奇数时，跨度=b
end
[mi_move,move_mi_num]=move(mi_fun,move_k);
[ai_move,move_ai_num]=move(ai_fun,move_k);
plot(t,mi_move,'g');
plot(t,ai_move,'c');
hold off;
subplot(3,1,3);plot(t,h);
hold on;plot(t,mi_move,'g');
plot(t,ai_move,'c');hold off;
 mi=mi_move;
ai=ai_move;
a=a.*ai;
si=(h-mi)./ai;
h=si;
ai_funmax=max(ai);
ai_funmin=min(ai);
if (ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1)
    tuichu=1
else
    tuichu=0
end

    
%------------------------------------------------
%测试emd.m求序列包络的效果
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
h=y;
t=1:length(h);
[envmin,envmax,envmoy,indmin,indmax,indzer] = envelope(t,h,'spline');
plot(t,h);hold on;plot(t,envmin,'g');plot(t,envmax,'c');plot(t,envmoy,'r');hold off;
%-----------------------------
%完整测试zhaochun.m函数的性能
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
h=y;
a=1;
wucha1=0.001;
[pf,a,si]=zhaochun(a,h,wucha1);
%----------------------------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
c=y;
k=1;
a=1;
wucha1=0.001;
%retry
k=k+1;
h=c;
[pf,a,si]=zhaochun(a,h,wucha1);
c=c-pf;
PF(k,:)=pf;
A(k,:)=a;
SI(k,:)=si;
c_pos=poss(c);
if length(c_pos)<=1
    tuichu=1
else
    tuichu=0
end
%-----------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
h=y;
a=1;
wucha1=0.001;
chun_num=0;
t=1:length(h);
%retry
chun_num=chun_num+1
subplot(2,1,1);plot(t,h);
[envmin,envmax,envmoy,indmin,indmax,indzer] = envelope(t,h,'spline');
minmax=cat(2,indmin,indmax);
pos=sort(minmax);
tpoint=t(pos);
hpoint=h(pos);
hold on;plot(tpoint,hpoint,'r+');hold off;
mi=(envmax+envmin)./2;
ai=abs(envmax-envmin)./2;
subplot(2,1,2);plot(t,h);
hold on;plot(t,mi,'g');
plot(t,ai,'c');hold off;
a=a.*ai;
si=(h-mi)./ai;
h=si;
ai_funmax=max(ai);
ai_funmin=min(ai);
if (ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1)
    tuichu=1
else
    tuichu=0
end
%----------------
%完整改进
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
h=y;
a=1;
wucha1=0.001;
[pf,a,si]=zhaochun1(a,h,wucha1);
%-------------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
c=y;
k=0;

wucha1=0.001;
%retry
k=k+1
h=c;
a=1;
[pf,a,si]=zhaochun1(a,h,wucha1);
c=c-pf;
PF(k,:)=pf;
A(k,:)=a;
SI(k,:)=si;
c_pos=poss(c);
if length(c_pos)<=1
    tuichu=1
else
    tuichu=0
end
plot(t,A(k,:));
%-------------

%------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);

n_l=nengliang(y);
c=y;
k=0;
wucha1=0.001;
%retry
k=k+1
h=c;
a=1;
[pf,a,si]=zhaochun1(a,h,wucha1);
c=c-pf;
    PF(k,:)=pf;
    A(k,:)=a;
    SI(k,:)=si;
c_pos=poss(c);
nn=nengliang(c);
if length(c_pos)<=3||nn<n_l/10
    PF(k+1,:)=c;
    tuichu=1
else
    tuichu=0
end
%--------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
t=0:0.001:1-0.001;
y=x(t);
[pf,a,si]=lmd(y);
%------------
x=@(t) (1+0.5*cos(9*pi*t)).*cos(200*pi*t+2*cos(10*pi*t))+sin(pi*t).*sin(30*pi*t);
fs=1000;%采样率
t=0:1/fs:1-0.001;
y=x(t);
[pf,a,si]=lmd2(y);
subplot(3,1,1);plot(t,pf(1,:));
subplot(3,1,2);plot(t,pf(2,:));
subplot(3,1,3);plot(t,pf(3,:));
ff=shunpin(si).*fs;%实际频率要乘以采样率
figure;
subplot(211);plot(t,ff(1,:));
subplot(212);plot(t,ff(2,:));
%subplot(3,1,3);plot(t,ff(3,:));
%-----------------------------
x=@(t) (2+cos(90*t).*cos(500*t+1800.*t.*t));
fs=5000;
t=0:1/fs:0.341;
y=x(t);
subplot(5,1,1);plot(t,y);xlabel('原始信号');
[pf,a,si]=lmd1(y);
subplot(5,1,2);plot(t,pf(1,:));xlabel('PF1');
subplot(5,1,3);plot(t,pf(2,:));xlabel('PF2');
subplot(5,1,4);plot(t,pf(3,:));xlabel('PF3');
subplot(5,1,5);plot(t,pf(4,:));xlabel('残量信号');