

files = [{'DT1_rLS_20150721_720um'},...
         {'DT1_rLS_20150721_1008um'},...
         {'DT1_rLS_20150722_1152um'},...
         {'DT1_rLS_20150722_1440um'},...
         {'DT1_rLS_20150723_1584um'},...
         {'DT1_rLS_20150723_1728um'},...
         {'DT1_rLS_20150724_1800um'},...
         {'DT1_rLS_20150724_1872um'},...
         {'DT1_rLS_20150725_1944um'},...         
         {'DT1_rLS_20150726_2016um'},... 
         {'DT1_rLS_20150727_2088um'},...
         {'DT1_rLS_20150727_2232um'},...
         {'DT1_rLS_20150728_2304um'},...
         {'DT1_rLS_20150729_2448um'},...
         {'DT1_rLS_20150729_2592um'},...
         {'DT1_rLS_20150730_2736um'}];

for f = 1:length(files)
    cd(files{f})
%     xml = dir('*xml');
    SetCurrentSession
    units = GetUnits();
    for n = 1:size(units,1);
        spikes{n} = GetSpikeTimes(units(n,:));
    end
%     if exist([files{f} '-states.mat'])
%         load([files{f} '-states.mat'])
%     end
%   if ~exist('meta.mat')
%         meta = generateMetaInfo('DT5');
%         save meta
%   else
%   load meta
%   end
%     
%     lfp = GetLFP([2]);
%     lfp = GetLFP([2 5 12 15 18 23 28 29]);
       amp_corr =[];
    amp_corr_flip = [];
    
    


%     h = spectrum.periodogram;  
%     hpsd = psd(h,lfp(:,2),'Fs',1250);
%     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
%     freqs = 0:.125:200;
%     save('psd.mat','freqs','ps');


    %% waveforms
%     if ~exist('waveforms.mat')
c=1;
for t = 1:units(end,1)
    u = find(units(:,1)==t);
    for i=1:length(u)
        temp = GetSpikeWaveforms(units(u(i),:));
        waveforms{t}(i,:,:)=squeeze(mean(temp));
    end
    if ~isempty(u)
    getfeatures = LoadSpikeFeatures([files{f} '.clu.' num2str(t)],20000);
    [L(c), Lratio(c), df(c)] = L_Ratio(getfeatures(:,4:end), find(getfeatures(:,3)==units(u(i),2)));
    IsolDist(c) = IsolationDistance(getfeatures(:,4:end), find(getfeatures(:,3)==units(u(i),2)));
    end
end
save('waveforms.mat','waveforms','L','Lratio','df','IsolDist','units','-v7.3')
clear ii count clusters kwx f u shank waveforms



%     end

%     
%     if ~exist('ripples_ls.rLS.evt')
% %        delete('ripples_ls.rLS.evt')
%  
%    [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,2);
%     save('ripples_ls.mat','ripples_ls','stdev_ls','noise_ls')
%     clear ripples stdev noise
%     end

