% load processedRecordings
% addpath(genpath('/zpool/Dropbox/Tools'))
% addpath(genpath('/zpool/Dropbox/code'))
files = dir('2016*')
% Get a logical vector that tells which is a directory.
dirFlags = [files.isdir]
% Extract only those that are directories.
subFolders = files(dirFlags)

ls_data =[];
hpc_data =[];
processedRecordings = subFolders ;

for f = 22:length(processedRecordings)
    cd(processedRecordings(f).name)
    disp(['working on ' processedRecordings(f).name '...'])
    if ~exist('meta.mat')
        meta = generateMetaInfo('DT5');
        save meta
    end
    load meta
    
    xml = dir('*xml');
    SetCurrentSession(xml.name)
    if ~isempty(dir('*clu*'))
        units = GetUnits();
        f = find(units(:,1)<5);
        units = units(f,:);
          for n = 1:size(units,1);
        spikes{n} = GetSpikeTimes(units(n,:));
  end

    end

% lfp = GetLFP([10 93 55]);
% [b a] = butter(3,[100/625 200/625],'bandpass');
% [ripples_ls,stdev_ls,noise_ls] = FindRipples([lfp(:,1) filtfilt(b,a,lfp(:,2))],...
% 'durations',[13 100],...
% 'noise',[lfp(:,1) filtfilt(b,a,lfp(:,4))]);
% [b a] = butter(3,[130/625 200/625],'bandpass');
% [ripples_hpc,stdev_hpc,noise_hpc] = FindRipples([lfp(:,1) filtfilt(b,a,lfp(:,3))],...
% 'durations',[20 120],...
% 'noise',[lfp(:,1) filtfilt(b,a,lfp(:,4))]);
% for i=1:3
% pow(:,i)=fastrms(filtfilt(b,a,lfp(:,i+1)),40);
% end
% hold on
% for i=1:size(ripples_hpc,1)
% [pk lo] = max(pow(round((ripples_hpc(i,2)-.02)*1250):round((ripples_hpc(i,2)+.02)*1250),1));
% [pkh loh] = max(pow(round((ripples_hpc(i,2)-.02)*1250):round((ripples_hpc(i,2)+.02)*1250),2));
% scatter(pk,pkh,'.k')
% hpc_data = [hpc_data;pk,pkh];
% end
% for i=1:size(ripples_ls,1)
% [pk lo] = max(pow(round((ripples_ls(i,2)-.02)*1250):round((ripples_ls(i,2)+.02)*1250),1));
% [pkh loh] = max(pow(round((ripples_ls(i,2)-.02)*1250):round((ripples_ls(i,2)+.02)*1250),2));
% scatter(pk,pkh,'.r')
% ls_data = [ls_data;pk,pkh];
% end
% pause(.2)
% clear pow lfp ripples*

    
%       if ~isempty(dir('*pos')) 
%          lfp = GetLFP(23);
%          b = dir('*pos');
%          pos = importdata(b.name);
%          pos(isnan(pos(:,2)),:)=[];
%          startBehav = pos(1,1);
%          stopBehav = pos(end,1);
%          lfp = lfp(round(startBehav*1250):round(stopBehav*1250),:);
%          h = spectrum.periodogram;  
%          hpsd = psd(h,lfp(:,2),'Fs',1250);
%          ps = zeros(length(0:.125:200),1);
%          c=1;r
%          for hz = 0:.125:200
%          t = hz-.125 <= hpsd.Frequencies & hpsd.Frequencies<hz;
%          ps(c) = mean(log(hpsd.Data(t)));
%          c=c+1;
%          end
%          freqs = 0:.125:200;
%         save('psd_behav.mat','freqs','ps');
%       end
%       
      
