rate_normalized_coherence.m

computes the coherence in hanning windows

  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)
  phaseopool = [];
  %****************
  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);
    phaseo = lfp_sta ./ 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];
    phaseopool = [phaseopool ; phaseo];
  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));
      phaso = 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 );
      phaso = mean(phaseopool);
    else
      coho = cohopool;
      phaso = phaseopool;
      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,0,1];[0,1,0];[1,1,0]];
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];[0,0.5,0]];
colormap(colo);

totspike = [];
for cubo = 1:size(spmat,2)
    totspike = [totspike; (cubo*spmat{1,cubo}) ];
end
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)

colo = 'rbgy';

totspike = [];
CNUM = size(spmat,2);
for ii = 1:size(spmat,2)
    totspike = [totspike ; spmat{1,ii}];
end
splito = size(totspike,1)/CNUM;

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));
   cubo = 1 + floor((ii-1)/splito);
   plot(1:size(totspike,2),yy,[colo(cubo),'-']); hold on;
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;