% EMD 計算經驗模式分解
%
%
%   語法
%
%
% IMF = EMD(X)
% IMF = EMD(X,...,'Option_name',Option_value,...)
% IMF = EMD(X,OPTS)
% [IMF,ORT,NB_ITERATIONS] = EMD(...)
%
%
%   描述
%
%
% IMF = EMD(X) X是一個實向量，計算方法參考[1]，計算結果包含在IMF矩陣中，每一行包含一個IMF分量，
% 最後一行是殘餘分量，預設的停止條件如下[2]：
%
%   在每一個點, mean_amplitude < THRESHOLD2*envelope_amplitude （注：平均幅度與包絡幅度的比值小於門限2）
%   &
%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE 
%  （注：平均幅度與包絡幅度比值大於門限的點數佔訊號總點數中的比例小於容限）
%   &
%   |#zeros-#extrema|<=1 （注：過零點和極值點個數相等或者相差1）
%
% 這裡 mean_amplitude = abs(envelope_max+envelope_min)/2 （注：平均幅度等於上下包絡相互抵消後殘差的一半的絕對值，理想情況等於0）
% 且 envelope_amplitude = abs(envelope_max-envelope_min)/2 （注：包絡幅度等於上下包絡相對距離的一半，理想情況等於上下包絡本身的絕對值）
% 
% IMF = EMD(X) X是一個實向量，計算方法參考[3]，計算結果包含在IMF矩陣中，每一行包含一個IMF分量，
% 最後一行是殘餘分量，預設的停止條件如下[2]：
%
%   在每一個點, mean_amplitude < THRESHOLD2*envelope_amplitude（注：平均幅度與包絡幅度的比值小於門限2）
%   &
%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE
%  （注：平均幅度與包絡幅度比值大於門限的點數佔訊號總點數中的比例小於容限）
%
% 這裡平均幅度和包絡幅度的定義與前面實數情況下類似
%
% IMF = EMD(X,...,'Option_name',Option_value,...) 設定特定引數（見選項）
%
% IMF = EMD(X,OPTS) 與前面等價，只是這裡OPTS是一個結構體，其中每一個域名與相應的選項名稱一致。
%
% [IMF,ORT,NB_ITERATIONS] = EMD(...) 返回正交指數
%                       ________
%         _  |IMF(i,:).*IMF(j,:)|
%   ORT = \ _____________________
%         /
%         -       || X ||^2        i~=j
%
% 和提取每一個IMF時進行的迭代次數。
%
%
%   選擇
%
%
%  停止條件選項:
%
% STOP: 停止引數 [THRESHOLD,THRESHOLD2,TOLERANCE]
% 如果輸入向量長度小於 3, 只有第一個引數有效，其他引數採用預設值
% 預設值: [0.05,0.5,0.05]
%
% FIX (int): 取消預設的停止條件，進行  指定次數的迭代
%
% FIX_H (int): 取消預設的停止條件，進行  指定次數的迭代，僅僅保留 |#zeros-#extrema|<=1 的停止條件，參考 [4]
%
%  復 EMD 選項:
%
% COMPLEX_VERSION: 選擇復 EMD 演算法(參考[3])
% COMPLEX_VERSION = 1: "algorithm 1"
% COMPLEX_VERSION = 2: "algorithm 2" (default)
% 
% NDIRS: 包絡計算的方向個數 (預設 4)
% rem: 實際方向個數 (根據 [3]) 是 2*NDIRS
% 
%  其他選項:
%
% T: 取樣時刻 (線性向量) (預設: 1:length(x))
%
% MAXITERATIONS: 提取每個IMF中，採用的最大迭代次數（預設：2000）
%
% MAXMODES: 提取IMFs的最大個數 (預設: Inf)
%
% DISPLAY: 如果等於1，每迭代一次自動暫停（pause）
% 如果等於2，迭代過程不暫停 (動畫模式)
% rem: 當輸入是複數的時候，演示過程自動取消
%
% INTERP: 插值方法 'linear', 'cubic', 'pchip' or 'spline' (預設)
% 詳情見 interp1 文件
%
% MASK: 採用 masking 訊號，參考 [5]
%
%
%   例子
%
%
% X = rand(1,512);
%
% IMF = emd(X);
%
% IMF = emd(X,'STOP',[0.1,0.5,0.05],'MAXITERATIONS',100);
%
% T = linspace(0,20,1e3);
% X = 2*exp(i*T)+exp(3*i*T)+.5*T;
% IMF = emd(X,'T',T);
%
% OPTIONS.DISLPAY = 1;
% OPTIONS.FIX = 10;
% OPTIONS.MAXMODES = 3;
% [IMF,ORT,NBITS] = emd(X,OPTIONS);
%
%
%   參考文獻
%
%
% [1] N. E. Huang et al., "The empirical mode decomposition and the
% Hilbert spectrum for non-linear and non stationary time series analysis",
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% [2] G. Rilling, P. Flandrin and P. Goncalves
% "On Empirical Mode Decomposition and its algorithms",
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% [3] G. Rilling, P. Flandrin, P. Goncalves and J. M. Lilly.,
% "Bivariate Empirical Mode Decomposition",
% Signal Processing Letters (submitted)
%
% [4] N. E. Huang et al., "A confidence limit for the Empirical Mode
% Decomposition and Hilbert spectral analysis",
% Proc. Royal Soc. London A, Vol. 459, pp. 2317-2345, 2003
%
% [5] R. Deering and J. F. Kaiser, "The use of a masking signal to improve 
% empirical mode decomposition", ICASSP 2005
%
%
% 也可以參考
%  emd_visu (visualization),
%  emdc, emdc_fix (fast implementations of EMD),
%  cemdc, cemdc_fix, cemdc2, cemdc2_fix (fast implementations of bivariate EMD),
%  hhspectrum (Hilbert-Huang spectrum)
%
%
% G. Rilling, 最後修改: 3.2007
% 
 