%     refs = meta.refChans;
%     refs(ismember(refs,meta.rippleChans))=[];
%     r=randperm(length(refs));
%     refs = refs(r(1));
%    lfp = GetLFP([meta.rippleCha+ns r ]); % add manual checking to this section or a
% % %    tag to the saved .mat file
% % 
%    if exist('ripples_ls.rLS.evt')
%        delete('ripples_ls.rLS.evt')
% %    end
%    [b a] = butter(3,[100/625 200/625],'bandpass');
%    [ripples_ls,stdev_ls,noise_ls] = FindRipples([lfp(:,1) filtfilt(b,a,lfp(:,2))],...
%        'durations',[13 100],...
%        'noise',[lfp(:,1) filtfilt(b,a,lfp(:,4))]);
% %    SaveRippleEvents('ripples_ls.rLS.evt',ripples_ls,meta.rippleChans(1));
% %     save('ripples_ls.mat','ripples_ls','stdev_ls','noise_ls')
% %     clear ripples stdev noise
%   if exist('ripples_hpc.rHPC.evt')
%    delete('ripples_hpc.rHPC.evt')
%    end
%    [b a] = butter(3,[130/625 200/625],'bandpass');
%    [ripples_hpc,stdev_hpc,noise_hpc] = FindRipples([lfp(:,1) filtfilt(b,a,lfp(:,3))],...
%        'durations',[20 120],...
%        'noise',[lfp(:,1) filtfilt(b,a,lfp(:,4))]);
%    SaveRippleEvents('ripples_hpc.rHPC.evt',ripples_hpc,meta.rippleChans(2));
%     save('ripples_hpc.mat','ripples_hpc','stdev_hpc','noise_hpc')
%     clear ripples stdev noise
    
%% hpc theta phaselocking 
% if ~isempty(dir('*clu*')) 
%     for n = 1:length(find(units(:,1)<6))
%         spikes{n} = GetSpikeTimes(units(n,:));
%     end
%     load meta 
%     lfp = GetLFP(93); %meta.thetaChans(2));
%     if ~isempty(dir('*pos'))
%        t = dir('*pos');
%        pos = importdata(t.name);
%     end
%     [b a] = butter(3,[4/625 10/625],'bandpass');
%     filt = filtfilt(b,a,lfp(:,2));
%     phase = angle(hilbert(filt));
%     thresh = mean(fastrms(filt,200)) + std(fastrms(filt,200));
%     phase(fastrms(filt,200)<thresh)=nan;
%     phase(isnan(pos(:,2)))=nan;
%     for c = 1:length(spikes)          
%         histo{c} = phase(ceil(spikes{(c)}*1250));
%         histo{c}(isnan(histo{c})) = [];
%     end
%     save('hpc_theta_phaselocking.mat')
%     clear histo thresh phase lfp spikes
%     end

