跪求一个emd 去噪的程序 matlab 代码 带中文解释的 方便理解

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 = [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_

供参考：

lev=5

[c,l]=wavedec(x,lev,wname)

sigma=wnoisest(c,l,1)

alpha=2

thr1=wbmpen(c,l,sigma,alpha)

[thr2,nkeep]=wdcbm(c,l,alpha)

xd1=wdencmp('gbl',c,l,wname,lev,thr1,'s',1)

[xd2,cxd,lxd,perf0,perfl2]=wdencmp('lvd',c,l,wname,lev,thr2,'h')

[thr,sorh,keepapp]=ddencmp('den','wv',x)

xd3=wdencmp('gbl',c,l,wname,lev,thr,'s',1)

subplot(411)plot(x)title('原始信号','fontsize',12)

subplot(412)plot(xd1)title('使用penalty阈值降噪后信号','fontsize',12)

subplot(413)plot(xd2)title('使用Birge-Massart阈值降噪后信号','fontsize',12)

subplot(414)plot(xd3)title('使用缺省阈值降噪后信号','fontsize',12)

s=[-1.58 0.42 0.46 0.78 -0.49 0.59 -1.3 -1.42 -0.16 -1.47 -1.35 0.36 -0.44 -0.14 1 -0.5 -0.2 -0.06 -0.6 0.42 -1.52 0.51 0.76 -1.5 0.16 -1.29 -0.65 -1.48 0.6 -1.65 -0.55]

[C,L]=wavedec(s,1,'db3')

ca1=wrcoef('a',C,L,'db3',1)

x1=ca1

[C,L]=wavedec(s,2,'db3')

ca2=wrcoef('a',C,L,'db3',2)

x2=ca2

[C,L]=wavedec(s,3,'db3')

ca3=wrcoef('a',C,L,'db3',3)

x3=ca3

[C,L]=wavedec(s,4,'db3')

ca4=wrcoef('a',C,L,'db3',4)

x4=ca4

cg = wrcoef('a',C,L,'sym5',1)

x5=cg

p=1:31

subplot(6,1,1)plot(p,s)ylabel('s')

subplot(6,1,2)plot(p,x1)ylabel('ca1')

subplot(6,1,3)plot(p,x2)ylabel('ca2')

subplot(6,1,4)plot(p,x3)ylabel('ca3')

subplot(6,1,5)plot(p,x4)ylabel('ca4')

subplot(6,1,6)plot(p,x5)ylabel('ca5') %加入的重构，是不是你要的？

下载的工具箱中的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环境中运行，算法的主要功能函数并不受影响。

---