function info = coherence_100ms_window(data,lfplabel,showplot)
%******** function info = coherence_100ms_window(data)
%    one example:  load 'msiv3_u8';
%                  info = coherence_100ms_window(data,'LFP1',1)
%
%                  load 'mnar4_u2';
%                  info = coherence_100ms_window(data,'LFP1',1);
%
%                  load 'msir27_u3';
%                  info = coherence_100ms_window(data,'LFP2',1);  
%
%***** general: computes the coherence in 100 ms hanning windows 
%*****          hopefully the same as in Womelsdorf and Fries
%*****  3/12/08: it runs coherence analysis off the isolated unit
%*****  in the data input 'Su1' with the labeled LFP channel (lpflabel) 
%****** and plots the computed coherence for attended and unattended trials
%*****  NOTE: data has been sorted of LFP1 is same lfp channel as the
%*****        isolated unit was recorded from, and same for MU1 mulit-unit
%*****  
%
%******* input:
%******
%****** data.label(nChannels) = {'SU1' ; 'MU1'; 'LFP1'; 'EYE1'};
%****** data.trials{ntrials} = [[Su1 data 1:T];[Mu1 data from 1:T];...];
%****** data.time{ntrials} = [1:T]
%****** data.attend(ntrials) = 1 if attended in RF
%****** data.fsample = 1000;
%****** data.sustained = XA:XB  % to mark the 800 ms sustained period
%****** data.waveform = [1:32];
%*****
%****** showplot - set this to 1 to see the results plotted out
%*****
%******  output:
%******      info.acoho  - attended coherence [1 x F]
%******      info.ascoho - std from permutes of estimate
%******      info.arcoho  - shuffled trials, attended coherence [1 x F]
%******      info.arscoho - shuffled trials, std from permutes of estimate
%******      info.afreq - frequencies [1 x F]
%******      info.ucoho  - attended coherence [1 x F]
%******      info.uscoho - std from jacknife of estimate
%******      info.urcoho  - shuffled trials, attended coherence [1 x F]
%******      info.usrcoho - shuffled trials, std from jacknife 
%******      info.ufreq - frequencies [1 x F]

%*********** locate su1
it_su = -1;
it_lfp = -1;
for zz = 1:size(data.label,2)
        if (strcmp(data.label{zz},'SU1'))
            it_su = zz;
        end
        if (strcmp(data.label{zz},lfplabel))
            it_lfp = zz;
        end
end
if ( (it_su == -1) | (it_lfp == -1) )
    disp(sprintf('Unable to identify SU1 or lfp label %s',lfplabel));
    info = [];
    return;
end

%********* get the single unit spikes for attended and unattended trials
attspikes = [];   % single unit attended spikes (in sustatined period)
uttspikes = [];
attlfp = [];
uttlfp = [];
for trial = 1:size(data.trials,2)
    if (data.attend(trial))
        attspikes = [attspikes ; data.trials{trial}(it_su,:)];
        attlfp = [attlfp ; data.trials{trial}(it_lfp,:)];
    else
        uttspikes = [uttspikes ; data.trials{trial}(it_su,:)];
        uttlfp = [uttlfp ; data.trials{trial}(it_lfp,:)];
    end
end

