%Analysis for Matt Kittelberger

%A 'b' in the file name means single unit
%--no b means whole nerve recording

%Data selection:
%--for vocal nerve, discard spikes before stim end
%--for single unit, subtract spikes before stim start
%--stim length is (Npulse-1)*3+1 mSec (col 3 on spread sheet)
%--time before stim is always 200 mSec

%For Vocal nerve:
% #spikes/sweep, average #/set, average #/group
% latency to onset (by sweep, set, stim group)
% ISI -> freq plot and mean ISI
% IBI and # of Bursts (threshold for burst start 50 mSec)
% #spikes/burst
%
%For single unit:
% #spikes/sweep, average #/set, average #/group
% subtract off the numbers from the zero stim sets
%separate PST section
%====================================================================

clear all

path = 'C:\Documents and Settings\bruce land\My Documents\Matlab\Kittelberger\Analysis\';
filename= 'pn04078.062304b.002timesCluster1.xls';
%filename = 'pn04165.092904b.002timesCluster1.xls';

%[filename,path] = uigetfile('*.txt','Choose a data file');
%check filename return
if (filename==0)
    error('invalid filename')
end

%cut off '.xls'
filebase=filename(1:strfind(filename,'.xls')-1);
fname  =  ([path,filename]);

allData = xlsread (fname)  ;
Ntrials = size(allData,1);

%vocal or single unit
isSingle=1;
%burst threshold for single unit
burstThreshold = 500; %mSec
%histogram for singleunit
histcenters = [0:49,50:10:500];
trialhist = zeros(Ntrials,length(histcenters));
PSTcenters = [0:5:1000];
PSTcentersAligned = [-300:5:700];
PSThist = zeros(Ntrials,length(PSTcenters));
PSThistAligned = zeros(Ntrials,length(PSTcentersAligned));
AllSpikes = cell(Ntrials,1);
GoodSpikes  = cell(Ntrials,1);

if (isempty(strfind(filename,'b')))
    %if the data is a vocal burst 
    %set the parameters 
    isSingle=0;
    %burst threshold for vocal 
    burstThreshold = 50; %mSec
    %histogram for vocal
    trialhist = zeros(Ntrials,50);
end
MaxBursts = 4; %arbitrary maximum for an array

if isSingle
    %now read -- JUST the latencies from the vocal response sweeps
    %This will be read in again in the sigle unit pass
    bloc = strfind(filename,'b');
    filebaseb = filebase([1:bloc-1,bloc+1:end]);
    filein = [path,filebaseb,'VocalLatency','.txt'];
    lat = load(filein);
    for trial=1:Ntrials
        VLatencyBegin(trial)=lat(trial,1);
        VLatencyEnd(trial)=lat(trial,2);
    end
end

