声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

声振论坛 展示 振动理论 信号处理 查看内容

送给搞EMD或者HHT但没有下载到完整程序的朋友

2015-11-5 03:44| 发布者: aspen| 查看: 1567| 评论: 95|原作者: eight

摘要: 下载完整程序的方法:(考虑到程序更新问题,建议使用第二、三种方法,flandrin就分别在2006年2月、2007年3月进行了更新) 1. 本论坛(自己搜索一下)——已经过时 2. (适用于2007年3月之前)登录以下网站:http: ...
下载完整程序的方法:(考虑到程序更新问题,建议使用第二、三种方法,flandrin就分别在2006年2月、2007年3月进行了更新)
1. 本论坛(自己搜索一下——已经过时
2. (适用于2007年3月之前)登录以下网站:http://perso.ens-lyon.fr/patrick.flandrin/emd.html,这是flandrin的一个关于EMD和HHT的网页,进入后就可以看见主要的几个m文件,对只搞EMD的朋友这些程序就够用了,如果要进行HHT后续分析,则需要点击Time-Frequency ToolBox链接下载工具箱,具体位置在Time-Frequency ToolBox页面的download处,下载后得到tftb-0.1.tar.gz压缩包,压缩包中有个mfiles文件夹,所有程序就齐全了。

注:instfreq 函数就在tftb-0.1.tar.gz里面;至于extr函数,如果阁下是在2006年2月之后从flandrin网站下载的程序,那么根本不需要这个m文件,因为该函数已集成到此版本的emd.m中,如果阁下是在2006年2月之前下载的,那么的确需要这个m文件,如果当时没有下载到,那么请下载此版本的emd.m,然后单独把extr函数提取出来;至于程序看不懂的问题,请按照emd.m文件中的帮助信息下载相关的两篇文章,就是Huang和Flandrin写的关于EMD和HHT的经典著作

Huang98年经典著作的下载地址:http://www.fuentek.com/technologies/hht.htm

3. (适用于2007年3月之后)登录上述网站下载压缩包,里面有详细说明。具体介绍可以参阅以下链接的帖子:EMD估计有新进展

注:程序看不懂的问题,请按照emd.m文件中的帮助信息下载相关的两篇文章(帮助信息列出的5篇文章中的前两篇),就是Huang和Flandrin写的关于EMD和HHT的经典著作;另外,instfreq函数的实现请参阅 refguide.pdf 和 tutorial.pdf 两个文件中的该函数解释,并根据给出的参考文献查找相关资料阅读,这两个pdf文件在 hhspectrum.m 帮助信息关于时频工具箱(Time-Frequency Toolbox)的下载地址中可以找到

[ 本帖最后由 eight 于 2007-11-17 22:34 编辑 ]
发表评论

最新评论

引用 MVH 2006-10-21 11:23
哈哈,eight有点火了
引用 hnlzx 2006-10-23 15:10
eigth说很全,instfreq 函数和extr函数就在那里。我以前都试通过!

但hhspectrum.m还需要调用hilbert.m,我是使用matlab7.0工具箱里自带的,不知有没有问题?
引用 eight 2006-10-23 15:31
原帖由 hnlzx 于 2006-10-23 15:10 发表
eigth说很全,instfreq 函数和extr函数就在那里。我以前都试通过!

但hhspectrum.m还需要调用hilbert.m,我是使用matlab7.0工具箱里自带的,不知有没有问题?



用的就是matlab自带的hilbert函数,没有问题。不过我对EMD比较熟悉一点
引用 hdwang 2006-10-24 21:57
感谢搂住!但是我调用emd.m这个程序进行分解出现了如下错误,这是怎么回事啊?
??? Undefined function or variable 'io'.

Error in ==> G:\study\EMD分解程序\emd.m
On line 185  ==> ort = io(x,imf);

Error in ==> G:\study\EMD分解程序\Untitled1.m
On line 3  ==> imf = emd(x);
我就使用它本身给出的例子程序:
%examples:
%
%>>x = rand(1,512);
%
%>>imf = emd(x);
%
%>>imf = emd(x,struct('stop',[0.1,0.5,0.05],'maxiterations',100));
%Remark: the following syntax is equivalent
%>>imf = emd(x,'stop',[0.1,0.5,0.05],'maxiterations',100);
%
%>>options.dislpay = 1;
%>>options.fix = 10;
%>>options.maxmodes = 3;
%>>[imf,ort,nbits] = emd(x,options);
引用 eight 2006-10-24 22:17
原帖由 hdwang 于 2006-10-24 21:57 发表
感谢搂住!但是我调用emd.m这个程序进行分解出现了如下错误,这是怎么回事啊?
??? Undefined function or variable 'io'.

Error in ==> G:\study\EMD分解程序\emd.m
On line 185  ==> ort = io(x,imf ...



flandrin网页就有io.m
引用 hunt_girl 2006-10-25 09:32
t=0:0.0004:1;
x=sin(2*pi*50*t)+0.3*sin(5.5*pi*50*t);
% x=sin(2*pi*50*t)+0.3*sin(6.5*pi*50*t);
tst=1;
stop = [0.01,0.1,0.01];
引用 hunt_girl 2006-10-25 09:36
我发现运行程序后得到的IMF1分量的右端点发散,怎么回事啊???
clc;
clear;
t=0:0.0004:1;
x=sin(2*pi*50*t)+0.3*sin(5.5*pi*50*t);
% x=sin(2*pi*50*t)+0.3*sin(6.5*pi*50*t);
tst=1;
stop = [0.01,0.1,0.01];
NBSYM = 2;
lx = length(x);
sd = stop(1);
sd2 = stop(2);
tol = stop(3);
MAXITERATIONS=2000;

sdt(lx) = 0;
sdt = sdt+sd;
sd2t(lx) = 0;
sd2t = sd2t+sd2;

% maximum number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;
r = x;
imf = [];
k = 1;
% iterations counter for extraction of 1 mode
nbit=0;

% total iterations counter
NbIt=0;

while ner > 2
        
  % current mode
  m = r;
  
  % mode at previous iteration
  mp = m;
  
  sx = sd+1;
  
  % tests if enough extrema to proceed
  test = 0;
  [indmin,indmax,indzer] = extr(m);
  lm=length(indmin);
  lM=length(indmax);
  nem=lm + lM;
  nzm=length(indzer);
  
  j=1;
  
  % sifting loop
  while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS
   
    if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)
      disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])
      disp(['stop parameter mean value : ',num2str(s)])
    end
   
    % boundary conditions for interpolations :
        lm=length(indmax);
        ln=length(indmin);
     if indmax(1) < indmin(1)
      if m(1) > m(indmin(1))
        lmax = fliplr(indmax(2:min(lm,NBSYM+1)));
        lmin = fliplr(indmin(1:min(ln,NBSYM)));
        lsym = indmax(1);
      else
        lmax = fliplr(indmax(1:min(lm,NBSYM)));
        lmin = [fliplr(indmin(1:min(ln,NBSYM-1))),1];
        lsym = 1;
      end
    else

      if m(1) < m(indmax(1))
        lmax = fliplr(indmax(1:min(lm,NBSYM)));
        lmin = fliplr(indmin(2:min(ln,NBSYM+1)));
        lsym = indmin(1);
      else
        lmax = [fliplr(indmax(1:min(lm,NBSYM-1))),1];
        lmin = fliplr(indmin(1:min(ln,NBSYM)));
        lsym = 1;
      end
    end
   
    if indmax(end) < indmin(end)
      if m(end) < m(indmax(end))
        rmax = fliplr(indmax(max(lm-NBSYM+1,1):end));
        rmin = fliplr(indmin(max(ln-NBSYM,1):end-1));
        rsym = indmin(end);
      else
        rmax = [lx,fliplr(indmax(max(lm-NBSYM+2,1):end))];
        rmin = fliplr(indmin(max(ln-NBSYM+1,1):end));
        rsym = lx;
      end
    else
      if m(end) > m(indmin(end))
        rmax = fliplr(indmax(max(lm-NBSYM,1):end-1));
        rmin = fliplr(indmin(max(ln-NBSYM+1,1):end));
        rsym = indmax(end);
      else
        rmax = fliplr(indmax(max(lm-NBSYM+1,1):end));
        rmin = [lx,fliplr(indmin(max(ln-NBSYM+2,1):end))];
        rsym = lx;
      end
    end
   
    tlmin = 2*t(lsym)-t(lmin);
    tlmax = 2*t(lsym)-t(lmax);
    trmin = 2*t(rsym)-t(rmin);
    trmax = 2*t(rsym)-t(rmax);
   
    % in case symmetrized parts do not extend enough
    if tlmin(1) > t(1) | tlmax(1) > t(1)
      if lsym == indmax(1)
        lmax = fliplr(indmax(1:min(lm,NBSYM)));
      else
        lmin = fliplr(indmin(1:min(ln,NBSYM)));
      end
      if lsym == 1
        error('bug')
      end
      lsym = 1;
      tlmin = 2*t(lsym)-t(lmin);
      tlmax = 2*t(lsym)-t(lmax);
    end   
   
    if trmin(end) < t(lx) | trmax(end) < t(lx)
      if rsym == indmax(end)
        rmax = fliplr(indmax(max(lm-NBSYM+1,1):end));
      else
        rmin = fliplr(indmin(max(ln-NBSYM+1,1):end));
      end
      if rsym == lx
        error('bug')
      end
      rsym = lx;
      trmin = 2*t(rsym)-t(rmin);
      trmax = 2*t(rsym)-t(rmax);
    end
         
    mlmax =m(lmax);
    mlmin =m(lmin);
    mrmax =m(rmax);
    mrmin =m(rmin);
     
    % definition of envelopes from interpolation
        
    envmax = interp1([tlmax t(indmax) trmax],[mlmax m(indmax) mrmax],t,'spline');
    envmin = interp1([tlmin t(indmin) trmin],[mlmin m(indmin) mrmin],t,'spline');

    envmoy = (envmax + envmin)/2;

    m = m - envmoy;
   
    [indmin,indmax,indzer] = extr(m);
    lm=length(indmin);
    lM=length(indmax);
    nem = lm + lM;
    nzm = length(indzer);

    % evaluation of mean zero
    sx=2*(abs(envmoy))./(abs(envmax-envmin));
    s = mean(sx);
     mp = m;
    nbit=nbit+1;
    NbIt=NbIt+1;

    if(nbit==(MAXITERATIONS-1))
      warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
    end
