⑴ matlab編的有關EMD去噪的程序，處理的是核磁共振測井信號，急求
function imf = emd(x,n);%%最好把函數名改為emd1之類的，以免和Grilling的emd沖突
%%n為你想得到的IMF的個數
c = x('; % 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 resiefor 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 resie %%判斷信號是不是已經符合殘余分量定義
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 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; % 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 = [imf; h]; % store the extracted IMF in the matrix imf
% if size(maxmin,2)<2, then h is the resie
% 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