%% load data
clc;clear;
tp_sh = load('d:\document\study\group_wu\code\index\tp_sh_stat_avg_1980-2014.data','-ascii');
%tp_sh_1 = load('d:\document\study\group_wu\code\index\tp_sh_stat.data','-ascii');
%tp_sh =tp_sh_1(:,2:end);
%% do ceemdan 
rmpath('c:/Users/rex18/Documents/MATLAB/eemd');
addpath('c:/Users/rex18/Documents/MATLAB/CEEMDAN_V00');

data = reshape(tp_sh',1,[]);
Nstd = 0.2;
NR = 500;
MaxIter = 5000;

imfs=ceemdan(data,Nstd,NR,MaxIter,1)';
%% period
%mean period
addpath('c:/Users/rex18/Documents/MATLAB/eemd');
rmpath('c:/Users/rex18/Documents/MATLAB/CEEMDAN_V00');

dimIMFs=size(imfs)
for j=1:dimIMFs(2)-1
    omega = ifndq(imfs(:,j),1 );
    T(j) = 2*pi/mean(omega);
end 
rmpath('c:/Users/rex18/Documents/MATLAB/eemd');
addpath('c:/Users/rex18/Documents/MATLAB/CEEMDAN_V00');

%zero method
for j=1:dimIMFs(2)-1
    kk=0;
    c1=imfs(:,j);
    for i=1:(length(c1)-1)
        if c1(i)*c1(i+1)<0
            kk=kk+1;
        end
    end
    T0(j) = length(c1)*2/kk;
end

%% imfs plot
[data_length num_imfs]=size(imfs);
t=1:data_length;

%plot origin data
figure;
subplot(num_imfs+1,1,1);
plot(t,data);
ylabel('ORIGIN')
set(gca,'xtick',[])
axis tight;
%plot IMFs
for i=1:num_imfs-1
    subplot(num_imfs+1,1,i+1);
    plot(t,imfs(:,i));
    ylabel (['IMF ' num2str(i)]);
    set(gca,'xtick',[])
    xlim([1 length(data)])
    axis tight;
    %每张图的右上角写字符串
    a_handle=get(gca);
    x=a_handle.XLim;%获取横坐标上下限
    y=a_handle.YLim;%获取纵坐标上下限
    k=[0.8 1.2];%给定text相对位置
    x0=x(1)+k(1)*(x(2)-x(1));%获取text横坐标
    y0=y(1)+k(2)*(y(2)-y(1));%获取text纵坐标
    text(x0,y0,['mean period:  ' num2str(T0(i),'%3.1f')]);
end;

%plot R
i=i+1;
if num_imfs>2
    subplot(num_imfs+1,1,i+1)
    plot(t,imfs(:,i))
    ylabel('R')
    axis tight;
end

%x轴的label改为年
set(gca,'XTick',[1:60:data_length]);
set(gca,'XTickLabel',cellstr(num2str([1980:5:2014]','%04g')));