中医现代化研究分享 http://blog.sciencenet.cn/u/baishp 当你们都还在想象时,我已经在路上了。这注定是一场一个人的战争吗?

博文

35、(五)切包络线_基于IMF自身的瞬时频率估计

已有 4353 次阅读 2012-11-30 15:32 |个人分类:斤斤计较|系统分类:论文交流| 包络线, 瞬时频率估计, 瞬时振幅, 切点, 极小值

各位博友:端午节快乐!粽子飘香,亲情满堂!

 

 

    言归正传。求IMF“切包络线”方法步骤如下:

 

1、将信号两端延拓;

 

2、求延拓信号极大、极小值点;

 

3、根据上述极大、极小值点求信号上下极点包络线;

 

4、将信号分别减去上下极点包络线,分别得上下“络基信号(我创立的名词,表示以包络线为0基线的信号)”;

 

5、分别求上下“络基信号”的极大值与极小值;

 

6、分别去除上述极大值与极小值点中,小于、大于上下包络线的极大值与极小值,剩下来的极大值与极小值,就是阶段性上下切点;

 

7、分别对上述阶段性上下切点作样条插值,所得插值函数为阶段性上下切包络线;

 

8、将上述阶段性上下切包络线代替上面步骤“4”中的上下极点包络线,重复至步骤“6”;

 

9、将步骤“6”所得阶段性上下切点,与上次所得阶段性上下切点进行比较。如果不完全相同,则重复步骤“7->8->4->5->6”;如果完全相同,则停止循环;

 

10、将步骤“7”最后所得插值函数,截去两端延拓部分,就是最后得到的信号的上下切包络线。

 

11、去除步骤“6”最后所得上下切点中,位于信号两端延拓部分中的上下切点,剩下来的就是原分析信号上下切包络线的上下切点。

 

    上面只是大概的步骤,有些细节问题就不列出了。如果分析信号不是IMF,而是任意信号,那就还有更多的细节问题需要处理。

 

    由于IMF切包络线在接近“0”时,其对称性不总是太好(见图35-1、图35-2),当时我还没有想到信号取绝对值的方法,就总是想到用切包络线代替EMD中的极点包络线进行迭代运算,以期得到更加对称的切包络线(真是一念之差,绕城十匝)。而且感觉用切包络线代替极点包络线似乎更加合情合理:包络线嘛,就应该尽量将信号包络起来才对,否则包络一部分,留一部分在外面,那就不是完全的包络线了。用切包络线代替极点包络线,最终所得的IMF可能更加合符实际情况,产生虚假IMF的可能性也许会更小,IMF的正交指数也可能会更小。

 

    所以我编了求任意信号上下切包络线的程序tangent_envelope(实际效果见图35-3),并且将它代替emd程序中的极点包络线进行了迭代运算。虽然现在好几列信号的运行都已经通过,但运行时间需增加几十倍,而且这些信号序列数据都较短(10000点以下)。所以我现在还无法评价它的性能较既有EMD程序的优劣性。

 

    绘图看看。 所编tangent_envelope函数格式:

 

[tenvb,tenva,tpoib,tpoia]=tangent_envelope(x)

 

输出tenvb,tenva,tpoib,tpoia分别是分析信号x的下、上切包络线;下、上切点。tangent_envelope内部未对x作延拓处理。

 

    下面是IMF上下切包络线某局部接近0值时的情况:

 

x=imfMBrd(4,:);

[tenvb,tenva,tpoib,tpoia]=tangent_envelope(x);

 

[indmin, indmax] = extr(x);%

 

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

 

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

set(gca,'position',[0.035,0.06,0.95,0.88]);

 

plot(x,'r-')
hold on%
plot(tenvb,'b--')%下切包络线

plot(tenva,'b--')%上切包络线

plot(zs,'k-.')%

plot(tpoib,x(tpoib),'k.')%下切点

plot(tpoia,x(tpoia),'k.')%上切点

 

plot(indmin,x(indmin),'ko')%极小值点
plot(indmax,x(indmax),'ko')%极大值点


xlim([15500,15700])
title('imfMBrd的4分量上下切包络线及其切点、极值点')

 

     运行,得:


    图35-1


    图35-2
 

    可见IMF切包络线在接近“0”时,其对称性不总是太好。图中是选了比较极端的情况,一般情况下其对称性当然没有这么差。

 

    将上面程序中的x=imfMBrd(4,:)换成任意信号x=MB2917_rd,修改一下绘图句柄、x轴显示范围,运行,得:



    图35-3

 

    从上图来看,在EMD中用“切包络线”代替“极点包络线”进行求取IMF的迭代运算,似乎更加“合情合理”。

 

    上面两图只是看看数据中部局部情况,故分析信号未作延拓。下面看看求整个IMF“切包络线综合振幅”的情况:

 

x=imfMBrd(13,:);

 

[tl,t,tr,xl,x,xr] = signal_extension(x);

 

lxr=[xl,x,xr];
ltr=[tl,t,tr];

 

abslxr=abs(lxr);
absxl=abs(xl);
absx=abs(x);
absxr=abs(xr);


[indmin, indmax] = extr(abslxr);

 

[tenvb,tenva,tpoib,tpoia]=tangent_envelope(abslxr);


tenval=tenva(1:length(tl));
tenvam=tenva(length(tl)+1:length(tl)+length(t));
tenvar=tenva(length(tl)+length(t)+1:length(tl)+length(t)+length(tr));

 

llxr = length(lxr);
zs=zeros(1,llxr);

 

fg3=figure(3);
set(fg3,'position',[50,250,1050,700])
set(gca,'position',[0.035,0.07,0.95,0.88]);


plot(t,absx)
hold on%
plot(tl,absxl,'r-')%
plot(tr,absxr,'r-')%

 

plot(ltr,tenva,'k--')


plot(indmax+tl(1),abslxr(indmax),'ko')%%
plot(tpoia+tl(1),abslxr(tpoia),'.k')

 

plot(ltr,zs,'k-.')


axis([tl(1)-fix(length(ltr)/50),tr(end)+fix(length(ltr)/50),-max(tenva)*0.04,max(tenva)*1.03])


title('imfMBrd之13分量延拓并取绝对值后上切包络线及其切点、极值点')

 

    运行,得:


    图35-4

 

    “切包络线综合振幅”情况良好。 续谈。

 

 

(本文首发于:http://blog.sina.com.cn/s/blog_6ad0d3de01016two.html

首发时间:2012-06-22 22:49:20)

 



https://blog.sciencenet.cn/blog-825323-637855.html

上一篇:34、(四)信号延拓_基于IMF自身的瞬时频率估计
下一篇:36、(六)四种包络_基于IMF自身的瞬时频率估计
收藏 IP: 14.153.189.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-6-17 15:26

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部