[email protected]
% 
% 翻譯：xray	11.2007

function [imf,ort,nbits] = emd(varargin)
% 採用可變引數輸入

% 處理輸入引數
[x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:});
% 引數說明：
% x 訊號
% t 時間向量
% sd 門限
% sd2 門限2
% tol 容限值
% MODE_COMPLEX 是否處理覆信號
% ndirs 方向個數
% display_sifting 是否演示迭代過程
% sdt 將門限擴充套件為跟訊號長度一樣的向量
% sd2t 將門限2擴充套件為跟訊號長度一樣的向量
% r 等於x
% imf 如果使用mask訊號，此時IMF已經得到了
% k 記錄已經提取的IMF個數
% nbit 記錄提取每一個IMF時迭代的次數
% NbIt 記錄迭代的總次數
% MAXITERATIONS 提取每個IMF時採用的最大迭代次數
% FIXE 進行指定次數的迭代
% FIXE_H 進行指定次數的迭代，且保留 |#zeros-#extrema|<=1 的停止條件
% MAXMODES 提取的最大IMF個數
% INTERP 插值方法
% mask mask訊號

% 如果要求演示迭代過程，用 fig_h 儲存當前圖形視窗控制代碼
if display_sifting
  fig_h = figure;
end

% 主迴圈 : 至少要求存在3個極值點，如果採用mask訊號，不進入主迴圈
while ~stop_EMD(r,MODE_COMPLEX,ndirs) && (k < MAXMODES+1 || MAXMODES == 0) && ~any(mask)

  % 當前模式
  m = r;

  % 前一次迭代的模式
  mp = m;

  % 計算均值和停止條件
  if FIXE % 如果設定了迭代次數
    [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);
  elseif FIXE_H % 如果設定了迭代次數，且保留停止條件|#zeros-#extrema|<=1
    stop_count = 0;
    [stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);
  else % 採用預設停止條件
    [stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);
  end

  % 當前模式幅度過小，機器精度就可能引起虛假極值點的出現
  if (max(abs(m))) < (1e-10)*(max(abs(x)))	% IMF的最大值小於訊號最大值的1e-10
    if ~stop_sift % 如果篩過程沒有停止
      warning('emd:warning','forced stop of EMD : too small amplitude')
    else
      disp('forced stop of EMD : too small amplitude')
    end
    break
  end


  % 篩迴圈
  while ~stop_sift && nbitMAXITERATIONS/5 && mod(nbit,floor(MAXITERATIONS/10))==0 && ~FIXE && nbit > 100)
      disp(['mode ',int2str(k),', iteration ',int2str(nbit)])
      if exist('s','var')
        disp(['stop parameter mean value : ',num2str(s)])
      end
      [im,iM] = extr(m);
      disp([int2str(sum(m(im) > 0)),' minima > 0; ',int2str(sum(m(iM) < 0)),' maxima < 0.'])
    end

    % 篩過程
    m = m - moyenne;

    % 計算均值和停止條件
    if FIXE
      [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);
    elseif FIXE_H
      [stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);
    else
      [stop_sift,moyenne,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);
    end

    % 演示過程
    if display_sifting && ~MODE_COMPLEX
      NBSYM = 2;
      [indmin,indmax] = extr(mp);
      [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,mp,mp,NBSYM);
      envminp = interp1(tmin,mmin,t,INTERP);
      envmaxp = interp1(tmax,mmax,t,INTERP);
      envmoyp = (envminp+envmaxp)/2;
      if FIXE || FIXE_H
        display_emd_fixe(t,m,mp,r,envminp,envmaxp,envmoyp,nbit,k,display_sifting)
      else
        sxp = 2*(abs(envmoyp))./(abs(envmaxp-envminp));
        sp = mean(sxp);
        display_emd(t,m,mp,r,envminp,envmaxp,envmoyp,s,sp,sxp,sdt,sd2t,nbit,k,display_sifting,stop_sift)
      end
    end

    mp = m;
    nbit = nbit+1;	% 單輪迭代計數
    NbIt = NbIt+1;	% 總體迭代計數

    if (nbit==(MAXITERATIONS-1) && ~FIXE && nbit > 100)
      if exist('s','var')
        warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
      else
        warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'.'])
      end
    end

  end % 篩迴圈
  
  imf(k,:) = m;
  if display_sifting
    disp(['mode ',int2str(k),' stored'])
  end
  nbits(k) = nbit;	% 記錄每個IMF的迭代次數
  k = k+1;		% IMF計數


  r = r - m;		% 從原訊號中減去提取的IMF
  nbit = 0;		% 單輪迭代次數清0


