% preparing high-f/low-f hemispherical land(CRU) or land+ocean(HAD) mean surface temperature
% series from Brohab et al's CRUTem3v/HadCRUT3v.


clear all; close all
%nnn=load('/holocene/s1/zuz10/work1/temann/newtemp/CRU_NH_reform')
%sss=load('/holocene/s1/zuz10/work1/temann/newtemp/CRU_SH_reform')
nnn=load('/holocene/s1/zuz10/work1/temann/newtemp/HAD_NH_reform')
sss=load('/holocene/s1/zuz10/work1/temann/newtemp/HAD_SH_reform')

nt=nnn(1:146,2);
st=sss(1:146,2);
gl=(nnn(1:146,2)+sss(1:146,2))/2;


save('/holocene/s1/zuz10/work1/temann/newtemp/tempn95infxxx','nt','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/temps95infxxx','st','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/tempg95infxxx','gl','-ascii')

%------------ get standardized NH/SM/GL mean (1850-1995)
respliced=nt;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtempn95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtempn95infxxx','stdtempn95','-ascii')

respliced=st;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtemps95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtemps95infxxx','stdtemps95','-ascii')

respliced=gl;
% standardize sst/sat based on interval 1850-1995 for the full series(1850-1998)
[nrows,ncols]=size(respliced)
means=nanmean(respliced(1:146,:));
stdevs=nanstd(respliced(1:146,:));
respliced=respliced-repmat(means,nrows,1);
respliced=respliced./repmat(stdevs,nrows,1);
stdtempg95=respliced;
save('/holocene/s1/zuz10/work1/temann/newtemp/stdtempg95infxxx','stdtempg95','-ascii')

%---------------------------------------------------------------

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtempn95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtempn95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.05;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
%%%%%

save('/holocene/s1/zuz10/work1/temann/newtemp/nhhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/nhlinfxxx','lowf','-ascii')

%%%%%%%

%---------------------------------------------------------------

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtemps95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtemps95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.05;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
%%%%%

save('/holocene/s1/zuz10/work1/temann/newtemp/shhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/shlinfxxx','lowf','-ascii')

%%%%%%%

%---------------------------------------------------------------

clear all

load /holocene/s1/zuz10/work1/temann/newtemp/stdtempg95infxxx

    cutoff=20; %years
	lowlim=1./(cutoff./2)

stdkeepers=stdtempg95infxxx; % because standardization have done in doannsst.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[nrows,ncols]=size(stdkeepers) 
gotstart=zeros(1,ncols);
gotstop=zeros(1,ncols);
start=zeros(1,ncols);
stop=zeros(1,ncols);
stop(:)=nrows;
[lowa,lowb]=butter(5,lowlim,'low');

%%%%%%%%%%%%%%  
%%%% This section finds the starting and stopping year for each column. It assumes that
%%%%  there are no missing values between start and stop, otherwise filtering chokes.

for i=1:nrows
	for j=1:ncols
		if ~isnan(stdkeepers(i,j)) & gotstart(j)~=1
			start(j)=i;
			gotstart(j)=1;
		elseif isnan(stdkeepers(i,j)) & gotstart(j)==1 & gotstop(j)~=1
			stop(j)=i-1;
			gotstop(j)=1;	
		end
	end
end
%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%%	Apply the zero-phase butterworth filter but only to the proxies
	highf=zeros(nrows,ncols);
	lowf=zeros(nrows,ncols); 
	lowf=lowf./lowf;
	highf=highf./highf;
	for k=1:ncols
		clear y y2
		data=stdkeepers(start(k):stop(k),k);
%		lowf(start(k):stop(k),k)=lowpass(data,0.05,2,2);
        frequency = 0.05;
        [smoothed0,icb,ice,mse0] = lowpassmin(data,frequency);
        lowf(start(k):stop(k),k)= smoothed0;
%		lowf(start(k):stop(k),k)=filtfilt(lowa,lowb,data);
		highf(start(k):stop(k),k)=data-lowf(start(k):stop(k),k);
	end
%%%%%

save('/holocene/s1/zuz10/work1/temann/newtemp/glhinfxxx','highf','-ascii')
save('/holocene/s1/zuz10/work1/temann/newtemp/gllinfxxx','lowf','-ascii')

%%%%%%%

