ThresholdSpikes

ThresholdSpikes - Post-hoc threshold correction for spike detection.

  USAGE

    spikes = ThresholdSpikes(amplitudes,factor,<options>)

    amplitudes     [time, electrode group, cluster, amplitude] <a href="matlab:help samples">samples</a>
                   e.g. obtained using <a href="matlab:help GetSpikeAmplitudes">GetSpikeAmplitudes</a>
    factor         the new threshold will be computed as the former threshold
                   multiplied by this factor
    <options>      optional list of property-value pairs (see table below)

    =========================================================================
     Properties    Values
    -------------------------------------------------------------------------
     'units'       existing units (this is used to generate new cluster IDs)
    =========================================================================

  OUTPUT

    The subset of spikes exceeding the new threshold. The cluster IDs are
    unchanged, unless the optional parameter 'units' is provided, in which
    case each unit is assigned a new cluster ID.

  SEE

    See also GetSpikeAmplitudes.

0001 function spikes = ThresholdSpikes(amplitudes,factor,varargin)
0002 
0003 %ThresholdSpikes - Post-hoc threshold correction for spike detection.
0004 %
0005 %  USAGE
0006 %
0007 %    spikes = ThresholdSpikes(amplitudes,factor,<options>)
0008 %
0009 %    amplitudes     [time, electrode group, cluster, amplitude] <a href="matlab:help samples">samples</a>
0010 %                   e.g. obtained using <a href="matlab:help GetSpikeAmplitudes">GetSpikeAmplitudes</a>
0011 %    factor         the new threshold will be computed as the former threshold
0012 %                   multiplied by this factor
0013 %    <options>      optional list of property-value pairs (see table below)
0014 %
0015 %    =========================================================================
0016 %     Properties    Values
0017 %    -------------------------------------------------------------------------
0018 %     'units'       existing units (this is used to generate new cluster IDs)
0019 %    =========================================================================
0020 %
0021 %  OUTPUT
0022 %
0023 %    The subset of spikes exceeding the new threshold. The cluster IDs are
0024 %    unchanged, unless the optional parameter 'units' is provided, in which
0025 %    case each unit is assigned a new cluster ID.
0026 %
0027 %  SEE
0028 %
0029 %    See also GetSpikeAmplitudes.
0030 
0031 % Copyright (C) 2009-2012 by Michaël Zugaro
0032 %
0033 % This program is free software; you can redistribute it and/or modify
0034 % it under the terms of the GNU General Public License as published by
0035 % the Free Software Foundation; either version 3 of the License, or
0036 % (at your option) any later version.
0037 
0038 % Default values
0039 allUnits = [];
0040 
0041 % Parse parameters
0042 if nargin < 2,
0043     error('Missing parameters (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).');
0044 end
0045 if ~isdmatrix(amplitudes),
0046     error('Incorrect amplitudes (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).');
0047 end
0048 if ~isdscalar(factor),
0049     error('Incorrect factor (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).');
0050 end
0051 
0052 if mod(length(varargin),2) ~= 0,
0053     error('Incorrect number of parameters (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).');
0054 end
0055 
0056 % Parse options
0057 for i = 1:2:length(varargin),
0058     if ~ischar(varargin{i}),
0059         error(['Parameter ' num2str(i+1) ' is not a property (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).']);
0060     end
0061     switch(lower(varargin{i})),
0062         case 'units',
0063             allUnits = varargin{i+1};
0064             if ~isdmatrix(allUnits,'>=0') || size(allUnits,2) ~= 2,
0065                 error('Incorrect value for property ''units'' (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).');
0066             end
0067         otherwise,
0068             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help ThresholdSpikes">ThresholdSpikes</a>'' for details).']);
0069     end
0070 end
0071 
0072 % Which units are we processing?
0073 units = unique(amplitudes(:,2:3),'rows');
0074 
0075 spikes = [];
0076 for i = 1:size(units,1),
0077     group = units(i,1);
0078     cluster = units(i,2);
0079     % Amplitudes for this unit
0080     this = amplitudes(amplitudes(:,2)==group&amplitudes(:,3)==cluster,:);
0081     % Guess former threshold
0082     peak = max(abs(this(:,4:end)),[],2);
0083     threshold = min(peak);
0084     % Update
0085     threshold = factor*threshold;
0086     % Spikes to keep
0087     keep = this(peak>=threshold,1:3);
0088     if ~isempty(allUnits),
0089         % New cluster ID = largest existing ID + 1
0090         newCluster = max(allUnits(allUnits(:,1)==group,2)) + 1;
0091         allUnits = [allUnits;group newCluster];
0092         keep(:,3) = newCluster;
0093         disp(['Unit (' int2str(group) ',' int2str(cluster) ') -> (' int2str(group) ',' int2str(newCluster) ')']);
0094     end
0095     spikes = [spikes;keep];
0096 end