|||
上一篇谈了RCADA的eemd函数存在的问题。这一篇谈谈我的eemd方法。我还是打算用G.Rilling编写的emd程序作eemd的基程序。但分析信号叠加各次白噪声再经emd分解后,所得imf分量个数不一样,怎么样把它们合成在一起呢?我的方法是分组:分量个数相同的imf放在同一组,总体计算只在组内进行。至于在相同的初始条件下为什么会分解出不同个数的imf分量来,这本身就是一个问题,应该另行研究,不应该回避。因为“真理只有一个”。
先交代一个小问题。我在《29、黄锷院士的一个疏忽》中提到,用Rilling的emd函数计算所得到的jmf,偶尔会出现大于0的极小值或小于0的极大值情况,我称其为非严格的imf,否则,若所有的极大值都大于0,所有的极小值都小于0,则称其为严格的imf。由于我以后使用的我自编的瞬时频率估计函数,是针对严格的imf,因此我对Rilling的emd函数作了一点点修改,使其计算所得imf都是严格的imf。修改方法是将其反映“过零点与极值点数目相差不超过1”的相关代码改为“所有的极大值都大于0,所有的极小值都小于0”代码。为了以示区别,改后函数以emdb表示。
将《29、黄锷院士的一个疏忽》中相关程序再运行一遍:
imf_MB2917_rd = emdb(MB2917_rd);
x1=imf_MB2917_rd(1,:);
[env,envmoy] = cenvelope(x1,1);%求上下包络线及其均值
env=env';
z=zeros(1,length(x1));
plot(x1,'r-')
hold on
plot(env)
plot(envmoy,'k--')
plot(z,'m-')
xlim([33900 34100])
得:
图51-1.
与图29-1比较,可知已解决了“非严格imf”问题。
取一段数据MB1387,为本人2003、6、17~2007、4、3之间的脉搏信号。先用emdb计算其imf, 根据所得imf分量个数,设定加噪信号分解得到的的imf分组条件。
程序如下:
x=MB1387;
nstd=0.1;
n=100;
[imf,o,v]=emdb(x);
STR0.imf0=imf;
STR0.ort0=o;
dv=length(v);%分解所得imf个数
lx=length(x);
sd=std(x);
i=1;
STR=[];
STR3=[];
STR2=[];
STR1=[];
STR_1=[];
STR_2=[];
STR_3=[];
m3=1;
m2=1;
m1=1;
m_1=1;
m_2=1;
m_3=1;
while i<=n
r=randn(1,lx);
ran=r*sd*nstd;
xr=x+ran;
[imf o v]=emdb(xr);
%分解所得imf个数等于原分析信号分解所得imf个数时
if length(v)==dv
STR(i).imf=imf;
STR(i).ort=o;
STR(i).ran=ran;
i=i+1;
%分解所得imf个数比原分析信号分解所得imf个数多出3个时
elseif length(v)==dv+3
STR3(m3).imf=imf;
STR3(m3).ort=o;
STR3(m3).ran=ran;
m3=m3+1;
%分解所得imf个数比原分析信号分解所得imf个数多出2个时
elseif length(v)==dv+2
STR2(m2).imf=imf;
STR2(m2).ort=o;
STR2(m2).ran=ran;
m2=m2+1;
%分解所得imf个数比原分析信号分解所得imf个数多出1个时
elseif length(v)==dv+1
STR1(m1).imf=imf;
STR1(m1).ort=o;
STR1(m1).ran=ran;
m1=m1+1;
%分解所得imf个数比原分析信号分解所得imf个数少出1个时
elseif length(v)==dv-1
STR_1(m_1).imf=imf;
STR_1(m_1).ort=o;
STR_1(m_1).ran=ran;
m_1=m_1+1;
%分解所得imf个数比原分析信号分解所得imf个数少出2个时
elseif length(v)==dv-2
STR_2(m_2).imf=imf;
STR_2(m_2).ort=o;
STR_2(m_2).ran=ran;
m_2=m_2+1;
%分解所得imf个数比原分析信号分解所得imf个数少出3个时
elseif length(v)==dv-3
STR_3(m_3).imf=imf;
STR_3(m_3).ort=o;
STR_3(m_3).ran=ran;
m_3=m_3+1;
end
end
end
看看运行结果:
图51-2
加噪信号与原分析信号分解的imf分量个数不一样的,只有STR1与STR_1有数据,其余都是空集。经过大量运行,该数量差异,极少有2个以上的。说明程序中设置最多差3个,足够适应实际情况。
将前面程序改成函数形式,调用格式为
[STR0,STR,STR3,STR2,STR1,STR_1,STR_2,STR_3]=imfstr(x,nstd,n);
根据变量STR,计算一下总体imf:
allmode=zeros(size(STR(1).imf));
for j=1:n
allmode=allmode+(STR(j).imf);
end
allmode=allmode/n;
%画图:
zs=zeros(1,lx);
W=1080;
H=720;
w=0.80;
h=0.70;
m=5;n=3;k=13;
fg=figure(3);
set(fg,'position',[50,50,W,H]);
for i=1:m
for j=1:n
p=j+(i-1)*n;
if p<=k
s(p)=subplot(m,n,p);
set(s(p),'position',[0.04+(j-1)/n,(1-0.95*i/m),w/n,h/m]);
x=allimf(p,:);
plot(x)
hold on
plot(zs,'k--')
mx=max(x);
mn=min(x);
xn=mx-mn;
axis([fix(-0.1*lx) fix(1.1*lx) mn-0.1*xn mx+0.1*xn])
end
end
end
title('allmode_MB1387_01_100')