前文说了我对学界既有的瞬时频率估计方法的失望，虽然我并没有将所有的那些方法都试遍，但我已经没有信心与耐心仔细去研究它们了，以后有时间、有心情时再去慢慢详细了解吧。现在我打算抛开一切“高深的”理论，从EMD当中所得到的基本物理量，根据瞬时频率最原始的含义，即各时间点上信号源转动的角度，来估计瞬时频率。

    众所周知，EMD之所以叫EMD（经验模式分解），是因为它最初并不是基于一种理论推导，而是基于这种分解过程所具有的物理意义被世人所公认，而每一组分解结果都能够经受理论的检验(各IMF之间的正交性)。因此上面我打算采用的瞬时频率估计法，是EMD的一以贯之、很自然的延伸。姑且可称之为“基本模式估计法”

    这种方法的关键，在于要明白、承认，IMF的包络线，其实就是信号的瞬时振幅(振幅函数)；而IMF是瞬时振幅在实轴上的投影；某时刻IMF值与IMF包络线值之比，就是该时刻瞬时相位的余弦值。有了瞬时相位，通过相位差分即得到瞬时频率；通过振幅、相位，还可以求出振幅在虚轴上的投影，从而得到IMF对应的复数函数。

    IMF的包络线有两条，哪一条是IMF的瞬时振幅呢？从理论上来说，IMF的两条包络线的均值为0，两条包络线是完全对称的，因此两条包络线的绝对值是相等的，都等于IMF的瞬时振幅。但在实际处理时，由于两条包络线不可能做到绝对对称，所以还是应该分别处理。在复平面上，当信号源运行在1、4象限时，IMF由上包络线在实轴上的投影形成；当信号源运行在2、3象限时，IMF由下包络线在实轴上的投影形成。因此，瞬时振幅应该根据IMF的情况由上下包络线分段组成。下面就根据这个原则来求IMF的瞬时振幅、瞬时频率等。

    下面imfMBrd含义见前两篇博文。

x=imfMBrd(11,:);%先任取imfMBrd的一个分量

lx=length(x);
zs=zeros(1,lx);

[indmin, indmax, indzer] = extr(x);%求极值点、
[env,envmoy] = cenvelope(x,1);%求包络线及其均值

envb=env(1,:);%将负包络线提取出来
enva=env(2,:);%将正包络线提取出来

fg1=figure(1);
set(fg1,'position',[50,5,700,1050]);

w=0.04;
W=0.94;

h=0.04;
d=0.003;
h1=h+0.75-3*d;
h2=h+0.50-2*d;
h3=h+0.25-d;
h4=h;

H=0.18;

sb1=subplot(4,1,1);
set(sb1,'position',[w,h1,W,H]);
plot(x,'r-')
hold on
plot(envb,'g--')
plot(enva,'b--')
plot(envmoy,'m--')
plot(zs,'k-')
grid
title('分析信号及其包络线')
xlim([0,lx+10])

c41=(x>0);
c23=(x<0);

dx=diff(x);
dx1=[0,dx];
dx2=[dx,0];
dx=(dx1+dx2)/2;

c12=(dx<0);
c34=(dx>0);

c1=c41.*c12;
c2=c23.*c12;
c3=c23.*c34;
c4=c41.*c34;

sb2=subplot(4,1,2);
set(sb2,'position',[w,h2,W,H]);
stem(c1,'.b-')
hold on
stem(c2,'.r-')
stem(c3,'.g-')
stem(c4,'.m-')
grid
title('四象限掩码')
xlim([0,lx+10])

x1=x.*c1;
x2=x.*c2;
x3=x.*c3;
x4=x.*c4;

sb3=subplot(4,1,3);
set(sb3,'position',[w,h3,W,H]);
stem(x1,'.b-')
hold on
stem(x2,'.r-')
stem(x3,'.g-')
stem(x4,'.m-')
grid
title('经掩码作用后的分析信号')
xlim([0,lx+10])

aenvb=abs(envb);
aenva=abs(enva);

tampb=c23.*aenvb;
tampa=c41.*aenva;
tamp=tampa+tampb;%瞬时振幅

sb4=subplot(4,1,4);
set(sb4,'position',[w,h4,W,H]);
stem(tampb,'.g-')
hold on
stem(tampa,'.b-')
plot(tamp,'k-')
grid
title('经掩码作用后的瞬时振幅')
xlim([0,lx+10])

    运行，得： 

    图31-1

    下面求瞬时相位、瞬时频率与复信号虚部。

    运行，得：

    图31-2

    从图31-2的第3个子图可以看到，所估计的“瞬时频率”没有前几篇博文中所提到的“跳变”。图中也有一些类似“线谱”一样，或者说象两根“门柱”一样高高耸立的频率成分，但它们是处在同一数量级下的变化，与以前所讲的由接近“0Hz”的负频变成接近“0.5Hz”的正频，性质截然不同的。为了与以前的“频率跳变”相区别，我就将此处的频率变化称之为“频率突变”吧。

    此“频率突变”是由于分析信号与它的包络线相割造成的。包络线以外部分，由于信号、包络线之比值绝对值大于1，反余弦无意义，所以被强制等于1与-1，反余弦值(瞬时相位)为0与0.5*2*pi，两割点之间部分，瞬时频率等于0；两割点之外靠近割点处，瞬时频率即形成“突变”。

    虽然“频率突变”是频率在同一数量级下的变化，但它毕竟还是人为造成的，不是信号本身的属性，因此也是不完全符合我预期的目标的。我希望的瞬时频率估计值，更加连续光滑。

    这一篇就暂到这里，相关问题后面继续分析。

（本文首发于：http://blog.sina.com.cn/s/blog_6ad0d3de01016f0a.html

首发时间：2012-06-10 19:35:44）