将脉搏序列基本模式函数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）