%%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 = [maxmin, i]%%找到极值点在信号中的坐标(不分极大值和极小值点)
end
elseif sign(d(i))~=sign(d(i+1)) % we are straddling a zero so%%如y=|x|在x=0处是极值点,但该点倒数不存在,所以不能用上面的判
断方法
maxmin = [maxmin, i+1] % 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 = [1 maxes N]
mins = [1 mins N]
%------------------------------------------------------------------------- % spline interpolate to get max and min envelopesform 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 = [imfh]% 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
参考资料
http://zhidao.baidu.com/link?url=Dv2ef87TOGx8zcbCT1UJsZ2kutWrm4FuT5kbMZY5mAAn5yv7APibQ1y8fSag5JvbF2fKlI5jhgpTXu95SDRgi_
下载的工具箱中的emd.m文件里的注释有详细的用法介绍。工具箱的安装
运行install_emd.m文件可以实现此工具箱的安装,uninstall_emd.m实现卸载。
安装中的问题
但是安装的时候,如果使用的是VS的编译器(mbuild –setup、mex –setup设宏族置),会报找不到complex.h的问题(用Linux下的Matlab不会出错),从而使cemdc2_fix.c等文件编译失败,这几个文件是为了快速实现计算EMD而用c编写的,所以嫌罩即使编译失败,也不影响直接使用emd.m实现EMD功能。如果想编译成功,可如下修改:(摘自http://www.chinavib.com/thread-79866-1-1.html)
G. Rilling 07年3月份的程序,运行作者的install_emd.m,出现找不到complex.h的问题,以下是个人的理解和解决过程:(个人的运行环境为matlab6.5)
complex.h的问题
产生原因:采用matlab的C编译函数mex时,定蔽者弊义了C99_OK的宏(EMDS/make_emdc.m (28行)),利用的是ANSI C99标准如果个人的电脑中没有相关的支持,就会出现这个问题。
解决方法:EMDS/make_emdc.m中第28行中mex(’-DC99_OK‘,args(:))语句中的 '-DC99_OK' 即可。
注意:
改完之后,运行install_emd,会出现M_PI没有定义的问题,缺少了常数PI的宏定义,导致一些.c文件编译失败。
产生原因:去掉C99_OK之后,程序中使用的是作者提供的 emd_complex.h和emd_complex.c两个文件来支持复数运算,这两个文件中,并没有定义M_PI这个宏。
解决方法:M_PI这个宏,只在两个文件中(clocal_mean.c和clocal_mean2.c)使用,个人的解决方法是,在相应的头文件(clocal_mean.h和clocal_mean2.h)中加入M_PI的宏定义即可。
在两个.h文件中分别加入一下语句:
#define CLOCAL_MEAN_H
#ifndef M_PI
#define M_PI 3.1415926
#endif
安装完成后,编译输出的.dll文件会出现,重复后缀名的问题,及 xxx.dll 变成了 xxx.dll.dll自己去掉多余的.dll即可
最后,关于版本问题:作者推荐使用7.1+版本,但只是针对个别的函数有影响,主要是作者提供的例子程序,无法在matlab6.5环境中运行,算法的主要功能函数并不受影响。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)