%% psd/csd 
% if ~exist('csd_psd.mat')
%     lfp = GetLFP([10 7 25 20 55 34 85 95 74 93 35 89 39]);
%     ps = NaN(length(0:.125:200),1)
%     if length(lfp(:,2)) > 1 & sum(~isnan(lfp(:,2)))~=0
%     h = spectrum.periodogram;  
%     hpsd = psd(h,lfp(:,2),'Fs',1000);
%     ps = zeros(length(0:.125:200),1);
%     c=1;
%     for hz = 0:.125:200
%     t = hz-.125 <= hpsd.Frequencies & hpsd.Frequencies<hz;
%     ps(c) = mean(log(hpsd.Data(t)));
%     c=c+1;
%     end
%     end
%     freqs = 0:.125:200;
%     
%     freqRanges = [4 9; 11 16; 35 52;60 80;85 160];
% %     for f = 1:length(freqRanges)
% %         [b a] = butter(3,[freqRanges(f,1)/625 freqRanges(f,2)/625],'bandpass');
% %         for ch = 1:size(lfp,2)-1
% %            lfp_filt(:,ch+1) = filtfilt(b,a,lfp(:,ch+1)); 
% %         end
% %         lfp_filt(:,1) = lfp(:,1);
% %         [pks locs] = findpeaks(lfp_filt(:,2),'MINPEAKHEIGHT',mean(lfp_filt(:,2))+3*std(lfp_filt(:,2)));
% %         window = round(1250/freqRanges(f,1));
% %         csd = zeros(f*2+1,size(lfp_filt,2)-2);
% %         for event = 1:length(locs)
% %             if locs(event) > window & locs(event) < length(lfp_filt)-window
% %            csd = csd + CSD(lfp_filt(locs(event)-window : locs(event)+window,:)); 
% %             end
% %         end
% %         subplot(3,2,f)
% %         csd(:,1) = 1:window*2+1;
% %         PlotCSD(csd)
% %         pause(.01)
% %         
% %         csd_all{f} = csd;
% %     end
% %     clear lfp*
% %     save('csd_psd.mat','csd_all','freqRanges','ps','freqs','-v7.3')
%     save('csd_psd.mat','freqRanges','ps','freqs','-v7.3')
% end
% % 
% % %% extract and save waveform shapes
% if ~exist('waveforms.mat') && ~isempty(dir('*clu*'))
%     count=1;
%     for shank=1:max(units(:,1))
%         cd(num2str(shank))
%         kwx{shank} = h5read([xml.name(1:end-4) '_sh' num2str(shank) ...
%             '.kwx'],['/channel_groups/' num2str(shank) '/waveforms_filtered']);
%     %     clusters = h5read([xml.name(1:end-4) '_sh'...
%     %         num2str(shank) '.kwik'],['/channel_groups/' num2str(shank) '/spikes/clusters/main']);
%     clusters = units(units(:,1)==shank,2);
%     u = unique(clusters);
%     for ii=1:length(u)
%         f = find(clusters==u(ii));
%         waveforms{shank}(ii,:,:) = squeeze(mean(kwx{shank}(:,:,f),3));
%         count=1+count;
%     end
%     cd ..
%     end
% save('waveforms.mat','waveforms','-v7.3')
% clear ii count clusters kwx f u shank waveforms
% end
% % %% phaselocking
%     if ~exist(['phaselocking_power_thresh.mat']) && ~isempty(dir('*clu*'))
%           for n = 1:size(units,1)
%         spikes{n} = GetSpikeTimes(units(n,:));
%           end
%     lfp = GetLFP([[10 7 25 20 55 34 85 95 74 93 35 89 39]]);
%         for ch = 1:size(lfp,2)
%             cells = find(units(:,1)==ch);
%             if ~isempty(cells)
%                 
%             for h = 1:100
%                [b a] = butter(3,[(h*2)/625 (h*2+4)/625],'bandpass'); 
%                lfp_filt = (filtfilt(b,a,lfp(:,ch+1)));
%                phase = angle((hilbert(lfp_filt)));
%                thresh = mean(fastrms(lfp_filt,round(1000/h))) + std(fastrms(lfp_filt,round(1000/h)));
%                phase(fastrms(lfp_filt,round(1000/h))<thresh)=nan;
%                
%                for c = 1:length(cells)          
%                        histo{ch}{c,h} = phase(ceil(spikes{cells(c)}*1250));
%                        histo{ch}{c,h}(isnan(histo{ch}{c,h})) = [];
%                        if length(histo{ch}{c,h}) > 50
%                        [phl_p{ch}(c,h) phl_z{ch}(c,h)] = circ_rtest(histo{ch}{c,h});
%                        end
%                 end
%               
% 
%             end 
%              end 
%             ch
%         end 
%         u=1;
%         for ch = 1:size(lfp,2)
%            cells = find(units(:,1)==ch);
%            if ~isempty(cells)
%            for c = 1:length(cells)    
%                for h=1:100
%                    if length(histo{ch}{c,h}) > 50
%                 spike_phase_histograms{u,h} = histo{ch}{c,h};
%                 spike_phase_pval(u,h) = phl_p{ch}(c,h);
%                 spike_phase_zscore(u,h) = phl_z{ch}(c,h);
%                    end
%                end
%                 u=1+u;
%            end
%            end
%         end
%         freqs = linspace(2,204,100);
%         save('phaselocking_power_thresh.mat','spike_phase*','freqs','-v7.3')
% 
%         clear amp* phl* histo spike_phase*
%     end
% % %     
% 
% %% all assemblies
% % if ~exist(['assemblies.mat'])
% %     lfp = GetLFP(23);
% %     if exist('behav.mat')
% %     load('behav.mat')
% %     nTrials = 0;
% %     
% %     for i=1:length(trials)
% %         for j=1:length(trials{i})
% %             nTrials = nTrials+1;
% %             all_trials{nTrials}=trials{i}{j};
% %         end
% %     end
% %     pair=[];
% %     for i=1:length(units)
% %         for j=1:length(units)
% %     %         if units(i,1) < 9 && units(j,1) > 8  % implement for cross
% %     %         regions only
% %     %             if s(i) > median(s)-median(s)/2  && s(j) > median(s)  % sets a threshold on firing rate
% %     if units(i,1) ~= units(j,1)
% %                     if i ~= j 
% %                     pair = [pair;i j];
% %                     end
% %     end
% %     %             end
% %     %         end
% %         end
% %     end
% % 
% %     sp = single(zeros(length(units),nTrials,6000));
% %     for i=1:length(units)
% %         for j=1:nTrials
% %         f = find(spikes{i}>all_trials{j}(1,5));
% %         ff = find(spikes{i}<all_trials{j}(1,5)+6);
% %         sp(i,j,ceil((spikes{i}(intersect(f,ff))-all_trials{j}(1,5))*1000))=1;
% %         end
% %     end
% %     u = unique(pair(:,1));
% %     c=1;
% %     for i=1:length(u)
% %         f = find(pair(:,1)==u(i));
% %         tic
% %         disp(['running assembly analysis on ' xml.name ' on ' ...
% %             num2str(length(f)) ' cell pairs, of ' num2str(length(pair)) ...
% %             ' total pairs'])
% %         l = size(sp,2);
% %         [d dc] = cell_assembly(sp,0:150,pair(f,:));
% % %         [dd ddc] = cell_assembly(sp,0:150,pair(f,:));
% % 
% %         for j =1:length(f)
% %             dev(:,c,:) = squeeze(mean([mean(d(:,j,:),3)],2));
% %             devControl(:,c,:) = squeeze(mean([mean(dc(:,j,:),3)],2));
% %             c=1+c;
% %         end
% %         save('assemblies.mat','dev*','pair','units')
% %         toc
% %     end
% %     clear sp
% %     save('assemblies.mat')
% %     clear dev* pair 
% %     end
% % end
%     %% head direction
% % % 
% if ~exist('hd.mat')
% if exist('behav.mat')
%     load('behav')
% % frames = GetLFP(66);
% % [pos] = alignPos(pos,frames);
% 
% angl = atan2(pos(:,4)-pos(:,2),pos(:,3)-pos(:,1));
% vel=fastrms(sum(abs(diff(pos(:,1:4))')),5);
% 
% for i=1:length(units)
%     hd{i} = [];
% end
% 
% for i=1:length(units)
%     for j=1:length(spikes{i})
%     [a b] = min(abs(spikes{i}(j)-pos(:,1)));
%     if a < .034
%         % velocity thresholdlfp(f-1250:f+1250,2)
%         if vel(b) > 1.5 && vel(b) < 60
%         hd{i} = [hd{i}; angl(b)];
%         end
%     end
%     end
% end
% 
% for i=1:length(units)
%     head_direction(i,:) = (hist([-pi; hd{i}; pi],36)./hist(angl,36));
% end
% save('hd.mat','-v7.3','head_direction','hd','angl','units')
% clear head_direction hd angl 
% end
% end
%  %% home cage vs track 
% % % 
% if ~exist(['state_rate.mat'])  && ~isempty(dir('*pos'))  && ~isempty(dir('*clu*'))
b = dir('*pos');
pos = importdata(b.name);
frames = GetLFP(1);
% [pos] = alignPos(pos,frames);
pos(isnan(pos(:,5)),:)=[];
for i =1:length(units)
   f = find(spikes{i}<pos(1,1));
   ff = find(spikes{i}>pos(end,1));
   
   rate_homecage(i) = (length(f)+length(ff))./(length(find(frames(:,1)<pos(1,1)))+...
       length(find(frames(:,1)>pos(end,1))))*1250;
   rate_track(i) = (length(spikes{i})-(length(f)+length(ff)))./...
       (length(intersect(find(frames(:,1)>pos(1,1)),...
       find(frames(:,1)<pos(end,1)))))*1250;
end

save('state_rate.mat','-v7.3','rate*')
clear rate*
% end

% %% power/rate correlations
% 
%     if ~exist(['amp_rate_corrs.mat'])
%         lfp = GetLFP([23 2 53 63 42 61 84 94 76 116 97 107 116]);
%         amp_corr_log = [];
%         amp_corr_flip_log=[];
%         amp_corr_lin = [];
%         amp_corr_flip_lin=[];
%         
%         for ch = 1:size(lfp,2)-1
%             amp_corr_lin_temp = [];
%             amp_corr_flip_lin_temp = [];
%             amp_corr_log_temp = [];
%             amp_corr_flip_log_temp = [];
%             cells = find(units(:,1)==ch);         
%                   if ~isempty(cells)
%                       prate = zeros(length(cells),length(lfp));
%                       for c =1:length(cells)
%                         prate(c,ceil(spikes{cells(c)}*1250))=1;
%                         prate(c,:) = fastrms(prate(c,:),40);
%                       end
%                    f=1;
%                    cr = 1:1000000:length(lfp)-1000001;
%                    cr = cr(randperm(length(cr)));
%                    for chunk =  1:1000000:length(lfp)-1000001
%                       [ freqs, tt, spec ] = WaveSpec(lfp(chunk:chunk+1000000,ch+1),[1 200],100,5,1/1250,'log');
%                       spec = single(abs(spec));
%                       [ freqs, tt, specControl ] = WaveSpec(lfp(cr(f):cr(f)+1000000,ch+1),[1 200],100,5,1/1250,'log');
%                       specControl = single(abs(specControl));
%                     for c = 1:length(cells)
%                       for h = 1:100
%                           temp = corrcoef(abs(spec(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_log_temp(c,h,f) = temp(2);
%                           
%                           temp = corrcoef(abs(specControl(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_flip_log_temp(c,h,f) = temp(2);
%                       end
%                     end
%                     
%                     [ freqs, tt, spec ] = WaveSpec(lfp(chunk:chunk+1000000,ch+1),[1 200],100,5,1/1250,'lin');
%                       spec = single(abs(spec));
%                       [ freqs, tt, specControl ] = WaveSpec(lfp(cr(f):cr(f)+1000000,ch+1),[1 200],100,5,1/1250,'lin');
%                       specControl = single(abs(specControl));
%                     for c = 1:length(cells)
%                       for h = 1:100
%                           temp = corrcoef(abs(spec(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_lin_temp(c,h,f) = temp(2);
%                           
%                           temp = corrcoef(abs(specControl(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_flip_lin_temp(c,h,f) = temp(2);
%                       end
%                     end
%                     f=1+f;
%                    end
%                    end 
%             amp_corr_log = [amp_corr_log; amp_corr_log_temp];
%             amp_corr_log_flip = [amp_corr_flip_log; amp_corr_flip_log_temp];
%              
%             amp_corr_lin = [amp_corr_lin; amp_corr_lin_temp];
%             amp_corr_flip_lin = [amp_corr_flip_lin; amp_corr_flip_lin_temp];
%             ch
%             
%         end
%         amp_corr_log = nanmean(amp_corr_log,3);
%         amp_corr_log_flip = nanmean(amp_corr_log_flip,3);
%         
%         amp_corr_lin = nanmean(amp_corr_lin,3);
%         amp_corr_flip_lin = nanmean(amp_corr_flip_lin,3);
% 
%         save('amp_rate_corrs.mat','amp_corr*','freqs','-v7.3')
% 
%         clear amp* 
%     end
% 
%     %% GTE
%     
% % %     
% if ~exist(['gte.mat'])
%     t = length(spikes);
%     sp = spikes;
%     sp{t+1} = [];
%     sp{t+2} = [t length(lfp)./1250];
%     [gte ci, all_te] = ASDFTE(sp,1:60);
%     clear sp t;
%     save('gte.mat','gte','ci','all_te','-v7.3')
% end
% 
% %% ripple modulation
% 
% if ~exist(['ripple_modulation.mat'])
    load('ripples_hpc.mat')
    ripple_mod = zeros(size(units,1),1251);
    [b a] = butter(3,[140/625 190/625],'bandpass'); 
    n = 1;
    for c = 1:length(units)
        for i = 1:length(ripples_hpc)
%             ll(n,:) = ll(n,:) + lfp_filt(locs_phase(i)-625:locs_phase(i)+625)';

        ff = find(spikes{n}>ripples_hpc(i,2)-.625);
        fff = find(spikes{n}<ripples_hpc(i,2)+.625);
%                 amp_corr = mean(amp_corr,3);
        ripple_mod(n,round((spikes{n}(intersect(ff,fff))-ripples_hpc(i,2))*1000)+626) = ...
            ripple_mod(n,round((spikes{n}(intersect(ff,fff))-ripples_hpc(i,2))*1000)+626)+1;
        end
        numRipples(n) = length(ripples_hpc);
        n=1+n;
    end
save('ripple_modulation.mat','ripple_mod','units','numRipples','-v7.3')
% clear avg* ll ripple* pks* locs* lfp_filt spec thresh numRipples
% end
% 
% %%     LFP phase and amplitude analysis
% clear lfp;
% lfp=GetLFP([23 97]);
%    
%  if ~exist(['lfp_phase_amplitude.mat'])
%      [b a] = butter(3,[4/625 9/625],'bandpass'); % filter bounds for theta
%      phase_amplitude = zeros(size(lfp,2)-1, 100,2501);
%      amplitude_amplitude = zeros(size(lfp,2)-1, 100,100);
%      phase_phase = zeros(size(lfp,2)-1, 99,361,361);
%      theta_wave = zeros(size(lfp,2)-1,2501);
%      
%      for ch = 1:size(lfp,2)-1       
%          for chunk = 1:1000000:length(lfp) - 999999
%              [ freqs, tt, spec ] = WaveSpec(lfp(chunk:chunk+1000000,ch+1),[1 200],100,5,1/1250,'lin');    
%              
%              % phase/amp
%              lfp_filt = (filtfilt(b,a,lfp(chunk:chunk+1000000,ch+1)));
%              [pks locs] = findpeaks(lfp_filt);
%              while locs(1) < 1251
%                  pks(1) = [];
%                  locs(1) = [];
%              end
%              while locs(end) + 1251 > length(spec)
%                  pks(end) = [];
%                  locs(end) = [];
%              end
%              for i=1:length(pks)
%                  phase_amplitude(ch,:,:) = squeeze(phase_amplitude(ch,:,:)) + abs(spec(:,locs(i)-1250:locs(i)+1250));
%                  theta_wave(ch,:) = theta_wave(ch,:) + lfp_filt(locs(i)-1250:locs(i)+1250)';
%              end
%  
%              % amp/amp
%              for f1 = 1:100
%                  for f2 = 1:100
%                      temp = corrcoef(abs(spec(f1,:)),abs(spec(f2,:)));    
%                      amplitude_amplitude(ch,f1,f2) = amplitude_amplitude(ch,f1,f2) + temp(2);
%                  end
%              end
% % 
%              % phase/phase
%              for f1 = 1:99
%                   [bb aa] = butter(3,[freqs(f1)/625 (freqs(f1+1)+1)/625],'bandpass');
%                   filt = (filtfilt(bb,aa,lfp(chunk:chunk+1000000,ch+1)));
%                   phase_theta = angle(hilbert(lfp_filt));
%                   phase_other = angle(hilbert(filt));
%                   for i = 1:length(phase_theta)
%                      t = round(rad2deg(phase_theta(i))+181);
%                      o = round(rad2deg(phase_other(i))+181);
%                      phase_phase(ch,f1,t,o) = phase_phase(ch,f1,t,o) + 1; 
%                   end
%                   ranges(f1,:) = [freqs(f1) freqs(f1+1)];
%              end
%          end            
%      end        
%   
%  save('lfp_phase_amplitude.mat','*phase*','ranges','*amplitude*','-v7.3')
%   
%  clear *phase* *amplitude* o t f1 f2 temp b a bb aa filt lfp_filt pks locs
%  end
%  
%  %% cross region analysis here
%  
% %  
% %%    event trees 
% if ~exist(['eventtree.mat'])
%     [tree counts] = EventTree(spikes,.02,3,20);
%     
%     
%     save('eventtree.mat','tree','counts','-v7.3')
%     
% end
%     
%     
%     %% ACG analysis
%     if ~exist('acg.mat')
% frames = GetLFP(66);
% 
% 
% if exist('behav.mat')
%     load('behav')
%     pos(isnan(pos(:,1)),:) = [];
%     
%         sp = single(zeros(size(units,1),ceil(length(pos)./1250*1000)));
%         for i =1:size(units,1)
%            sp(i,ceil(spikes{i}*1000)) =1;
%         end
%         for i=1:size(units,1)
%             acg(i,:) = crosscorr(sp(i,:),sp(i,:),1000);
%             disp([num2str(((i-1)*size(units,1))./size(units,1)) ' percent done...'])
%         end
%         sp(:,1:round(pos(1,1)*1000))=0;
%         sp(:,round(pos(end,1)*1000):end)=0;
%         
% 
%         for i=1:size(units,1)
%             acg_behav(i,:) = crosscorr(sp(i,:),sp(i,:),1000);
%             disp([num2str(((i-1)*size(units,1))./(size(units,1))*2) ' percent done...'])
%         end
%      
%         
%         save('acg.mat','acg*','-v7.3')
%         clear sp acg*
% end
%        end
%     %%ISI analysis
%         if ~exist('isis.mat')
% 
%     for i=1:length(units)
%         d = diff(spikes{i});
%         d(d>.5) = nan;
%         isi(i,:) = hist([d' 0],500);
%     end
%     save('isis.mat','isi');
% end
    % CCG analysis
%     if ~exist('ccg.mat')
% % lfp = GetLFP([23]);
%         sp = uint8(zeros(size(units,1),ceil(length(lfp)./1250*1000)));
%         for i =1:size(units,1)
%            sp(i,ceil(spikes{i}*1000)) =1;
%         end
%         for i=1:size(units,1)
%             for j = i:size(units,1)
%                 ccg{i,j} = [];
%             end
%         end
%         for i=1:size(units,1)
%             for j = i:size(units,1)
%                 if units(i,1) ~= units(j,1)
% %                   if i==j
%                     ccg{i,j} = crosscorr(sp(i,:),sp(j,:),250);
%                 end
%             end
%             disp([num2str(i./size(units,1)) ' percent done...'])
%         end
%         time = -250:250;
% % [times groups blah] = spikes2sorted(spikes);
% % [ccg,time] = CCG(times,groups,'binSize',.001);
%         save('ccg.mat','ccg','time','units','-v7.3')
%         clear ccg time times groups blah spikes
%     end
% %     
    cd ..
end
