% 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> 指定次數的迭代%% FIX_H (int): 取消默認的停止條件，進行 <FIX_H> 指定次數的迭代，僅僅保留 |#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% gabriel.rilling@ens-lyon.fr% % 翻譯：xray	11.2007function [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 && nbit<MAXITERATIONS    if(~MODE_COMPLEX && nbit>MAXITERATIONS/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;		% 單輪迭代次數清0end % 主循環% 計入殘余信號if any(r) && ~any(mask)  imf(k,:) = r;end% 計數正交指數ort = io(x,imf);% 關閉圖形if display_sifting  closeendend%---------------------------------------------------------------------------------------------------% 測試是否存在足夠的極值點（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;endend%---------------------------------------------------------------------------------------------------% 計數包絡均值和模式幅度估計值，返回包絡均值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  endelse	% 實信號情況  [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;  	% 計算包絡幅度  endendend%-------------------------------------------------------------------------------% 默認停止條件，這是單輪迭代停止條件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  endcatch  stop = 1;  envmoy = zeros(1,length(m));  s = NaN;endend%-------------------------------------------------------------------------------% 針對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;endend%-------------------------------------------------------------------------------% 針對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);  endcatch  moyenne = zeros(1,length(m));  stop = 1;endend%-------------------------------------------------------------------------------% 演示分解過程（默認準則）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  offsubplot(4,1,2)plot(t,sx)hold onplot(t,sdt,'--r')plot(t,sd2t,':k')title('stop parameter')set(gca,'XTick',[])hold offsubplot(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')endif display_sifting == 2  pause(0.01)else