end % 主迴圈

% 計入殘餘訊號
if any(r) && ~any(mask)
  imf(k,:) = r;
end

% 計數正交指數
ort = io(x,imf);

% 關閉圖形
if display_sifting
  close
end

end

%---------------------------------------------------------------------------------------------------
% 測試是否存在足夠的極值點（3個）進行分解，極值點個數小於3個則返回1，這是整體停止條件
function stop = stop_EMD(r,MODE_COMPLEX,ndirs)
if MODE_COMPLEX  % 覆信號情況
  for k = 1:ndirs
    phi = (k-1)*pi/ndirs;
    [indmin,indmax] = extr(real(exp(i*phi)*r));
    ner(k) = length(indmin) + length(indmax);
  end
  stop = any(ner < 3);
else % 實訊號情況
  [indmin,indmax] = extr(r);
  ner = length(indmin) + length(indmax);
  stop = ner < 3;
end
end

%---------------------------------------------------------------------------------------------------
% 計數包絡均值和模式幅度估計值，返回包絡均值
function [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs)
NBSYM = 2;	% 邊界延拓點數
if MODE_COMPLEX		% 覆信號情況
  switch MODE_COMPLEX
    case 1
      for k = 1:ndirs
        phi = (k-1)*pi/ndirs;
        y = real(exp(-i*phi)*m);
        [indmin,indmax,indzer] = extr(y);
        nem(k) = length(indmin)+length(indmax);
        nzm(k) = length(indzer);
        [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,m,NBSYM);
        envmin(k,:) = interp1(tmin,zmin,t,INTERP);
        envmax(k,:) = interp1(tmax,zmax,t,INTERP);
      end
      envmoy = mean((envmin+envmax)/2,1);
      if nargout > 3
        amp = mean(abs(envmax-envmin),1)/2;
      end
    case 2
      for k = 1:ndirs
        phi = (k-1)*pi/ndirs;
        y = real(exp(-i*phi)*m);
        [indmin,indmax,indzer] = extr(y);
        nem(k) = length(indmin)+length(indmax);
        nzm(k) = length(indzer);
        [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,y,NBSYM);
        envmin(k,:) = exp(i*phi)*interp1(tmin,zmin,t,INTERP);
        envmax(k,:) = exp(i*phi)*interp1(tmax,zmax,t,INTERP);
      end
      envmoy = mean((envmin+envmax),1);
      if nargout > 3
        amp = mean(abs(envmax-envmin),1)/2;
      end
  end