%first pass analysis -- same for single or vocal
for trial=1:Ntrials
    %find end of stimulus
    if allData(trial,3)>0
        StimEnd = 200 + (allData(trial,3)-1)*3+1;
    else %there was no stimulus
        StimEnd = 0;
    end
    TrialDuration(trial) = 1000-StimEnd; %mSec
    %get all the spikes
    AllSpikes{trial}= ...
        allData(trial,5:end);
    %find spikes after end of stim
    GoodSpikes{trial} = ...
        allData(trial,4+(find(allData(trial,5:end)>StimEnd)));
    
    %count the good spikes
    Nspikes(trial) = length(GoodSpikes{trial});
    
    %get the PST 
    if (length(AllSpikes{trial})>0) 
        PSThist(trial,:) = ...
            histc(AllSpikes{trial},PSTcenters);
    end
    
    %get the latency to the first spike
    %need at least one spike
    if (length(GoodSpikes{trial})>0)
        %latency from end of stim
        LatencyEnd(trial) = GoodSpikes{trial}(1)-StimEnd;
        %latency from beginning of stim
        LatencyBegin(trial) = GoodSpikes{trial}(1)-200;
    else
        LatencyEnd(trial) = nan;
        LatencyBegin(trial) = nan;
    end
    
    if isSingle
        %get the ALIGNED PST
        if (length(AllSpikes{trial})>0) & (~isnan(VLatencyBegin(trial)))
            PSThistAligned(trial,:) = ...
                histc(AllSpikes{trial}-VLatencyBegin(trial)-200,...
                PSTcentersAligned);
        else
            PSThistAligned(trial,:) = 0;
        end
        %get latency difference (su-vocal)
        if (length(AllSpikes{trial})>0) & (~isnan(LatencyBegin(trial)))
            LatencyDiff(trial) = ...
                LatencyBegin(trial)-VLatencyBegin(trial);
        else
            LatencyDiff(trial) = nan;
        end
    end
    
    % mean ISI < burstThreshold
    %need at least two spikes
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        diffspike(diffspike>burstThreshold) = nan;
        trialISI{trial} = diffspike ; %list for later averaging
        meanISI(trial) = nanmean(diffspike);
        stdISI(trial) = nanstd(diffspike);
    else
        meanISI(trial) = nan;
        stdISI(trial) = nan ;
        trialISI{trial} = nan; 
    end

    % average burst interval--need at least 2 spikes
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        %eliminate all intevals less than burst threshold
        diffspike(diffspike<=burstThreshold) = nan;
        meanIBI(trial) = nanmean(diffspike);
        stdIBI(trial) = nanstd(diffspike);
    else
        meanIBI(trial) = nan;
        stdIBI(trial) = nan ;
    end

    % number of bursts--need at least 2 spikes
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        %eliminate all intevals less than burst threshold
        diffspike(diffspike<=burstThreshold) = nan;
        numberBursts(trial) = sum(~isnan(diffspike)) + 1 ;
    else
        numberBursts(trial) = 0;
    end

    % # spikes/burst--need at least 2 spikes
    spikeburst(trial,1:MaxBursts)=0;
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        %find all intevals greater than burst threshold
        diffspike(diffspike>burstThreshold) = nan;
        diffspike = [0 diffspike];
        burstcounter=1;
        for i=1:length(GoodSpikes{trial})
            if (isnan(diffspike(i)))
                burstcounter = burstcounter + 1;
                %spikeburst{trial}(burstcounter)=0;
            end
            %ignore multiple bursts above maxBursts
            if burstcounter<=MaxBursts
                spikeburst(trial,burstcounter) = ...
                    spikeburst(trial,burstcounter) + 1;
            end
        end
    else
        %spikeburst{trial} = nan;
    end

    % burst duration--need at least 2 spikes
    burstdur(trial)=nan;
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        %find all intevals greater than burst threshold
        diffspike(diffspike>burstThreshold) = nan;
        diffspike = [0 diffspike];
        if sum(isnan(diffspike))==0 %more than one burst?
            burstdur(trial) = ...
                GoodSpikes{trial}(end)-GoodSpikes{trial}(1);
        else %if yes
            lastone = find(isnan(diffspike)==1);
            lastone = lastone(1)-1;
            burstdur(trial) = ...
                GoodSpikes{trial}(lastone)-GoodSpikes{trial}(1);
        end
    end

    %histogram of Goodspikes
    if (length(GoodSpikes{trial})>1)
        diffspike = diff(GoodSpikes{trial});
        %eliminate all intevals less than burst threshold
        diffspike(diffspike>burstThreshold) = nan;
        if isSingle
            trialhist(trial,:) = histc(diffspike,histcenters);
        else
            trialhist(trial,:) = histc(diffspike,0:49);
        end
        
    end
end

if ~isSingle
    %now make an results file -- JUST the latencies
    %This will be read in again in the sigle unit pass
    fileout = [filebase,'VocalLatency','.txt'];
    fileout = ([path,fileout]);
    fid = fopen(fileout,'w');
    for trial=1:Ntrials
        fprintf(fid,'%8.2f\t%8.2f\r', ...
            LatencyBegin(trial), LatencyEnd(trial))
    end
    fprintf(fid,'\r');
    fclose(fid)
