|
楼主 |
发表于 2011-11-16 21:23
|
显示全部楼层
本帖最后由 huazi071783 于 2011-11-16 21:28 编辑
这是我的EMD代码
clear all;clc
acce0_2=load('fi0137.txt');
Y=acce_2(:,3);
% NOTE:
% Nstd: ratio of the standard deviation of the added noise and that
% of Y;(加性噪声的标准差与Y的比)
% NE: Ensemble number for the EEMD
% It should be noted that when Nstd is set to zero and NE is set to 1, the
% program degenerates to a EMD program.(如果Nstd=0,NE=1,该函数退化成普通的EMD函数)
%
Nstd=0;NE=1;
allmode=eemd(Y,Nstd,NE);
figure(1)
subplot(4,1,1)
plot(allmode(:,1))
axis([ 0 2500 -0.2 0.2])
subplot(4,1,2)
plot(allmode(:,3))
subplot(4,1,3)
plot(allmode(:,4))
subplot(4,1,4)
plot(allmode(:,5))
%——————function eemd.m
function allmode=eemd(Y,Nstd,NE)
xsize=length(Y);
dd=1:1:xsize;
Ystd=std(Y); %计算Y的标准差
Y=Y/Ystd;
TNM=fix(log2(xsize))-1; %初始化用于存储原始信号和IMF及趋势项分量的矩阵
TNM2=TNM+2;
for kk=1:1:TNM2,
for ii=1:1:xsize,
allmode(ii,kk)=0.0;
end
end
% 初始化输入信号,可以加加性噪声,当Nstd=0时,不加加性噪声,原始信号
%直接参与EMD分解
for iii=1:1:NE,
for i=1:xsize,
temp=randn(1,1)*Nstd;
X1(i)=Y(i)+temp;
end
%mode矩阵初始化为原始数据
for jj=1:1:xsize,
mode(jj,1) = Y(jj);
end
xorigin = X1;
xend = xorigin;
nmode = 1;
while nmode <= TNM,
xstart = xend;
iter = 1;
while iter<=10,
[spmax, spmin, flag]=extrema(xstart);
upper= spline(spmax(:,1),spmax(:,2),dd);
lower= spline(spmin(:,1),spmin(:,2),dd);
mean_ul = (upper + lower)/2;
xstart = xstart - mean_ul;
iter = iter +1;
end
xend = xend - xstart;
nmode=nmode+1;
for jj=1:1:xsize,
mode(jj,nmode) = xstart(jj);
end
end
for jj=1:1:xsize,
mode(jj,nmode+1)=xend(jj);
end
allmode=allmode+mode;
end
allmode=allmode/NE;
allmode=allmode*Ystd;
%-------——————function extrema.m
function [spmax, spmin, flag]= extrema(in_data)
flag=1;
dsize=length(in_data);
spmax(1,1) = 1;
spmax(1,2) = in_data(1);
jj=2;
kk=2;
while jj<dsize,
if ( in_data(jj-1)<=in_data(jj) & in_data(jj)>=in_data(jj+1) )
spmax(kk,1) = jj;
spmax(kk,2) = in_data (jj);
kk = kk+1;
end
jj=jj+1;
end
spmax(kk,1)=dsize;
spmax(kk,2)=in_data(dsize);
if kk>=4
slope1=(spmax(2,2)-spmax(3,2))/(spmax(2,1)-spmax(3,1));
tmp1=slope1*(spmax(1,1)-spmax(2,1))+spmax(2,2);
if tmp1>spmax(1,2)
spmax(1,2)=tmp1;
end
slope2=(spmax(kk-1,2)-spmax(kk-2,2))/(spmax(kk-1,1)-spmax(kk-2,1));
tmp2=slope2*(spmax(kk,1)-spmax(kk-1,1))+spmax(kk-1,2);
if tmp2>spmax(kk,2)
spmax(kk,2)=tmp2;
end
else
flag=-1;
end
msize=size(in_data);
dsize=max(msize);
xsize=dsize/3;
xsize2=2*xsize;
spmin(1,1) = 1;
spmin(1,2) = in_data(1);
jj=2;
kk=2;
while jj<dsize,
if ( in_data(jj-1)>=in_data(jj) & in_data(jj)<=in_data(jj+1))
spmin(kk,1) = jj;
spmin(kk,2) = in_data (jj);
kk = kk+1;
end
jj=jj+1;
end
spmin(kk,1)=dsize;
spmin(kk,2)=in_data(dsize);
if kk>=4
slope1=(spmin(2,2)-spmin(3,2))/(spmin(2,1)-spmin(3,1));
tmp1=slope1*(spmin(1,1)-spmin(2,1))+spmin(2,2);
if tmp1<spmin(1,2)
spmin(1,2)=tmp1;
end
slope2=(spmin(kk-1,2)-spmin(kk-2,2))/(spmin(kk-1,1)-spmin(kk-2,1));
tmp2=slope2*(spmin(kk,1)-spmin(kk-1,1))+spmin(kk-1,2);
if tmp2<spmin(kk,2)
spmin(kk,2)=tmp2;
end
else
flag=-1;
end
flag=1;
|
|