% 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 = 3; % peak must be rmsfac*RMS(sig) to be considered a relavant peak
minsep = 1; % minimum peak separation in ms (compound PSCs within 1 ms will be considered one synaptic event)
smoothwind = 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);
    
    % 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);
        reclen = length(samp);
        
        % 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 = samp(1+(k-1)*rmswindsamp:k*rmswindsamp);
            else
                sampseg = samp(1+(k-1)*rmswindsamp:end);
            end
            minamp = rmsfac*sqrt(mean(sampseg.^2)); % adaptive minimum amplitude
            [pospks poslocs] =  findpeaks(sampseg,'minpeakheight',...
                minamp,'minpeakdistance',minsepsamp);
            [negpks neglocs] =  findpeaks(-sampseg,'minpeakheight',...
                minamp,'minpeakdistance',minsepsamp);
            locs = [poslocs neglocs];
            locs = sort(locs);
            locs(locs<negsnip) = [];
            locs(locs>(length(curr(:,i))-possnip)) = [];
            locstemp = [locstemp locs+(k-1)*rmswindsamp];
            pkstemp = [pkstemp curr(locs,i)];
        end
        
        snipmat = zeros(size(locstemp,2),possnip+negsnip+1);
        for k = 1:size(locstemp,2)
            %             if (locstemp(k) - negsnip) > 1 && (locstemp(k)+possnip) <
            %             length(curr(:,i))
            snipmat(k,:) = curr(locstemp(k) - negsnip: locstemp(k)+possnip,i);
            %             else
            %                 snipmat(k,:) = NaN(1,negsnip+possnip+1);
            %             end
            
        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