end
%=================================
%Vocal nerve response
%=================================
if isSingle==0
    %now make an results file
    fileout = [filebase,'Output','.txt'];
    fileout = ([path,fileout]);
    fid = fopen(fileout,'w');

    %=================================
    %data for individual sweeps
    fprintf(fid,'Data by SWEEP: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tSwp\tStim\tTime\t#Spk\tLat(End)\tLat(Begin)\t',...
        'meanISI\tstdISI\tmeanIBI\tstdIBI\t',...
        '#Brst \tDur1\t#Spk/b1\t#Spk/b2\t#Spk/b3\t#Spk/b4',...
        '\tISIhist\n\r'])
    for trial=1:Ntrials
        %print the set,sweep and stim
        fprintf(fid,'%d\t%d\t%d\t%d\t', ...
            allData(trial,1),allData(trial,2),allData(trial,3),allData(trial,4));
        %print the Nspikes, Latency, meanISI
        fprintf(fid,'%d\t%6.1f\t%6.1f\t%6.1f\t%6.1f\t', ...
            Nspikes(trial), LatencyEnd(trial),LatencyBegin(trial),...
            meanISI(trial),stdISI(trial));
        %print the meanIBI, numberBursts
        fprintf(fid,'%6.1f\t%6.1f\t%4d\t', ...
            meanIBI(trial),stdIBI(trial), numberBursts(trial));
        %burstduration
        fprintf(fid,'%6.1f\t',burstdur(trial));
        %#spike/burst
        for i=1:MaxBursts
            fprintf(fid,'%6.2f\t', spikeburst(trial,i));
        end
        %histogram of GoodSpikes
        for i=1:50
            fprintf(fid,'%3.0f\t', trialhist(trial,i))
        end
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');
    %====================================
    %data clumped by Set
    MaxSet = max(allData(:,1));
    fprintf(fid,'Data by SET: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tStim\tTime\t#Spike\t#S-std\t',...
        ' Lat(E)\tL(E)std\t',...
        ' Lat(B)\tL(B)std\t',...
        ' ISI  \tISIstd\t',...
        '  IBI  \tIBIstd\t',...
        '  #Burst\t#Bstd \tDur1\tsDur1\t'...
        '#Spk/b1\t#Spk/b2\t#Spk/b3\t#Spk/b4',...
        '\tISIhist\r'])

    for set=1:MaxSet
        trial = find(allData(:,1)==set);
        fprintf(fid,'%d\t%d\t%d\t', ...
            set,allData(trial(1),3),allData(trial(1),4));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(Nspikes(trial)), nanstd(Nspikes(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
        %make a vector of all the ISIs
        allISI=[];
        for i=trial
            allISI = [allISI, trialISI{i}];
        end
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(allISI),nanstd(allISI));       
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(meanIBI(trial)),nanstd(meanIBI(trial)));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(numberBursts(trial)),nanstd(numberBursts(trial)));
        %burstduration
        fprintf(fid,'%6.1f\t%6.1f\t',...
            nanmean(burstdur(trial)),nanstd(burstdur(trial)));
        %#spike/burst
        for i=1:MaxBursts
            ss = spikeburst(trial,i);
            ss(ss==0)=nan;
            fprintf(fid,'%6.2f\t',nanmean(ss));
        end
        %histogram of GoodSpikes
        for i=1:50
            fprintf(fid,'%3.0f\t', sum(trialhist(trial,i)))
        end
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');
    %================================
    %data clumped by Stimulus

    stimuli = unique(allData(:,3));
    fprintf(fid,'Data by STIMULUS: %s\n\r', filebase)
    fprintf(fid,...
        ['Stim\t#Spike\t#S-std\t',...
        ' Lat(E)\tL(E)std\t',...
        ' Lat(B)\tL(B)std\t',...
        ' ISI  \tstd\t',...
        '  IBI  \tstd\t',...
        '  #Burst\t#Bstd\tDur1\tsDur1\t',...
        '#Spk/b1\t#Spk/b2\t#Spk/b3\t#Spk/b4r',...
        '\tISIhist\r'])

    for stim=1:length(stimuli)
        trial = find(allData(:,3)==stimuli(stim));
        fprintf(fid,'%d\t', ...
            stimuli(stim));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(Nspikes(trial)), nanstd(Nspikes(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
        allISI=[];
        for i=trial
            allISI = [allISI, trialISI{i}];
        end
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(allISI),nanstd(allISI));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(meanIBI(trial)),nanstd(meanIBI(trial)));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(numberBursts(trial)),nanstd(numberBursts(trial)));
        %burstduration
        fprintf(fid,'%6.1f\t%6.1f\t',...
            nanmean(burstdur(trial)),nanstd(burstdur(trial)));
        for i=1:MaxBursts
            ss = spikeburst(trial,i);
            ss(ss==0)=nan;
            fprintf(fid,'%6.2f\t',nanmean(ss));
        end
        %histogram of GoodSpikes
        for i=1:50
            fprintf(fid,'%3.0f\t', sum(trialhist(trial,i)))
        end
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by TIME
    times = unique(allData(:,4));
    times(find(isnan(times)))=[];
    if sum(~isnan(times))>0
        fprintf(fid,'Data by TIME: %s\n\r', filebase)
        fprintf(fid,...
            ['Time\t#Spike\t#S-std\t',...
            ' Lat(E)\tL(E)std\t',...
            ' Lat(B)\tL(B)std\t',...
            ' ISI  \tISIstd\t',...
            '  IBI  \tIBIstd\t',...
            '  #Burst\t#Bstd\tDur1\tsDur1\t'...
            '#Spk/b1\t#Spk/b2\t#Spk/b3\t#Spk/b4',...
            '\tISIhist\r'])

        for time=1:length(times)
            trial = find(allData(:,4)==times(time));
            fprintf(fid,'%d\t', ...
                times(time));
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(Nspikes(trial)), nanstd(Nspikes(trial)))
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
            allISI=[];
            for i=trial
                allISI = [allISI, trialISI{i}];
            end
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(allISI),nanstd(allISI));
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(meanIBI(trial)),nanstd(meanIBI(trial)));
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(numberBursts(trial)),nanstd(numberBursts(trial)));
            %burstduration
            fprintf(fid,'%6.1f\t%6.1f\t',...
                nanmean(burstdur(trial)),nanstd(burstdur(trial)));
            for i=1:MaxBursts
                ss = spikeburst(trial,i);
                ss(ss==0)=nan;
                fprintf(fid,'%6.2f\t',nanmean(ss));
            end
            %histogram of GoodSpikes
            for i=1:50
                fprintf(fid,'%3.0f\t', ...
                    sum(trialhist(trial,i)));
            end
            fprintf(fid,'\r');
        end
    end %of nonempty times
    %======================================
    %VOCAL -- PST only
    %======================================
    
    %data for individual sweeps -- VOCAL PST
    fprintf(fid,'PST by SWEEP: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tSwp\tStim\tTime\t',...
        'PST\n\r'])
    for trial=1:Ntrials
        %print the set,sweep and stim
        fprintf(fid,'%d\t%d\t%d\t%d\t', ...
            allData(trial,1),allData(trial,2),...
            allData(trial,3),allData(trial,4));
        
        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', (PSThist(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %====================================
    %data clumped by Set -- VOCAL PST
    MaxSet = max(allData(:,1));
    fprintf(fid,'PST by SET: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tStim\tTime\t',...
        'PST\r'])

    for set=1:MaxSet
        trial = find(allData(:,1)==set);
        fprintf(fid,'%d\t%d\t%d\t', ...
            set,allData(trial(1),3),allData(trial(1),4));
        
        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by Stimulus -- VOCAL PST

    stimuli = unique(allData(:,3));
    fprintf(fid,'PST by STIMULUS: %s\n\r', filebase)

    fprintf(fid,...
        ['tStim\t',...
        'PST\r'])

    for stim=1:length(stimuli)
        trial = find(allData(:,3)==stimuli(stim));
        fprintf(fid,'%d\t', ...
            allData(trial(1),3));

        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
        end

        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by TIME -- VOCAL PST
    times = unique(allData(:,4));
    times(find(isnan(times)))=[];
    if sum(~isnan(times))>0
        fprintf(fid,'PST by TIME: %s\n\r', filebase)
        fprintf(fid,...
            ['Time\t',...
            'PST\r'])

        for time=1:length(times)
            trial = find(allData(:,4)==times(time));
            fprintf(fid,'%d\t', ...
                times(time));
         
            %PST of AllSpikes
            for i=1:length(PSTcenters)
                fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
            end
            
            fprintf(fid,'\r');
        end %each time
    end
    fprintf(fid,'\r');
    %======================================
    %finish the file
    fclose(fid);
