%%%%%% this file is for splitting proxy series into
%%%%%% high-frequency/low-frequency bands
clear all
close all

% do high/low pass tree ring data
load /holocene/s1/zuz10/work1/temann/data/itrdbmatrix % proxy dataset
%trtestt=coldslp(:,1:10);
%save('d:\database\tempupdated\trtestt','trtestt','-ascii')
trtest=itrdbmatrix;

[n,m]=size(trtest)
%%
llow(1:n,1:m)=NaN;
hhigh(1:n,1:m)=NaN;
llow(:,1)=trtest(:,1);
hhigh(:,1)=trtest(:,1);
llow(1:3,:)=trtest(1:3,:);
hhigh(1:3,:)=trtest(1:3,:);
%%
%%---- for screening
llow1(1:n,1:m)=NaN;
hhigh1(1:n,1:m)=NaN;
llow1(:,1)=trtest(:,1);
hhigh1(:,1)=trtest(:,1);
llow1(1:3,:)=trtest(1:3,:);
hhigh1(1:3,:)=trtest(1:3,:);
%%
%% 
llow_1(1:n,1:m)=NaN;
hhigh_1(1:n,1:m)=NaN;
llow_1(:,1)=trtest(:,1);
hhigh_1(:,1)=trtest(:,1);
llow_1(1:3,:)=trtest(1:3,:);
hhigh_1(1:3,:)=trtest(1:3,:);
%%----



keepers=trtest(4:n,2:m);


% upper I add to do proxies====================

% standardize proxy based on interval 1900-1960 for the full series
[nrows,ncols]=size(keepers)




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

%load ./prepped   %load the input file. This need only have the matrix to be imputed.
%%%%%%%%%%%%%%
means=nanmean(keepers(nrows-146:nrows-1,:));
stdevs=nanstd(keepers(nrows-146:nrows-1,:));
stdkeepers=keepers-repmat(means,nrows,1);
stdkeepers=stdkeepers./repmat(stdevs,nrows,1);
%stdkeepers=stdkeepers(1:561,:); % end the reconstruction at 1960



% start with some weights
%  tricky here, a higher number is a lower weight. The actual weight in RegEM
%  is the reciprocal of the weight here. i.e. RegEM devides by the weights.

%	weight(1:length(stdkeepers))=1;  %thats the pcproxies
%	weight(1:ncols)=1;  %thats the pcproxies
%	weight=diag(weight);


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

[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
    %%
%% Now scale lowf and highf proxies only by the percentage of variance each contains
%  relative to the unfiltered series (ignore NaNs for this).
%%
	lowfpct=(nanstd(lowf).^2)./(nanstd(stdkeepers).^2);
	lowfpct=1-lowfpct;
	lowfpct(find(lowfpct<0.010))=0.01;  %can't have a zero weight
	highfpct=(nanstd(highf).^2)./nanstd(stdkeepers).^2;
	highfpct=1-highfpct;
	highfpct(find(highfpct<0.01))=0.01;  %can't have a zero weight
%%	highfpct(1:2008)=1; % keep instrumental weights=1
%%	lowfpct(1:2008)=1; % keep instrumental weights 1
%%  The weights need to be a diagonal matrix
%	highweight=diag(highfpct);
%	highweight=highweight*weight; %combine the weights
%	lowweight=diag(lowfpct);
%	lowweight=lowweight*weight;  %combine the weights
%%	save highweight highweight;  % save the highfrequency input file
%%	save lowweight lowweight;    % save the low frequency input file
%%
%%%%%%%%

%  Just use to 1980 
		%lowf=lowf(10:580,:);
		%highf=highf(10:580,:);
%		save highfinput highf
%		save lowfinput lowf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%%%% Ok, so that stuff above splits the data matrix into a high frequency
%%%% and a low frequency data matrix.  Now I've got to call two Matlabs and
%%%% send one the highf and one the lowf.
%!hybrid.pl

%=== below i add for my file

% because we need do PCs in the propick_3.m, here we normalize data
% for both high and low band series

%----------------------------------------
% standardize low band proxy based on interval 1856-1995 for the full series
%[nrows,ncols]=size(lowf)
%means=nanmean(lowf(nrows-140:nrows-1,:)); % because last yr is 1996
%stdevs=nanstd(lowf(nrows-140:nrows-1,:));
%lowf=lowf-repmat(means,nrows,1);
%lowf=lowf./repmat(stdevs,nrows,1);

% standardize high band proxy based on interval 1856-1995 for the full series
%[nrows,ncols]=size(highf)
%means=nanmean(highf(nrows-140:nrows-1,:)); % because last yr is 1996
%stdevs=nanstd(highf(nrows-140:nrows-1,:));
%highf=highf-repmat(means,nrows,1);
%highf=highf./repmat(stdevs,nrows,1);
%----------------------------------------

llow(4:n,2:m)=lowf;
hhigh(4:n,2:m)=highf;
trtest(4:n,2:m)=stdkeepers; % get normalized tree ring dataset(my form)

   
   save('/holocene/s1/zuz10/work1/temann/data/itrdbmatrixl','llow','-ascii')
   save('/holocene/s1/zuz10/work1/temann/data/itrdbmatrixh','hhigh','-ascii')

