Description of FitCCG

0001 function [dt1,dt2,travel,index] = FitCCG(t,ccg,varargin)
0002 
0003 %FitCCG - Fit dampened sinewave to the cross-correlogram of a pair of theta-modulated cells.
0004 %
0005 %  USAGE
0006 %
0007 %    [dt1,dt2,travel,index] = FitCCG(t,ccg,<options>)
0008 %
0009 %    t              time bins
0010 %    ccg            cross-correlogram
0011 %    <options>      optional list of property-value pairs (see table below)
0012 %
0013 %    =========================================================================
0014 %     Properties    Values
0015 %    -------------------------------------------------------------------------
0016 %     'show'        display and plot results (default = 'off')
0017 %    =========================================================================
0018 %
0019 %  OUTPUT
0020 %
0021 %    dt1            non-parametric estimate of first theta-modulated peak
0022 %    dt2            parametric estimate of first theta-modulated peak
0023 %    travel         travel time between field peaks (see Geisler et al. 2011)
0024 %    index          theta modulation index (see Royer et al. 2010)
0025 
0026 % Copyright (C) 2012-2015 by Michaël Zugaro
0027 %
0028 % This program is free software; you can redistribute it and/or modify
0029 % it under the terms of the GNU General Public License as published by
0030 % the Free Software Foundation; either version 3 of the License, or
0031 % (at your option) any later version.
0032 
0033 % Default values
0034 show = 'off';
0035 
0036 % Check parameters
0037 if nargin < 2,
0038   error('Incorrect number of parameters (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).');
0039 end
0040 
0041 if ~isdvector(t) || ~isdvector(ccg),
0042   error('Variables must be vectors (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).');
0043 end
0044 if length(t) ~= length(ccg),
0045   error('Variables are of different lengths (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).');
0046 end
0047 
0048 % Parse parameter list
0049 for i = 1:2:length(varargin),
0050     if ~ischar(varargin{i}),
0051         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).']);
0052     end
0053     switch(lower(varargin{i})),
0054         case 'show',
0055             show = varargin{i+1};
0056             if ~isastring(show,'off','on'),
0057                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).');
0058             end
0059         otherwise,
0060             error(['Unknown property ''' num2str(varargin{j}) ''' (type ''help <a href="matlab:help FitCCG">FitCCG</a>'' for details).']);
0061     end
0062 end
0063 
0064 global FitCCG_t FitCCG_ccg;
0065 FitCCG_t = t;
0066 FitCCG_ccg = ccg;
0067 
0068 ccg = ccg(:);
0069 
0070 % First estimate of dt: smooth the CCG and find the first peak for t>0
0071 smoothed = Smooth(ccg,3);
0072 maxima = IsExtremum([t smoothed]);
0073 peaksPos = t(t>0&maxima);
0074 peaksNeg = t(t<0&maxima);
0075 
0076 if isempty(peaksPos),
0077     if isempty(peaksNeg),
0078         dt1 = NaN;
0079     else
0080         dt1 = peaksNeg(size(peaksNeg,1));
0081     end
0082 else
0083     if isempty(peaksNeg) || peaksPos(1) < abs(peaksNeg(size(peaksNeg,1))),
0084         dt1 = peaksPos(1);
0085     else
0086         dt1 = peaksNeg(size(peaksNeg,1));
0087     end
0088 end
0089 
0090 
0091 
0092 
0093 
0094 
0095 % Try to find decent initial values
0096 tau = 0.5;
0097 mu = 0;
0098 omega = 7; % near theta
0099 b = mean(ccg);
0100 phi = 2*pi*omega*dt1;
0101 a = b;
0102 
0103 beta = [a tau mu omega phi b];
0104 
0105 % Compute best fit
0106 [beta,~,code] = fminsearch(@Fit,beta);
0107 a = abs(beta(1));
0108 tau = -abs(beta(2));
0109 mu = beta(3);
0110 omega = beta(4);
0111 phi = beta(5);
0112 b = beta(6);
0113 
0114 % Get the fit corresponding to the estimated parameters
0115 [~,fit] = Fit(beta);
0116 
0117 % Second estimate of dt: find the first peak of the model for t>0 or t<0
0118 % (the closest to zero)
0119 maxima = IsExtremum([t fit]);
0120 peaksPos = t(t>0&maxima);
0121 peaksNeg = t(t<0&maxima);
0122 
0123 if isempty(peaksPos),
0124     if isempty(peaksNeg),
0125         dt2 = NaN;
0126     else
0127         dt2 = peaksNeg(size(peaksNeg,1));
0128     end
0129 else
0130     if isempty(peaksNeg) || peaksPos(1) < abs(peaksNeg(size(peaksNeg,1))),
0131         dt2 = peaksPos(1);
0132     else
0133         dt2 = peaksNeg(size(peaksNeg,1));
0134     end
0135 end
0136 
0137 % Theta index
0138 index = a/b;
0139 
0140 % Determine travel time (see Geisler et al., 2011)
0141 carrier = Smooth(fit,5);
0142 [~,i] = max(carrier);
0143 travel = t(i);
0144 
0145 % Plot and display results
0146 if strcmp(show,'on'),
0147     printf('a       = %.3f',a);
0148     printf('tau     = %.3f',tau);
0149     printf('mu      = %.3f',mu);
0150     printf('omega   = %.3f',omega);
0151     printf('phi     = %.3f',phi);
0152     printf('b       = %.3f',b);
0153 
0154     figure;
0155     hold on;
0156     bar(t,ccg);
0157     plot(t,fit,'r','linewidth',5);
0158 end
0159 
0160 % ------------------------------------------------------------------------------------------------
0161 % Compute estimation error using a dampened sinewave model
0162 
0163 function [error,model] = Fit(beta)
0164 
0165 global FitCCG_t FitCCG_ccg;
0166 t = FitCCG_t;
0167 ccg = FitCCG_ccg;
0168 
0169 a = beta(1);
0170 tau = beta(2);
0171 mu = beta(3);
0172 omega = beta(4);
0173 phi = beta(5);
0174 b = beta(6);
0175 
0176 % Fit a dampened sinewave to the data
0177 model = (abs(a)*(sin(2*pi*omega*t + phi)+1) + b) .* exp(-abs(tau)*abs(t-mu));
0178 % Compute summed square error
0179 error = sum((ccg-model).^2);
0180 
0181 
0182 
0183 
0184 
0185 
0186 
0187 %  options = optimset('Display','iter','MaxFunEvals',50,'MaxIter',50,'TolX',1e-6,'TolFun',1e-6,'PlotFcns',@optimplotfval);%'FunValCheck');
0188 %  [beta,~,code] = fminsearch(@Fit,beta);%,options);
FitCCG

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