else	% 實訊號情況
  [indmin,indmax,indzer] = extr(m);	% 計數最小值、最大值和過零點位置
  nem = length(indmin)+length(indmax);
  nzm = length(indzer);
  [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,m,m,NBSYM);	% 邊界延拓
  envmin = interp1(tmin,mmin,t,INTERP);
  envmax = interp1(tmax,mmax,t,INTERP);
  envmoy = (envmin+envmax)/2;
  if nargout > 3
    amp = mean(abs(envmax-envmin),1)/2;  	% 計算包絡幅度
  end
end
end

%-------------------------------------------------------------------------------
% 預設停止條件，這是單輪迭代停止條件
function [stop,envmoy,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs)
try
  [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);
  sx = abs(envmoy)./amp;
  s = mean(sx);
  stop = ~((mean(sx > sd) > tol | any(sx > sd2)) & (all(nem > 2)));  % 停止準則（增加了極值點個數大於2）
  if ~MODE_COMPLEX
    stop = stop && ~(abs(nzm-nem)>1);	% 對於實訊號，要求極值點和過零點的個數相差1
  end
catch
  stop = 1;
  envmoy = zeros(1,length(m));
  s = NaN;
end
end

%-------------------------------------------------------------------------------
% 針對FIX選項的停止條件
function [stop,moyenne]= stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs)
try
  moyenne = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);	% 正常情況下不會導致停止
  stop = 0;
catch
  moyenne = zeros(1,length(m));
  stop = 1;
end
end

%-------------------------------------------------------------------------------
% 針對FIX_H選項的停止條件
function [stop,moyenne,stop_count]= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs)
try
  [moyenne,nem,nzm] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);
  if (all(abs(nzm-nem)>1))
    stop = 0;
    stop_count = 0;
  else	% 極值點與過零點個數相差1後，還要達到指定次數才停止
    stop_count = stop_count+1;
    stop = (stop_count == FIXE_H);
  end
catch
  moyenne = zeros(1,length(m));
  stop = 1;
end
end

%-------------------------------------------------------------------------------
% 演示分解過程（默認準則）
function display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift)

subplot(4,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold  off
subplot(4,1,2)
plot(t,sx)
hold on
plot(t,sdt,'--r')
plot(t,sd2t,':k')
title('stop parameter')
set(gca,'XTick',[])
hold off
subplot(4,1,3)
plot(t,m)
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(4,1,4);
plot(t,r-m)
title('residue');
disp(['stop parameter mean value : ',num2str(sb),' before sifting and ',num2str(s),' after'])
if stop_sift
  disp('last iteration for this mode')
end
if display_sifting == 2
  pause(0.01)
else
  pause
end
end

%---------------------------------------------------------------------------------------------------
% 演示分解過程（FIX和FIX_H停止準則）
function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)
subplot(3,1,1)
plot(t,mp);hold on;
plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
set(gca,'XTick',[])
hold  off
subplot(3,1,2)
plot(t,m)
title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(3,1,3);
plot(t,r-m)
title('residue');
if display_sifting == 2
  pause(0.01)
