|||
关于IMF原信号与其包络线相割的问题,解决之道当然就是求取与与原信号相切的包络线(以后可称之为“切包络线”),这是能够做到的。为了同时解决瞬时振幅在信号过零点处的阶梯突变,可以将原信号取绝对值之后,再求其切包络线(只求信号上部一条包络线)。
为了解决端点效应,求极点包络线时需要先将信号极值点向两端外面延拓。求切包络线时则需先将整个信号数据向两端外面延拓。这个工作比极点延拓稍微麻烦一点,因为延拓部分既要包含延拓极值点,又要与原信号的端点连接是连续、大致光滑的(至少不要变成新极值点了)。
步骤如下:
1、将极值点延拓函数
[tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
的输出tmin,tmax,zmin,zmax改为:
[ tlmin,tlmax,trmin,trmax,zlmin,zlmax,zrmin,zrmax ] = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym)。
上面输出tlmin,tlmax是左延极小值、极大值点时间序号,zlmin,zlmax是相应函数值;trmin,trmax是右延极小值、极大值点时间序号,zrmin,zrmax是相应函数值。将原函数名改一下(boundary_conditions_lr)保存。
2、调用上述函数,并将待分析的IMF原函数左、右两个端点分别“混入”上面左、右极值点中;
3、将包括上述“混入”点在内的左延极小值、极大值点时间序号,统一按从小到大顺序排序,并保持对应的函数值跟随排序。同样原则处理右延部分;
4、将上述左延点作样条插值,输出的插值函数时间区段,为左延点最小时间序号点更左一个采样点,至原信号左端点时间序号更左一个采样点(即“0”)。同样原则处理右延部分;
5、令上述所得左延插值函数最左边一点的函数值,等于左边第3点函数值。此举是为了保证左边第2点(原函数极值点延拓而来)仍为极值点。同样原则处理右延部分;
信号延拓完成。特别说明一下,本延拓方法适用于任意信号,不局限于IMF。
根据上述步骤所编程序,函数文件名存为signal_extension,调用格式为
[tl,t,tr,xl,x,xr] = signal_extension(x)。其
tl,t,tr分别为左延信号、原信号、右延信号时间序号。
xl,x,xr分别为左延信号、原信号、右延信号函数值。
绘图:
x=imfMBrd(11,:);
[tl,t,tr,xl,x,xr] = signal_extension(x);
lxxx = length(xl)+length(x)+length(xr);
zs=zeros(1,lxxx);
%下面是为了绘图而重复运行
[indmin, indmax, indzer] = extr(x);
t=1:length(x);
z=x;
nbsym=2;
[ tlmin,tlmax,trmin,trmax,zlmin,zlmax,zrmin,zrmax ] = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym);
%上面是为了绘图而重复运行
fg1=figure(1);
set(fg1,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.06,0.95,0.90]);
plot(t,x)
hold on%
plot(tl,xl,'r-')
plot(tr,xr,'r-')
plot([tl,t,tr],zs,'k--')
plot(indmin,x(indmin),'ko')
plot(indmax,x(indmax),'ko')
plot(tlmin,zlmin,'ko')%
plot(tlmax,zlmax,'ko')%
plot(trmin,zrmin,'ko')%
plot(trmax,zrmax,'ko')%
grid
title('imfMBrd之11分量左右延拓及其极值点 ')
运行,得:
图34-1
将x换成一个IMF高频分量再看看:
x=imfMBrd(3,:);
[tl,t,tr,xl,x,xr] = signal_extension(x);
lxxx = length(xl)+length(x)+length(xr);
zs=zeros(1,lxxx);
%下面是为了绘图而重复运行
[indmin, indmax, indzer] = extr(x);
t=1:length(x);
z=x;
nbsym=2;
[ tlmin,tlmax,trmin,trmax,zlmin,zlmax,zrmin,zrmax ] = boundary_conditions_lr(indmin,indmax,t,x,z,nbsym);
%上面是为了绘图而重复运行
fg2=figure(2);
set(fg2,'position',[50,250,1400,525])
li=1;
lj=2;
for i=1:li
for j=1:lj
n=(i-1)*lj+j;
s(n)=subplot(li,lj,n);
set(s(n),'position',[0.03+(j-1)*(0.997^j)/lj,0.06+(li-0.98*(0.999^i)*i)/li,0.92/lj,0.84/li]);
plot(t,x)
hold on%
plot(tl,xl,'r-')%
plot(tr,xr,'r-')%
plot([tl,t,tr],zs,'k--')%
plot(indmin,x(indmin),'ko')%
plot(indmax,x(indmax),'ko')%
plot(tlmin,zlmin,'ko')%
plot(tlmax,zlmax,'ko')%
plot(trmin,zrmin,'ko')%
plot(trmax,zrmax,'ko')%
grid
if j==1
xlim([tl(1)-10 t(100)])
title(['imfMBrd之3分量左延部分及其极值点'])%
else j==2
xlim([t(end-99) tr(end)+10])
title(['imfMBrd之3分量右延部分及其极值点'])%
end
end
end
运行,得:
图34-2
将x换成一个非IMF分量MB2917_rd看看。只将上图程序绘图句柄及时轴显示区限xlim改一下就可以了。
运行,得:
图34-3
信号边界延拓完成。后篇续谈。
(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de01016qx7.html
首发时间:2012-06-19 22:35:46)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 07:05
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社