emd程序，EMD分解算法

emd程序，可用于信号处理，信号去躁，emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁emd程序，可用于信号处理，信号去躁

% EMD.M 
%
% G. Rilling July 2002
%
% 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
%
% 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）
%          - t （optional） : sampling times （line vector） （default : 1:length（x））
%          - stop （optional） : threshold threshold2 and tolerance （optional）
%                              for sifting stopping criterion 
%                              default : [0.050.50.05]
%          - tst （optional） : if equals to 1 shows sifting steps with pause
%                             if equals to 2 no pause
%
% outputs: - imf : intrinsic mode functions （last line = residual）
%          - ort : index of orthogonality
%          - nbits : number of iterations for each mode
%
% calls:   - extr （finds extrema and zero-crossings）
%          - io : computes the index of orthogonality

function [imfortnbits] = emd（xtstoptst）

% default for stopping
defstop = [0.050.50.05];

if（nargin==1）
  t = 1:length（x）;
  stop = defstop;
  tst = 0;
end

if（nargin==2）
  stop = defstop;
  tst = 0;
end

if （nargin==3）
  tst=0;
end

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（‘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 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（‘stop must have only one row or one column of max three elements‘）
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

sd = stop（1）;
sd2 = stop（2）;
tol = stop（3）;

if tst
  figure
end

% maximum number of iterations
MAXITERATIONS=2000;

% maximum number of symmetrized points for interpolations
NBSYM = 2;

lx = length（x）;

sdt（lx） = 0;
sdt = sdt+sd;
sd2t（lx） = 0;
sd2t = sd2t+sd2;

% maximum number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;

r = x;
imf = [];
k = 1;

% iterations counter for extraction of 1 mode
nbit=0;

% total iterations counter
NbIt=0;

while ner > 2
        
  % current mode
  m = r;
  
  % mode at previous iteration
  mp = m;
  
  sx = sd+1;
  
  % tests if enough extrema to proc