else
  pause
end
end

%---------------------------------------------------------------------------------------
% 處理邊界條件（映象法）
function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
% 實數情況下，x = z

lx = length(x);

% 判斷極值點個數
if (length(indmin) + length(indmax) < 3)
  error('not enough extrema')
end

% 插值的邊界條件
if indmax(1) < indmin(1)	% 第一個極值點是極大值
  if x(1) > x(indmin(1))	% 以第一個極大值為對稱中心
    lmax = fliplr(indmax(2:min(end,nbsym+1)));
    lmin = fliplr(indmin(1:min(end,nbsym)));
    lsym = indmax(1);
  else	% 如果第一個取樣值小於第一個極小值，則將認為該值是一個極小值，以該點為對稱中心
    lmax = fliplr(indmax(1:min(end,nbsym)));
    lmin = [fliplr(indmin(1:min(end,nbsym-1))),1];
    lsym = 1;
  end
else
  if x(1) < x(indmax(1))	% 以第一個極小值為對稱中心
    lmax = fliplr(indmax(1:min(end,nbsym)));
    lmin = fliplr(indmin(2:min(end,nbsym+1)));
    lsym = indmin(1);
  else  % 如果第一個取樣值大於第一個極大值，則將認為該值是一個極大值，以該點為對稱中心
    lmax = [fliplr(indmax(1:min(end,nbsym-1))),1];
    lmin = fliplr(indmin(1:min(end,nbsym)));
    lsym = 1;
  end
end

% 序列末尾情況與序列開頭類似
if indmax(end) < indmin(end)
  if x(end) < x(indmax(end))
    rmax = fliplr(indmax(max(end-nbsym+1,1):end));
    rmin = fliplr(indmin(max(end-nbsym,1):end-1));
    rsym = indmin(end);
  else
    rmax = [lx,fliplr(indmax(max(end-nbsym+2,1):end))];
    rmin = fliplr(indmin(max(end-nbsym+1,1):end));
    rsym = lx;
  end
else
  if x(end) > x(indmin(end))
    rmax = fliplr(indmax(max(end-nbsym,1):end-1));
    rmin = fliplr(indmin(max(end-nbsym+1,1):end));
    rsym = indmax(end);
  else
    rmax = fliplr(indmax(max(end-nbsym+1,1):end));
    rmin = [lx,fliplr(indmin(max(end-nbsym+2,1):end))];
    rsym = lx;
  end
end
    
% 將序列根據對稱中心，映象到兩邊
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
    
% 如果對稱的部分沒有足夠的極值點
if tlmin(1) > t(1) || tlmax(1) > t(1)	% 對摺後的序列沒有超出原序列的範圍
  if lsym == indmax(1)
    lmax = fliplr(indmax(1:min(end,nbsym)));
  else
    lmin = fliplr(indmin(1:min(end,nbsym)));
  end
  if lsym == 1	% 這種情況不應該出現，程式直接中止
    error('bug')
  end
  lsym = 1;	% 直接關於第一取樣點取映象
  tlmin = 2*t(lsym)-t(lmin);
  tlmax = 2*t(lsym)-t(lmax);
end   
    
% 序列末尾情況與序列開頭類似
if trmin(end) < t(lx) || trmax(end) < t(lx)
  if rsym == indmax(end)
    rmax = fliplr(indmax(max(end-nbsym+1,1):end));
  else
    rmin = fliplr(indmin(max(end-nbsym+1,1):end));
  end
  if rsym == lx
    error('bug')
  end
  rsym = lx;
  trmin = 2*t(rsym)-t(rmin);
  trmax = 2*t(rsym)-t(rmax);
end 

