% Jon Newman
% Ga Tech

% 15 Apr 2010

% Finds PSC times in data. Creates PSC snips. Appends this information to
% the data file (as a new field) and resaves

% Parameters defining PSC
% *************************************************************************
dirID = 'C:\TempData\20100409\20100409-1\';
dirdat = dir(strcat([dirID '*.mat']));
rmsfac = [4 1.5]; % 4*RMS for initial slope and 2*RMS for decay slope
minsep = 1; % minimum peak separation in ms (compound PSCs within 1 ms will be considered one synaptic event)
ppscwind = 30; % ms
smoothwind = 0.5; % milliseconds to smooth over to before doing peak detection
rmswind = 500; % adaptive RMS caclulation every rmswind ms
possnip = 200; %milliseconds after PSC peak to save in snip
negsnip = 75; %milliseconds before PSC peak to save in snip
% *************************************************************************

for j = 1:size(dirdat,1)
    
    fid = dirdat(j).name;
    load(strcat([dirID fid]));
    
    if isfield(data,'arasid') && ~isempty(data.arasid.iart)
        curr =  data.arasid.iart;
    else
        disp('Probably should do some preprosing on the data if you want decrease false positive rate')
        curr =  data.intrarec.imon;
    end
    
    % Set up some data storage structures
    psccoords = cell(1,size(curr,2));
    pscwave = cell(1,size(curr,2));
    pscplot = [];
    
    % Convert parameters to sample time
    Fs = data.sampfreq;
    minsepsamp = floor((minsep/1000)*Fs);
    rmswindsamp = floor((rmswind/1000)*Fs);
    smoothwindsamp = floor((smoothwind/1000)*Fs);
    ppscwindsamp = floor((ppscwind/1000)*Fs);
    cutsamp = 100;
    
    % Find PSC times and amplitudes
    for i = 1:size(curr,2)
        
        % Smooth recording to get rid of high freq noise peaks
        samp = smooth(curr(:,i),smoothwindsamp);
        sampdiff = diff(samp(cutsamp:end-cutsamp));
        reclen = length(sampdiff);
        
        % Set up adaptive rms calculation
        if rmswindsamp < reclen
            numwind = floor(reclen/rmswindsamp);
        else
            numwind = 1;
        end
        
        % Temp peak storage
        locstemp = [];
        pkstemp = [];
        
        % Peak detection algorithm
        for k = 1:numwind
            if k ~=numwind
                sampseg = sampdiff(1+(k-1)*rmswindsamp:k*rmswindsamp);
            else
                sampseg = sampdiff(1+(k-1)*rmswindsamp:end);
            end
            minrise = rmsfac(1)*sqrt(mean(sampseg.^2)); % adaptive minimum amplitude
            mindecay = rmsfac(2)*sqrt(mean(sampseg.^2)); % adaptive minimum amplitude
            
            
            % Step 1: Search for positive and negative peaks in
            % differential waveform that are greater than a set cuttoff
            % ampltitude (some multiple of the RMS)
            
            % IPSC search
            [pospks posinflocs] =  findpeaks(sampseg,'minpeakheight',...
                minrise,'minpeakdistance',minsepsamp);
            killind = [];
            for l = 1:length(posinflocs)
                if posinflocs(l)+ppscwindsamp < length(sampseg)
                    minvals = sort(sampseg(posinflocs(l):posinflocs(l)+ppscwindsamp));
                    if mean(minvals(1:100)) > -mindecay
                        killind  = [killind l];
                    end
                else
                    killind  = [killind l];
                end
            end
            posinflocs(killind) = [];
            pospks(killind) = [];
            
            
            % EPSC search
            [negpks neginflocs] =  findpeaks(-sampseg,'minpeakheight',...
                minrise,'minpeakdistance',minsepsamp);
            negpks = -negpks;
            killind = [];
            for l = 1:length(neginflocs)
                if neginflocs(l)+ppscwindsamp < length(sampseg)
                    minvals = sort(sampseg(neginflocs(l):neginflocs(l)+ppscwindsamp));
                    if mean(minvals(end-100:end)) < mindecay
                        killind  = [killind l];
                    end
                else
                    killind  = [killind l];
                end
            end
        end
        neginflocs(killind) = [];
        negpks(killind) = [];
        
        locs = zeros(size([posinflocs neginflocs]));
        h = 1;
        for l = posinflocs
            locs(h) =  l+find(sampdiff(l:end)<0,1)-1;
            h = h+1;
        end
        for l = neginflocs
            locs(h) =  l+find(sampdiff(l:end)>0,1)-1;
            h = h+1;
        end
        
        locs = cutsamp+sort(locs);
        
        % remove detected psc's that you can't get a snip for
        locs(locs<negsnip) = [];
        locs(locs>(length(curr(:,i))-possnip)) = [];
        locstemp = [locstemp locs+(k-1)*rmswindsamp];
        pkstemp = [pkstemp curr(locs+(k-1)*rmswindsamp,i)];
        
        
        
        snipmat = zeros(size(locstemp,2),possnip+negsnip+1);
        for k = 1:size(locstemp,2)
            snipmat(k,:) = curr(locstemp(k) - negsnip: locstemp(k)+possnip,i);
        end
        
        psccoords{i} = [locstemp' pkstemp];
        pscwave{i} = snipmat;
        pscplot = [pscplot;psccoords{i}];
        
    end
    
    figure()
    subplot(211)
    hold on
    plot(curr)
    plot(pscplot(:,1),pscplot(:,2),'k*')
    for k = 1:size(pscwave,2)
        subplot(212)
        hold all
        plot(pscwave{k}')
    end
    findfigs
    
    data.psc.coords = psccoords;
    data.psc.snips = pscwave;
    disp(strcat(['Found PSCs. Re-saving file as: ' fid]));
    save(strcat([dirID fid]),'data');
end