end %vocal nerve file

%======================================
%single unit
%======================================
if isSingle==1
    %get the null-stimulus background
    %find all zero stimuli
    zeroTrials = find(allData(:,3)==0);
    zeroMeanNumber = 0;
    %get the total spikes in these null sweeps
    for i=1:length(zeroTrials)
        zeroSpikes{i} = ...
            allData(zeroTrials(i),4+find(allData(zeroTrials(i),5:end)>0));
        zeroNumber(i) = length(zeroSpikes{i});
    end
    %average number per sweep
    %Since each sweep is 1 sec, this is also the rate/sec
    zeroMean = mean(zeroNumber);
    zeroStd = std(zeroNumber);

    %=================================
    %now make an results file -- single unit
    fileout = [filebase,'Output','.txt'];
    fileout = ([path,fileout]);
    fid = fopen(fileout,'w');

    %=================================
    %data for individual sweeps -- SINGLE UNIT
    fprintf(fid,'Data by SWEEP: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tSwp\tStim\tTime\t',...
        '#Spk\t#Spon \t#net  \tstd#net\t',...
        'Lat(End)\tLat(Begin)\tLag \t',...
        'meanISI\tstdISI\',...
        '\tISIhist\n\r'])
    for trial=1:Ntrials
        %print the set,sweep and stim
        fprintf(fid,'%d\t%d\t%d\t%d\t', ...
            allData(trial,1),allData(trial,2),...
            allData(trial,3),allData(trial,4));
        %print the Nspikes, expected spond #, difference
        fprintf(fid,'%d\t%6.2f\t%6.2f\t%6.2f\t', ...
            Nspikes(trial), ...
            zeroMean*TrialDuration(trial)*0.001, ...
            Nspikes(trial)-zeroMean*TrialDuration(trial)*0.001,...
            zeroStd*TrialDuration(trial)*0.001);
        %print the ISI
        fprintf(fid,'%6.1f\t%6.1f\t%6.1f\t%6.1f\t',...
            LatencyEnd(trial),LatencyBegin(trial),...
            LatencyDiff(trial),...
            meanISI(trial),stdISI(trial));
        %histogram of GoodSpikes
        for i=1:length(histcenters)
            fprintf(fid,'%3.0f\t',...
                trialhist(trial,i));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %====================================
    %data clumped by Set -- SINGLE UNIT
    MaxSet = max(allData(:,1));
    fprintf(fid,'Data by SET: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tStim\tTime\t',...
        '#Spk  \ts#Spk \t#Spon \t#net  \tstd#net\t',...
        ' Lat(E)\tL(E)std\t',...
        ' Lat(B)\tL(B)std\t',...
        'Lag \tLagstd\t',...
        ' ISI  \tISIstd\t',...
        '\tISIhist\r'])

    for set=1:MaxSet
        trial = find(allData(:,1)==set);
        fprintf(fid,'%d\t%d\t%d\t', ...
            set,allData(trial(1),3),allData(trial(1),4));
        %fprintf(fid,'%6.2f\t%6.2f\t', ...
        %nanmean(Nspikes(trial)), nanstd(Nspikes(trial)))
        %print the Nspikes, expected spond #, difference
        fprintf(fid,'%6.2f\t%6.2f\t%6.2f\t%6.2f\t', ...
            nanmean(Nspikes(trial)),nanstd(Nspikes(trial)), ...
            zeroMean*TrialDuration(trial(1))*0.001, ...
            nanmean(Nspikes(trial))-zeroMean*TrialDuration(trial(1))*0.001,...
            sqrt((zeroStd*TrialDuration(trial(1))*0.001)^2+...
            nanstd(Nspikes(trial))^2));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyDiff(trial)),nanstd(LatencyDiff(trial)))
        allISI=[];
        for i=trial
            allISI = [allISI, trialISI{i}];
        end
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(allISI),nanstd(allISI));

        %histogram of GoodSpikes
        for i=1:length(histcenters)
            fprintf(fid,'%3.0f\t', sum(trialhist(trial,i)))
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by Stimulus -- SINGLE UNIT

    stimuli = unique(allData(:,3));
    fprintf(fid,'Data by STIMULUS: %s\n\r', filebase)

    fprintf(fid,...
        ['tStim\t',...
        '#Spk  \ts#Spk \t#Spon \t#net  \tstd#net\t',...
        ' Lat(E)\tL(E)std\t',...
        ' Lat(B)\tL(B)std\t',...
        'Lag \tLagstd\t',...
        ' ISI  \tISIstd\t',...
        '\tISIhist\r'])

    for stim=1:length(stimuli)
        trial = find(allData(:,3)==stimuli(stim));
        fprintf(fid,'%d\t', ...
            allData(trial(1),3));
        %fprintf(fid,'%6.2f\t%6.2f\t', ...
        %nanmean(Nspikes(trial)), nanstd(Nspikes(trial)))
        %print the Nspikes, expected spond #, difference
        fprintf(fid,'%6.2f\t%6.2f\t%6.2f\t%6.2f\t', ...
            nanmean(Nspikes(trial)),nanstd(Nspikes(trial)), ...
            zeroMean*TrialDuration(trial(1))*0.001, ...
            nanmean(Nspikes(trial))-zeroMean*TrialDuration(trial(1))*0.001,...
            sqrt((zeroStd*TrialDuration(trial(1))*0.001)^2+...
            nanstd(Nspikes(trial))^2));
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyDiff(trial)),nanstd(LatencyDiff(trial)))
        allISI=[];
        for i=trial
            allISI = [allISI, trialISI{i}];
        end
        fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(allISI),nanstd(allISI));

        %histogram of GoodSpikes
        for i=1:length(histcenters)
            fprintf(fid,'%3.0f\t', sum(trialhist(trial,i)))
        end

        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by TIME -- SINGLE UNIT
    times = unique(allData(:,4));
    times(find(isnan(times)))=[];
    if sum(~isnan(times))>0
        fprintf(fid,'Data by TIME: %s\n\r', filebase)
        fprintf(fid,...
            ['Time\t',...
            '#Spk  \ts#Spk \t#Spon \t#net  \tstd#net\t',...
            ' Lat(E)\tL(E)std\t',...
            ' Lat(B)\tL(B)std\t',...
            'Lag \tLagstd\t',...
            ' ISI  \tISIstd\t',...
            '\tISIhist\r'])

        for time=1:length(times)
            trial = find(allData(:,4)==times(time));
            fprintf(fid,'%d\t', ...
                times(time));
            fprintf(fid,'%6.2f\t%6.2f\t%6.2f\t%6.2f\t', ...
                nanmean(Nspikes(trial)),nanstd(Nspikes(trial)), ...
                zeroMean*TrialDuration(trial(1))*0.001, ...
                nanmean(Nspikes(trial))-zeroMean*TrialDuration(trial(1))*0.001,...
                sqrt((zeroStd*TrialDuration(trial(1))*0.001)^2+...
                nanstd(Nspikes(trial))^2));
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(LatencyEnd(trial)),nanstd(LatencyEnd(trial)))
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(LatencyBegin(trial)),nanstd(LatencyBegin(trial)))
            fprintf(fid,'%6.2f\t%6.2f\t', ...
            nanmean(LatencyDiff(trial)),nanstd(LatencyDiff(trial)))
            allISI=[];
            for i=trial
                allISI = [allISI, trialISI{i}];
            end
            fprintf(fid,'%6.2f\t%6.2f\t', ...
                nanmean(allISI),nanstd(allISI));

            %ISI histogram -- GoodSpikes
            for i=1:length(histcenters)
                fprintf(fid,'%3.0f\t', sum(trialhist(trial,i)));
            end
           
            fprintf(fid,'\r');
        end %each time
    end
    fprintf(fid,'\r');
    
    %======================================
    %single unit -- PST only -- PSTcenters = [0:5:1000]
    %======================================
    
    %data for individual sweeps -- SINGLE UNIT PST
    fprintf(fid,'PST by SWEEP: PSTcenters = [0:5:1000] %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tSwp\tStim\tTime\t',...
        'PST\n\r'])
    for trial=1:Ntrials
        %print the set,sweep and stim
        fprintf(fid,'%d\t%d\t%d\t%d\t', ...
            allData(trial,1),allData(trial,2),...
            allData(trial,3),allData(trial,4));
        
        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', (PSThist(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %====================================
    %data clumped by Set -- SINGLE UNIT PST
    MaxSet = max(allData(:,1));
    fprintf(fid,'PST by SET: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tStim\tTime\t',...
        'PST\r'])

    for set=1:MaxSet
        trial = find(allData(:,1)==set);
        fprintf(fid,'%d\t%d\t%d\t', ...
            set,allData(trial(1),3),allData(trial(1),4));
        
        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by Stimulus -- SINGLE UNIT PST

    stimuli = unique(allData(:,3));
    fprintf(fid,'PST by STIMULUS: %s\n\r', filebase)

    fprintf(fid,...
        ['tStim\t',...
        'PST\r'])

    for stim=1:length(stimuli)
        trial = find(allData(:,3)==stimuli(stim));
        fprintf(fid,'%d\t', ...
            allData(trial(1),3));

        %PST of AllSpikes
        for i=1:length(PSTcenters)
            fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
        end

        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by TIME -- SINGLE UNIT PST
    times = unique(allData(:,4));
    times(find(isnan(times)))=[];
    if sum(~isnan(times))>0
        fprintf(fid,'PST by TIME: %s\n\r', filebase)
        fprintf(fid,...
            ['Time\t',...
            'PST\r'])

        for time=1:length(times)
            trial = find(allData(:,4)==times(time));
            fprintf(fid,'%d\t', ...
                times(time));
         
            %PST of AllSpikes
            for i=1:length(PSTcenters)
                fprintf(fid,'%3.0f\t', sum(PSThist(trial,i)));
            end
            
            fprintf(fid,'\r');
        end %each time
    end
    fprintf(fid,'\r');
    
    %======================================
    %single unit -- ALLIGNED PST only
    %======================================
    
    %data for individual sweeps --  ALLIGNED SINGLE UNIT PST
    fprintf(fid,'Aligned PST by SWEEP: PSTcentersAligned = [-300:5:700]; %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tSwp\tStim\tTime\t',...
        'PST\n\r'])
    for trial=1:Ntrials
        %print the set,sweep and stim
        fprintf(fid,'%d\t%d\t%d\t%d\t', ...
            allData(trial,1),allData(trial,2),...
            allData(trial,3),allData(trial,4));
        
        %PST of AllSpikes
        for i=1:length(PSTcentersAligned)
            fprintf(fid,'%3.0f\t', (PSThistAligned(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %====================================
    %data clumped by Set --  ALLIGNED SINGLE UNIT PST
    MaxSet = max(allData(:,1));
    fprintf(fid,'Aligned PST by SET: %s\n\r', filebase)
    fprintf(fid,...
        ['Set\tStim\tTime\t',...
        'PST\r'])

    for set=1:MaxSet
        trial = find(allData(:,1)==set);
        fprintf(fid,'%d\t%d\t%d\t', ...
            set,allData(trial(1),3),allData(trial(1),4));
        
        %PST of AllSpikes
        for i=1:length(PSTcentersAligned)
            fprintf(fid,'%3.0f\t', sum(PSThistAligned(trial,i)));
        end
        
        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by Stimulus --  ALLIGNED SINGLE UNIT PST

    stimuli = unique(allData(:,3));
    fprintf(fid,'Aligned PST by STIMULUS: %s\n\r', filebase)

    fprintf(fid,...
        ['tStim\t',...
        'PST\r'])

    for stim=1:length(stimuli)
        trial = find(allData(:,3)==stimuli(stim));
        fprintf(fid,'%d\t', ...
            allData(trial(1),3));

        %PST of AllSpikes
        for i=1:length(PSTcentersAligned)
            fprintf(fid,'%3.0f\t', sum(PSThistAligned(trial,i)));
        end

        fprintf(fid,'\r');
    end
    fprintf(fid,'\r');

    %================================
    %data clumped by TIME --  ALLIGNED SINGLE UNIT PST
    times = unique(allData(:,4));
    times(find(isnan(times)))=[];
    if sum(~isnan(times))>0
        fprintf(fid,'Aligned PST by TIME: %s\n\r', filebase)
        fprintf(fid,...
            ['Time\t',...
            'PST\r'])

        for time=1:length(times)
            trial = find(allData(:,4)==times(time));
            fprintf(fid,'%d\t', ...
                times(time));
         
            %PST of AllSpikes
            for i=1:length(PSTcentersAligned)
                fprintf(fid,'%3.0f\t', sum(PSThistAligned(trial,i)));
            end
            
            fprintf(fid,'\r');
        end %each time
    end
    fprintf(fid,'\r');
    %===============================================================
    %finish the file
    fclose(fid);
end %end single unit file