% 延拓點上的取值       
zlmax = z(lmax); 
zlmin = z(lmin);
zrmax = z(rmax); 
zrmin = z(rmin);
     
% 完成延拓
tmin = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
zmin = [zlmin z(indmin) zrmin];
zmax = [zlmax z(indmax) zrmax];

end
    
%---------------------------------------------------------------------------------------------------
% 極值點和過零點位置提取
function [indmin, indmax, indzer] = extr(x,t)

if(nargin==1)
  t = 1:length(x);
end

m = length(x);

if nargout > 2
  x1 = x(1:m-1);
  x2 = x(2:m);
  indzer = find(x1.*x2<0);	% 尋找訊號符號發生變化的位置

  if any(x == 0)	% 考慮訊號取樣點恰好為0的位置
    iz = find( x==0 );  % 訊號取樣點恰好為0的位置
    indz = [];
    if any(diff(iz)==1) % 出現連0的情況
      zer = x == 0;	% x=0處為1，其它地方為0
      dz = diff([0 zer 0]);	% 尋找0與非0的過渡點
      debz = find(dz == 1);	% 0值起點
      finz = find(dz == -1)-1;  % 0值終點
      indz = round((debz+finz)/2);	% 選擇中間點作為過零點
    else
      indz = iz;	% 若沒有連0的情況，該點本身就是過零點
    end
    indzer = sort([indzer indz]);	% 全體過零點排序
  end
end

% 提取極值點
d = diff(x);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;	% 最小值
indmax = find(d1.*d2<0 & d1>0)+1;	% 最大值


% 當連續多個取樣值相同時，把最中間的一個值作為極值點，處理方式與連0類似
if any(d==0)

  imax = [];
  imin = [];

  bad = (d==0);
  dd = diff([0 bad 0]);
  debs = find(dd == 1);
  fins = find(dd == -1);
  if debs(1) == 1	% 連續值出現在序列開頭
    if length(debs) > 1
      debs = debs(2:end);
      fins = fins(2:end);
    else
      debs = [];
      fins = [];
    end
  end
  if length(debs) > 0
    if fins(end) == m	% 連續值出現在序列末尾
      if length(debs) > 1
        debs = debs(1:(end-1));
        fins = fins(1:(end-1));

      else
        debs = [];
        fins = [];
      end
    end
  end
  lc = length(debs);
  if lc > 0
    for k = 1:lc
      if d(debs(k)-1) > 0	% 取中間值
        if d(fins(k)) < 0
          imax = [imax round((fins(k)+debs(k))/2)];
        end
      else
        if d(fins(k)) > 0
          imin = [imin round((fins(k)+debs(k))/2)];
        end
      end
    end
  end

  if length(imax) > 0
    indmax = sort([indmax imax]);
  end

  if length(imin) > 0
    indmin = sort([indmin imin]);
  end

end
end

%---------------------------------------------------------------------------------------------------

function ort = io(x,imf)
% ort = IO(x,imf) 計算正交指數
%
% 輸入 : - x    : 分析訊號
%        - imf  : IMF訊號

n = size(imf,1);

s = 0;
% 根據公式計算
for i = 1:n
  for j = 1:n
    if i ~= j
      s = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));
    end
  end
end

ort = 0.5*s;
end

%---------------------------------------------------------------------------------------------------
% 函式引數解析
function [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin)

x = varargin{1};
if nargin == 2
  if isstruct(varargin{2})
    inopts = varargin{2};
  else
    error('when using 2 arguments the first one is the analyzed signal X and the second one is a struct object describing the options')
  end
elseif nargin > 2
  try
    inopts = struct(varargin{2:end});
  catch
    error('bad argument syntax')
  end
end

% 預設停止條件
defstop = [0.05,0.5,0.05];

opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask','ndirs','complex_version'};
% 時間序列，停止引數，是否演示，最大迭代次數，每一輪迭代次數，IMF個數，插值方法，每一輪迭代次數（帶條件），mask訊號，方向數，是否採用複數模式

