matlab_混沌信号经过EMD分解，得到的IMF序列，求每个IMF序列与原始序列的互相关系数

% 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 % [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 && 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 warnin