声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1270|回复: 0

[HHT] 次端点镜像延拓问题

[复制链接]
发表于 2011-12-4 22:15 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
这里有个次端点镜像延拓的程序,不知道怎么用
  1. function imf = TryEmd(ip)
  2. %--------------------------------------------------------------------------
  3. %  
  4. N=length(ip);
  5. x=1:N;
  6. [max,max_No,min,min_No]=MaxMin(ip);

  7. %--------------------------------------------------------------------------
  8. %----函数说明(152):
  9. %    输入ip--待分解数据序列,一维列向量;
  10. %    输出imf--分解后得到的IMF分量,二维数组表示;
  11. %--------------------------------------------------------------------------
  12. %----调用函数:
  13. %    Enovelope--次端点镜像法求解上、下包络线;(39)
  14. %    MaxMin--求局部极值点及其位置;(24)
  15. %    MaxMax--求数组极大值;(6)
  16. %    ZeroNum--求过零点数目;(17)
  17. %--------------------------------------------------------------------------
  18. %----参数设置:
  19.      PRE_yuzhi=0.005; %--偏差阈值;
  20.      NUM1=5; %--分解得到IMF分量的最多数目;
  21.      NUM2=800; %--提取IMF分量的最多筛分次数;
  22. %--------------------------------------------------------------------------
  23. flag=0;
  24. I=0;
  25. num1=0;
  26. num2=0;
  27. while(max_No(1)~=0 && min_No(1)~=0)
  28.     if(num1>=NUM1) break; end
  29.     while(1)
  30.         num2=num2+1;
  31.         [yy1,yy2]=enovelope(x,max,max_No,min,min_No);
  32.         mean=(yy1+yy2)/2;
  33.         h1=ip-mean;
  34.         [max,max_No,min,min_No]=MaxMin(h1);
  35.         if(max_No(1)==0 || min_No(1)==0)
  36.             flag=1;
  37.             break;
  38.         end
  39.         [yy1,yy2]=enovelope(x,max,max_No,min,min_No);
  40.         mean=(yy1+yy2)/2;
  41.         hmax=MaxMax(mean);
  42.         N1=ZeroNum(h1);
  43.         tiaojian1=hmax-PRE_yuzhi;
  44.         tiaojian2=abs(length(max_No)+length(min_No)-N1);
  45.         if((tiaojian1<0 && tiaojian2<=1) || num2<=NUM2)
  46.             I=I+1;
  47.             imf(I,:)=h1;
  48.             break;
  49.         else
  50.             ip1=h1;
  51.             [max,max_No,min,min_No]=MaxMin(ip1);
  52.         end
  53.     end
  54.    
  55.     if(flag==1)
  56.         flag=0;
  57.         break;
  58.     else
  59.         ip=ip-imf(I,:);
  60.         num1=I;
  61.         [max,max_No,min,min_No]=MaxMin(ip);
  62.     end
  63. end
  64. imf(I+1,:)=ip;
  65. end
  66. %%%%%%%%%%%%%%

  67. function [y_up,y_low]= Enovelope(x,max,max_No,min,min_No)
  68. %--次端点镜像延拓法获得的上下包络线;

  69. l=length(x);
  70. l1=length(max);
  71. N1=max_No(1)-1;
  72. N2=length(x)-max_No(end);

  73. l2=length(min);
  74. N3=min_No(1)-1;
  75. N4=length(x)-min_No(end);
  76. max1_No(1)=1;
  77. for i=2:(l1+1)
  78.     max1_No(i)=max_No(i-1)+N1;
  79. end
  80. max1_No(l1+2)=l+N1+N2;
  81. max1(1)=max(1);
  82. for i=2:(l1+1)
  83.     max1(i)=max(i-1);
  84. end
  85. max1(l1+2)=max(end);
  86. min1_No(1)=1;
  87. for i=2:(l2+1)
  88.     min1_No(i)=min_No(i-1)+N3;
  89. end
  90. min1_No(l2+2)=l+N3+N4;
  91. min1(1)=min(1);

  92. for i=2:(l2+1)
  93.     min1(i)=min(i-1);
  94. end
  95. min1(l2+2)=min(end);

  96. x1=1:(l+N1+N2);
  97. cs1=spline(max1_No,max1);%三次样条差值
  98. y1=ppval(cs1,x1);%求出插值曲线
  99. for i=1:l
  100.     y_up(i)=y1(i+N1);
  101. end

  102. x2=1:(l+N3+N4);
  103. cs2=spline(min1_No,min1);%三次样条差值
  104. y2=ppval(cs2,x2);%求出插值曲线
  105. for
  106.     i=1:l
  107.     y_low(i)=y2(i+N3);
  108. end
  109. end
复制代码
回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 18:27 , Processed in 0.056426 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表