% EMD 計(jì)算經(jīng)驗(yàn)?zāi)Ｊ椒纸?/span>%%%   語(yǔ)法%%% IMF = EMD(X)% IMF = EMD(X,...,'Option_name',Option_value,...)% IMF = EMD(X,OPTS)% [IMF,ORT,NB_ITERATIONS] = EMD(...)%%%   描述%%% IMF = EMD(X) X是一個(gè)實(shí)矢量，計(jì)算方法參考[1]，計(jì)算結(jié)果包含在IMF矩陣中，每一行包含一個(gè)IMF分量，% 最后一行是殘余分量，默認(rèn)的停止條件如下[2]：%%   在每一個(gè)點(diǎn), mean_amplitude < THRESHOLD2*envelope_amplitude （注：平均幅度與包絡(luò)幅度的比值小于門(mén)限2）%   &%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE %  （注：平均幅度與包絡(luò)幅度比值大于門(mén)限的點(diǎn)數(shù)占信號(hào)總點(diǎn)數(shù)中的比例小于容限）%   &%   |#zeros-#extrema|<=1 （注：過(guò)零點(diǎn)和極值點(diǎn)個(gè)數(shù)相等或者相差1）%% 這里 mean_amplitude = abs(envelope_max+envelope_min)/2 （注：平均幅度等于上下包絡(luò)相互抵消后殘差的一半的絕對(duì)值，理想情況等于0）% 且 envelope_amplitude = abs(envelope_max-envelope_min)/2 （注：包絡(luò)幅度等于上下包絡(luò)相對(duì)距離的一半，理想情況等于上下包絡(luò)本身的絕對(duì)值）% % IMF = EMD(X) X是一個(gè)實(shí)矢量，計(jì)算方法參考[3]，計(jì)算結(jié)果包含在IMF矩陣中，每一行包含一個(gè)IMF分量，% 最后一行是殘余分量，默認(rèn)的停止條件如下[2]：%%   在每一個(gè)點(diǎn), mean_amplitude < THRESHOLD2*envelope_amplitude（注：平均幅度與包絡(luò)幅度的比值小于門(mén)限2）%   &%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE%  （注：平均幅度與包絡(luò)幅度比值大于門(mén)限的點(diǎn)數(shù)占信號(hào)總點(diǎn)數(shù)中的比例小于容限）%% 這里平均幅度和包絡(luò)幅度的定義與前面實(shí)數(shù)情況下類似%% IMF = EMD(X,...,'Option_name',Option_value,...) 設(shè)置特定參數(shù)（見(jiàn)選項(xiàng)）%% IMF = EMD(X,OPTS) 與前面等價(jià)，只是這里OPTS是一個(gè)結(jié)構(gòu)體，其中每一個(gè)域名與相應(yīng)的選項(xiàng)名稱一致。%% [IMF,ORT,NB_ITERATIONS] = EMD(...) 返回正交指數(shù)%                       ________%         _  |IMF(i,:).*IMF(j,:)|%   ORT = \ _____________________%         /%         -       || X ||^2        i~=j%% 和提取每一個(gè)IMF時(shí)進(jìn)行的迭代次數(shù)。%%%   選擇%%%  停止條件選項(xiàng):%% STOP: 停止參數(shù) [THRESHOLD,THRESHOLD2,TOLERANCE]% 如果輸入矢量長(zhǎng)度小于 3, 只有第一個(gè)參數(shù)有效，其他參數(shù)采用默認(rèn)值% 默認(rèn)值: [0.05,0.5,0.05]%% FIX (int): 取消默認(rèn)的停止條件，進(jìn)行 <FIX> 指定次數(shù)的迭代%% FIX_H (int): 取消默認(rèn)的停止條件，進(jìn)行 <FIX_H> 指定次數(shù)的迭代，僅僅保留 |#zeros-#extrema|<=1 的停止條件，參考 [4]%%  復(fù) EMD 選項(xiàng):%% COMPLEX_VERSION: 選擇復(fù) EMD 算法(參考[3])% COMPLEX_VERSION = 1: "algorithm 1"% COMPLEX_VERSION = 2: "algorithm 2" (default)% % NDIRS: 包絡(luò)計(jì)算的方向個(gè)數(shù) (默認(rèn) 4)% rem: 實(shí)際方向個(gè)數(shù) (根據(jù) [3]) 是 2*NDIRS% %  其他選項(xiàng):%% T: 采樣時(shí)刻 (線性矢量) (默認(rèn): 1:length(x))%% MAXITERATIONS: 提取每個(gè)IMF中，采用的最大迭代次數(shù)（默認(rèn)：2000）%% MAXMODES: 提取IMFs的最大個(gè)數(shù) (默認(rèn): Inf)%% DISPLAY: 如果等于1，每迭代一次自動(dòng)暫停（pause）% 如果等于2，迭代過(guò)程不暫停 (動(dòng)畫(huà)模式)% rem: 當(dāng)輸入是復(fù)數(shù)的時(shí)候，演示過(guò)程自動(dòng)取消%% INTERP: 插值方法 'linear', 'cubic', 'pchip' or 'spline' (默認(rèn))% 詳情見(jiàn) interp1 文檔%% MASK: 采用 masking 信號(hào)，參考 [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);%%%   參考文獻(xiàn)%%% [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)% 采用可變參數(shù)輸入% 處理輸入?yún)?shù)[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{:});% 參數(shù)說(shuō)明：% x 信號(hào)% t 時(shí)間矢量% sd 門(mén)限% sd2 門(mén)限2% tol 容限值% MODE_COMPLEX 是否處理復(fù)信號(hào)% ndirs 方向個(gè)數(shù)% display_sifting 是否演示迭代過(guò)程% sdt 將門(mén)限擴(kuò)展為跟信號(hào)長(zhǎng)度一樣的矢量% sd2t 將門(mén)限2擴(kuò)展為跟信號(hào)長(zhǎng)度一樣的矢量% r 等于x% imf 如果使用mask信號(hào)，此時(shí)IMF已經(jīng)得到了% k 記錄已經(jīng)提取的IMF個(gè)數(shù)% nbit 記錄提取每一個(gè)IMF時(shí)迭代的次數(shù)% NbIt 記錄迭代的總次數(shù)% MAXITERATIONS 提取每個(gè)IMF時(shí)采用的最大迭代次數(shù)% FIXE 進(jìn)行指定次數(shù)的迭代% FIXE_H 進(jìn)行指定次數(shù)的迭代，且保留 |#zeros-#extrema|<=1 的停止條件% MAXMODES 提取的最大IMF個(gè)數(shù)% INTERP 插值方法% mask mask信號(hào)% 如果要求演示迭代過(guò)程，用 fig_h 保存當(dāng)前圖形窗口句柄if display_sifting  fig_h = figure;end% 主循環(huán) : 至少要求存在3個(gè)極值點(diǎn)，如果采用mask信號(hào)，不進(jìn)入主循環(huán)while ~stop_EMD(r,MODE_COMPLEX,ndirs) && (k < MAXMODES+1 || MAXMODES == 0) && ~any(mask)  % 當(dāng)前模式  m = r;  % 前一次迭代的模式  mp = m;  % 計(jì)算均值和停止條件  if FIXE % 如果設(shè)定了迭代次數(shù)    [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);  elseif FIXE_H % 如果設(shè)定了迭代次數(shù)，且保留停止條件|#zeros-#extrema|<=1    stop_count = 0;    [stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);  else % 采用默認(rèn)停止條件    [stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);  end  % 當(dāng)前模式幅度過(guò)小，機(jī)器精度就可能引起虛假極值點(diǎn)的出現(xiàn)  if (max(abs(m))) < (1e-10)*(max(abs(x)))	% IMF的最大值小于信號(hào)最大值的1e-10    if ~stop_sift % 如果篩過(guò)程沒(méi)有停止      warning('emd:warning','forced stop of EMD : too small amplitude')    else      disp('forced stop of EMD : too small amplitude')    end    break  end  % 篩循環(huán)  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    % 篩過(guò)程    m = m - moyenne;    % 計(jì)算均值和停止條件    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    % 演示過(guò)程    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;	% 單輪迭代計(jì)數(shù)    NbIt = NbIt+1;	% 總體迭代計(jì)數(shù)    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 % 篩循環(huán)    imf(k,:) = m;  if display_sifting    disp(['mode ',int2str(k),' stored'])  end  nbits(k) = nbit;	% 記錄每個(gè)IMF的迭代次數(shù)  k = k+1;		% IMF計(jì)數(shù)  r = r - m;		% 從原信號(hào)中減去提取的IMF  nbit = 0;		% 單輪迭代次數(shù)清0end % 主循環(huán)% 計(jì)入殘余信號(hào)if any(r) && ~any(mask)  imf(k,:) = r;end% 計(jì)數(shù)正交指數(shù)ort = io(x,imf);% 關(guān)閉圖形if display_sifting  closeendend%---------------------------------------------------------------------------------------------------% 測(cè)試是否存在足夠的極值點(diǎn)（3個(gè)）進(jìn)行分解，極值點(diǎn)個(gè)數(shù)小于3個(gè)則返回1，這是整體停止條件function stop = stop_EMD(r,MODE_COMPLEX,ndirs)if MODE_COMPLEX  % 復(fù)信號(hào)情況  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 % 實(shí)信號(hào)情況  [indmin,indmax] = extr(r);  ner = length(indmin) + length(indmax);  stop = ner < 3;endend%---------------------------------------------------------------------------------------------------% 計(jì)數(shù)包絡(luò)均值和模式幅度估計(jì)值，返回包絡(luò)均值function [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs)NBSYM = 2;	% 邊界延拓點(diǎn)數(shù)if MODE_COMPLEX		% 復(fù)信號(hào)情況  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	% 實(shí)信號(hào)情況  [indmin,indmax,indzer] = extr(m);	% 計(jì)數(shù)最小值、最大值和過(guò)零點(diǎn)位置  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;  	% 計(jì)算包絡(luò)幅度  endendend%-------------------------------------------------------------------------------% 默認(rèn)停止條件，這是單輪迭代停止條件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)));  % 停止準(zhǔn)則（增加了極值點(diǎn)個(gè)數(shù)大于2）  if ~MODE_COMPLEX    stop = stop && ~(abs(nzm-nem)>1);	% 對(duì)于實(shí)信號(hào)，要求極值點(diǎn)和過(guò)零點(diǎn)的個(gè)數(shù)相差1  endcatch  stop = 1;  envmoy = zeros(1,length(m));  s = NaN;endend%-------------------------------------------------------------------------------% 針對(duì)FIX選項(xiàng)的停止條件function [stop,moyenne]= stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs)try  moyenne = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);	% 正常情況下不會(huì)導(dǎo)致停止  stop = 0;catch  moyenne = zeros(1,length(m));  stop = 1;endend%-------------------------------------------------------------------------------% 針對(duì)FIX_H選項(xiàng)的停止條件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	% 極值點(diǎn)與過(guò)零點(diǎn)個(gè)數(shù)相差1后，還要達(dá)到指定次數(shù)才停止    stop_count = stop_count+1;    stop = (stop_count == FIXE_H);  endcatch  moyenne = zeros(1,length(m));  stop = 1;endend%-------------------------------------------------------------------------------% 演示分解過(guò)程（默認(rèn)準(zhǔn)則）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