|||
前文说了我对学界既有的瞬时频率估计方法的失望,虽然我并没有将所有的那些方法都试遍,但我已经没有信心与耐心仔细去研究它们了,以后有时间、有心情时再去慢慢详细了解吧。现在我打算抛开一切“高深的”理论,从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)
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 08:22
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社