matlab编的有关EMD去噪的程序，处理的是核磁共振测井信号，急求

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

基础背景自行搜索

希尔伯特-黄变换包括两个步骤：

    1.使用 EMD 算法获得本质模态函数 (IMF)。

    2.通过对上一步获得的结果应用希尔伯特-黄变换得到初始序列的瞬时频率谱。HHT 能够获得非线性和非静态序列的清物郑瞬时频率谱。之后，可以使用经验模态分解处理这些序列。

本文重点讲解第一部分

从 EMD 生成的本质模态函数（ IMF） 应满足以下要求：

1.IMF 极值的数量（最大值和最小值的数量答颂之和）与零穿越的数量必须相等或最多相差 1；

2.在 IMF 的任意点，局部最大值定义的包络线的平均值和局部最小值定义的包络线的平均值应等于零。

1.求极值点

2.拟合包络函数

用三次样条插值法拟合包络

3.均值包络线  

将两条极值曲线平均获得平均包络线

4.获得中间信号

原始信号减均值包络线，得到中间信号。若中间信号中还存在负的局部极大值和正的局部极小值（此判据等效于是否满足上述经验模态分解条件），说明这还不是一个本征模函数IMF，需要继续进行“筛选”。筛选的过程就是以该中间信号为新的输入信号继续重复1~4的步骤。筛选过程通常在蚂凳残数只包含不超过两个极值时停止。

5.迭代

用上述方法得到第一个IMF后，用原始信号减IMF1，作为新的原始信号，再通过上述的1~4步骤，可以得到IMF2，以此类推，完成EMD分解。

6.分解完需要的IMF后，至此，EMD分解完成。

这篇文章能让你明白经验模态分解（EMD）——基础理论篇

经验模态分解法介绍

经验模态分解 (Empirical Mode Decomposition)

EMD算法PPT演示迭代过程

---