%band pass filter the time series using the given cutoff periods
function [ts_out]=low_filter(ts_in, fs, cut_per); 

%ts_in time series in
%fs data points per uniot of time
% cut_per cutoff period

n=4; %filter order

Wc = (1/cut_per)/(fs/2); 
[b,a] = butter(n,Wc,'low'); %defined digitial filter for the bandpass filter- this is the opposite of the filter we want
%plot the filter response
[h,w]=freqz(b,a,128);
r=abs(h); %imaginary part has 
rhigh=1-r.*r; %r r conjugate for double pass
fr=fs*w/pi/2.0;  %convert from interval 0 to pi (angular frequency) to 0 to .5
FSz=20; FN='Arial';
% figure(69); clf;
% plot(fr,r, 'b', 'LineWidth', 2.5)
% 
% xlabel('Frequency (cycle per year)', 'FontSize', FSz, 'FontName', FN)
% ylabel('Filter Response Function', 'FontSize', FSz, 'FontName', FN)
% title('Filter Response Function', 'FontSize', FSz, 'FontName', FN);
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth',2); box on;
% 
% 
 ts_out=filtfilt(b,a,ts_in);
% 
% 
% figure(2); clf;
% plot(1/fs:1/fs:length(ts_in)/fs,ts_in,'r',1/fs:1/fs:length(ts_in)/fs , ts_out, 'g')
% legend('original', 'low pass')
% xlabel('Time', 'FontSize', FSz, 'FontName', FN)
% ylabel('Value', 'FontSize', FSz, 'FontName', FN)
% title('Original and filtered time series', 'FontSize', FSz, 'FontName', FN);
% set(gca, 'FontSize', FSz, 'FontName', FN, 'LineWidth',2); box on;