%   %% phaselocking
%     if ~exist(['phaselocking_power_thresh.mat'])
%         
%         n=1;
%         for ch = 1:size(lfp,2)-1
%             cells = find(units(:,1)==ch);
%             for c = 1:length(cells)
%                 
%             for h = 1:100
%                [b a] = butter(3,[h/625 (h*2+4)/625],'bandpass'); 
%                lfp_filt = (filtfilt(b,a,lfp(:,ch+1)));
%                phase = angle(abs(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;
%                 
% 
%                    if ~isempty(cells)
%     %                 prate = zeros(length(lfp),1);
%     %                 prate(round(spikes{cells(c)}*1000))=1;
%     %                 prate = fastrms(prate,40);
% 
%                        histo{n,h} = phase(round(spikes{n}*1250));
%                        histo{n,h}(isnan(histo{n,h})) = [];
%                        if length(histo{n,h}) > 50
%                        [phl_p(n,h) phl_z(n,h)] = circ_rtest(histo{n,h});
%                        end
%     %                f=1;
%     %                for chunk = 1:1000000:length(lfp)-9999999
%     %                   cc = cwt(lfp(chunk:chunk+100000,ch+1),freq,'cgau4');
%     %                   S = abs(cc.*cc); cc=[];
%     %                   temp = corrcoef((S(h,:)),prate(chunk:chunk+100000));
%     %                   amp_corr(n,h,f) = temp(2);
%     %                   f=1+f;
%     %                end
%                    if h == 100
%                        n=1+n;
%                    end
%                    end
%                    freqRange(h,:)=[h h*2+4];
%                 end 
%             end 
%             ch
%         end
%         amp_corr = mean(amp_corr,3);
% 
%         save('phaselocking_power_thresh.mat','histo','phl*','freqRange','-v7.3')
% 
%         clear amp* phl* histo
%     end
%% power/rate correlations
%     if ~exist(['amp_rate_corrs.mat'])
%         for ch = 1:size(lfp,2)-1
%             amp_corr_temp = [];
%             amp_corr_flip_temp = [];
%             cells = find(units(:,1)==ch);         
%                   if ~isempty(cells)
%                       prate = zeros(length(cells),length(lfp));
%                       for c =1:length(cells)
%                         prate(c,round(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_temp(c,h,f) = temp(2);
%                           
%                           temp = corrcoef(abs(specControl(h,200:1000000-200)),prate(c,chunk+200:chunk+1000000-200));
%                           amp_corr_flip_temp(c,h,f) = temp(2);
%                           
%                       end
%                     end
%                     f=1+f;
%                    end
%                    end 
%             amp_corr = [amp_corr; amp_corr_temp];
%             amp_corr_flip = [amp_corr_flip; amp_corr_flip_temp];
%             ch
%             
%         end
%         amp_corr = nanmean(amp_corr,3);
%         amp_corr_flip = nanmean(amp_corr_flip,3);
%         amp_corr_flip_std = nanstd(amp_corr_flip_temp,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:40);
%     clear sp t;
%     save('gte.mat','gte','ci','all_te','-v7.3')
% end
% 
% %% ripple modulation

% if ~exist(['ripple_modulation.mat'])
%     
%     histo_phase = zeros(size(units,1),1251);
%     ll = zeros(size(units,1),1251);
%     [b a] = butter(3,[140/625 190/625],'bandpass'); 
%     n = 1;
%     for ch = 1:size(lfp,2)-1
%         cells = find(units(:,1)==ch);
%         if ~isempty(cells)
%         [ freqs, tt, spec ] = WaveSpec(lfp(:,ch+1),[140 190],10,5,1/1250,'lin');
%         spec = mean(abs(spec));
%         thresh = mean(abs(spec))+2*std(abs(spec));
%         
%         lfp_filt = (filtfilt(b,a,lfp(:,ch+1)));
%         lfp_filt(spec<thresh) = 0;
%         [pks_phase locs_phase] = findpeaks(lfp_filt,'MINPEAKHEIGHT',1);
%         while locs_phase(1) < 1250
%             pks_phase(1) = [];
%             locs_phase(1) = [];
%         end
%         while locs_phase(end) + 1250 > length(lfp)
%             pks_phase(end) = [];
%             locs_phase(end) = [];
%         end
% 
%         
%         for c = 1:length(cells)
%            if ~isempty(cells)
%             for i = 1:length(pks_phase)
%             ll(n,:) = ll(n,:) + lfp_filt(locs_phase(i)-625:locs_phase(i)+625)';
%             
%             ff = find(spikes{n}*1250>locs_phase(i)-625);
%             fff = find(spikes{n}*1250<locs_phase(i)+625);
%             
%             histo_phase(n,round(spikes{n}(intersect(ff,fff))*1250)-locs_phase(i)+626) = ...
%                 histo_phase(n,round(spikes{n}(intersect(ff,fff))*1250)-locs_phase(i)+626)+1;
%             end
%            end 
%             n=1+n;
%         end
%         end
%     end
% 
% avg_filtered_lfp = ll;
% save('ripple_modulation.mat','histo_phase','avg_filtered_lfp','-v7.3')
% clear avg* ll histo* pks* locs* lfp_filt spec thresh 
% end

%     % LFP phase and amplitude analysis
% if ~exist(['lfp_phase_amplitude.mat'])
%     [b a] = butter(3,[4/625 9/625],'bandpass'); % filter bounds for theta
%     phase_amplitude = zeros(100,2501);
%     amplitude_amplitude = zeros(100,100);
%     phase_phase = zeros(99,361,361);
%     theta_wave = zeros(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 = phase_amplitude + abs(spec(:,locs(i)-1250:locs(i)+1250));
%                 theta_wave = theta_wave + 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(f1,f2) = amplitude_amplitude(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(f1,t,o) = phase_phase(f1,t,o) + 1; 
%                  end
%                  ranges(f1,:) = [freqs(f1) freqs(f1+1)];
%             end
%         end            
%     end        
%  
% save('lfp_phase_amplitude.mat','*phase*','ranges','*amplitude*','theta_wave','-v7.3')
% end  
% clear *phase* *amplitude* o t f1 f2 temp b a bb aa filt lfp_filt pks locs

    %% event trees 
% if ~exist(['eventtree.mat'])
%     [tree counts] = EventTree(spikes,.03,6,6);
% %     
% %     
%     save('eventtree.mat','tree','counts','-v7.3')
% %     
% end
    

%% CCG 
%     if ~exist('ccg.mat')
%     sp = single(zeros(size(units,1),length(lfp)./1250*1000));
%     for i =1:size(units,1)
%        sp(i,round(spikes{i}*1000)) =1;
%     end
%     for i=1:size(units,1)
%         for j = i:size(units,1)
%             ccg{i,j} = [];
% %             if units(i,1) ~= units(j,1) || i==j
%                 ccg{i,j} = crosscorr(sp(i,:),sp(j,:),2000);
% %             end
%         end
%         disp([num2str(((i-1)*size(units,1))./size(units,1)) ' percent done...'])
%     end
%     save('ccg.mat','ccg','-v7.3')
%     clear sp ccg
%     end
    cd ..
end