%Copy These Codes to a New MATLAB Function File and Named As TREC.M

%copy these codes to a new MATLAB function file and named as “TREC.m”

% This is a program used to develop chronologies based on EEMD

% USEAGE EXAMPLE: y=TREC('cana325.rwl_raw')

% NOTE:

% eemd program should be added to current path or the path of MABLAB toolbox (C:\Program Files\MATLAB\R2012a)

% eemd program could be downloaded from the website http://rcada.ncu.edu.tw/

% reference: Wu and Huang, 2009, Ensemble empirical mode decomposition:

% a noise-assisted data analysis method. Advances in adaptive data

% analysis. 1(1): 1-41

% function chronology=TREC(file)

% INPUT:

% file: raw measurement file in format rwl_raw

% FILE NOTE:

% Other raw measurement format such as tuson could be transformed to rwl_raw format by

% ARSTAN program

% The file format is listed as:

% year site1 site2 site3

% 1901 0.31 0.34 -9.99

% 1902 0.45 0.56 -9.99

% (-9.99 was the missing value)

%OUTPUT:

% Chronology:tree ring chronology

% Trend.tiff: trend for every tree-ring series%

%Reference: Xianliang Zhang, Zhenju Chen, 2015, A new method to

% remove the tree-growth trend based on Ensemble Empirical Mode

% Decomposition. Trees-Structure and Function

% The code is prepared by Xianliang Zhang. For questions, please contact

%

function chronology=TREC(file)

a1=importdata(file,' ',1);

data=a1.data;

clear a1

ss=size(data);

initial=data(1,1); % inial year

tt1=ss(1,1); % number of years

tt2=ss(1,2)-1; % number of cores

year1=data(:,1);

trin=zeros(tt1,tt2);

tmp1=data(:,1);

for i=2:tt2+1

tmp2=data(:,i);

tmp3=tmp1(tmp2>0);

year(i-1,:)=tmp3(1,1);

end

for i=1:tt2

data1=data(:,i+1);

data2=data1(data1>0);

yy=eemd(data2,0.2,500);

fi=['trend',num2str(i),'.tiff'];

clf

figure

plot(yy(:,1))

hold on

plot(yy(:,end))

saveas(gcf,fi)

clear data2

tree=yy(:,1)./yy(:,end);

t1=year(i,1);

t2=size(tree);

t3=t2(1,1);

m=t1-initial+1;

n=m+t3;

trin(m:n-1,i)=tree;

clear tree

end

%% average mean of tree ring indice

%trin(end,:)=[];

% sample size

%for i=1:tt1

% a1=trin(i,:);

% ll(i,1)=length(find(a1>0));

%end

% average mean

%ind_m=sum(trin,2)./ll;

%% biweight results

for j=1:tt1

th=1e-18;

a1=trin(j,:);

a1=a1';

a2=a1(a1>0);

it1=mean(a2);

st=median(a2);

ss1=size(a2);

t=1;

while t>=th

for i=1:ss1

tmp1=((a2(i,1)-it1)/6/st)^2;

if tmp1<1

wt(i,:)=(1-((a2(i,1)-it1)/6/st).*((a2(i,1)-it1)/6/st))^2;

else wt=0;

end

end

it2=wt.*a2;

clear wt a2a1

it3=mean(it2);

t=it3-it1;

it1=it3;

end

ind_bi(j,1)=it1;

end

chronology=[year1,ind_bi];

std_crn=trin;

save chronology.mat chronology

save std_crn.mat std_crn

fname1=file(1:end-8);

fname=[fname1, 'chronology.txt'];

fid = fopen(fname,'wt');

fprintf(fid,'%4.0f %6.4f \r ', chronology');

fclose(fid);