EMD信号处理的程序
附件不能发压缩包，好奇怪的，，只好复制粘贴过来了， 是以前我用过的EMD信号处理程序 ，看对大家能不能有用 。% EMD.M
%
% computes EMD (Empirical Mode Decomposition) according to:
%
% 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
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gon峚lv弒
% "On Empirical Mode Decomposition and its algorithms",
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% default stopping criterion for sifting :
% at each point : mean amplitude < threshold2*envelope amplitude
% &
% mean of boolean array ((mean amplitude)/(envelope amplitude) > threshold) < tolerance
% &
% |#zeros-#extrema|<=1
%
% inputs:
% - x: analysed signal (line vector)
% - opts (optional): struct object with (optional) fields:
% - t: sampling times (line vector) (default: 1:length(x))
% - stop: threshold, threshold2 and tolerance (optional)
% for sifting stopping criterion
% default:
% - display: if equals to 1 shows sifting steps with pause
% if equals to 2 shows sifting steps without pause (movie style)
% - maxiterations: maximum number of sifting steps for the computation of each mode
% - fix (int): disable the stopping criterion and do exactly
% the value of the field number of sifting steps for each mode
% - maxmodes: maximum number of imfs extracted
% - interp: interpolation scheme: 'linear', 'cubic' or 'spline' (default)
% - fix_h (int): do <fix_h> sifting iterations with |#zeros-#extrema|<=1 to stop
% according to 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
% - mask: masking signal used to improve the decomposition
% according to R. Deering and J. F. Kaiser, "The use of a masking signal to
% improve empirical mode decomposition",
% ICASSP 2005
%
% outputs:
% - imf: intrinsic mode functions (last line = residual)
% - ort: index of orthogonality
% - nbits: number of iterations for each mode
%
% calls:
% - io: computes the index of orthogonality
%
%examples:
%
%>>x = rand(1,512);
%
%>>imf = emd(x);
%
%>>imf = emd(x,struct('stop',,'maxiterations',100));
%Remark: the following syntax is equivalent
%>>imf = emd(x,'stop',,'maxiterations',100);
%
%>>options.dislpay = 1;
%>>options.fix = 10;
%>>options.maxmodes = 3;
%>> = emd(x,options);
function = emd(varargin);
= init(varargin{:});
if display_sifting
figure
end
% maximum number of iterations
% MAXITERATIONS=2000;
%main loop : requires at least 3 extrema to proceed
while ~stop_EMD(r) & (k < MAXMODES+1 | MAXMODES == 0) & ~any(mask)
% current mode
m = r;
% mode at previous iteration
mp = m;
if FIXE
= stop_sifting_fixe(t,m,INTERP);
elseif FIXE_H
stop_count = 0;
= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);
stop_count = 0;
else
= stop_sifting(m,t,sd,sd2,tol,INTERP);
end
if (max(m) - min(m)) < (1e-10)*(max(x) - min(x))
if ~stop_sift
warning('forced stop of EMD : too small amplitude')
else
disp('forced stop of EMD : too small amplitude')
end
break
end
% sifting loop
while ~stop_sift & nbit<MAXITERATIONS
if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0 & ~FIXE & nbit > 100)
disp(['mode ',int2str(k),', iteration ',int2str(nbit)])
if exist('s')
disp(['stop parameter mean value : ',num2str(s)])
end
= extr(m);
disp()
end
%sifting
m = m - moyenne;
%computation of mean and stopping criterion
if FIXE
= stop_sifting_fixe(t,m,INTERP);
elseif FIXE_H
= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);
else
= stop_sifting(m,t,sd,sd2,tol,INTERP);
end
% display
if display_sifting
= envelope(t,mp,INTERP);
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')
warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
else
warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'.'])
end
end
end % sifting loop
imf(k,:) = m;
if display_sifting
disp(['mode ',int2str(k),' stored'])
end
nbits(k) = nbit;
k = k+1;
r = r - m;
nbit=0;
end %main loop
if sum(r.^2) & ~any(mask)
imf(k,:) = r;
end
ort = io(x,imf);
if display_sifting
close
end
%---------------------------------------------------------------------------------------------------
function stop = stop_EMD(r)
= extr(r);
ner = length(indmin) + length(indmax);
stop = ner <3;
%---------------------------------------------------------------------------------------------------
function = stop_sifting(m,t,sd,sd2,tol,INTERP)
try
= envelope(t,m,INTERP);
nem = length(indmin) + length(indmax);
nzm = length(indzer);
% evaluation of mean zero
sx=2*(abs(envmoy))./(abs(envmax-envmin));
s = mean(sx);
stop = ~((mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (nem > 2));
catch
stop = 1;
envmoy = zeros(1,length(m));
s = NaN;
end
%---------------------------------------------------------------------------------------------------
function = stop_sifting_fixe(t,m,INTERP)
try
= envelope(t,m,INTERP);
stop = 0;
catch
moyenne = zeros(1,length(m));
stop = 1;
end
%---------------------------------------------------------------------------------------------------
function = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H)
try
= envelope(t,m,INTERP);
nem = length(indmin) + length(indmax);
nzm = length(indzer);
if (abs(nzm-nem)>1)
stop = 0;
stop_count = 0;
else
stop_count = stop_count+1;
stop = (stop_count == FIXE_H);
end
catch
moyenne = zeros(1,length(m));
stop = 1;
end
%---------------------------------------------------------------------------------------------------
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',[])
holdoff
subplot(4,1,2)
plot(t,sx)
hold on
plot(t,sdt,'--r')
plot(t,sd2t,':k')
title('stop parameter')
set(gca,'XTick',[])
hold off
subplot(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')
end
if display_sifting == 2
pause(0.01)
else
pause
end
%---------------------------------------------------------------------------------------------------
function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)
subplot(3,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',[])
holdoff
subplot(3,1,2)
plot(t,m)
title(['IMF ',int2str(k),'; iteration ',int2str(nbit),' after sifting']);
set(gca,'XTick',[])
subplot(3,1,3);
plot(t,r-m)
title('residue');
if display_sifting == 2
pause(0.01)
else
pause
end
%---------------------------------------------------------------------------------------------------
function = envelope(t,x,INTERP)
%computes envelopes and mean with various interpolations
NBSYM = 2;
DEF_INTERP = 'spline';
if nargin < 2
x = t;
t = 1:length(x);
INTERP = DEF_INTERP;
end
if nargin == 2
if ischar(x)
INTERP = x;
x = t;
t = 1:length(x);
end
end
if ~ischar(INTERP)
error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end
if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end
if min() > 1
error('x and t must be vectors')
end
s = size(x);
if s(1) > 1
x = x';
end
s = size(t);
if s(1) > 1
t = t';
end
if length(t) ~= length(x)
error('x and t must have the same length')
end
lx = length(x);
= extr(x,t);
%boundary conditions for interpolation
= boundary_conditions(indmin,indmax,t,x,NBSYM);
% definition of envelopes from interpolation
envmax = interp1(tmax,xmax,t,INTERP);
envmin = interp1(tmin,xmin,t,INTERP);
if nargout > 2
envmoy = (envmax + envmin)/2;
end
%---------------------------------------------------------------------------------------
function = boundary_conditions(indmin,indmax,t,x,nbsym)
% computes the boundary conditions for interpolation (mainly mirror symmetry)
lx = length(x);
if (length(indmin) + length(indmax) < 3)
error('not enough extrema')
end
if indmax(1) < indmin(1)
if x(1) > x(indmin(1))
lmax = fliplr(indmax(2:min(end,nbsym+1)));
lmin = fliplr(indmin(1:min(end,nbsym)));
lsym = indmax(1);
else
lmax = fliplr(indmax(1:min(end,nbsym)));
lmin = ;
lsym = 1;
end
else
if x(1) < x(indmax(1))
lmax = fliplr(indmax(1:min(end,nbsym)));
lmin = fliplr(indmin(2:min(end,nbsym+1)));
lsym = indmin(1);
else
lmax = ;
lmin = fliplr(indmin(1:min(end,nbsym)));
lsym = 1;
end
end
if indmax(end) < indmin(end)
if x(end) < x(indmax(end))
rmax = fliplr(indmax(max(end-nbsym+1,1):end));
rmin = fliplr(indmin(max(end-nbsym,1):end-1));
rsym = indmin(end);
else
rmax = ;
rmin = fliplr(indmin(max(end-nbsym+1,1):end));
rsym = lx;
end
else
if x(end) > x(indmin(end))
rmax = fliplr(indmax(max(end-nbsym,1):end-1));
rmin = fliplr(indmin(max(end-nbsym+1,1):end));
rsym = indmax(end);
else
rmax = fliplr(indmax(max(end-nbsym+1,1):end));
rmin = ;
rsym = lx;
end
end
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
% in case symmetrized parts do not extend enough
if tlmin(1) > t(1) | tlmax(1) > t(1)
if lsym == indmax(1)
lmax = fliplr(indmax(1:min(end,nbsym)));
else
lmin = fliplr(indmin(1:min(end,nbsym)));
end
if lsym == 1
error('bug')
end
lsym = 1;
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
end
if trmin(end) < t(lx) | trmax(end) < t(lx)
if rsym == indmax(end)
rmax = fliplr(indmax(max(end-nbsym+1,1):end));
else
rmin = fliplr(indmin(max(end-nbsym+1,1):end));
end
if rsym == lx
error('bug')
end
rsym = lx;
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
end
xlmax =x(lmax);
xlmin =x(lmin);
xrmax =x(rmax);
xrmin =x(rmin);
tmin = ;
tmax = ;
xmin = ;
xmax = ;
%---------------------------------------------------------------------------------------------------
function = extr(x,t);
%extracts the indices corresponding to extrema
if(nargin==1)
t=1:length(x);
end
m = length(x);
if nargout > 2
x1=x(1:m-1);
x2=x(2:m);
indzer = find(x1.*x2<0);
if any(x == 0)
iz = find( x==0 );
indz = [];
if any(diff(iz)==1)
zer = x == 0;
dz = diff();
debz = find(dz == 1);
finz = find(dz == -1)-1;
indz = round((debz+finz)/2);
else
indz = iz;
end
indzer = sort();
end
end
d = diff(x);
n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;
indmax = find(d1.*d2<0 & d1>0)+1;
% when two or more consecutive points have the same value we consider only one extremum in the middle of the constant area
if any(d==0)
imax = [];
imin = [];
bad = (d==0);
dd = diff();
debs = find(dd == 1);
fins = find(dd == -1);
if debs(1) == 1
if length(debs) > 1
debs = debs(2:end);
fins = fins(2:end);
else
debs = [];
fins = [];
end
end
if length(debs) > 0
if fins(end) == m
if length(debs) > 1
debs = debs(1:(end-1));
fins = fins(1:(end-1));
else
debs = [];
fins = [];
end
end
end
lc = length(debs);
if lc > 0
for k = 1:lc
if d(debs(k)-1) > 0
if d(fins(k)) < 0
imax = ;
end
else
if d(fins(k)) > 0
imin = ;
end
end
end
end
if length(imax) > 0
indmax = sort();
end
if length(imin) > 0
indmin = sort();
end
end
%---------------------------------------------------------------------------------------------------
function = init(varargin)
x = varargin{1};
if nargin == 2
if strcmp(class(varargin{2}),'struct')
inopts = varargin{2};
else
error('when using 2 arguments the first one is the analysed signal x and the second one is a struct object describing the options')
end
elseif nargin > 2
try
inopts = struct(varargin{2:end});
catch
error('bad argument syntax')
end
end
% default for stopping
defstop = ;
opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask'};
defopts.stop = defstop;
defopts.display = 0;
defopts.t = 1:max(size(x));
defopts.maxiterations = 2000;
defopts.fix = 0;
defopts.maxmodes = 0;
defopts.interp = 'spline';
defopts.fix_h = 0;
defopts.mask = 0;
opts = defopts;
if(nargin==1)
inopts = defopts;
elseif nargin == 0
error('not enough arguments')
end
names = fieldnames(inopts);
for nom = names'
if length(strmatch(char(nom), opt_fields)) == 0
error(['bad option field name: ',char(nom)])
end
eval(['opts.',char(nom),' = inopts.',char(nom),';'])
end
t = opts.t;
stop = opts.stop;
display_sifting = opts.display;
MAXITERATIONS = opts.maxiterations;
FIXE = opts.fix;
MAXMODES = opts.maxmodes;
INTERP = opts.interp;
FIXE_H = opts.fix_h;
mask = opts.mask;
S = size(x);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('x must have only one row or one column')
end
if S(1) > 1
x = x';
end
S = size(t);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
error('option field t must have only one row or one column')
end
if S(1) > 1
t = t';
end
if (length(t)~=length(x))
error('x and option field t must have the same length')
end
S = size(stop);
if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 3) | (S(2) > 3) | (length(S) > 2)
error('option field stop must have only one row or one column of max three elements')
end
if ~all(isfinite(x))
error('data elements must be finite')
end
if S(1) > 1
stop = stop';
S = size(stop);
end
if S(2) < 3
stop(3)=defstop(3);
end
if S(2) < 2
stop(2)=defstop(2);
end
if ~ischar(INTERP)
error('interp field must be ''linear'', ''cubic'' or ''spline''')
end
if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
error('interp field must be ''linear'', ''cubic'' or ''spline''')
end
%special procedure when a masking signal is specified
if any(mask)
S = size(mask);
if min(S) > 1
error('masking signal must have the same dimension as the analyzed signal x')
end
if S(1) > 1
mask = mask';
end
if max(S) ~= max(size(x))
error('masking signal must have the same dimension as the analyzed signal x')
end
opts.mask = 0;
imf1 = emd(x+mask,opts);
imf2 = emd(x-mask,opts);
if size(imf1,1) ~= size(imf2,1)
warning(['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])
end
S1 = size(imf1,1);
S2 = size(imf2,1);
if S1 ~= S2
if S1 < S2
tmp = imf1;
imf1 = imf2;
imf2 = tmp;
end
imf2(max(S1,S2),1) = 0;
end
imf = (imf1+imf2)/2;
end
sd = stop(1);
sd2 = stop(2);
tol = stop(3);
lx = length(x);
sdt = sd*ones(1,lx);
sd2t = sd2*ones(1,lx);
if FIXE
MAXITERATIONS = FIXE;
if FIXE_H
error('cannot use both ''fix'' and ''fix_h'' modes')
end
end
% number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;
r = x;
if ~any(mask) % if a masking signal is specified "imf" already exists at this stage
imf = [];
end
k = 1;
% iterations counter for extraction of 1 mode
nbit=0;
% total iterations counter
NbIt=0;
感谢您的分享您是大学教师吧代码放在这里更规范
页:
[1]