%************ plot the spike rasters and mean firing rates to sanity check
%************ and also show the local fields and mean local fields 
if (showplot == 1)
   
   spiker{1} = attspikes;
   spiker{2} = uttspikes;
   %****************
   disp('Plotting rastergrams (slow) ...');
   subplot('position',[0.1 0.4 0.4 0.55]);
   make_nice_spike_raster(spiker);
   V = axis;
   axis([0 size(attspikes,2) V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title(sprintf('Unit %s rasters',data.name));
   %*************
   subplot('position',[0.1 0.1 0.4 0.3]);
   smooth_window = 25;  % give sigma of 12.5ms
   make_nice_mean_raster(spiker,smooth_window);
   V = axis;
   axis([0 size(attspikes,2) V(3) V(4)]);
   plot([data.sustained(1),data.sustained(1)],[V(3),V(4)],'k-'); hold on;
   plot([data.sustained(end),data.sustained(end)],[V(3),V(4)],'k-'); hold on;
   ylabel('Firing Rate');
   xlabel('Time (ms)');
   
   spiker{1} = attlfp;
   spiker{2} = uttlfp;
   %****************
   disp('Plotting lfps per trial (slow) ...');
   subplot('position',[0.57 0.4 0.4 0.55]);
   make_nice_lfp_raster(spiker);
   V = axis;
   axis([0 size(attspikes,2) V(3) V(4)]);
   grid on;
   ylabel('Trial Number');
   title(sprintf('Unit %s lfp %s',data.name,lfplabel));
   %*************
   subplot('position',[0.57 0.1 0.4 0.3]);
   smooth_window = 25;  % give sigma of 12.5ms
   make_nice_mean_raster(spiker,smooth_window);
   V = axis;
   axis([0 size(attspikes,2) V(3) V(4)]);
   plot([data.sustained(1),data.sustained(1)],[V(3),V(4)],'k-'); hold on;
   plot([data.sustained(end),data.sustained(end)],[V(3),V(4)],'k-'); hold on;
   ylabel('LFP');
   xlabel('Time (ms)');
   
end

%**************** now call for the coherence estimates ***********
Win = 100;  % use 100 ms hanning tapered windows
MinCount = 200;  % use minimum 200 spikes per permuted estimate
%************************************************
disp('Computing coherence on attended trials');
[acoho,ascoho,arcoho,asrcoho,afreq] = fries_coherence_match(attlfp,attspikes,Win,MinCount);
disp('Computing coherence on non-attended trials');
[ucoho,uscoho,urcoho,usrcoho,ufreq] = fries_coherence_match(uttlfp,uttspikes,Win,MinCount);
%*********** store results *********
info.acoho = acoho;
info.ascoho = ascoho;
info.arcoho = arcoho;
info.asrcoho = asrcoho;
info.afreq = afreq;
info.ucoho = ucoho;
info.uscoho = uscoho;
info.urcoho = urcoho;
info.usrcoho = usrcoho;
info.ufreq = ufreq;
%***********************************

if (showplot == 1)
   figure;
   H = plot(afreq,acoho,'r-'); hold on;
   set(H,'Linewidth',2);
   H = plot(afreq,acoho+ascoho,'r-',afreq,acoho-ascoho,'r-'); hold on;
   set(H,'Linewidth',1);
   set(H,'Color',[1,0.5,0.5]);
   
   H = plot(afreq,arcoho,'r--'); hold on;
   set(H,'Linewidth',1.5);
   H = plot(afreq,arcoho+asrcoho,'r--',afreq,arcoho-asrcoho,'r--'); hold on;
   set(H,'Linewidth',1);
   set(H,'Color',[1,0.5,0.5]);
   
   H = plot(ufreq,ucoho,'b-'); hold on;
   set(H,'Linewidth',2);
   H = plot(ufreq,ucoho+uscoho,'b-',ufreq,ucoho-uscoho,'b-'); hold on;
   set(H,'Linewidth',1);
   set(H,'Color',[0.5,0.5,1]);
   
   H = plot(ufreq,urcoho,'b--'); hold on;
   set(H,'Linewidth',1.5);
   H = plot(ufreq,urcoho+usrcoho,'b--',ufreq,urcoho-usrcoho,'b--'); hold on;
   set(H,'Linewidth',1);
   set(H,'Color',[0.5,0.5,1]);
   
   ylabel('Coherence');
   xlabel('Frequency (Hz)');
   title(sprintf('%s with %s (dashed shuffled)',data.name,lfplabel));
   
end

return;


function [coho,scoho,urcoho,srcoho,ifcoher] = fries_coherence_match(speco,spiko,Win,MinCount)
%******* [coho,scoho,urcoho,srcoho,ifcoher] = fries_coherence(speco,spiko,Win)
%****** computes the coherence value using spike-triggered LFP similar to 
%****** that of Fries et al, 2001. in 100ms Hanning windowed segments 
%****** ALSO, to rate normalized, total set of spikes is divided into
%****** independent permutations of size MinCount spikes
%******
%****** inputs:
%******    speco - a MxT array where M is trials and T is time in ms
%******    spiko - a MxT array with same dimensions, 0 or 1 for spikes
%******    Win - a window around which to compute LFP segments 
%******    MinCount - number of spikes to use in Fries comp per permutation
%******             - if set to 0, don't rate normalize, run entire set in
%******               of spikes.
%****** outputs:
%******    coho - 1xF array of coherence per frequencies F
%******    scoho - std of the independent permutes in analysis
%******    urcoho - mean coherence of randomly shuffled trials
%******    srcoho - std of coherence from random shuffles
%******    ifcoher - list of frequency values, 1xF

  %**** compute number of spikes, must be at least 200 to run
  if (MinCount == 0)
      yspikes = find( spiko > 0);
      PN = size(yspikes,1);
      PRR = 1:size(yspikes,1);
      NPERM = 1;
      MinCount = PN;   % include every single spike
  else
      sumo = sum(sum(spiko));
      yspikes = find( spiko > 0 );
      NPERM = floor( sumo / MinCount );   % number of random permutations
      PN = size(yspikes,1);  % returns as a column of numbers
      PRR = randperm(PN);
  end
  
  %************* random permutations for shuffle corrected trials
  RB = 20;
  trialperm = [];
  for i = 1:RB
      trialperm = [trialperm ; randperm(size(speco,1))];  %random permutes of trial order
  end
  
  %*********** prep multi-taper fft ***************
  TT = size(speco,2);
  N = Win;
  if (N>(TT/4))
      N = floor(TT/4);
      disp(sprintf('Window size reduced to 1/4 total interval (%d to %d ms)',Win,N));
  end
  Trials = size(speco,1);
  %*******************
  TTa = 1+floor(Win/2);   %tighten up interval of analysis so Win does not fall outside
  TTb = TT-floor(Win/2);
  taper = hanning(Win,'symmetric')';
  %********** compute range of frequencies for FFT analysis
  pad = 0;
  fpass = [0.003 0.088];   % our LFP has a restricted range due to hardware filter
  Fs = 1;
  nfft=2^(nextpow2(N));
  df=Fs/nfft;
  f=0:df:Fs; % all possible frequencies
  f=f(1:nfft);
  findx=find(f>=fpass(1) & f<=fpass(end));
  f=f(findx);
  %**********************
  
  cohopool = [];     % LFP coherences on spike subsets
  urcohopool = [];   % LFP coherences over random permutes of trials (shuffle corrections)
  %****************
  for uber = 1:NPERM
      
    %********** do analysis many times on different permutes of same data **    
    spika = zeros(size(spiko));
    aa = 1 + ((uber-1)*MinCount);
    bb = (uber*MinCount);
    spika( yspikes(PRR(aa:bb)) ) = 1;   % select random subset of MinCount spikes
    disp(sprintf('Computing Subset %d of %d (contains %d spikes)',uber,NPERM,MinCount));
    
    %******************
    spikecount = 0;
    lfp_pow = [];
    lfp_sta = [];   % time rep of LFP
    rlfp_sta = cell(1,RB);  % random permutations (shuffle correction)
  
    %********* run through all trials, compute for each spike in the
    %********* subset of spikes the LFP segment, its FFT, power etc...
    for ii = 1:Trials
      for tt = TTa:TTb
         if (spika(ii,tt)>0)
            spikecount = spikecount + 1;
            toa = (tt-floor(Win/2));
            tob = toa + Win - 1;
            lfper = speco(ii,toa:tob) .* taper;
            %*************************
            J1=fft(lfper,nfft);
            J1=J1(1,findx);
            %*************************
            if (spikecount == 1)
                lfp_pow = conj(J1) .* J1;  % power spectrum 
                lfp_sta = J1;  
                %****************
                for rr = 1:RB
                   lfper = speco(trialperm(rr,ii),toa:tob);
                   lfper = lfper .* taper;
                   jj1 = fft(lfper,nfft);
                   jj1 = jj1(1,findx);
                   rlfp_sta{1,rr} = jj1;
                end
            else
                lfp_pow = lfp_pow + conj(J1) .* J1;   %Fries 2001 take average of power spectra
                lfp_sta = lfp_sta + J1;
                %****************
                for rr = 1:RB
                   lfper = speco(trialperm(rr,ii),toa:tob);
                   lfper = lfper .* taper;
                   jj1 = fft(lfper,nfft); 
                   jj1 = jj1(1,findx);
                   rlfp_sta{1,rr} = rlfp_sta{1,rr} + jj1;
                end
            end    
                    
         end
      end
    end  % for all trials    
    %**********************************
  
    %********** compute the coherence value
    lfp_sta = (lfp_sta / spikecount);
    S1 = sqrt(lfp_pow / spikecount);
    S12 = abs(lfp_sta);
    coho = (S12 ./ S1);
    
    %********* compute null distribution coherences ****
    %****** as average of several randomly permuted trials
    rcoho = [];
    for rr = 1:RB
      rlfp_sta{rr} = (rlfp_sta{rr}/spikecount);
      jja = abs(rlfp_sta{rr}); 
      cobo = (jja ./ S1);
      rcoho = [rcoho ; cobo];
    end
    urcoho = mean(rcoho);
  
    %********* pool the different permutes for later *********
    cohopool = [cohopool ; coho];
    urcohopool = [urcohopool ; urcoho];
  end
  
  ifcoher = f * 1000;
  if (NPERM < 1)  % require min set of samples to estimate cohere
      disp(sprintf('Ret Nan only %d spikes',spikecount));
      coho = NaN * ones(1,size(f,2));
      scoho = NaN * ones(1,size(f,2));
      urcoho = NaN * ones(1,size(f,2));
      srcoho = NaN * ones(1,size(f,2));
  else
    if (NPERM > 1)
      coho = mean(cohopool);
      scoho = std(cohopool) / sqrt( NPERM );
      urcoho = mean(urcohopool);
      srcoho = std(urcohopool) / sqrt( NPERM );
    else
      coho = cohopool;
      urcoho = urcohopool;
      %******** punt on std if too few spikes - could do a real jacknife
      %******** but leave it for later
      scoho = NaN * ones(1,size(f,2));
      srcoho = NaN * ones(1,size(f,2));
    end
  end
  
return;





%************** below here are some boring graphics routines to plot things

%****************** function to make a nice raster plot ****************
function make_nice_mean_raster(spmat,smooth_window)
%*********** spmat1 and spmat2 are spike matrices of two conditions you wish to compare
%*********** smooth_window ... gaussian smoothing in millisecs
numconds = size(spmat,2);
if (numconds==2)
    colo = [[1,0,0];[0,0,1]];
else
    colo = [[1,0,0];[0,1,0];[0,0,1]];
end

for k = 1:numconds
  spud = spmat{1,k};  
  numtrials(k) = size(spud,1);
  smorate = gauss_smooth(sum( spud(1:numtrials(k),:))/....
                      numtrials(k),smooth_window)*1000;
  H = plot(smorate,'k-'); hold on;
  set(H,'Color',colo(k,:));
end

return;


%******************* make a rastergram of the actual spikes on each trial
function make_nice_spike_raster(spmat)

colo = [[1,1,1];[0.5,0,0];[0,0,0.5]];
colormap(colo);

spmat{1,2} = spmat{1,2} * 2;  
totspike = [spmat{1,1} ; spmat{1,2}];
totspike = totspike + 1;
imagesc(totspike);
V = axis;
axis([V(1) V(2) 0 size(totspike,1)]);

return;

%******************* make a rastergram of the actual spikes on each trial
function make_nice_lfp_raster(spmat)

totspike = [spmat{1,1} ; spmat{1,2}];
mino = min(min(totspike)) * 0.5;  %allow some partial overlap
maxo = max(max(totspike)) * 0.5;
for ii = 1:size(totspike,1)
   offset = (size(totspike,1)-ii);
   yy = offset + ((totspike(ii,:)-mino)/(maxo-mino));
   if (ii <= (size(totspike,1)/2))
       plot(1:size(totspike,2),yy,'r-'); hold on;
   else
       plot(1:size(totspike,2),yy,'b-'); hold on;
   end
end
V = axis;
axis([V(1) V(2) 0 size(totspike,1)]);

return;

%**************************************************************
function output = gauss_smooth(input, window)
% Smoothing function: 
% output = smooth(input, window)
% "Window" is the total kernel width.
% Input array must be one-dimensional.

input_dims = ndims(input);
input_size = size(input);
if input_dims > 2 | min(input_size) > 1,
   disp('Input array is too large.');
   return
end

if input_size(2) > input_size(1),
   input = input';
   toggle_dims = 1;
else
   toggle_dims = 0;
end

if window/2 ~= round(window/2),
   window = window + 1;
end
halfwin = window/2;

input_length = length(input);
%********* gauss window +/- 1 sigma
x = -halfwin:1:halfwin;
kernel = exp(-x.^2/(window/2)^2);
kernel = kernel/sum(kernel);

padded(halfwin+1:input_length+halfwin) = input;
padded(1:halfwin) = ones(halfwin, 1)*input(1);
padded(length(padded)+1:length(padded)+halfwin) = ones(halfwin, 1)*input(input_length);

output = conv(padded, kernel);
output = output(window:input_length+window-1);

if toggle_dims == 1,
   output = output';
end

return;