end
  imf(k,:) = m;
  if tst
    disp(['mode ',int2str(k),' enregistre'])
  end
  nbits(k) = nbit;
  k = k+1;
  r = r - m;
  [indmin,indmax,indzer] = extr(r);
  ner = length(indmin) + length(indmax);
  nzr = length(indzer);
  nbit=1;

  if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
    if ner > 2
      warning('forced stop of EMD : too small amplitude')
    else

      disp('forced stop of EMD : too small amplitude')
    end
    break
  end
  
end

imf(k,:) = r;

if tst
  close
end
引用 无水1324 2006-10-25 09:41
"实在受不了","检查人品"

帮助解决一个问题
需要提升到这个高度吗?解决了,心里还是感觉有点凉!
引用 hunt_girl 2006-10-25 09:43
附件是我运行得到的图像。imf1,大家看看我的程序哪儿出错了??
引用 eight 2006-10-25 11:02
原帖由 无水1324 于 2006-10-25 09:41 发表
"实在受不了","检查人品"

帮助解决一个问题
需要提升到这个高度吗?解决了,心里还是感觉有点凉!



呵呵,我可没有这个意思。只是不明白大家在找资料的时候为何缺少耐性,更何况那个页面就摆在面前
引用 eight 2006-10-25 11:03
原帖由 hunt_girl 于 2006-10-25 09:43 发表
附件是我运行得到的图像。imf1,大家看看我的程序哪儿出错了??



