|||
将脉搏序列基本模式函数imfMBrd第10分量的瞬时振幅、瞬时频率、复信号虚部再作几张放大、对比图:
(1)IMF原信号及其复信号虚部:
x=imfMBrd(10,:);
[indmin, indmax, indzer] = extr(x);%求极值点、
[env,envmoy] = cenvelope(x,1);%求包络线
envb=env(1,:);%将负包络线提取出来
enva=env(2,:);%将正包络线提取出来
lx=length(x);
zs=zeros(1,lx);
[ a f y] = afyimf( x );
fg1=figure(1);
set(fg1,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.05,0.95,0.90]);
plot(x,'r-')
hold on
plot(envb,'k--')
plot(enva,'k--')
plot(envmoy,'k:')
plot(y,'b-')
plot(zs,'k-')
grid
title('imfMBrd 10分量及其复信号之虚部 ')
xlim([5000,8200])
运行,得:
图33-1
fg2=figure(2);
set(fg2,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.05,0.95,0.90]);
plot(x,'r-')
hold on
plot(envb,'k--')
plot(enva,'k--')
plot(f*150,'b-')%f*150是为了放大好对比观察
plot(zs,'k-')
grid
title('imfMBrd 10分量瞬时频率突变位置分析')
xlim([5000,8200])
运行,得:
图33-2
fg3=figure(3);
set(fg3,'position',[50,250,1050,700])
set(gca,'position',[0.05,0.05,0.93,0.90]);
tt=6500:7500;
stem3(tt,f(tt),a(tt));%hold on
view(-15,30)
grid on
title('imfMBrd 10分量局部三维时频图')
运行,得:
图33-3
从上面几张图,可以更加清楚地看到,虽然没有了以前所说的瞬时频率跳变,但所求得的瞬时振幅、瞬时频率、复信号虚部,都还存在着突变,没有完全达到我希望的结果。其中瞬时振幅的求取方式,是造成问题的根源。因此我打算针对其过零点与割点问题,采取两个改进措施。下面先看第一个措施。
此前的瞬时振幅,是由IMF上下包络线取绝对值后,在信号过零点处分段取值、连接而成。除了造成瞬时振幅自身在信号过零点处的突变,也造成瞬时频率与复信号虚部在过零点处突变。我已经想到了一个很好的方法解决这个问题,就是先将正负极值点统一取绝对值后,再作样条插值(只得到一个插值函数)。具体如下:
lx=length(x);
t=1:lx;
z=x;
nbsym=2;
INTERP='spline';
[tmin,tmax,zmin,zmax] = boundary_conditions_emd;(indmin,indmax,t,x,z,nbsym);
zmin=abs(zmin);
TZ=[tmin,tmax;zmin,zmax];
TZ=TZ';
TZTZ=sortrows(TZ,1);
TZTZ=TZTZ';
tmm=TZTZ(1,:);
zmm=TZTZ(2,:);
amp = interp1(tmm,zmm,t,INTERP);%振幅
此振幅综合由上下极值点插值而成,因此可以叫它“综合振幅”。它也是原imf信号的上包络线,取负值后又是原信号的下包络线,可以统称为“综合包络线”。
将上述amp的计算方法与计算结果代替函数[ a f y] = afyimf( x )中a的计算方法与计算结果,再重复本文前面的计算过程:
x=imfMBrd(10,:);
[indmin, indmax, indzer] = extr(x);%求极值点、
[env,envmoy] = cenvelope(x,1);%求包络线
envb=env(1,:);%将负包络线提取出来
enva=env(2,:);%将正包络线提取出来
lx=length(x);
zs=zeros(1,lx);
[ a f y] = afyimf( x );
fg4=figure(4);
set(fg4,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.05,0.95,0.90]);
plot(x,'r-')
hold on
plot(envb,'k--')
plot(-a,'k:')%综合包络线
plot(enva,'k--')
plot(a,'k:')%综合包络线
plot(envmoy,'k--')
plot(y,'b-')
plot(zs,'k-')
grid
title('综合包络线 imfMBrd 10分量及其复信号之虚部')
xlim([5000,8200])
运行,得:
图33-4
fg5=figure(5);
set(fg5,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.05,0.95,0.90]);
plot(x,'r-')
hold on
plot(envb,'k--')
plot(-a,'k:')%综合包络线
plot(enva,'k--')
plot(a,'k:')%综合包络线
plot(f*150,'b-')
plot(zs,'k-')
grid
title('综合包络线 imfMBrd 10分量瞬时频率突变位置分析')
xlim([5000,8200])
运行,得:
图33-5
fg6=figure(6);
set(fg6,'position',[50,250,1050,700])
set(gca,'position',[0.05,0.05,0.93,0.90]);
tt=6500:7500;
stem3(tt,f(tt),a(tt));%hold on
view(-15,30)
grid on
title('综合包络线 imfMBrd 10分量局部三维时频图')
运行,得:
图33-6
综合振幅自身在过零点处、对瞬时频率、复信号虚部在过零点处的突变的消除作用,图中已经表达的很清楚了,不再赘述。
关于信号与包络线相割的问题,后面继续研究。
(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de01016odk.html
首发时间:2012-06-17 12:08:03)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 09:26
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社