chao0922 发表于 2008-12-2 22:03

最原始的EMD程序(论坛里的老帖子有点小问题)

我是大地测量专业做数据处理的:
最近在看些EMD方面的文献,首先不谈该方法的利弊,大概扫了下论坛,现在论坛里做EMD貌似少了不少;
因为是新手,所以先学着写最基本的EMD程序,有个老帖子如下:

转网上的一个简单直观的EMD程序,有助初学理解                                                                                                                                                                                                                                                                                                                        呵呵,可能大家都有了这个程序了,不过我觉得这个程序对我现在刚学EMD很有帮助,
因为程序简单,太长了初学的时候没有耐心看那么多行的程序
我根据我的理解对程序作了一些简单的注释,EMD原理基本理解,但是对于程序具体实现时有些地方不太懂,欢迎大家指正,共同探讨。

这个程序是一个最基本的程序,没有包括任何改进算法,如端点延拓等。

% EMD:Emprical mode decomposition
%
% imf = emd(x,n)
%
% x   - input signal (must be a column vector)
% n   - number of intrinsic mode functions
%
% imf - Matrix of intrinsic mode functions (each as a row)
%
% See:Huang et al, Royal Society Proceedings on Math, Physical,
%       and Engineering Sciences, vol. 454, no. 1971, pp. 903-995,
%       8 March 1998
%
% Author: Ivan Magrin-Chagnolleau<ivan@ieee.org>
%
function imf = emd(x,n);%%最好把函数名改为emd1之类的,以免和Grilling的emd冲突                        
                                  %%n为你想得到的IMF的个数
c = x(:)'; % copy of the input signal (as a row vector)
N = length(x);-
% loop to decompose the input signal into n successive IMFs
imf = []; % Matrix which will contain the successive IMF, and the residuefor t=1:n

while(1)%最外层循环

% loop on successive IMFs   
%-------------------------------------------------------------------------   
% inner loop to find each imf   
h = c; % at the beginning of the sifting process, h is the signal
SD = 1; % Standard deviation which will be used to stop the sifting process      
while SD > 0.3 % while the standard deviation is higher than 0.3 (typical value) %%筛选停止准则            
% find local max/min points   
d = diff(h); % approximate derivative %%求各点导数   
maxmin = []; % to store the optima (min and max without distinction so far)   
for i=1:N-2      
if d(i)==0                        % we are on a zero %         
      if sign(d(i-1))~=sign(d(i+1))% it is a maximum %         
      maxmin = ;                         %%
      end      
elseif sign(d(i))~=sign(d(i+1))   % we are straddling a zero so%
                                                                                                   
   maxmin = ;      % define zero as at i+1 (not i) %%         
end      
end         
if size(maxmin,2) < 2 % then it is the residue %%判断信号是不是已经符合残余分量定义      
break      
end         
% divide maxmin into maxes and mins%% 分离极大值点和极小值点      
if maxmin(1)>maxmin(2)            % first one is a max not a min      
maxes = maxmin(1:2:length(maxmin));      
mins= maxmin(2:2:length(maxmin));      
else                              % is the other way around      
maxes = maxmin(2:2:length(maxmin));         
mins= maxmin(1:2:length(maxmin));      
end            % make endpoints both maxes and mins   
maxes = ;   
mins= ;               
%-------------------------------------------------------------------------   % spline interpolate to get max and min envelopes; form imf   
maxenv = spline(maxes,h(maxes),1:N);    %%用样条函数插值拟合所有的极大值点         
minenv = spline(mins, h(mins),1:N);   %%用样条函数插值拟合所有的极小值点   
m = (maxenv + minenv)/2; % mean of max and min enveloppes%%求上下包络的均值   
prevh = h; % copy of the previous value of h before modifying it %%h为分解前的信号   
h = h - m; % substract mean to h                         %% 减去包络均值            
% calculate standard deviation   
eps = 0.0000001; % to avoid zero values   
SD = sum ( ((prevh - h).^2) ./ (prevh.^2 + eps) );         %% 计算停止准则      
end   
imf = ; % store the extracted IMF in the matrix imf   
% if size(maxmin,2)<2, then h is the residue      
% stop criterion of the algo. if we reach the end before n
if size(maxmin,2) < 2   
break
end   
c = c - h; % substract the extracted IMF from the signal
end
return

这里贴的代码和下面一个筒子贴的代码只是缺了:while(1)

其实这段代码里面之所以只分解出一个IMF,和函数形参中n设定没关系,而是因为缺乏while这个外层循环控制,
导致只算出一个IMF。加上while(1)才是正确的(应该说是最基本的EMD)。当然也可以用while(n的表达式)来控制分解层数,当然这没必要。
F的文章还没怎么读,所以先看看这个最基本的,虽然这个很多问题解决不了。

baobao1982 发表于 2008-12-3 09:14

谢谢分享
共同学习

anycall 发表于 2008-12-24 11:49

做的时间序列EMD分解

对非线性数据做的EMD,大家帮忙分析下...

anycall 发表于 2008-12-24 11:51

分解后的数据和原始数据之和有差异,特别是在最后余项部分,数值偏大,有谁能帮忙解决下分解数据偏大的问题。

zuozemin2000 发表于 2011-3-8 22:16

好帖,学习有帮助

lxy2646 发表于 2011-3-15 11:10

很好的学习资料,谢谢。

tru_balla 发表于 2011-3-28 13:00

试试!!!

gefei1989 发表于 2011-4-7 20:33

最近正在看emd的程序,有些地方不是太懂,这个简单的程序给了我不少的启示,多谢楼主了

lulei71 发表于 2011-4-13 11:26

分析的不错……

tomato_zfour 发表于 2011-4-13 19:58

回复 1 # chao0922 的帖子

多谢多谢

whapril 发表于 2011-4-27 11:51

这个程序,我调试时,有错误

??? Error using ==> minus
Matrix dimensions must agree.

Error in ==> emd1fun at 63
      h = h - m; % substract 减去mean to h   %减去包络均值

Error in ==> test3 at 16
imf=emd1fun(data,3);

哪里的问题?

憨憨帝 发表于 2011-4-27 15:26

这个程序为什么我不能运行出来啊?

憨憨帝 发表于 2011-4-27 15:26

出现错误
??? Error: File: D:\Documents and Settings\Administrator\桌面\emd123.m Line: 19 Column: 16
Incomplete or misformed expression or statement.

sunsun1212 发表于 2011-5-22 17:44

好,解释详细,学习了

ben1984 发表于 2011-6-18 11:17

maxmin = []; % to store the optima (min and max without distinction so far)貌似把极值点存入
页: [1] 2
查看完整版本: 最原始的EMD程序(论坛里的老帖子有点小问题)