应该是正常的,这是EMD算法的端点效应导致的
引用 hunt_girl 2006-10-25 11:20
楼上你能不能把我所给的信号,用你的程序运行一下把结果发上来给我看看,我现在心理没什么底儿.因为我上次用极值延拓的右边界很好,但是左边界误差大,我也做过AR模型预测端点极值的,效果对有些信号好,有些不行.(等会我把那篇文章上传上来.这次看见镜像延拓的.本来文章上说的效果应该不错,不过运行后感觉不尽如人意.
引用 hunt_girl 2006-10-25 11:26
参考《应用自回归模型处理EMD方法中的边界问题》
引用 eight 2006-10-25 12:29
原帖由 hunt_girl 于 2006-10-25 11:20 发表
楼上你能不能把我所给的信号,用你的程序运行一下把结果发上来给我看看,我现在心理没什么底儿.因为我上次用极值延拓的右边界很好,但是左边界误差大,我也做过AR模型预测端点极值的,效果对有些信号好,有些不 ...



我运行了一下,跟你贴的结果一样。这是正常的,不同信号效果不同,EMD算法本身就是一个近似的算法,端点延拓是一个没有完全解决的问题
引用 computer 2006-10-26 16:02
To 老八:
        不好意思,我找了很久还是没有找到你说的如果看不懂,就在emd.m下面下载两篇著作,请具体说一下,行吗?我回去一定会好好检查我的人品的,呵呵.不甚感激!
引用 computer 2006-10-26 16:43
我从google下面下载了那些程序,可运行emd.m的文件时还是出现错误了啊,请帮我指点一下,错误如下:
??? Index exceeds matrix dimensions.

Error in ==> emd>init at 575
x = varargin{1};

Error in ==> emd at 75
[x,t,sd,sd2,tol,display_sifting,sdt,sd2t,ner,nzr,lx,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:});       

>>
引用 computer 2006-10-26 16:46
To:老八,刚才问的那两篇著作已经找到,不用帮我解答了,谢谢1能不能帮我看看哪个程序怎么出问题了吗?不甚感激!
引用 eight 2006-10-26 18:43
原帖由 computer 于 2006-10-26 16:46 发表
To:老八,刚才问的那两篇著作已经找到,不用帮我解答了,谢谢1能不能帮我看看哪个程序怎么出问题了吗?不甚感激!



x是一个向量,我忘记是列向量还是行向量了,这些小问题自己调试跟踪一下就ok了
引用 hdwang 2006-10-28 19:58
八楼的程序直接运行后出现了下面结果,没有看到图?
mode 1 enregistre
mode 2 enregistre
mode 3 enregistre
mode 4 enregistre
mode 5 enregistre
mode 6 enregistre

查看全部评论(95)

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

GMT+8, 2024-5-16 09:48 , Processed in 0.062738 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

返回顶部