% 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
% gabriel.rilling@ens-lyon.fr
%
% 翻译：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
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;
stop = (stop_count == FIXE_H);
end
catch
moyenne = zeros(1,length(m));
stop = 1;
end
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',[])
hold off
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
end
%---------------------------------------------------------------------------------------------------
% 演示分解过程（FIX和FIX_H停止准则）
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',[])
hold off
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
end
%---------------------------------------------------------------------------------------
% 处理边界条件（镜像法）
function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)
% 实数情况下，x = z
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 = [fliplr(indmin(1:min(end,nbsym-1))),1];
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 = [fliplr(indmax(1:min(end,nbsym-1))),1];
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 = [lx,fliplr(indmax(max(end-nbsym+2,1):end))];
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 = [lx,fliplr(indmin(max(end-nbsym+2,1):end))];
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);
% 如果对称的部分没有足够的极值点
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
% 延拓点上的取值
zlmax = z(lmax);
zlmin = z(lmin);
zrmax = z(rmax);
zrmin = z(rmin);
% 完成延拓
tmin = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
zmin = [zlmin z(indmin) zrmin];
zmax = [zlmax z(indmax) zrmax];
end
%---------------------------------------------------------------------------------------------------
% 极值点和过零点位置提取
function [indmin, indmax, indzer] = extr(x,t)
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) % 考虑信号采样点恰好为0的位置
iz = find( x==0 ); % 信号采样点恰好为0的位置
indz = [];
if any(diff(iz)==1) % 出现连0的情况
zer = x == 0; % x=0处为1，其它地方为0
dz = diff([0 zer 0]); % 寻找0与非0的过渡点
debz = find(dz == 1); % 0值起点
finz = find(dz == -1)-1; % 0值终点
indz = round((debz+finz)/2); % 选择中间点作为过零点
else
indz = iz; % 若没有连0的情况，该点本身就是过零点
end
indzer = sort([indzer indz]); % 全体过零点排序
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; % 最大值
% 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连0类似
if any(d==0)
imax = [];
imin = [];
bad = (d==0);
dd = diff([0 bad 0]);
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 = [imax round((fins(k)+debs(k))/2)];
end
else
if d(fins(k)) > 0
imin = [imin round((fins(k)+debs(k))/2)];
end
end
end
end
if length(imax) > 0
indmax = sort([indmax imax]);
end
if length(imin) > 0
indmin = sort([indmin imin]);
end
end
end
%---------------------------------------------------------------------------------------------------
function ort = io(x,imf)
% ort = IO(x,imf) 计算正交指数
%
% 输入 : - x : 分析信号
% - imf : IMF信号
n = size(imf,1);
s = 0;
% 根据公式计算
for i = 1:n
for j = 1:n
if i ~= j
s = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));
end
end
end
ort = 0.5*s;
end
%---------------------------------------------------------------------------------------------------
% 函数参数解析
function [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=[-178.45,-86.45,-69.45,-105.45,-90.45,-117.45,-102.45,-86.45,-60.45,-31.45,-41.45,-56.45,-57.45,-24.45,-44.45,-39.45,10.55,-2.45,5.55,-17.45,27.55,-60.45,-21.45,39.55,37.55,29.55,12.55,111.55,84.55,29.55,85.55,80.55,111.55,13.55,94.55,117.55,59.55,46.55,114.55,181.55];
%x=[3920.5,231.5,2394.5,-3303.5,327.5,1497.5,824.5,-3470.5,-1777.5,-2780.5,3451.5,253.5,1143.5,3139.5,1436.5,757.5,-1711.5,-3628.5,-1015.5,838.5,2588.5,-2705.5,-506.5,2419.5,-289.5,-699.5,685.5,-562.5,2143.5,-579.5,4674.5,-71.5,-4004.5,-1778.5,-1394.5,-769.5,-504.5,886.5,352.5,-2413.5];
%x=[4815,5570,6576,9609,9314,11152,4868,6914,14046,5264,5995,3669,7755,8172,2614,5267,5924,3865,9129,5970,5110,3740,6981,4745,3912,6828,7789,6643,7182,3807,3932,5444,4899,4888,7210,6653,6839,6733,4422,6973,7479,5415,5067,8132,5725,7009,4309,7317,2669,3711,3389,3704,4449,4835,4107,3180,4839,6263,4806,5225];
%x=[1336,1312,1323,1326,1413,1256,1305,1427,1447,1443,1532,1482,1474,1385,1407,1324,1332,1393,1260,1309,1345,1367,1379,1332,1466,1301,1387,1393,1330,1325,1480,1536,1558,1428,1384,1455,1470,1524,1583,1523,1501,1541,1554,1643,1600,1529,1566,1573,1576,1540,1546,1578,1550,1624,1583,1610,1681,1605,1590,1507];
%50353
%x=[-11.33,-27.75,-16.83,-27.08,-1.42,-22.5,-23.08,-11.92,-19.25,-17.5,-15.17,-8.08,-7.5,-17.5,-14.83,-6.08,-24.5,1.17,1.58,8.58,-8.83,-5.5,-1.5,-5.75,4.92,-6.17,-0.58,-2.17,-6.67,-9.92,-2,7.25,-1.92,2.42,6.33,-10.33,8.5,12.5,-7.17,-0.67];
%x=[291.75,495.58,347.75,362.25,342.42,319.5,380.25,350.17,344.25,347,373.42,510.5,479.5,430.42,474.92,419.67,395.5,394.92,421.17,477.92,497.92,345.08,449.42,327.08,396.92,436.58,334.83,331,315.67,387.33,307.92,328.58,475.25,400.58,245.5,343.92,254,373.33,486.25,420.17];
%50434
%x=[-46.17,-54.75,-48.25,-54.33,-37.92,-55,-50.92,-50.83,-48.83,-53.08,-51.92, -40.5,-41.08,-51.5,-48.25,-43.33,-56.08,-38.58,-34.33,-27, -43.58,-37.58,-34.83,-34.58,-31.42,-46.42,-39.42,-32.83,-42.83,-44.5,-45.83,-33.83,-41.08,-37.83,-40.08,-44.75,-28.83,-32.83,-41.5,-45.17];
%x=[339.08,442.25,281.5,321.67,367.25,411.42,377.58,418.75,317.42,410.42,299.33,464.42,446.58,519.58,443.92,281.42,312.42,543.08,358.17,452.5,442.5,319.92,446.33,335.75,330.75,479.75,429.5,370.08,354,281.17,305,242.75,376.42,288.58,322,342.67,236.92,313.75,394.58,278.33];
%50527
%x=[-19.08,-22.42,-16.33,-19.67,1.25,-18.33,-21.83,-17.5,-15,-21.75,-20.5,-5.92,-12.83,-24.5,-18.33,-6.25,-16.17,-8.58,-2.08,5,-7.25,1.17,-1,4.17,8.67,-3.92,-1.42,6.83,-4.92,-8.92,0.83,4.75,0.75,4.5,-3.17,-4.92,12.58,4.75,-10.25,-12.33];
%x=[228.17 260.5 303.25 268.42 260.75 338.5 400 321.75 285.58 285.58 240.67 368.92 275.75 446.25 296.58 103.75 214.83 363.42 367 451.42 281.58 292.83 370.5 272 221.67 355.5 307.58 415.5 230.33 350.5 205 293.25 265.33 295.08 259.67 207.5 193.75 302.17 307.08 242.42];
%50557
%x=[5.25 -4.08 -2 -7.58 16.08 -1.92 -9.5 -1.5 -0.58 -3.58 1.25 9.58 -0.42 -2.58 -0.67 0.92 -11.25 9.92 12.67 17.67 4.75 5.67 8.58 8.33 19.67 6.42 12.58 13.08 7.75 1.42 8.25 16.83 16.08 15.67 12.25 6.67 25.83 21.92 3.17 6.42];
%x=[396.92 427.92 292.5 383.67 374.42 338.83 334.75 408.08 329.42 445.5 589.58 335.83 425.42 502.33 326.75 361.67 449.17 482.75 467.17 455.75 452.5 387.17 499.5 488.17 426.67 357.42 346.67 510.08 367.42 334.67 257.33 325.33 444.42 363 294.67 320.42 242.67 398.17 387.17 454.83];
%50564x=[-11.08 -18.08 -12.75 -21.75 1.83 -18.92 -21.42 -12.83 -17.75 -16.67 -12.42 -2.17 -6.5 -11.42 -9.58 -4.08 -17.17 -2.08 7.33 11.42 -4 0.75 2.92 -0.17 8.25 -1.75 8 5.33 1 -1.58 2.92 10.75 10.83 11.92 8.17 2.17 20 17.25 -4.42 2.83];
%x=[324.92 585.17 360.58 391.83 443.5 343.25 327.67 389.25 410.83 399.33 473.17 502.92 575.5 594.25 377 419.92 485.33 460 311.17 604.17 506.83 380.58 572.92 479.08 436.08 423 347.67 520.58 616 381.92 285.25 441.58 639.58 383.5 342.67 392.5 280.83 338.33 553.08 451.67];
%50632
%x=[-7.75 -13 -3.5 -13.42 4.83 -14.25 -13.17 -6.83 -8 -15.08 -8 0.67 -0.92 -14.92 -12.92 -3.33 -14.83 -3.08 4.5 8.67 -3.5 3 1.5 4.33 6.92 -7 1.92 4.42 -4.92 -5.08 -3.17 3.58 5.42 4.42 -4.42 -5.67 15.58 8.08 -4 -10.08];
%x=[346 406.92 402.25 314.83 381.25 384.42 424.92 349.5 322.17 530.25 432.58 438.75 414.67 426.17 479.67 416.83 356.75 452.92 500.75 470.17 485.92 277.92 533.5 367.25 283.92 388.58 469 625.33 274.42 276.67 323.42 373.67 412.33 275.58 340.5 360.75 186.83 382.42 440.42 344.17];
%50658
%x=[17.33 8.67 13.92 10.67 31.92 12.83 7.58 13.83 15.5 10.17 17.42 27 15 10.33 10.75 16 4.08 21 26.58 30.75 21.58 20.58 22.25 21.83 30.58 23.58 28.58 30 20.58 20.58 21.5 27.17 33.25 30.58 22.92 19.58 39.42 33.08 15.67 17.92];
%x=[328.25 571.5 435.67 433.75 360.33 210.42 439.75 345.5 374.58 440.92 446.5 331.17 440.17 564.5 462.92 428.17 482.67 424.08 335.83 404.83 499.5 371.33 491.83 463.25 254.5 485.25 497 750.83 359.42 311.75 300.67 402.58 656.17 338.08 409.92 433.33 275.92 299.33 593.42 356.83];
%50727
%x=[-31.67 -32.58 -28 -36.5 -18.25 -36.92 -33 -31.92 -26.5 -37.17 -35.75 -20.67 -26.58 -39.25 -36.5 -31.5 -33.75 -27.17 -18.33 -14.42 -26.83 -19.25 -22.83 -17.58 -19.17 -26 -19.67 -11 -23.58 -27.75 -28.5 -19.83 -16.83 -15.75 -26 -20.92 -4.17 -11.83 -22.25 -28];
%x=[376.25 321.42 328.08 324.58 388.42 375.42 348.42 294.83 340.92 354.75 409.42 437.33 432.33 390.83 436.08 348.42 304.75 489.83 368.92 515.25 380.17 364.83 370.08 506.67 332.42 401 326.33 534.33 226.5 293.58 283.25 339 365.5 236.33 376.83 317.42 259.92 340.5 349.67 297.58];
%50745
%x=[36.08 28.42 34 24.67 47.08 27.5 29 35.92 35.42 30.58 37.42 47.17 36.92 29.67 31.83 37.75 29.58 43.08 48.42 48.75 42.58 41.25 45.17 46.33 52.83 43.58 51.25 54.17 43.75 43.08 44.08 45.17 51.58 49.5 40.5 40.08 58.83 55.83 38.08 34.42];
%x=[293.42 283.42 271.67 316.83 236.92 274.17 303.58 275.17 248.67 358.33 481.83 293.33 432.92 384.83 371.83 413.08 408.5 543.08 327.08 370.25 440.75 295.33 391.08 384.08 286.58 301.17 270.17 490 345.25 288.33 208.58 309.75 500.17 245.5 409.58 387.58 332.08 313.17 420.5 426.67];
%50756
%x=[19.92 10.75 13.25 10.17 30.83 13.33 9.08 14.5 18.83 12.33 17 29.33 18.08 11.92 10.67 19.58 7.5 20.58 28.25 33.58 23.58 22 23.5 24.42 31.08 26.25 30.33 33.92 21.58 22.42 23.83 30.08 34.83 32.5 23.83 22.75 42.5 36 20.75 20];
%x=[398.33 609.17 466.17 505.33 492.75 312 416.92 402.25 398.17 289.42 414.83 397.08 504.58 523 508.67 375.08 558.83 506.83 381.92 457.92 636.08 482.5 469.08 396.67 381.33 480 574.75 480.67 373.5 422 249.75 531.83 540 297.75 502.17 511.75 374.67 501.67 494.33 412.58];
%50788
%x=[26.25 21.5 25.67 24.58 40.08 24.42 19.92 26.92 29.08 24.5 27.33 38.5 27.17 21.75 28.5 32.5 24 37.42 42.83 45.42 37 37.5 36.67 36.92 40.75 33 37 43.67 29.5 28.17 29.58 31.25 36.25 34.17 31.08 27.67 42.92 40.92 24.33 26.33];
%x=[561.75 515.75 487.75 443.67 307.33 341.83 282.08 343.75 309.08 420.58 553.83 413 541.17 565.67 477.58 322.42 589.17 395.83 333 483.83 490.58 374.67 466.08 626.92 350.5 436.5 432.17 455.08 257.08 365.83 301.25 451 287.92 371.58 327.92 507.17 390.5 334.75 392.33 512.83];
%50854
%x=[39.75 26.75 28 22.67 47.58 27.08 26.25 34 34.42 25.92 33.33 46 35.75 30.25 28.42 34.83 28.17 38.25 45 49.5 41.17 40.92 40.92 42.58 51.5 44.17 50.83 49.58 41.17 39.25 40.25 45.75 51.25 51.25 38.17 40.75 57.17 52.25 37.08 32.58];
%x=[365.67 370.58 350.33 351.67 345.17 432.67 338 306.92 291.67 304 384.5 281.08 439.42 350.17 440 333.67 411.25 500.42 289.92 334.92 409.92 264 446.75 320.67 246.75 420.83 384.33 456.67 304.33 223.92 207.08 379.83 433.58 267.83 471.83 350.83 325.42 327.08 429.83 376.25];
%50915
%x=[9 5.58 14.75 1.67 27.08 6.75 3.83 5.25 13.33 10.5 3.67 14.92 17.58 6.67 3.83 9.42 5.5 16.75 25 22.83 7.92 14.33 11.92 19.42 24.92 12.17 24.92 30.92 27.67 19.08 21.42 31.33 22.83 29.92 22.58 25.33 41.5 32.08 22.42 15.75];
%x=[204 148.83 293.08 196.25 132.83 206.25 205.42 222.25 182.5 199.58 264 163.33 163.75 205.58 181.58 159.75 253.5 183.92 208.67 347.75 269.25 310.17 291.75 224.58 148.5 228.33 158.67 381.08 150.33 181.92 166.83 243.25 203.83 161.75 139.58 117 124.75 252.92 213.25 174.08];
%50949
%x=[49.67 46.08 48.17 41.58 62.42 42.75 45.42 49.67 51.5 44.67 48.42 60.67 54.67 47.75 44.17 50.75 45.42 55.92 61.83 64.5 56.25 57.25 58.83 61.33 66.08 60.42 68.75 70 60.08 48.42 51.92 60.5 64.83 65.25 50.83 58.42 72.5 68.42 54.5 46.75];
%x=[322.42 385.33 352.08 412 383.83 273.5 385 310.75 345.75 354.75 434.42 202.67 374.75 322.25 473.67 367.33 449.25 384.42 262.5 493.83 329.17 296.67 346.58 391.75 269.5 255.75 335.58 441.17 298.92 302.67 227.42 442.67 364.42 246.92 415 282 228 518.17 328.5 306.83];
%50953
%x=[38.33 31.75 35.25 29.17 52.75 32.92 32.67 37.33 41.08 30.25 37.58 52 41.75 34.17 33.33 39.42 33.67 41.75 50.58 54.33 44.33 42.83 43.83 47 52.83 50.33 56.58 60.5 48.83 46.17 48.5 54.83 60 58.58 47.42 52.67 66.42 66 49.5 45.08];
%x=[504.08 413.17 335.92 349.5 294.92 302.58 454.58 352.58 348.5 500.83 518.75 443.75 445.42 515.92 621.58 426.92 570 507.58 287.92 408.83 497.75 385 460.58 682.33 353.92 400.58 401.58 549.33 365.5 406.67 321 500.58 428.08 438.17 423.25 406.58 370.08 365.83 445.08 492.75];
%50963
%x=[24.92 18.08 19.25 15.17 35.42 17.17 16.08 20.5 25.25 14.75 19.58 34.67 24.75 18.5 18.58 23.58 17.75 27.08 35 39.17 28.25 27.08 28.08 30.83 36.25 28.5 33.25 39.5 26.17 26.17 26.83 31.17 36.5 34.83 28.33 29.33 44.58 41.25 24.67 21.67];
%x=[601.58 480.58 512.67 520.08 384.25 384.67 553.17 547.42 323.58 521 562.83 420.33 443.83 493.83 585.25 383.75 592.17 582 535.83 457.33 505.08 511.92 552.58 727 521.67 459.25 552.75 516.5 280.42 565 284.17 523.5 431.67 423.25 409.5 483.75 265.67 378.17 490.58 451.17];
%50968
%x=[27.42 20.42 23.58 19.75 37.58 17.25 21 24.92 27.75 19.92 24.92 37.5 29.67 21.92 21.58 24.25 20.58 29.5 36.42 40.92 30.92 31.25 31.25 35.17 40.25 33.58 41.92 47.58 33.08 29.33 30.67 37.25 43.17 41.83 33.67 38.92 51.33 45.83 31.92 27.92];
%x=[671.58 451.5 536.83 517.33 505.75 374.42 551.08 397.92 376.17 578.17 656.67 391.17 596.25 676 685.25 559.17 655.58 660.42 467.75 478.67 727.08 464.67 508.08 853.83 601.67 437.67 518.92 607.08 389.42 615.67 451.67 535.25 585.83 532.25 569.42 478.08 336.17 440.25 461.92 469.83];
%50978
%x=[37.5 34.25 39.08 31 49.92 35 36 40.75 39.75 34.5 37.25 50.92 40.92 33.42 39.17 39.5 36 44.67 51.75 54.25 44.58 44.75 43.33 46.83 51.92 42.58 47.17 53.25 46 40.25 45.92 47.33 54.67 49.08 39.42 42.25 54.08 54 38.5 39.92];
%x=[583.83 450.08 511.08 494.33 261.42 359.25 284.25 429.42 321.92 546.67 682.33 393.33 457.58 482 451.33 491.25 581.08 417.33 427.67 525.08 571 380.75 459.25 558.5 411.67 421.17 412.83 366 316.5 495.75 313.67 517.75 262.92 370 398.75 481.58 404.33 536.08 498.08 462.08];
%51076
%x=[42.33 31.58 45.08 34.25 43 36 46.92 53.83 38.33 43.25 45.08 60.67 56.92 20.67 33.58 43.92 38.92 40 56.67 55.92 53 44.42 36.42 50 55.17 41.92 62.17 49.17 52.33 46 45 57.75 42.08 50.92 42 52.17 62.75 56.42 46.75 33.25];
%x=[170.33 153.5 147.83 102.58 99.33 128.75 123.92 129.92 162.67 144.92 165.58 78.75 142.42 240 110.67 111.25 219.58 180.83 123.17 190.83 124.25 216.17 258.67 175.83 159.5 229.08 121.58 208 130 232.83 213.83 168.83 162.67 183.17 208.5 195.67 179.08 121.42 179.42 282.08];
%51087
%x=[22.42 11.75 18.75 8.75 13 12 21.25 26.58 26.33 30.5 31.83 43.5 42.17 9.25 16.83 30.5 26.92 26.25 43.58 42.92 41.92 32.33 24.75 41.83 51.92 34 49.58 41.67 41.67 36.25 35.33 50.08 33.42 43.75 33.42 43.58 60.17 54.58 41.42 28.08];
%x=[157.58 132.25 144.58 90 93 133.75 124.42 92.17 173.5 129.17 126.08 76.75 138.75 216.25 110.42 168.75 228.92 166.08 142.08 158.5 147.25 245.42 236.25 195.92 137.08 218 164.83 156.42 177.92 257.75 177.83 148.83 201.58 164.92 196.42 157.25 119.17 135.33 178 292];
%51156
%x=[36 27 36.67 31.25 30.25 28.25 39.58 40.33 33.75 37.58 34.58 46.5 43.75 18.75 30.75 36.08 36.5 32.67 40.67 46.08 45.08 37.67 30 38.33 48.42 35.75 57.08 48.33 47.25 40.5 43.17 49 39.25 47.42 37.42 51.08 50.75 51.75 43.92 36.67];
%x=[98.83 174.42 132.25 63 126.5 125 57.92 85.42 93 109.5 124.08 70.83 149.83 112.58 75.25 127.92 98.33 124.08 147.17 112.25 76.25 173.75 214.17 115.92 113.33 89.17 59.92 111.92 127.5 117.83 146.25 157.08 123 137.42 140.17 66.75 237.83 77.17 93.42 174.92];
%51243
%x=[85 77.25 87.83 77.42 84 73.25 83.75 87.42 82.75 84.75 81.08 100.58 102.25 69.08 79.42 80.42 83.83 79.33 92.5 100.08 94.33 83.67 76.42 84.67 92 81.42 101.67 92.33 94.75 89.25 91.75 91.83 85.5 95 83.5 94.33 95.25 99.5 90.5 83.58];
%x=[66.58 84.17 123.67 53.17 53.92 49 70.42 76.5 80.08 112.42 92.25 80.58 61 93.33 58.83 53.58 141.58 152.08 66.25 156.83 80.33 128.08 111 125.5 98.58 68.25 31.25 102.92 103.08 67.58 123.75 144.92 112.75 148.92 125.17 89.92 104.33 53.42 92.75 105.75];
%51334
%x=[77.25 67 83 68.67 72.08 67.75 76.92 78.33 73.42 81.17 75.75 90.75 93.08 65.25 70.42 68.83 74.75 73.42 80.42 88.75 83 74.17 68.08 75.83 80 74.75 93.08 84.42 87.08 85.83 87.67 87 80.33 88.75 78.25 94.83 93.75 94.42 91.67 83.75];
%x=[70.33 74.08 84.75 67.58 45.17 69.17 67.83 84.83 68.42 93.75 94.5 86.08 130.33 71.17 75.92 95.75 136.58 116.67 64.92 82.08 58.08 94.5 124.25 108.33 62.17 85.83 50.67 100.5 101.75 84.33 137 118.83 114.83 101.58 81.58 83.58 114 56.75 80.58 88];
%51379
%x=[49.92 42.92 51.92 49.83 48.42 46 53 56.75 45.42 50.75 52.25 63.08 63.58 34.08 45.42 52.92 49.58 49.92 59 60.5 59.67 51.83 45.58 51.75 54.58 48.67 60.17 57.5 55.33 50.5 57.17 60.92 45.75 55.5 53.67 61.83 59.83 61.42 55.17 52.75];
%x=[157.58 184.83 131.67 96.25 128.5 130.5 128.42 149.5 151.17 130.92 140.33 114.92 137.83 224.5 110 108.33 271.25 157.25 124.33 198.92 188.25 161.25 155.58 229.33 141 171 114.08 263.83 219.25 178.67 105.67 168.42 160 155 157.08 113.5 243 134.25 195 182.25];
%51431
%x=[89 74.17 86.83 75.08 88 79.67 87.92 94.17 91.08 97.08 93.33 90.83 101.58 72.42 81.92 88.5 95 91.08 90.33 97.17 92.08 93.33 85.75 86.67 94.92 86.67 102.5 90.83 101.33 98.58 100.92 104.17 92.33 105.58 99 104.67 104.5 105 101.33 103.33];
%x=[220.25 206.5 228.67 201 135.5 261 188.83 168.25 206.25 229.92 233.92 160.58 151.58 221.75 254.83 214.92 354.83 282.08 194.5 235.83 137 228 307.17 259.33 128.67 232.5 191.5 374.25 224.92 289.33 204.33 348 294.67 413.58 237.5 247.17 243.83 141.33 274.58 387.83];
%51463
%x=[74.42 68.75 77 73.58 74 55.25 68 70.5 61 64.75 66.08 76.25 77.92 50.25 60.25 66.08 64.08 63.25 71.58 76 76.5 66.83 61.33 68.25 77.08 67 88.33 78.25 79.75 72.58 77.08 81 67.67 79.83 75.33 85.92 84.83 87 79.75 73.92];
%x=[189.75 202.83 130.17 109.42 187.25 206.17 133.75 227.83 271.17 199.83 202.17 202.58 230.83 287.42 202.42 210.92 302.83 311.25 218.5 243 148.17 290.33 250 259.5 200.58 325.75 133.17 349.33 280.33 276.92 231.42 285.75 308.33 261.83 230.25 196.58 349.58 143.17 294.25 235.33];
%51573
%x=[140.58 137 141.67 137.75 130.17 126.25 135.58 137.92 134.42 144.67 142.58 150.33 146.42 134.67 140.42 148.5 148.5 148.08 152.33 152.42 151.17 143.25 147.5 149.08 146.42 142.25 152.67 151.75 157.33 153.25 158.75 161.58 148.25 154.33 153.92 158.08 163 161 160.08 156.25];
%x=[19.92 9.92 15.25 16.58 14 3.58 25.08 6.67 19.17 4.17 11.58 3.17 4.08 25.17 6.75 7.17 22.25 22.5 17.42 13.58 7.08 19.33 6 17.75 9.5 8.67 4.58 27.83 8.25 13.67 13.92 21.33 25.75 8.75 7.5 6.83 10.25 19.33 5.75 5.83];
%51644
%x=[118.25 112.67 118.92 107.83 109.08 105 113.75 111.33 114 121.17 109.5 116.83 113.83 107 111 115.75 113.92 112 113.5 117.17 112.42 109.42 110.75 115.42 108.92 105 117.75 114.17 118.33 113.67 115.92 117.83 107.17 113.17 107.67 108.33 116 109.17 113.67 108.83];
%x=[70.58 41.42 43.33 59.08 36.25 57 68.67 78.58 82.33 32.83 66.58 54.17 36.33 29 61.17 72.67 106.75 57.67 82.75 64.08 63.17 83.92 69.42 36.33 85.58 58.17 82 68.42 76.17 38.92 63.83 51.67 106.42 63.75 79.17 56.08 29.92 61.25 40.17 77.17];
%51709
%x=[124.92 114.67 125.58 104.58 111.08 110.42 122.92 119.08 124.75 123.83 116.58 119.75 119.58 112.08 120.17 118.25 118.33 119.42 114 121.08 112.25 110.75 113.5 120.5 110 109.83 126.42 126.92 130 129.17 129.58 129.92 127.92 134.42 126.92 131.25 139.5 129.33 135.83 130.5];
%x=[63.83 51.92 25.25 121.83 36.33 81.83 52.42 51.25 14.42 44.92 85.67 66.42 47.5 58 15.58 41.75 100.83 35.42 47.92 53.83 49.75 56.08 83.75 13.5 56.92 132.17 17.67 55.75 17 21 45.33 65.08 102.08 87.33 86.58 49.92 32.83 52.83 46.08 159.67];
%51716
%x=[120.17 114.08 123.33 110.58 112.67 113.42 122.5 120.25 120.08 125 118.5 122.5 120.08 116.58 121.83 122.5 122.42 122.5 118.83 129.58 117.58 120 121.92 125.92 117.58 113.42 128.33 126.08 128.83 126.83 128.83 127.83 123.83 133.17 123.75 126.75 133.42 124.08 131.08 124.58];
%x=[65.25 79.67 16.33 85.08 10.33 44.17 62.58 11.83 23.25 11.75 103.33 71.83 39.83 31.75 7.17 35.58 73.17 73.08 52.42 42.83 52.75 65.75 114.17 29.75 79.5 97.17 15.67 56 24 50.92 24.17 77.83 101.08 33.83 76.17 57.75 25.17 38.58 26.83 124.33];
%51765
%x=[111.08 106.08 111 104.75 105 102.5 109.92 109.83 104.67 112.33 104.83 110.33 105.17 97.67 106.42 108.67 110.83 108.67 111.17 115.42 113.5 107.58 111.67 114.08 110.25 107.67 115.92 118.33 118.33 113.08 117.25 119.5 112 119.08 114.75 118.92 120.58 116.83 118.5 113.58];
%x=[37.5 34.17 53.33 63.08 34.5 52 20.67 12.5 44 29 62.5 20.75 13.08 22.83 23.58 10.17 49.75 60.67 10.25 49.83 49.58 20.5 53.83 14.33 19.58 15.58 15.17 36.17 33.42 24.08 2.83 15.42 26 23.75 35 35.42 37.75 42.08 9.42 29.33];
%51777
%x=[118.92 112.25 119.33 110.42 112.75 108.92 118.08 116 115.08 118.83 113.5 119.25 115.25 107.83 116.83 117.33 120.33 118.42 117.17 121.75 117.83 113.92 116.67 123.08 117.5 112.67 123.58 128.17 125.17 120.5 124.5 124.83 118.42 128.75 122.25 127 130.17 121.75 124.33 120.42];
%x=[10 17.33 16.92 27.92 9.25 25.08 15.5 7.83 27.42 13.42 89.08 25.58 15.58 34 12.75 9.08 20.33 53.58 15.92 15.75 18.75 34 38.25 7.75 24.5 19.58 11.58 50.42 24.75 31.83 4.25 33.42 47.42 16.75 97.5 33.58 56.92 20.58 19.25 48.17];
%52495
%x=[65.58 72.75 70.25 63.67 72.08 61.92 67.92 75.67 70 70.42 69.17 78.83 71.67 62.25 68.08 69.67 81.83 72.58 78.42 80.83 78.67 74.42 69.67 79.25 69.83 68.33 78.17 89.17 86 77.5 81.25 82.25 76.33 79.75 74.17 83.25 83.42 74.92 80.67 80.83];
%x=[51.58 34.08 137.67 95.25 81.25 85.08 165.33 55.08 153.83 94.67 70.83 31.25 68.08 84.25 41.58 53 50.5 77.25 66 42.17 63.33 50.42 101 112.5 113.58 139.58 62.08 145.92 64.42 82.25 55.92 125.17 123.25 94.25 84.75 101.17 120.83 100.5 99.33 121.33];
%52681
%x=[78 77.92 81.33 75.58 79.08 73.25 78.17 82.42 80.25 81.25 79.75 86.42 78.5 73.33 79.83 80.67 91.42 82.17 84.83 90.75 87.58 83.5 81.5 90.83 82.58 82.67 89.33 99.75 98 90.5 92.17 94.58 89.92 93 90.33 96.75 94.42 86.5 96.83 95];
%x=[114.5 68.75 154 63.42 101.83 60.42 125.5 140.92 113.17 84.92 76.08 43.33 112.58 80.17 120.67 72.58 96.33 92.33 64.42 67.42 66.25 96.25 89.75 168.33 83.42 115 77.17 95.75 78.67 101.75 80.75 129.83 124.17 83.5 81 105.67 124.08 120.42 88.92 93.17];
%53068
%x=[31.5 34.67 38.42 31.42 46.92 34.5 32.33 36.33 41.92 36.92 27.58 44.5 40.42 32.42 23 25.92 35.42 38.58 46.42 45.42 42.08 38.5 36.17 49.75 50.58 44.25 53.58 59.08 56.5 40.83 55.5 58.83 43.67 60.5 50.5 54.92 64.92 60 57.08 49.75];
%x=[95 100.25 203.58 79.17 149.67 150.17 130 78.83 120.17 78.08 145.5 75 100.67 110.83 116.17 171.67 81.08 150.67 55.33 198.5 95.25 83.58 103.25 126.42 125.83 214.08 103 141.42 84.83 89.17 33.08 106.92 162.58 94.25 52.58 99 107.33 112.92 125.67 108.58];
%54161
%x=[50.67 49.25 54 47.58 64.75 44.83 50.42 49.92 54.92 44 47.83 64.5 58.92 49.17 47.83 50.25 50.83 57.92 66.83 68.58 56.58 60.08 58.42 65.42 65 58.92 66.92 73.5 60.67 56.5 60.58 68.58 70.5 71.42 56.42 65.75 76.83 72.08 61.17 51.5];
%x=[419.67 365.42 624 476.08 400.42 496 515.67 386.33 456 525.17 447.33 274.75 490.67 544.58 684.92 653.5 494.17 391.42 519.92 482.75 583.42 432.75 391.75 574.5 466.67 408.5 479 519.42 408 346.92 324.92 408.83 431.42 396.83 567.5 527.17 445.17 597.33 400.83 731.92];
%54337
%x=[90.92 90.5 94.83 85.92 99.17 84 88.25 90.92 92.5 86.25 92.75 98.75 98.75 89.08 83 88.42 88.17 96.5 104 99.17 95.17 100.75 97.33 104.25 100.42 95.67 103.92 104.83 104.83 97.92 96.5 106.42 104.58 108.33 98.42 101.42 110.08 109.42 105.25 93.83];
%x=[414.75 297.75 479.83 681.08 531.08 438.33 642.08 376.83 442.58 324.17 365.5 290.42 405.58 446.5 617.92 488.5 562.33 428.42 300.58 595.08 650.67 303.42 450.08 697.92 524.75 551.08 423.75 765.25 325.08 371.33 353.58 349.83 322.25 609.58 548.17 419.17 469.92 531 370.17 674.42]
%53614
%x=[81.83 85.58 89.17 83.33 87 81.17 87.42 89.33 90.67 88.33 90.17 92 87.5 78.17 84.67 84.67 97.25 87.92 89.75 96.33 95.58 87.92 84.33 96.08 90.25 89.08 97.58 105.33 102.58 96.17 101.5 101.08 97.58 103.5 100.75 108.67 103.92 98.58 104.75 102.67];
%x=[108.83 107.58 250 110.08 156.42 204.25 203.08 224.83 195.75 81.83 83 104 149.75 185.83 153.08 130.67 134.58 167.42 177.75 211.08 162.42 239.25 102.83 134.83 194.25 127.67 130.33 175.83 137.58 111.5 136 253 162.33 120 62.42 163.17 178.92 162.17 150 171.92];
%53845
%x=[94.67 95.25 100.83 96.75 95.5 89.17 97.17 99.92 99.58 95.5 95.58 99.33 92.92 87.75 93.08 89.25 104.33 94.75 98.58 101.58 101.75 96.92 95.5 104.75 103.83 97.42 110.33 114.67 113.42 109.83 111.33 111.83 106.33 109.67 107.58 114.42 112.08 103.83 112.58 110.67];
%x=[399.58 307.33 494.33 275 557.83 428.33 393.42 479.58 459.42 399.5 645 385.42 567.25 425.58 507.08 374.33 407.92 619.58 372.83 464.75 376.42 385.42 472.83 505.75 300.58 387.33 310.08 473.17 286.92 306.08 477.92 448.75 548.33 330.25 367.67 395 447.5 359.83 524.17 388.17];
%53772
%x=[92.92 95.5 100.25 95.17 100.83 90.25 99.67 99.83 96.58 95.67 96.67 102.5 96.75 89.42 94 95.5 101.25 98.25 100.42 102.08 99.58 98.25 94 105.33 101.75 97.58 108.33 114.92 116.17 107.25 110 110.5 102 108.92 109.33 117.58 114.17 109 111.08 112.67];
%x=[397.25 180.08 520.83 298.25 359.92 410.92 454.5 377.58 433.5 277.58 284.75 372 415.83 288 452.08 214.75 338.17 483.92 349.17 391 368.25 323.75 321.75 357.83 407.67 543.33 206.5 310.75 290.25 349.42 248.33 350.08 438.58 314.33 228.92 354 446.17 296.08 520.92 313.83];
%54511
%x=[112.08 113.92 114.92 111 122.17 108.42 115.58 116.08 110.83 110.42 123.33 128 129.83 119 115.33 121.25 122.5 127 131.92 126.92 125.08 128.42 129.5 136.92 133.33 127.42 130.5 131.08 131.33 128.33 128.83 131.5 129.17 135.33 131.92 134.17 140.08 133.75 132.5 125.58];
%x=[425.83 311.67 581.75 395.42 326 569 649.08 553.58 598.5 317.25 327.67 453.67 408.25 407.33 600.83 554.42 569.92 561.08 368.5 581.08 623.25 451.25 422.25 677.67 477.08 584.08 359.08 609.75 222.42 309.25 282.42 308.67 370.75 402.92 342.25 265 403.25 521.92 400.5 435.42];
%54539
%x=[99 100.08 101.67 97.58 109.42 96.08 103.5 104.42 101.92 100.08 103.25 108.83 109.33 97.75 94.83 100.5 102.5 106.75 111.58 109.42 106.33 108.33 110.75 113.83 107.42 104.67 113.08 118.25 116.25 114.42 114.42 119.92 115.83 120.42 115.25 117.25 122.58 118.5 120.58 110.33];
%x=[661.33 272.67 589.75 653.33 514.08 592.83 631.92 386.17 569.92 421.92 509.58 373.67 349.75 427.17 726.25 680.08 634.83 673.83 361.58 568.75 471 464.67 258.08 620.83 617.42 403.75 457.17 568.58 309.17 402.08 258.58 253.25 641.92 519.5 497.25 367.5 561 469.5 358.5 530.58];
%55591
%x=[72.25 78.67 77.75 79.92 80.33 78.42 73.33 72 77.33 74.67 73.5 75.83 74.33 81.67 81.25 78.5 79.58 83.5 84.75 76.08 80.58 77.42 81.92 84.67 88.5 85.42 75.25 90.08 92 83.83 87.92 85.92 87.58 85.08 92.75 97.08 97.58 89.75 102.75 99.75];
%x=[346.08 347 345.25 402.67 270.25 288.17 418.58 375.92 384.92 406.33 275 252.33 191.33 344.92 441.5 240.5 322.08 348.33 303.58 511.5 396.92 243.08 447.67 309.83 431.5 374.08 267.92 484.08 448.42 441.42 410.25 449.42 458.08 462.67 413.17 282.75 397.75 444.83 286.67 299.83]
%56004
%x=[-39.67 -40.08 -44.42 -33.17 -38.75 -41.58 -40.83 -41.67 -40.33 -40.33 -36.58 -40.58 -48.08 -38.08 -72.83 -60.67 -40 -35.5 -41.83 -38.92 -35.58 -43.42 -39.5 -35.08 -37.25 -37.83 -50.17 -29.83 -42.67 -41.92 -34.75 -31.83 -29.58 -31.25 -28.42 -20.33 -29.25 -33 -21.92 -20.58];
%x=[275.92 288.33 208.58 324.08 249.33 217.5 244.08 208.25 150.33 224.25 274.58 297.92 251.5 136.83 382.83 168.67 192.33 181.33 258.25 171.67 208.5 180.92 251.92 135.58 186.33 224.75 204.08 241.33 279.42 268.67 278.92 351.83 209.58 233.17 314.08 200.08 293.25 275.42 419.17 276.42];
%56029
%x=[25.5 32.42 27 35.92 32.25 30.25 23.42 26.75 33 29.5 30.92 30.08 23.58 39.25 30.42 30.58 33.42 40.25 34.5 30.75 33.33 28.75 32.25 36.08 33.75 36.17 25.67 40.58 45.92 32.25 36.33 41 44.75 43.67 45.5 50.83 48.67 42.75 52.42 49.17];
%x=[354 388.33 401.08 435.42 409.25 438.33 318.33 367.58 446.17 484 512.67 457.58 462.33 268.08 473.58 313.58 375.75 359.67 531.92 404 409.75 389.08 436.17 354 373.25 350 355.5 457.75 396.25 424.25 409.33 291.33 481.67 364.33 456.25 332.25 315.92 427.75 503.33 353.42];
%56964
%x=[171.08 178.75 180.17 177.17 178.42 176.17 177.42 178.75 180.92 183.33 182.5 181.25 181.17 183.92 182.25 182.67 186.83 185.75 186.42 184.5 188.25 180.92 184.17 189.25 189.08 188.42 187 195.25 193 188.42 193.08 192.67 195.25 187.08 193.83 195.75 190.92 189.83 195.5 201.08];
%x=[1560.33 1252.92 1433.17 958.33 1260.92 1358.42 1302.33 944.42 1085.5 1001.92 1521.25 1254.75 1328.92 1495.25 1353.33 1296.75 1091 931.25 1149 1303.5 1449.33 1008.17 1191.42 1435.25 1209.5 1175.33 1290.75 1186.75 1412.25 1185.33 1623.17 1227.67 899.92 1085.42 1117.42 1169.5 1191.58 1307.5 1263 1032.5];
%x=[146 149.5 154.25 148.42 149.33 145.58 146 150.33 152.67 150.75 152.58 146.67 145.67 148.08 147.17 148.5 152.25 156.08 151.5 150.67 152.75 147.5 149.92 155.25 154.75 153.67 148.33 157 162.17 154.67 158.92 158.25 157.75 155.25 161.08 159.92 156.25 156.5 161.25 160.5];
%x=[1159.5 941 1470.67 1554.08 1330.33 1186 1306.42 1173.42 1224 1173.17 1102.83 1128.42 1561.08 1461.08 1388.5 1126 1247.08 1394.5 1329 1168.83 1470.08 1172.5 1381.58 1077 1491.25 1176.67 1193.83 1034.17 1496.08 1259.33 1449.92 1157.58 1122.75 1556.5 1166.42 979.08 1487.33 1190.42 987.58 1551.5];
%59293
%x=[210 212.83 214.33 210.08 211.58 207.42 215.17 213.33 212.5 215.67 214.33 214.17 211.42 208.17 210.75 215 221.08 210.67 217.25 219.75 222.25 213.33 214 219 213.25 215.92 216.67 225.75 222.5 221.58 223 227 224.5 222.75 221.67 218.83 216.75 209.83 217.75 212.17];
%x=[1267 1584.17 2137.67 1564.33 1934.67 1523.5 1426.92 1351.17 1626.5 1684.58 2102 1576.08 2259.17 1919.33 1637.92 1489.42 1628.75 1638.42 1399.17 1409.92 1049.92 2192.17 2125.75 1804.67 1554.92 1410.25 2338.5 1757 1127.83 1628 1922.42 1401.17 1142.67 1123.08 1703.67 1861 1498.75 1749.75 988.33 1675.33];
%59663
%x=[219.25 220.42 227 219.83 222.17 216.25 224.42 222.08 221.92 225.5 225.5 222.67 219.5 217.58 220.08 222.92 230.42 221.92 224.83 227.17 231 221.83 225.17 230.75 223.17 225.42 229.25 237.58 234.17 233.25 232.17 234.42 236.75 226.08 222.58 228.83 228 220 225.92 224.25];
%x=[1267.25 2781.67 2753.67 1683.08 2257.5 1928 997.58 2476.75 1977.25 1752.25 2581.5 1758.08 2053.83 1969.92 1931.25 2107.83 1995 2222.92 1209.17 2173.08 1335.58 2326.33 2309.33 2183.75 2583.42 2168.83 2554.67 2589.42 1714.5 1424.17 3009.42 2472.75 1409.17 1512.5 1370.75 1623 1210.58 2406.75 2348.08 2318.67];
%59758
%x=[231.83 238 244.17 232.33 237.25 234 238.42 239.08 241.75 244.17 243.5 240 237.58 233.25 235.33 239.17 248.17 239.67 239.83 243.33 249.25 241.42 243.42 247.58 241.92 243.08 246 254.08 246.42 245.25 248.67 251.25 252.58 247.42 250.58 253.75 241 234.33 242.75 245.83];
%x=[1261.67 1888.25 1570.5 1470.83 1327.17 1670.33 726.25 1521.42 1257.42 1235.33 1120.92 1952.25 1119.92 1143.83 1599 1152.42 856.33 1431.67 1798.5 1185.58 1082.08 1303.75 1269.08 1621.42 1494.17 1171.67 1827.42 1306.25 1129.92 1803.08 1872.5 1485.17 1185.5 819.67 989.17 1289.75 1182.75 1992.67 2190.17 2037.58];
%56096
%x=[147.42 147.25 149 144.42 144.42 136.33 145.42 147.92 147.58 144.25 143.5 140.75 142.17 136 142.67 145.42 151.5 142.08 142.42 148 148 141.08 141.67 151.25 148.5 147.58 155.17 159.33 153.33 150.67 156.42 157.67 152.83 153.42 151.92 162.75 156.17 151.58 153.92 153.67];
%x=[298.08 347 406 326.17 443.75 340.58 373.42 496.42 358.67 477.5 412.42 412.67 398.83 574.42 414.5 313.58 335.08 509 391.83 507.75 336.25 488 521.75 335.25 364.08 283.25 225.42 399.83 331 375 314.25 346.42 433.17 343 341.75 314.92 381.08 407.83 423.33 281.92];
%57127
%x=[139.25 142.08 147.92 141.75 144 136.42 146.92 146.67 145.58 139.83 140.67 140.42 138.75 138.25 140.17 139.92 144.92 141.92 138.33 146 145.08 140.42 139.75 147.67 149.17 141.25 150.83 152.5 150.67 149.83 150.25 154.25 146.58 152.67 151.17 157.5 154.75 151.58 152.58 151.17];
%x=[616.42 649.75 904.33 713.58 737.58 654.75 588.75 577.17 608.08 908.33 998.75 848.17 1219 832.33 713.58 533.17 777 569.42 942.17 777.42 500.08 644 661.92 612.25 441.92 628 432.58 867.25 608.75 748.75 488.58 525.42 755.33 604.67 685.17 598 663.92 811.42 706.33 751.33];
%57494
%x=[158 158.67 164.67 162.08 165.33 160.42 162.67 169 167.17 159.5 163.33 163.25 165 158.92 161.17 164.25 165.83 166.5 161.92 171.17 164.33 167.75 162.17 173.33 173.67 167.92 175.42 181.58 176.08 177 180.25 180.33 174.92 182.67 178 183.17 185.5 176.58 178.5 165.75];
%x=[666.83 896 1024.58 804.58 1100.17 742.25 995.83 676.42 834.92 1353 961.67 1360.33 1579.08 1007.5 858.08 875 1207.83 1110.25 1379.08 1129.17 1496 930.33 1320.5 871.25 1080.25 1099.58 788.83 1441 1150.5 983.17 749.83 1263.42 1155.08 1310.17 930.5 872.58 852.67 1055.67 965 1114.92];
%58238
%x=[151.25 147.08 154.17 151.5 155.83 149.5 151.75 160.42 156.08 145.75 150.5 153.33 153.17 147.58 150.67 150.33 152.67 153.5 153.67 160.92 152.17 154.08 151.17 166.58 157.58 154.42 161.75 167 158.5 164.5 165.92 166.83 160.75 169.58 163.67 169.5 173.5 161 163.5 162.17];
%x=[841.92 1066.33 773.42 1073.67 1115.75 743.08 966.83 445 797.42 852.58 874.92 922.08 936.67 799.92 844.42 602.08 1148.25 770.75 1051.08 792 1521.5 737.67 1034.58 539.92 642.25 1011.25 752.33 1032.5 1012.08 858 614.42 895.5 1381.92 812.58 826.92 922.33 892.42 812.5 1136.25 1082];
%58633
%x=[173.5 170.58 173.25 171.08 173.42 168.33 170.83 175.33 177.33 168.67 167.58 171 173 166.25 171.17 171.42 173 171.33 170.75 177.25 175 170.58 171.67 180.25 171.75 171.75 175.33 185 177.33 176.5 176.5 179.67 180.75 177.67 176.33 182.67 185.25 178.58 179.75 175.25];
%x=[1074.67 1319.42 1883.42 1379.25 1839.83 1478.42 1688.67 937 920 1395.08 1322 1227.92 2053.75 1203.33 1162.83 987.58 1573.58 1083 1811.67 1286.08 1144.5 1649.75 1776.42 1488.33 1687 1043.92 1515.42 1874.83 1448.08 1369.33 1208.67 1918.92 1249.08 980.75 1105 1111.92 1084.42 1105.17 1145.67 2010.42];
%57993
%x=[191.92 192.75 193.75 192.83 192.08 187.42 194.83 197.5 198.25 192 192.42 193 195.92 185.58 190 195.17 199 190.67 194.58 196.67 196.5 193.17 191.33 194.58 192.83 191.42 194.5 205.83 199.67 194 196.42 199.25 200 197.5 194 197.75 202.5 195.67 201.92 199.08];
%x=[854.58 1151.92 1547.08 988.58 1582.17 1320.92 1095.58 1054.58 1036.25 1223.42 1304.42 1414.5 1660.58 1239.58 1295.33 808 1033.67 1098.5 1022.17 1210.75 879.75 1600.92 1305.33 1544.17 1134.83 1005 1374.25 1352.67 1066.33 1324.75 1281.08 1705.92 819.17 1079.33 1104.42 1347.92 991.92 1081.58 964.17 1105.5];
%57957
%x=[185.42 186.75 191.5 187.58 185.83 181.75 187.42 190.08 191.42 187.33 186.17 187 189.17 179.42 185.58 189.75 191.83 185.42 188.08 192.25 189.83 191.67 185.83 189 190.67 187.25 189.08 198.92 195.5 188.5 192.67 193.92 195.92 196.5 191.58 196.17 201.25 193 199.83 195.58];
%x=[1505 1489.67 1726.75 1420.42 1713.17 1818.75 1784.83 1713.75 1405.33 1703.75 1591.25 1925.33 1619.5 1213.42 1397.75 1414.08 1581.25 1249.42 1366.67 1532.33 1320.42 1686.25 2215.58 1888.58 1186.33 1755.33 1624.58 1786.17 1681.67 1713.42 1182.58 2339.17 1295.58 1535.75 1564.5 1473.5 1164.42 1783.83 1555.33 1546.42];
%53068
%x=[31.50 34.67 38.42 31.42 46.92 34.50 32.33 36.33 41.92 36.92 27.58 44.50 40.42 32.42 23.00 25.92 35.42 38.58 46.42 45.42 42.08 38.50 36.17 49.75 50.58 44.25 53.58 59.08 56.50 40.83 55.50 58.83 43.67 60.50 50.50 54.92 64.92 60.00 57.08 49.75];
%51156
%x=[25.00 36.00 27.00 37.00 31.00 30.00 28.00 40.00 40.00 34.00 38.00 35.00 47.00 44.00 19.00 31.00 36.00 37.00 33.00 41.00 46.00 45.00 38.00 30.00 38.00 48.00 36.00 57.00 48.00 47.00 41.00 43.00 49.00 39.00 47.00 37.00 51.00 51.00 52.00 44.00];
%50968
%x=[23 27 20 24 20 38 17 21 25 28 20 25 38 30 22 22 24 21 30 36 41 31 31 31 35 40 34 42 48 33 29 31 37 43 42 34 39 51 46 32];
%53463
%x=[53 54 61 64 56 71 56 63 66 62 62 60 69 68 56 61 63 74 64 75 74 72 69 65 77 68 65 78 83 85 75 82 81 72 80 77 86 90 74 80];
%56182
%x=[56 55 58 57 58 58 52 51 56 59 59 61 53 54 59 58 56 64 63 58 60 61 54 59 65 60 61 57 70 69 60 65 65 69 63 69 73 68 65 71];
%56021
%x=[-28 -25 -18 -23 -18 -23 -27 -29 -26 -23 -24 -19 -23 -33 -18 -36 -29 -22 -16 -20 -22 -15 -25 -20 -16 -17 -14 -32 -11 -17 -19 -16 -18 -11 -14 -9 -5 -11 -14 -2];
%54405
%x=[83 87 91 93 85 101 86 94 96 89 90 92 99 98 89 87 91 97 97 101 99 96 97 100 107 99 95 105 111 110 100 106 109 103 106 102 107 112 99 102];
%58477
%x=[161 163 159 165 162 166 157 165 166 168 159 160 164 164 160 163 161 163 162 163 171 165 162 161 173 163 164 169 176 169 171 172 175 172 174 167 176 179 169 172];
%57713
%x=[148 151 152 158 148 153 146 150 154 156 151 154 150 154 147 153 152 158 154 150 160 154 153 151 157 154 151 155 162 159 154 160 160 162 157 158 166 164 158 163];
%降水
%57290
%x=[6584 10930 7687 8234 9201 12092 7473 11413 7109 11977 9782 9287 17916 11420 14965 8256 5869 11376 5664 11908 10632 10377 4763 7381 7021 7489 12393 5440 13503 6443 15750 6220 8848 13857 7122 13943 10762 11532 9891 8547];
%59082
%x=[17079 12805 16590 19559 14996 21208 16309 14432 14775 11312 14594 17909 16114 20381 14006 13602 11693 13999 12830 12065 14366 11743 18887 20279 21321 15069 16331 20453 18623 13143 15658 16898 18149 13882 12518 17722 17828 15023 15531 12755];
%50727
%x=[4552 4515 3857 3937 3895 4661 4505 4181 3538 4091 4257 4913 5248 5188 4690 5233 4181 3657 5878 4427 6183 4562 4378 4441 6080 3989 4812 3916 6412 2718 3523 3399 4068 4386 2836 4522 3809 3119 4086 4196];
%54337
%x=[5495 4977 3573 5758 8173 6373 5260 7705 4522 5311 3890 4386 3485 4867 5358 7415 5862 6748 5141 3607 7141 7808 3641 5401 8375 6297 6613 5085 9183 3901 4456 4243 4198 3867 7315 6578 5030 5639 6372 4442];
%56021
%x=[4356 4295 3846 3571 4735 4969 4392 4007 3267 2878 4187 4320 4575 5077 3553 4248 3803 4987 4426 5223 3068 4131 3850 4792 3414 3763 3600 3232 3631 4364 3455 4001 3284 3853 4859 5458 3967 4463 4953 5180];
%53772
%x=[2965 4767 2161 6250 3579 4319 4931 5454 4531 5202 3331 3417 4464 4990 3456 5425 2577 4058 5807 4190 4692 4419 3885 3861 4294 4892 6520 2478 3729 3483 4193 2980 4201 5263 3772 2747 4248 5354 3553 6251];
%51828
%x=[267 183 896 262 391 124 188 435 259 276 54 451 366 317 239 34 250 1009 816 223 261 393 617 690 320 287 629 106 326 323 192 260 638 861 403 884 417 238 343 192];
%56571
%x=[9586 10105 7288 10886 12635 8454 9591 8685 9109 8122 10516 9386 8416 8326 9591 11268 10438 10913 9938 8422 11120 12791 7326 12348 10587 10311 7577 9782 15492 13234 11402 11765 10033 8528 10190 8257 9856 9630 10649 9081];
%58921
%x=[17686 12418 15647 18695 14743 23348 15739 14229 17747 13806 14589 16635 16726 17881 14070 16081 14845 15833 16568 13095 15657 10648 17328 13588 16586 14692 11991 18948 18179 12257 16557 19068 16055 9737 12758 18272 18182 13298 12738 14489];
%53391
%x=[19 19 26 29 19 33 18 24 28 26 22 20 29 31 16 15 23 34 24 36 34 29 28 27 38 30 26 40 46 42 30 38 39 31 39 31 40 47 33 38];
%54662
x=[101 102 101 107 100 113 98 104 109 109 99 104 113 114 104 97 103 103 112 118 113 112 113 112 118 114 109 118 118 121 115 115 119 113 122 109 114 123 114 115];
if nargin == 2
if isstruct(varargin{2})
inopts = varargin{2};
else
error('when using 2 arguments the first one is the analyzed 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
% 默认停止条件
defstop = [0.05,0.5,0.05];
opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask','ndirs','complex_version'};
% 时间序列，停止参数，是否演示，最大迭代次数，每一轮迭代次数，IMF个数，插值方法，每一轮迭代次数（带条件），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;
defopts.ndirs = 4;
defopts.complex_version = 2;
opts = defopts;
if(nargin==1)
inopts = defopts;
elseif nargin == 0
error('not enough arguments')
end
names = fieldnames(inopts);
for nom = names'
if ~any(strcmpi(char(nom), opt_fields))
error(['bad option field name: ',char(nom)])
end
if ~isempty(eval(['inopts.',char(nom)]))
eval(['opts.',lower(char(nom)),' = inopts.',char(nom),';'])
end
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;
ndirs = opts.ndirs;
complex_version = opts.complex_version;
if ~isvector(x)
error('X must have only one row or one column')
end
if size(x,1) > 1
x = x.';
end
if ~isvector(t)
error('option field T must have only one row or one column')
end
if ~isreal(t)
error('time instants T must be a real vector')
end
if size(t,1) > 1
t = t';
end
if (length(t)~=length(x))
error('X and option field T must have the same length')
end
if ~isvector(stop) || length(stop) > 3
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 size(stop,1) > 1
stop = stop';
end
L = length(stop);
if L < 3
stop(3) = defstop(3);
end
if L < 2
stop(2) = defstop(2);
end
if ~ischar(INTERP) || ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
error('INTERP field must be ''linear'', ''cubic'', ''pchip'' or ''spline''')
end
% 使用mask信号时的特殊处理
if any(mask)
if ~isvector(mask) || length(mask) ~= length(x)
error('masking signal must have the same dimension as the analyzed signal X')
end
if size(mask,1) > 1
mask = mask.';
end
opts.mask = 0;
imf1 = emd(x+mask, opts);
imf2 = emd(x-mask, opts);
if size(imf1,1) ~= size(imf2,1)
warning('emd: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 % 如果两个信号分解得到的IMF个数不一致，调整顺序
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
MODE_COMPLEX = ~isreal(x)*complex_version;
if MODE_COMPLEX && complex_version ~= 1 && complex_version ~= 2
error('COMPLEX_VERSION parameter must equal 1 or 2')
end
% 极值点和过零点的个数
ner = lx;
nzr = lx;
r = x;
if ~any(mask) % 如果使用了mask信号，此时imf就已经计算得到了
imf = [];
end
k = 1;
% 提取每个模式时迭代的次数
nbit = 0;
% 总体迭代次数
NbIt = 0;
end
%---------------------------------------------------------------------------------------------------
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
