function [smoothtimeseries3] = platereadersmooth(timeseries,method)
%Barry Canton
%Data smoothing.  

%create empty arrays to make the loop execute faster.  Stupid smooth
%function won't operate on multiple columns.

%timeseries is a cell array of timecourses for each of the labels used in
%the protocol
if ~isempty(method) & (method == 'fb' | method == 'ff')
    cutofffreq = 0.2;
    neighbors = [];
elseif ~isempty(method) & method == 'co'
    cutofffreq = 0.15;
    neighbors = 9;
elseif ~isempty(method) & method == 'ma'
    cutofffreq = [];
    neighbors = 5;
end


smoothtimeseries = cell(size(timeseries));
for i = 1:size(timeseries,2)
    smoothtimeseries{i} = zeros(size(timeseries{i}));
end
smoothtimeseries2 = smoothtimeseries;
smoothtimeseries3 = smoothtimeseries2;
if method == 'ma'
    for i = 1:size(timeseries{1},2)
        for j = 1:size(timeseries{1},3)
            for k = 1:length(timeseries)
                smoothtimeseries3{k}(:,i,j) = smooth(timeseries{k}(:,i,j),neighbors,'rlowess');
            end
        end
    end
elseif method == 'ff'
    [b2,a2] = butter(2,cutofffreq);

    for i = 1:size(timeseries{1},2)
        for j = 1:size(timeseries{1},3)
            for k = 1:length(timeseries)
                smoothtimeseries3{k}(:,i,j) = filter(b2,a2,timeseries{k}(:,i,j)-timeseries{k}(1,i,j));
                smoothtimeseries3{k}(:,i,j) = smoothtimeseries3{k}(:,i,j)+timeseries{k}(1,i,j);
            end
        end
    end
elseif method == 'co'
    [b2,a2] = butter(2,cutofffreq);

    for i = 1:size(timeseries{1},2)
        for j = 1:size(timeseries{1},3)
            for k = 1:length(timeseries)
                smoothtimeseries{k}(:,i,j) = smooth(timeseries{k}(:,i,j),neighbors,'moving');
                smoothtimeseries3{k}(:,i,j) = filter(b2,a2,smoothtimeseries{k}(:,i,j)-smoothtimeseries{k}(1,i,j));
                smoothtimeseries3{k}(:,i,j) = smoothtimeseries3{k}(:,i,j)+smoothtimeseries{k}(1,i,j);
            end
        end
    end
% elseif method == 'fb'
%     [b2,a2] = butter(2,0.08);
% 
%     for i = 1:size(timeseries{1},2)
%         for j = 1:size(timeseries{1},3)
%             for k = 1:length(timeseries)
%                 smoothtimeseries2{k}(end:-1:1,i,j) = filter(b2,a2,timeseries{k}(:,i,j)-timeseries{k}(1,i,j));
%                 smoothtimeseries2{k}(:,i,j) = smoothtimeseries2{k}(:,i,j)+timeseries{k}(1,i,j);
%             end
%         end
%     end
%     for i = 1:size(timeseries{1},2)
%         for j = 1:size(timeseries{1},3)
%             for k = 1:length(timeseries)
%                 smoothtimeseries3{k}(end:-1:1,i,j) = filter(b2,a2,smoothtimeseries2{k}(:,i,j)-smoothtimeseries2{k}(1,i,j));
%                 smoothtimeseries3{k}(:,i,j) = smoothtimeseries3{k}(:,i,j)+smoothtimeseries2{k}(1,i,j);
%             end
%         end
%     end
elseif method == 'fb'
    [b2,a2] = butter(2,cutofffreq);
%     [b2,a2] = cheby1(2,0.5,cutofffreq,'low')
    for k = 1:length(timeseries)
        firstpoint = repmat(timeseries{k}(1,:,:),[size(timeseries{k},1),1,1]);
        smoothtimeseries2{k}(end:-1:1,:,:) = filter(b2,a2,timeseries{k}-firstpoint);
        smoothtimeseries2{k} = smoothtimeseries2{k}+firstpoint;
   
        firstpoint = repmat(smoothtimeseries2{k}(1,:,:),[size(smoothtimeseries2{k},1),1,1]);
        smoothtimeseries3{k}(end:-1:1,:,:) = filter(b2,a2,smoothtimeseries2{k}-firstpoint);
        smoothtimeseries3{k} = smoothtimeseries3{k}+firstpoint;
    end
    
end