% Jon Newman
% Ga Tech

% 13 Apr 2010
% 16 Apr 2010 - changed detrender so the detrend window size would not be
% affected by the length of the recording.

% Takes data that has already been sorted and pulled from TDT tank by
% stimulus parameters and ARASID's it. Then it appends the arasid data to
% to the file and resaves.

% Parameters regulating ARASID and data smoothing
% *************************************************************************
blockID = 'C:\TempData\20100409\20100409-1\';
dirdat = dir(strcat([blockID '*.mat']));
dtwind = 50; % size of the detrending window in ms [This is essentially a first order version of SALPA]
smoothwind = 50; % number of samples to smooth data over after ARASID
% *************************************************************************

for j = 1:size(dirdat,1)
    
    fid = dirdat(j).name;
    load(strcat([blockID fid]));
    
    % data parameters
    Fs = data.sampfreq;
    Ts = 1/Fs;
    
    % Convert detrend window to samples
    n = floor((dtwind/1000)*data.sampfreq);
    
    % Retrieve stimulus parameters
    amp = data.stimpara.amp;
    tau = data.stimpara.tau;
    
    if tau~=0
        
        % EC stim and IC Imon
        Vstim = data.stimulus.vmon;
        Imem = data.intrarec.imon;
        time = Ts*(1:size(Imem,1));
        
        % ARASID
        u = Vstim(:,1);
        y = Imem(:,1);
        arasid_data = iddata(y,u,Ts);
        artiSys = impulse(arasid_data,[0 0.02]);
        
        disp(strcat(['Running ARASID on data from: ' fid]));
        
        % Correct for direct activation/cap. artifact and detrend
        ImemSansDA = zeros(size(Imem));
        EstArt = zeros(size(Imem));
        modFits = zeros(size(Imem,2),1);
        for i = 1:size(Imem,2);
            Y = Imem(:,i);
            YH = sim(artiSys,Vstim(:,i));
            EstArt(:,i) = YH;
            breakpoints = n:n:length(Y);
            ImemSansDA(:,i) = detrend(Y-YH,'linear',breakpoints);
            modFits(i) = 100*(1-(norm(Y-YH)/norm(Y-mean(Y)))); % in %
        end
        
        % Append corrected data to data struct and resave
        data.arasid.model = artiSys;
        data.arasid.iart = ImemSansDA;
        data.arasid.estifact = EstArt;
        data.arasid.fits = modFits;
        
    else
        
        % Detrend data
        Imem = data.intrarec.imon;
        Imemdtrend =  zeros(size(Imem));
        for i = 1:size(Imem,2);
            breakpoints = n:n:length(Imem(:,i));
            samp = detrend(Imem(:,i),'linear',breakpoints);
            Imemdtrend(:,i) = samp;
        end
        
        data.arasid.model = 'NA';
        data.arasid.iart = Imemdtrend;
        data.arasid.estifact = 'NA';
        data.arasid.fits = 'NA';
    end
    
    disp(strcat(['Finshed. Re-saving file as: ' fid]));
    save(strcat([blockID fid]),'data');
    
end