function [x_butter_ann, x_butter_bia, x_butter_mos, x_lp_ann, x_lp_mos, x_prolate_ann, x_prolate_bia, x_prolate_mos, x_rmn_ann, x_rmn_bia, x_rmn_mos] = make_filtered_monthly_series (x, x_date)
%
% Assumes input is a monthly time series.
% Returns a "seasonal", annual and biannual filter of the data using running 
% means, prolate filters, butterworth filters etc.
%
% Best if there are no NaNs in the data.
%
% Names are assuming the data input is daily values.
%
% Assumes that you are passing a column vector of data

%
% Make sure you have columns of data
%
x = x(:);
x_date = x_date(:);


%
% Make a temporary variable that flips on the edges of that series
% so that no points are lost.
%
temp=[flipud(x); x; flipud(x);];


%
% Butterworth filtering
%
% Uses routine:  [B,A,FLT]=butter_filter(xx,nf,n,wn)
% where xx is the input series and FLT is the output filtered series
%
% nf =number of frequencies in filter-controls smoothness of
%                   frequency response plots
%   = best values on the order of 25*nyq/wn
% n = order of filter
% wn=low pass frequency cutoff
%   =0.04 is the half power point of the filter
%   = units should be cycles/user_units where user unit is time
%         interval of cable_t_day
% e.g. 40 hr half-pt using hourly data
% wn=1/40=0.025, nyq=0.5, nf=25*.5/.025=500
% our case:
% daily data, nyq=0.5 days,
% for monthly smoothing half power at about 25 days wn=1/25=0.04
%     nf=312.5 (I used twice that because I don't understand why
%     the nyquist isn't 1 day, not 0.5 days as Dennis insisted.
% for annual smoothing with a half power point at 300 days
%%%%%%%     wn=1/300=0.003333, nf=25*.5/.0033=3750
%     wn=1/330=0.0030303, nf=25*.5/.0033=4125
% for biannual smoothing (half power at 600 days)
%%%%%%%     wn=1/600=0.001666, nf=7500
%     wn=1/630=0.0015873, nf=7875
%
%
[m,n]=size(x);
%
% First a monthly filter
%
% Monthly data use:  
[a,b,temp_filt]=butter_filter(temp,100,3,0.5);
% Daily Data use:
%       [a,b,temp_filt]=butter_filter(temp,625,3,0.04);
x_butter_mos=temp_filt(m+1:m+m);
%
% Then an annual filter
%
% Monthly data use:
%[a,b,temp_filt]=butter_filter(temp,350,3,0.08333);
[a,b,temp_filt]=butter_filter(temp,63,3,0.2);
% Daily data use:
%       [a,b,temp_filt]=butter_filter(temp,3750,3,0.003333);
x_butter_ann=temp_filt(m+1:m+m);
%
% Then a biannual filter
%
% Monthly data use:
%[a,b,temp_filt]=butter_filter(temp,700,3,0.04166);
[a,b,temp_filt]=butter_filter(temp,113,3,1.1111e-01);
% Daily data use:
%        [a,b,temp_filt]=butter_filter(temp,7500,3,0.001666);
x_butter_bia=temp_filt(m+1:m+m);

%
% Check filtered versions of the array
%
%check_corr(x_butter_mos, x,['Butter Monthly vs. original']);
%check_corr(x_butter_ann, x,['Butter Annual vs. original']);
%check_corr(x_butter_bia, x,['Butter Biannual vs. original']);
check_corr4(x_butter_mos, x,['Butter Monthly vs. original'], 'Months', 30, 1);
check_corr4(x_butter_ann, x,['Butter Annual vs. original'], 'Months', 30, 1);
check_corr4(x_butter_bia, x,['Butter Biannual vs. original'], 'Months', 30, 1);


%
% Plot all the Butterworth Filtered series together
%
figure;
  plot(x_date, x);
  hold on;
  datetick;
  h=plot(x_date, x_butter_mos,'r');
  set(h,'LineWidth',3);
  h=plot(x_date, x_butter_ann,'g');
  set(h,'LineWidth',3);
  h=plot(x_date, x_butter_bia,'m');
  set(h,'LineWidth',3);
  xlabel('Time, years');
  ylabel('Original Data units')  ;
  title('Butterworth filtering');
  legend('Orig','Mos','Ann','Bi-Ann');


%
% Now make some running means
%  Use routine:
%       SMTH = RMN(dat,len,rep)
%  where dat is the input and SMTH is the filtered output.
%    Creates a rectangular filter of length <len>, then smooths
%    <dat> by convolution <rep> times. If the total equivalent
%    filter length is n, then the (n-1)/2 points at the start
%    and end of the series are lost in the filtrand (smth) and
%    are replaced by NaN's.
%
% Assmuning you have daily data, the filter widths should be
%    monthly = 31
%    annual  = 365
%    biannual = 731
% All values are odd so that the filter is symetrical
%
% Assuming you have monthly data, filters something like
%    monthly = 3
%    annual  = 13
%    biannual = 25
% What impact is there that these are *not* odd?
%
[m,n]=size(x);
temp_filt = rmn(temp,3,1);
x_rmn_mos = temp_filt(m+1:m+m);

temp_filt = rmn(temp,13,1);
x_rmn_ann = temp_filt(m+1:m+m);

temp_filt = rmn(temp,25,1);
x_rmn_bia = temp_filt(m+1:m+m);

%
% Check filtered versions of the array
%
%check_corr(x_rmn_mos,x,['Running Mean Monthly vs. original']);
%check_corr(x_rmn_ann,x,['Running Mean Annual vs. original']);
%check_corr(x_rmn_bia,x,['Running Mean Biannual vs. original']);
check_corr4(x_rmn_mos,x,['Running Mean Monthly vs. original'], 'Months', 30, 1);
check_corr4(x_rmn_ann,x,['Running Mean Annual vs. original'], 'Months', 30, 1);
check_corr4(x_rmn_bia,x,['Running Mean Biannual vs. original'], 'Months', 30, 1);

%
% Plot all the running mean series together
%
figure;
  plot(x_date, x);
  hold on;
  datetick;
  h=plot(x_date,x_rmn_mos,'r');
  set(h,'LineWidth',3);
  h=plot(x_date,x_rmn_ann,'g');
  set(h,'LineWidth',3);
  h=plot(x_date,x_rmn_bia,'m');
  set(h,'LineWidth',3);
  xlabel('Time, years');
  ylabel('Original Data units')  ;
  title('Running Mean filtering');
  legend('Orig','Mos','Ann','Bi-Ann');


%
% Now use a prolate filter
%
% biannual prolate filter *Very* much like annual
% butterworth filter except amplitude of prolate is smaller.
%
% Use routine:
%
% For daily data filters should be:
%   monthly = 31
%   annual  = 365
%   biannual = 731
% For monthly values
%   monthly = 5
%   annual  = 13
%   biannual = 25
%
[m,n]=size(x);
temp_filt = prolate_filt(temp,5,1);
x_prolate_mos = temp_filt(m+1:m+m);

temp_filt = prolate_filt(temp,13,1);
x_prolate_ann = temp_filt(m+1:m+m);

temp_filt = prolate_filt(temp,25,1);
x_prolate_bia = temp_filt(m+1:m+m);

%
% Check filtered versions of the array
%
%check_corr(x_prolate_mos,x,['Prolate Monthly vs. original']);
%check_corr(x_prolate_ann,x,['Prolate Annual vs. original']);
%check_corr(x_prolate_bia,x,['Prolate Biannual vs. original']);
check_corr4(x_prolate_mos,x,['Prolate Monthly vs. original'], 'Months', 30, 1);
check_corr4(x_prolate_ann,x,['Prolate Annual vs. original'], 'Months', 30, 1);
check_corr4(x_prolate_bia,x,['Prolate Biannual vs. original'], 'Months', 30, 1);

%
% Plot all the prolate series together
%
figure;
  plot(x_date,x);
  hold on;
  datetick;
  h=plot(x_date,x_prolate_mos,'r');
  set(h,'LineWidth',3);
  h=plot(x_date,x_prolate_ann,'g');
  set(h,'LineWidth',3);
  h=plot(x_date,x_prolate_bia,'m');
  set(h,'LineWidth',3);
  xlabel('Time, years');
  ylabel('Original Data units')  ;
  title('Prolate filtering');
  legend('Orig','Mos','Ann','Bi-Ann');



%
% Try Dennis's Lancho filter
%
% By definition Dennis's Lancho routine assume that the series is monthly.
%

[m,n]=size(x);

temp_filt = lpfilt(temp,5,1);

x_lp_mos = temp_filt(m+1:m+m);

temp_filt = lpfilt(temp,18,1);
x_lp_ann = temp_filt(m+1:m+m);

%
% Check filtered versions of the array
%
%check_corr(x_lp_mos, x,['Lancho Monthly vs. original']);
%check_corr(x_lp_ann, x,['Lancho Annual vs. original']);
check_corr4(x_lp_mos, x,['Lancho Monthly vs. original'], 'Months', 30, 1);
check_corr4(x_lp_ann, x,['Lancho Annual vs. original'], 'Months', 30, 1);

%
% Plot all the Lancho series together
%
figure;
  plot(x_date, x);
  hold on;
  datetick;
  h=plot(x_date, x_lp_mos,'r');
  set(h,'LineWidth',3);
  h=plot(x_date, x_lp_ann,'g');
  set(h,'LineWidth',3);
  xlabel('Time, years');
  ylabel('Original Data units')  ;
  title('Lancho filtering');
  legend('Orig','Mos','Ann');

return;


