Frequency

Frequency - Compute instantaneous frequency for a point process (e.g. a spike train).

 Compute instantaneous frequency either by smoothing the point process using a fixed
 or adaptive Gaussian kernel, or by simply taking the inverse of the interevent intervals.
 Adaptive kernel smoothing is described in Richmond et al. (1990).

  USAGE

    frequency = Frequency(timestamps,<options>)

    timestamps     list of timestamps
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'method'      convolute with a fixed ('fixed', default) or adaptive
                   ('adaptive') Gaussian kernel, or compute inverse inter-
                   event intervals ('inverse')
     'limits'      [start stop] in seconds (default = approx. first and last
                   timestamps)
     'binSize'     bin size in seconds (default = 0.050)
     'smooth'      Gaussian kernel width in number of samples (default = 2)
     'show'        plot results (default = 'off')
    =========================================================================

  NOTE

    While the convolution method returns a Mx2 matrix (time bins in column 1,
    frequencies in column 2), the inverse interspike intervals method returns
    a Nx2 matrix (original spike timestamps in column 1, frequencies in column 2)

0001 function frequency = Frequency(timestamps,varargin)
0002 
0003 %Frequency - Compute instantaneous frequency for a point process (e.g. a spike train).
0004 %
0005 % Compute instantaneous frequency either by smoothing the point process using a fixed
0006 % or adaptive Gaussian kernel, or by simply taking the inverse of the interevent intervals.
0007 % Adaptive kernel smoothing is described in Richmond et al. (1990).
0008 %
0009 %  USAGE
0010 %
0011 %    frequency = Frequency(timestamps,<options>)
0012 %
0013 %    timestamps     list of timestamps
0014 %    <options>      optional list of property-value pairs (see table below)
0015 %
0016 %    =========================================================================
0017 %     Properties    Values
0018 %    -------------------------------------------------------------------------
0019 %     'method'      convolute with a fixed ('fixed', default) or adaptive
0020 %                   ('adaptive') Gaussian kernel, or compute inverse inter-
0021 %                   event intervals ('inverse')
0022 %     'limits'      [start stop] in seconds (default = approx. first and last
0023 %                   timestamps)
0024 %     'binSize'     bin size in seconds (default = 0.050)
0025 %     'smooth'      Gaussian kernel width in number of samples (default = 2)
0026 %     'show'        plot results (default = 'off')
0027 %    =========================================================================
0028 %
0029 %  NOTE
0030 %
0031 %    While the convolution method returns a Mx2 matrix (time bins in column 1,
0032 %    frequencies in column 2), the inverse interspike intervals method returns
0033 %    a Nx2 matrix (original spike timestamps in column 1, frequencies in column 2)
0034 
0035 % Copyright (C) 2004-2011 by Michaël Zugaro
0036 %
0037 % This program is free software; you can redistribute it and/or modify
0038 % it under the terms of the GNU General Public License as published by
0039 % the Free Software Foundation; either version 3 of the License, or
0040 % (at your option) any later version.
0041 
0042 % Defaults
0043 method = 'fixed';
0044 binSize = 0.05;
0045 smooth = 2;
0046 limits = [];
0047 show = 'off';
0048 
0049 % Check number of parameters
0050 if nargin < 1,
0051   error('Incorrect number of parameters (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0052 end
0053 
0054 % Check parameter sizes
0055 if size(timestamps,2) ~= 1,
0056     error('Parameter ''timestamps'' is not a vector (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0057 end
0058 
0059 % Parse parameter list
0060 for i = 1:2:length(varargin),
0061     if ~ischar(varargin{i}),
0062         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).']);
0063     end
0064     switch(lower(varargin{i})),
0065 
0066         case 'method',
0067             method = lower(varargin{i+1});
0068             if ~isastring(method,'fixed','adaptive','inverse','iisi'),
0069                 error('Incorrect value for property ''method'' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0070             end
0071 
0072         case 'binsize',
0073             binSize = varargin{i+1};
0074             if ~isdscalar(binSize,'>=0'),
0075                 error('Incorrect value for property ''binSize'' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0076             end
0077 
0078         case 'limits',
0079             limits = varargin{i+1};
0080             if ~isdvector(limits,'#2'),
0081                 error('Incorrect value for property ''limits'' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0082             end
0083 
0084         case 'smooth',
0085             smooth = varargin{i+1};
0086             if ~isdscalar(smooth,'>=0'),
0087                 error('Incorrect value for property ''smooth'' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0088             end
0089 
0090         case 'show',
0091             show = varargin{i+1};
0092             if ~isastring(show,'on','off'),
0093                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).');
0094             end
0095 
0096         otherwise,
0097             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help Frequency">Frequency</a>'' for details).']);
0098 
0099   end
0100 end
0101 
0102 if isempty(limits),
0103     limits = [timestamps(1)-10*binSize timestamps(end)+10*binSize];
0104 end
0105 
0106 if strcmp(method,'inverse') | strcmp(method,'iisi'),
0107     % Inter-spike intervals
0108     ds = diff(timestamps);
0109     f1 = 1./ds;
0110     t = timestamps(1:end-1) + ds/2;
0111     f2 = Interpolate([t f1],timestamps(2:end-1));
0112     f2 = [timestamps(1) 0;f2;timestamps(end) 0];
0113     frequency = f2;
0114 else
0115     % Smoothing
0116     % 1) Fixed-kernel
0117     t = (limits(1):binSize:limits(2))';
0118     T = Bin(timestamps,t);
0119     binned = Accumulate(T,1,size(t));
0120     f = Smooth(binned/binSize,smooth);
0121     frequency = [t f];
0122     if strcmp(method,'adaptive'),
0123         % 2) Variable-kernel (requires the above 'pilot' fixed-kernel estimate)
0124         % Compute variable-kernel sigma
0125         N = length(t);
0126         mu = exp(1/N*sum(log(f(f~=0)))); % Correction: geometric mean is computed using only non-zero values
0127         lambda = sqrt(f/mu);
0128         sigma = (smooth*binSize)./lambda;
0129         sigma = Clip(sigma,0,N*binSize/3);
0130         % Perform variable-kernel smoothing
0131         binned = [flipud(binned);binned;flipud(binned)];
0132         for i = 1:N,
0133             % Gaussian
0134             x = (0:binSize:3*sigma(i))';
0135             x = [flipud(-x);x(2:end)];
0136             kernel = exp(-x.^2/sigma(i)^2);
0137             kernel = kernel/sum(kernel);
0138             % 'Pointwise' convolution
0139             n = (length(kernel)-1)/2;
0140             bins = N + i + (-n:n);
0141             S2(i) = sum(binned(bins).*kernel)/binSize;
0142         end
0143         frequency = [t S2'];
0144     end
0145 end
0146 
0147 if strcmp(lower(show),'on'),
0148     PlotXY(frequency);
0149     hold on;
0150     PlotTicks(timestamps,'size',10,'k');
0151 end