等模分解matlab代码,带端点延拓的经验模态分解matlab程序和例子代码

function [imf,ort,nbits] = emd1(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{:});

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

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')%查找是否存在变量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;