信号处理方法 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

zhlong 发表于 2007-5-16 11:54

转网上的一个简单直观的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
% 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 %%导数为0的点，即”驻点“，但驻点不一定都是极值点，如y=x^3的x=0处
   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%%如y=|x|在x=0处是极值点，但该点倒数不存在，所以不能用上面的判
                                                                                          断方法
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

[ 本帖最后由 zhlong 于 2007-5-16 12:01 编辑 ]

zhlong 发表于 2007-5-16 12:08

t=;
a=16;
x=sin(2*pi*(200*t))+sin(2*pi*(500*t));

分别利用上面的程序（令n=2)和Flandrin的程序得到的结果如图11和22所示

ytlzd 发表于 2007-5-16 20:04

转网上的一个简单直观的EMD程序，有助初学理解

多谢zhong:handshake ，不过还得请教一下，Flandrin的程序是什么，
能说详细点么，小弟刚开始学，多谢啦！！！！！！！！！！！！！

zhangnan3509 发表于 2007-5-16 20:06

回复 #3 ytlzd 的帖子

你为什么对这么多置顶贴视若无睹，我们的努力白费了！:@Q

qinle 发表于 2007-5-28 07:39

补一个和楼主出处在一起的EMD程序另一种写法:这两个程序的源应该是美国rice大学.
注解我就不加了,和楼主的有很多是一样的.大家看看吧
% EMD:Emprical mode decomposition
%
% imf = emd(x)
%
% x - input signal (must be a column vector)
%
% This version will calculate all the imf's (longer)
%
% imf - Matrix of intrinsic mode functions (each as a row)
%    with residual in last 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);

c = x(:)'; % copy of the input signal (as a row vector)
N = length(x);

%-------------------------------------------------------------------------
% loop to decompose the input signal into successive IMF

imf = []; % Matrix which will contain the successive IMF, and the residue

while (1) % the stop criterion is tested at the end of the loop

%-------------------------------------------------------------------------
% 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
         maxmin = ;
      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 - 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 size(maxmin,2) < 2
   break
end

c = c - h; % substract the extracted IMF from the signal

end

return

shmily0904 发表于 2007-7-30 15:02

比较了zhlong给的程序，qinle的多一个while(1)语句，我运行后者，发现matlab死循环，不知道是否这个原因，这个while语句到底如何解释呢？:@Q

kevin19821 发表于 2007-9-2 17:28

能否解释以下if size(maxmin,2) < 2中括弧中的2是什么意思?

kevin19821 发表于 2007-9-2 19:16

看了程序,以下问题有点不解:
(1) c = x(:)'这样等于c=x,并不是c=x',不知道编程的者的意图是什么?
(2)imf=emd(x,n)中imf的个数n在程序中那里有体现?
(3)size(maxmin,2)还有这个2是什么意思呀?

[ 本帖最后由 kevin19821 于 2007-9-2 20:05 编辑 ]

jingrenzhi 发表于 2007-9-2 20:00

不错,对初学者很有帮助!我有个问题?
对于下面这条语句中的"n为你想得到的IMF的个数"这句话中n是在运行程序前设定好的吗?根据设想随便定的吗?不是emd自己分解出几个imf就是几个吗?
function imf = emd(x,n);%%最好把函数名改为emd1之类的，以免和Grilling的emd冲突
%%n为你想得到的IMF的个数

kevin19821 发表于 2007-9-2 20:07

同意楼上的,

eight 发表于 2007-9-2 21:08

本帖最后由 wdhd 于 2016-9-6 14:37 编辑

原帖由 jingrenzhi 于 2007-9-2 20:00 发表
不错,对初学者很有帮助!我有个问题?
对于下面这条语句中的"n为你想得到的IMF的个数"这句话中n是在运行程序前设定好的吗?根据设想随便定的吗?不是emd自己分解出几个imf就是几个吗?
function imf = emd(x,n);%% ...

可以让它自己分解，最终得到的个数不定；也可以提供用户自己设定，分解到满足用户指定的个数后就终止，剩下的全部作为余量

jingrenzhi 发表于 2007-9-2 22:06

回复 #11 eight 的帖子

哦,是这样!谢谢!那么"如果分解到满足用户指定的个数后就终止，剩下的全部作为余量"的话,楼主的这句程序是不是应该写成
function imf = emd(x);?
如5楼的程序!

eight 发表于 2007-9-2 23:13

本帖最后由 wdhd 于 2016-9-6 14:37 编辑

原帖由 jingrenzhi 于 2007-9-2 22:06 发表
哦,是这样!谢谢!那么"如果分解到满足用户指定的个数后就终止，剩下的全部作为余量"的话,楼主的这句程序是不是应该写成
function imf = emd(x);?
如5楼的程序!
function imf = emd(x) 这是默认的分解方式，分解的个数不定
function imf = emd(x,n) 这是提供用户自定义的方式

kevin19821 发表于 2007-9-3 16:35

我自己用楼主程序加曲线延拓后所得的IMF和RES图如下:
晕!图片太大,发不上来

[ 本帖最后由 kevin19821 于 2007-9-3 16:44 编辑 ]

jingrenzhi 发表于 2007-9-3 16:48

回复 #13 eight 的帖子

哦!明白!谢谢!

页: [1] 2

声振论坛's Archiver

转网上的一个简单直观的EMD程序，有助初学理解

转网上的一个简单直观的EMD程序，有助初学理解

回复 #3 ytlzd 的帖子

回复 #11 eight 的帖子

回复 #13 eight 的帖子