CEEMDAN CEEMD EMDD EMD算法附matlab代码.zip

function [varargout] = SAM_EMD(x,IMF_num,varargin) %较原版emd，只增加了输出IMF数量的控制 %如果不使用varargin，可以不用输入IMF_num，否则，必须输入IMF_num，当不理睬时，输入0 %检测输入参数数量（1-无穷），输出参数数量（1-3） ssc=SAM_SSD(x,1024); narginchk(1,inf); nargoutchk(0,3); if nargin==1 IMF_num=0; end %验证和读取输入参数 [x,t,td,isTT,opt] = parseAndValidateInputs(x, varargin{:}); %输出 IMF 残差 信息 [IMF, residual, info] = localEMD(x,t,opt,IMF_num); if(isTT) t = td; end if (nargout == 0 && coder.target('MATLAB')) signalwavelet.internal.guis.plot.emdPlot(x, IMF, residual, t); end if(isTT && coder.target('MATLAB')) IMF = array2timetable(IMF,'RowTimes',t); residual = array2timetable(residual,'RowTimes',t); end IMF=SPT_ST(IMF); if nargout > 0 varargout{1} = IMF; end if nargout > 1 varargout{2} = residual; end % if nargout > 2 % varargout{3} = info; % end if nargout > 2 varargout{3} = ssc; end end %输入参数处理， function [x,t,td,isTT,opt] = parseAndValidateInputs(x, varargin) % 输入参数类型检查 isTT = isa(x,'timetable'); if ~isTT validateattributes(x, {'single','double','timetable'},{'vector'},'emd','X'); end % handle timetable and single if(isTT) signalwavelet.internal.util.utilValidateattributesTimetable(x, {'sorted','singlechannel'},'emd','X'); [x, t, td] = signalwavelet.internal.util.utilParseTimetable(x); else isSingle = isa(x,'single'); td = []; if(isSingle) t = single(1:length(x))'; else t = (1:length(x))'; end end % data integrity checking validateattributes(x, {'single','double'},{'nonnan','finite','real','nonsparse'},'emd','X'); validateattributes(t, {'single','double'},{'nonnan','finite','real'},'emd','T'); % turn x into column vector if isrow(x) x = x(:); end % parse and validate name-value pairs if(isempty(varargin)) opt = signalwavelet.internal.emd.emdOptions(); else opt = signalwavelet.internal.emd.emdOptions(varargin{:}); end validatestring(opt.Interpolation,{'spline', 'pchip'}, 'emd', 'Interpolation'); validateattributes(opt.SiftStopCriterion.SiftMaxIterations,... {'numeric'},{'nonnan','finite','scalar','>',0,'integer'}, 'emd', 'SiftMaxIterations'); validateattributes(opt.SiftStopCriterion.SiftRelativeTolerance,... {'numeric'},{'nonnan','finite','scalar','>=',0,'<',1},'emd', 'SiftRelativeTolerance'); validateattributes(opt.DecompStopCriterion.MaxEnergyRatio,... {'numeric'},{'nonnan','finite','scalar'}, 'emd', 'MaxEnergyRatio'); validateattributes(opt.DecompStopCriterion.MaxNumExtrema,... {'numeric'},{'nonnan','finite','scalar','>=',0,'integer'},'emd','MaxNumExtrema'); validateattributes(opt.DecompStopCriterion.MaxNumIMF,... {'numeric'},{'nonnan','finite','scalar','>',0,'integer'},'emd', 'MaxNumIMF'); end %%%%EMD分解函数 增加IMF_num：控制分解次数 function [IMFs, rsig, info] = localEMD(x, t, opt,IMF_num) isInMATLAB = coder.target('MATLAB'); isSingle = isa(x,'single'); % get name-value pairs Interpolation = opt.Interpolation; MaxEnergyRatio = opt.DecompStopCriterion.MaxEnergyRatio; MaxNumExtrema = opt.DecompStopCriterion.MaxNumExtrema; MaxNumIMF = opt.DecompStopCriterion.MaxNumIMF; SiftMaxIterations = opt.SiftStopCriterion.SiftMaxIterations; SiftRelativeTolerance = opt.SiftStopCriterion.SiftRelativeTolerance; Display = opt.Display; % initialization rsig = x; N = length(x); if(isSingle) ArrayType = 'single'; else ArrayType = 'double'; end IMFs = zeros(N, MaxNumIMF, ArrayType); info.NumIMF = zeros(MaxNumIMF, 1, ArrayType); info.NumExtrema = zeros(MaxNumIMF, 1, ArrayType); info.NumZerocrossing = zeros(MaxNumIMF, 1, ArrayType); info.NumSifting = zeros(MaxNumIMF, 1, ArrayType); info.MeanEnvelopeEnergy = zeros(MaxNumIMF, 1, ArrayType); info.RelativeTolerance = zeros(MaxNumIMF, 1, ArrayType); % preallocate memory rsigPrev = zeros(N, 1, ArrayType); mVal = zeros(N, 1, ArrayType); upperEnvelope = zeros(N, 1, ArrayType); lowerEnvelope = zeros(N, 1, ArrayType); % % Define intermediate print formats % if(isInMATLAB && opt.Display) % fprintf('Current IMF | #Sift Iter | Relative Tol | Stop Criterion Hit \n'); % formatstr = ' %5.0f | %5.0f | %12.5g | %s\n'; % end % use different functions under different environment if(isInMATLAB) if(~isSingle) localFindExtramaIdx = @(x) signalwavelet.internal.emd.cg_utilFindExtremaIdxmex_double(x); else localFindExtramaIdx = @(x) signalwavelet.internal.emd.cg_utilFindExtremaIdxmex_single(x); end else localFindExtramaIdx = @(x) signalwavelet.internal.emd.utilFindExtremaIdx(x); end % extract IMFs i = 0; outerLoopExitFlag = 0; if IMF_num~=0 MaxNumIMF=IMF_num; end while(i<MaxNumIMF) % convergence checking [peaksIdx, bottomsIdx] = localFindExtramaIdx(rsig); numResidExtrema = length(peaksIdx) + length(bottomsIdx); energyRatio = 10*log10(norm(x,2)/norm(rsig,2)); if energyRatio > MaxEnergyRatio outerLoopExitFlag = 1; break end if numResidExtrema < MaxNumExtrema outerLoopExitFlag = 2; break end % SIFTING process initialization rsigL = rsig; rtol = ones(1, ArrayType); k = 0; SiftStopCriterionHit = 'SiftMaxIteration'; % Sifting process while (k<SiftMaxIterations) % check convergence if(rtol<SiftRelativeTolerance) SiftStopCriterionHit = 'SiftMaxRelativeTolerance'; break; end % store previous residual rsigPrev(1:N) = rsigL; % finding peaks [peaksIdx, bottomsIdx] = localFindExtramaIdx(rsigL); if((length(peaksIdx) + length(bottomsIdx))>0) % compute upper and lower envelope using extremas [uLoc, uVal, bLoc, bVal] = computeSupport(t, rsigL, peaksIdx, bottomsIdx); upperEnvelope(:) = interp1(uLoc, uVal, t, Interpolation); lowerEnvelope(:) = interp1(bLoc, bVal, t, Interpolation); % subtract mean envelope from residual mVal(1:N) = (upperEnvelope + lowerEnvelope)/2; else mVal(1:N) = 0; end rsigL = rsigL - mVal; % residual tolerance rtol = (norm(rsigPrev-rsigL,2)/norm(rsigPrev,2))^2; k = k + 1; end % if(isInMATLAB && Display) % fprintf(formatstr, i+1, k, rtol, SiftStopCriterionHit); % end % record information [peaksIdx, bottomsIdx] = localFindExtramaIdx(rsigL); numZerocrossing = sum(diff(sign(rsigL))~=0); info.NumIMF(i+1) = i+1; info.NumExtrema(i+1) = length(peaksIdx) + length(bottomsIdx); info.NumZerocrossing(i+1) = numZerocrossing; info.MeanEnvelopeEnergy(i+1) = mean(mVal.^2); info.NumSifting(i+1) = k; info.RelativeTolerance(i+1) = rtol; % extract new IMF and subtract the IMF from residual signal IMFs(:,i+1) = rsigL; rsig = rsig - IMFs(:,i+1); i = i + 1; end % if(isInMATLAB && Display) % switch outerLoopExitFlag % case 0 % disp(getString(message('shared_signalwavelet:emd:general:MaxNumIMFHit', 'MaxNumIMF'))); % case 1 % disp(getString(message('shared_signalwavelet:emd:general:MaxEnergyRatioHit', 'MaxEnergyRatio'))); % case 2 % disp(getString(message('shared_signalwavelet:emd:general:MaxNumExtremaHit', 'MaxNumExtrema'))); % end % end % remove extra portion IMFs = IMFs(:,1:i); info.NumIMF = info.NumIMF(1:i); info.NumExtrema = info.NumExtrema(1:i); info.NumZerocrossing = info.NumZerocrossing(1:i); info.NumSifting = info.NumSifting(1:i); info.MeanEnvelopeEnergy = info.MeanEnvelopeEnergy(1:i); info.RelativeTolerance = info