% Jon Newman
% Ga Tech

% 13 Apr 2010

% 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\20100420\20100420-35\';
dirdat = dir(strcat([blockID '*.mat']));
smoothwind = 50; % number of samples to smooth data over after ARASID
percentcut = 1; % percent of record to cut off begining and end (gets rid of possible non-perfectly synced stim transitions)
% *************************************************************************

for j = 1:size(dirdat,1)
    
    fid = dirdat(j).name;
    load(strcat([blockID fid]));
    
    % data parameters
    Fs = data.sampfreq;
    Ts = 1/Fs;
    
    % 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
        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 = floor(length(Y)/10*(1:10));
            ImemSansDA(:,i) = detrend(Y-YH,'linear',breakpoints);
            modFits(i) = 100*(1-(norm(Y-YH)/norm(Y-mean(Y)))); % in %
        end
        
        % Smooth and trim data
        ImemSansDAsmooth = [];
        for i = 1:size(Imem,2);
            samp = smooth(ImemSansDA(:,i),smoothwind);
            reclen = length(samp);
            numcut = floor((percentcut/100)*reclen);
            samp = samp(numcut:end-numcut);
            ImemSansDAsmooth = [ImemSansDAsmooth, samp] ;
        end
        
        % Append corrected data to data struct and resave
        data.arasid.model = artiSys;
        data.arasid.iart = ImemSansDAsmooth;
        data.arasid.estifact = EstArt;
        data.arasid.fits = modFits;
        data.arasid.smoothwind = smoothwind;
        data.arasid.percentcut = percentcut;
        
    else
        
        % Smooth and trim data
        Imem = data.intrarec.imon;
        Imemsmooth = [];
        for i = 1:size(Imem,2);
            samp = smooth(Imem(:,i),smoothwind);
            reclen = length(samp);
            numcut = floor((percentcut/100)*reclen);
            samp = samp(numcut:end-numcut);
            Imemsmooth = [Imemsmooth, samp] ;
        end
        
        data.arasid.model = 'NA';
        data.arasid.iart = Imemsmooth;
        data.arasid.estifact = 'NA';
        data.arasid.fits = 'NA';
        data.arasid.smoothwind = smoothwind;
        data.arasid.percentcut = percentcut;
    end
    
    disp(strcat(['Finshed. Re-saving file as: ' fid]));
    save(strcat([blockID fid]),'data');
    
end

beep; beep; beep;