defopts.stop = defstop;
defopts.display = 0;
defopts.t = 1:max(size(x));
defopts.maxiterations = 2000;
defopts.fix = 0;
defopts.maxmodes = 0;
defopts.interp = 'spline';
defopts.fix_h = 0;
defopts.mask = 0;
defopts.ndirs = 4;
defopts.complex_version = 2;

opts = defopts;

if(nargin==1)
  inopts = defopts;
elseif nargin == 0
  error('not enough arguments')
end

names = fieldnames(inopts);
for nom = names'
  if ~any(strcmpi(char(nom), opt_fields))
    error(['bad option field name: ',char(nom)])
  end
  if ~isempty(eval(['inopts.',char(nom)])) 
    eval(['opts.',lower(char(nom)),' = inopts.',char(nom),';'])
  end
end

t = opts.t;
stop = opts.stop;
display_sifting = opts.display;
MAXITERATIONS = opts.maxiterations;
FIXE = opts.fix;
MAXMODES = opts.maxmodes;
INTERP = opts.interp;
FIXE_H = opts.fix_h;
mask = opts.mask;
ndirs = opts.ndirs;
complex_version = opts.complex_version;

if ~isvector(x)
  error('X must have only one row or one column')
end

if size(x,1) > 1
  x = x.';
end

if ~isvector(t)
  error('option field T must have only one row or one column')
end

if ~isreal(t)
  error('time instants T must be a real vector')
end

if size(t,1) > 1
  t = t';
end

if (length(t)~=length(x))
  error('X and option field T must have the same length')
end

if ~isvector(stop) || length(stop) > 3
  error('option field STOP must have only one row or one column of max three elements')
end

if ~all(isfinite(x))
  error('data elements must be finite')
end

if size(stop,1) > 1
  stop = stop';
end

L = length(stop);
if L < 3
  stop(3) = defstop(3);
end

if L < 2
  stop(2) = defstop(2);
end

if ~ischar(INTERP) || ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
  error('INTERP field must be ''linear'', ''cubic'', ''pchip'' or ''spline''')
end

% 使用mask訊號時的特殊處理
if any(mask)
  if ~isvector(mask) || length(mask) ~= length(x)
    error('masking signal must have the same dimension as the analyzed signal X')
  end

  if size(mask,1) > 1
    mask = mask.';
  end
  opts.mask = 0;
  imf1 = emd(x+mask, opts);
  imf2 = emd(x-mask, opts);
  if size(imf1,1) ~= size(imf2,1)
    warning('emd:warning',['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])
  end
  S1 = size(imf1,1);
  S2 = size(imf2,1);
  if S1 ~= S2	% 如果兩個訊號分解得到的IMF個數不一致，調整順序
    if S1 < S2
      tmp = imf1;
      imf1 = imf2;
      imf2 = tmp;
    end
    imf2(max(S1,S2),1) = 0;	% 將短的那個補零，達到長度一致
  end
  imf = (imf1+imf2)/2;

end


sd = stop(1);
sd2 = stop(2);
tol = stop(3);

lx = length(x);

sdt = sd*ones(1,lx);
sd2t = sd2*ones(1,lx);

if FIXE
  MAXITERATIONS = FIXE;
  if FIXE_H
    error('cannot use both ''FIX'' and ''FIX_H'' modes')
  end
end

MODE_COMPLEX = ~isreal(x)*complex_version;
if MODE_COMPLEX && complex_version ~= 1 && complex_version ~= 2
  error('COMPLEX_VERSION parameter must equal 1 or 2')
end


% 極值點和過零點的個數
ner = lx;
nzr = lx;

r = x;

if ~any(mask) % 如果使用了mask訊號，此時imf就已經計算得到了
  imf = [];
end
k = 1;

% 提取每個模式時迭代的次數
nbit = 0;

% 總體迭代次數
NbIt = 0;
end
%---------------------------------------------------------------------------------------------------
 

宋老師的版本，留著備用，後續實現c++移植