EMD的matlab程序

clear;
clc;
tic;
N=512;
fs=512;
t=0:1/fs:(N-1)/fs;
xinhao=sin(2*pi*50*t)+0.3*sin(5.5*pi*50*t);
plot(xinhao);
title('');
xlabel('');
ylabel('');
imf=sigmatching_emd(xinhao,[0.02,0.5,0.05],0);
plot_emd_result(imf);
toc
function imf=sigmatching_emd(x,stop_condition,display)
if nargin<3
display=0;
end
if nargin<2;
stop_condition=[0.05,0.5,0.05];
display=0;
end
r=x;
k=1;
while ~stop_EMD(r);
m=r;
[stop_sift,moyenne]=stop_sifting(m,stop_condition,display);
if (max(abs(m)))<(1e-16)*(max(abs(x)))
break
end
while ~stop_sift
m=m-moyenne;
[stop_sift,moyenne]=stop_sifting(m,stop_condition,display);
end
imf(k,:)=m;
k=k+1;
r=r-m;
end
imf(k,:)=r;
end
%--------------------------------------------------------------------------function stop=stop_EMD(r)
[indmin indmax]=extr(r);
ner=length(indmin)+length(indmax);
stop=ner<3;
end
%--------------------------------------------------------------------------function [envmoy,nem,nzm,amp]=compute_mean(m,t,display) NBSYM=2;
[indmin,indmax,indzer]=extr(m);
nzm=length(indzer);
nem=length(indmin)+length(indmax);
if display==1;
figure
subplot(3,1,1)
plot(m);
subplot(3,1,2)
box on
hold on
plot(m);
plot(indmin,m(indmin),'r.');
plot(indmax,m(indmax),'c.');
hold off
end
[lz,rz,sig]=ppyt(m);
[indmin,indmax]=extr(sig);
[tmin,tmax,mmin,mmax]=boundary_conditions_emd(indmin,indmax,1:length(sig),sig,sig,NBSYM); envmin=interp1(tmin,mmin,[(lz+1):(length(t)+lz)],'spline');
envmax=interp1(tmax,mmax,[(lz+1):(length(t)+lz)],'spline');
envmoy=(envmin+envmax)/2;
amp=mean(abs(envmax-envmin),1)/2;
if display==1;
subplot(3,1,3)
box on
hold on
plot(m);
plot(envmin,'r-');
plot(envmax,'c-');
plot(envmoy,'y-');
hold off
end
end
%--------------------------------------------------------------------------
function [stop,envmoy]=stop_sifting(m,stop_condition,display)
sd=stop_condition(1);
sd2=stop_condition(2);
tol=stop_condition(3);
t=1:length(m);
try
[envmoy,nem,nzm,amp]=compute_mean(m,t,display);
sx=abs(envmoy)./amp;
stop=~((mean(sx>sd)>tol|any(sx>sd2))&(all(nem>2)));
stop=stop&&~(abs(nzm-nem)>1);
catch
stop=1;
envmoy=zeros(1,length(m));
end
end
function [zjsy,begin]=sig_match(sig)
[indmax indmin]=extr(sig);
l=length(sig);
lmax=length(indmax);
lmin=length(indmin);
if lmax>=lmin
k=lmin;
else
k=lmax;
end
for i=1:(k-1)
tx(i)=(indmax(1)*indmin(i+1)-indmin(1)*indmax(i+1))/(indmax(1)-indmin(1));
end
n=1:l;
y=interp1(n,sig,tx,'linear');
for i=1:length(y)
e(i)=abs(sig(indmax(i+1))-sig(indmax(1)))+abs(sig(indmin(i+1))-sig(indmin(1)))+abs(y(i)-sig(1)); end
[ea,indea]=sort(e);
inde=indea(1);
zjsy=tx(inde)-1;
zjsy=floor(zjsy);
if inde<=1
begin=1;
else
begin=tx(inde-1);
end
begin=floor(begin);
end
function [lz,rz,sigu]=ppyt(sig)
lz=sig_match(sig);
l=length(sig);
sig1 = fliplr(sig);
rz=sig_match(sig1);
sig2=[sig1(1:rz),sig1];
sig3 = fliplr(sig2);
sigu=[sig(1:lz),sig3];
end
function plot_emd_result(imf)
[m n]=size(imf);
figure
t=1:n;
for k=1:m-1
subplot(m,1,k);
biaoqian(k)=gca;
plot(t,imf(k,:))
xlabel('t')
ylabel(['imf',int2str(k)])
set(gca,'XLim',[0,n])
end
title(biaoqian(1),'EMD·Ö½â½á¹û')
subplot(m,1,m); plot(t,imf(m,:),'r') set(gca,'XLim',[0,n]) ylabel('r')
xlabel('t')。