

%Correlation between vocalization times and spike trains
%Fig 1 is raster
%Fig 2 is individual vos
%Fig 3 is scan average
%Fig 4 is grand average

clear all

filename= '08-13-03_8_spiketimes.txt'; %original calibration file
%filename= '08-07-03_6_spiketimes.txt'; %null results
%filename= '08-26-03_1_spiketimes.txt';
%filename= '08-20-03_2_spiketimes.txt'; %small peak
%filename = '08-07-03u1r1.txt' ;
%filename = '08-07-03u2r7.txt' ;
fid = (['C:\My Documents\BigData\Weeg\SpikeTimesSept2003\',filename]);
%fid = (['e:\spike times text\',filename]);

allData = textread (fid,'%s','delimiter','\n','bufsize',2^14)  ;
%FILE FORMAT
%2 descriptor lines
% numberSweeps eff data lines (some may be blank)
% one blank line
%2 descriptor lines
% numberSweeps voc data lines (some may be blank)

%extract actual sweeps from file
numberSweeps=20;
%align on rising phase if 1, otherwise align on falling phase
rising=1;
%separate loger voc fro shorter ones
%define arbitrary cutoff
ShortCutoff = 1; %mSec
LongCutoff = inf;

%tokenize and convert file to a cell-array of spike times 
%the strange limits on the loops are a result of line offsets
%in the original file.
count=1;
for i=3:numberSweeps+2
    Spikelist=[];
    [T,R]=strtok(allData{i}); %get first token
    if (~isempty(T))          %if no data then skip it
        Spikelist = str2num(T);
        while(~isempty(R))
            [T,R]=strtok(R);
            %if (~isempty(T))
            Spikelist = [Spikelist str2num(T)];
            %end
        end
    end
    effSpike{count}=Spikelist ;
    count=count+1 ;
end

count=1;
for i=numberSweeps+6:2*numberSweeps+5
    Spikelist=[];
    [T,R]=strtok(allData{i}); %get first token
    if (~isempty(T))          %if no data then skip it
        Spikelist = str2num(T);
        while(~isempty(R))
            [T,R]=strtok(R); %get next token
            %if (~isempty(T))
            Spikelist = [Spikelist str2num(T)]; %and add it to list
            %end
        end
    end
    vocSpike{count}= Spikelist ;
    count=count+1 ;  
end

%======================================================================
%plot raw spikes
%and convolve spikes with gaussian

%set voc width and make normalized gaussian
vocwidth = 5;
voctk = -3*vocwidth:3*vocwidth ;
vockernel = exp(-(voctk/vocwidth).^2/2)/(vocwidth*sqrt(2*pi));

%set efferent averaging width and make gaussian
effwidth = 2;
efftk = -3*effwidth:3*effwidth ;
effkernel = exp(-(efftk/effwidth).^2/2)/(effwidth*sqrt(2*pi));

%make spike rasters and draw the raster plot
%then convolve the spike rasters with the gaussians
%produced above
figure(1)
clf

dt=1; %mSec
t=0:dt:5000;

count=1;
for i=1:numberSweeps
    spiketrainV = zeros(size(t));
    spiketrainV(fix((vocSpike{i}+1)/dt)) = 1;
    plot(t,0.9*spiketrainV+count)
    hold on
    count=count+1;
    vocConv(i,:) = conv(spiketrainV,vockernel);
end

count=1;
for i=1:numberSweeps
    spiketrainE = zeros(size(t));
    spiketrainE(fix((effSpike{i}+1)/dt)) = 1;
    plot(t,0.9*spiketrainE+count+.1,'r')
    hold on
    count=count+1;
    effConv(i,:) = conv(spiketrainE,effkernel);
end

%======================================================================
%plot individual voc-synched traces
figure(2)
clf
map = jet(numberSweeps);
%signal range +/- from selected edge of vocalization
signalRange=200;

%init ALL sweep averages
effMean = zeros(numberSweeps,2*signalRange+1);
vocMean = zeros(numberSweeps,2*signalRange+1);
    
%average over all traces -- init the avg
effGrandMean = zeros(1,2*signalRange+1);
vocGrandMean = zeros(1,2*signalRange+1);
vocGrandNum = 0;

for i=1:numberSweeps
    col = map(i,:);
    
    %find the 1/2 rise/fall point on the voc
    maxConv = max(vocConv(i,:)) ;
    %start point: difference of thresholded voc > 0
    vocStart = diff(vocConv(i,:)>maxConv/2) > 0 ; 
    %end point: difference of thresholded voc < 0
    vocEnd = diff(vocConv(i,:)>maxConv/2) < 0 ;
    
    %now get the actual starting index list,
    %the ending index list.
    vocStartindex = find(vocStart==1);
    vocEndindex = find(vocEnd==1);
    
    %init the current sweep count
    vocNum = 0;

    %test for align on rise/fall
    if (rising==0)
        vocindex = vocEndindex;
    else
        vocindex = vocStartindex;
    end
    
    %now recalculate voc conv with narrow kernel
    spiketrainV(fix((vocSpike{i}+1)/dt)) = 1;
    vocConv(i,:) = conv(spiketrainV,effkernel);
    
    %test for actual voc in the sweep
    if ~isempty(vocindex)
        %and then once thru for each voc
        for j =1:length(vocindex)
            
            %test to make sure whole range is in trace
            if (-signalRange+vocindex(j))>0 & (signalRange+vocindex(j))<5000
                
                %test for voc length for inclusion
                if ((vocEndindex(j)-vocStartindex(j))>=ShortCutoff) & ...
                        ((vocEndindex(j)-vocStartindex(j))<LongCutoff)
                    
                    vocNum = vocNum +1;
                    
                    effMean(i,:) = effMean(i,:) +...
                        effConv(i,[-signalRange:signalRange]+vocindex(j));
                    effGrandMean = effGrandMean +...
                        effConv(i,[-signalRange:signalRange]+vocindex(j));
                    
                    %show each single eff
                    subplot(2,1,1) 
                    plot(effConv(i,[-signalRange:signalRange]+vocindex(j)),'color',col,'linewidth',2)
                    hold on
                    
                    vocMean(i,:) = vocMean(i,:) +...
                        vocConv(i,[-signalRange:signalRange]+vocindex(j));
                    vocGrandMean = vocGrandMean +...
                        vocConv(i,[-signalRange:signalRange]+vocindex(j));
                    
                    %show each single voc
                    subplot(2,1,2)
                    plot(vocConv(i,[-signalRange:signalRange]+vocindex(j)),'color',col,'linewidth',2)
                    hold on
                end %end of test for voc length for inclusion
            end %end of test to make sure whole range is in sweep
        end %end of once thru for each voc
        
        %compute these only if vocNum>0
        effMean(i,:) = effMean(i,:)./vocNum; 
        vocMean(i,:) = vocMean(i,:)./vocNum;
        vocGrandNum = vocGrandNum + vocNum;
        
    end %end of check for empty start/end list 
end %end of sweep loop

vocGrandMean = vocGrandMean/vocGrandNum ;
effGrandMean = effGrandMean/vocGrandNum ;
%===============================================================
%plot voc-synched, sweep-average trace
figure(3)
clf
for i=1:numberSweeps
    col = map(i,:);
    subplot(2,1,1)
    plot(effMean(i,:),'color',col,'linewidth',2)
    title(['File=',filename,'  VocWidth=',num2str(vocwidth),...
            '  EffWidth=',num2str(effwidth),...
            '  VocDur=[',num2str(ShortCutoff),',',num2str(LongCutoff),']', ...
            '  N=', num2str(vocGrandNum)])
    hold on
    subplot(2,1,2)
    plot(vocMean(i,:),'color',col,'linewidth',2)
    hold on
end

%==================================================================
%plot grand average, voc-synched trace
figure(4)
clf
plot(effGrandMean,'r')
title(['File=',filename,'   VocWidth=',num2str(vocwidth),...
        '   EffWidth=',num2str(effwidth), ...
        '  VocDur=[',num2str(ShortCutoff),',',num2str(LongCutoff),']', ...
        '  N=', num2str(vocGrandNum)])
hold on
plot(vocGrandMean,'b')
[max,maxpos] = max(effGrandMean);
ypos = get(gca,'ylim');
text(20,ypos(2)/2,['Eff peak pos=',num2str(maxpos)])
