声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 6920|回复: 23

[HHT] 镜像延拓的emd程序

  [复制链接]
发表于 2009-3-8 15:19 | 显示全部楼层 |阅读模式

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

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

x
function imf = TryEmd(ip)
%--------------------------------------------------------------------------
%    Designer:万学功
N=length(ip);
x=1:N;
[max,max_No,min,min_No]=MaxMin(ip);
%    Released Time:2007-6-21
%    Email:bigbomb825826@yahoo.com
%    Department:江苏大学机械工程学院测控技术研究所
%--------------------------------------------------------------------------
%----函数说明(152):
%    输入ip--待分解数据序列,一维列向量;
%    输出imf--分解后得到的IMF分量,二维数组表示;
%--------------------------------------------------------------------------
%----调用函数:
%    Enovelope--次端点镜像法求解上、下包络线;(39)
%    MaxMin--求局部极值点及其位置;(24)
%    MaxMax--求数组极大值;(6)
%    ZeroNum--求过零点数目;(17)
%--------------------------------------------------------------------------
%----参数设置:
     PRE_yuzhi=0.005; %--偏差阈值;
     NUM1=5; %--分解得到IMF分量的最多数目;
     NUM2=800; %--提取IMF分量的最多筛分次数;
%--------------------------------------------------------------------------
flag=0;
I=0;
num1=0;
num2=0;
while(max_No(1)~=0 && min_No(1)~=0)
    if(num1>=NUM1) break; end
    while(1)
        num2=num2+1;
        [yy1,yy2]=enovelope(x,max,max_No,min,min_No);
        mean=(yy1+yy2)/2;
        h1=ip-mean;
        [max,max_No,min,min_No]=MaxMin(h1);
        if(max_No(1)==0 || min_No(1)==0)
            flag=1;
            break;
        end
        [yy1,yy2]=enovelope(x,max,max_No,min,min_No);
        mean=(yy1+yy2)/2;
        hmax=MaxMax(mean);
        N1=ZeroNum(h1);
        tiaojian1=hmax-PRE_yuzhi;
        tiaojian2=abs(length(max_No)+length(min_No)-N1);
        if((tiaojian1<0 && tiaojian2<=1) || num2<=NUM2)
            I=I+1;
            imf(I,:)=h1;
            break;
        else
            ip1=h1;
            [max,max_No,min,min_No]=MaxMin(ip1);
        end
    end
   
    if(flag==1)
        flag=0;
        break;
    else
        ip=ip-imf(I,:);
        num1=I;
        [max,max_No,min,min_No]=MaxMin(ip);
    end
end
imf(I+1,:)=ip;

我看过之后有几处不大明白
1  在程序中是怎么体现中止条件的也就是SD(0.2和0.3之间)。
2  在程序中用的控制循环的条件是max_No(1)~=0 && min_No(1)~=0 请问这个式子的含义是什么? 为什么要用它做循环条件
max_No(1)有没有更深刻的含义?

      我十分期待与您交流!!
                                                                   谢谢!!

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2009-3-13 11:08 | 显示全部楼层

回复 楼主 qqvirile 的帖子

已经解决了!!有同样问题的朋友可以交流!
发表于 2009-3-20 16:37 | 显示全部楼层

回复 沙发 qqvirile 的帖子

万学功的那个程序好像是次端点极值法的吧?楼主有没有次端点极值法的资料啊?
 楼主| 发表于 2009-3-20 16:57 | 显示全部楼层

回复 板凳 00yangmin 的帖子

你是要Enovelope函数么?
function [y_up,y_low]= Enovelope(x,max,max_No,min,min_No)
%--次端点镜像延拓法获得的上下包络线;


l=length(x);
l1=length(max);
N1=max_No(1)-1;
N2=length(x)-max_No(end);

l2=length(min);
N3=min_No(1)-1;
N4=length(x)-min_No(end);

max1_No(1)=1;
for i=2:(l1+1) max1_No(i)=max_No(i-1)+N1; end
max1_No(l1+2)=l+N1+N2;

max1(1)=max(1);
for i=2:(l1+1) max1(i)=max(i-1); end
max1(l1+2)=max(end);

min1_No(1)=1;
for i=2:(l2+1) min1_No(i)=min_No(i-1)+N3; end
min1_No(l2+2)=l+N3+N4;

min1(1)=min(1);
for i=2:(l2+1) min1(i)=min(i-1); end
min1(l2+2)=min(end);

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

x2=1:(l+N3+N4);
cs2 = spline(min1_No,min1);%三次样条差值
y2=ppval(cs2,x2);%求出插值曲线
for i=1:l y_low(i)=y2(i+N3); end
 楼主| 发表于 2009-3-20 17:07 | 显示全部楼层
这个效果一般 我用了这个程序处理了论坛里一个帖子里的数据 效果很不好 不过那个帖子最后也没有讨论出怎么解决问   你可以试一下 我用小波做了也不能提出频率!!至于小波包我到没有用过,你试一下吧

数据说明:
OR014_0是0马力负载下的故障数据。转速频率为:30Hz;故障频率为:107.3;
OR014_2是2马力负载下的故障数据,转速频率为:29.15HZ,故障频率为:104.5Hz。

OR014_0.mat

75.86 KB, 下载次数: 47

OR014_2.mat

75.37 KB, 下载次数: 25

发表于 2009-4-14 17:20 | 显示全部楼层

解决EMD端点效应问题

本人到目前为止研究解决EMD端点问题的方法有次端点镜像延拓法、包络极值延拓法和端点极值法(即直接在端点处取极值),但是感觉这三个方法分解出来的IMF效果都不怎么好,有人建议考虑曲线拟合问题,我用的是全区三次样条拟合问题,也尝试过分段三次样条拟合和其他的曲线拟合方法,得到的IMF与黄大吉等人发表的《希尔伯特_黄变换的端点延拓》中的图形相差较远,请各位大侠指条明路啊!答辩在即,而我到目前为止还没解决一维EMD分解问题,二维的HHT更是难上加难!
回复 支持 1 反对 0

使用道具 举报

发表于 2009-4-14 21:16 | 显示全部楼层
楼主,能不能解释一下你已经解决了的那两个问题啊?
我也对这两个问题感到不解。
发表于 2009-4-16 15:47 | 显示全部楼层

困惑

各位大侠,我是新手,请问从论坛里下的那个emd.m那个800多行的程序,看起来好复杂,请问怎么来调用对含噪信号进行EMD分解呢,我运行总是出错,说emd 那个function不能定义什么的,请高手给点指点,小妹感激不尽啊
 楼主| 发表于 2009-4-20 15:03 | 显示全部楼层

回复 7楼 news 的帖子

这个问题我很长时间之前看的 具体的是有点忘了
但是我记得我看一个文章的时候 在实际的编程过程中只要能分解到单调的时候 就可以认为是imf分解结束
具体当时怎么想明白的 我现在也有点忘了 再交流吧
发表于 2010-7-15 15:55 | 显示全部楼层

h1=ip-mean;
怎么理解,好像不对
发表于 2010-9-27 12:47 | 显示全部楼层
新手 看看
发表于 2010-11-10 13:40 | 显示全部楼层
发表于 2011-10-22 18:53 | 显示全部楼层
想请问一下镜像延拓的函数怎么调用,也就是我要用的话,怎么编写程序呢,谢谢各位大虾
发表于 2011-12-1 09:36 | 显示全部楼层
感谢分享,学习了
发表于 2012-2-20 17:27 | 显示全部楼层
是不是没有maxmin 函数啊
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-5 22:01 , Processed in 0.091307 second(s), 27 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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