Description of FindSpindles

0001 function spindles = FindSpindles(filtered,varargin)
0002 
0003 %FindSpindles - Find thalamo-cortical spindles (9-17Hz oscillations).
0004 %
0005 %  USAGE
0006 %
0007 %    spindles = FindSpindles(filtered,<options>)
0008 %
0009 %    filtered       spindle-band filtered LFP <a href="matlab:help samples">samples</a> (one channel). This must
0010 %                   be restricted to slow wave sleep periods for the algorithm
0011 %                   to perform best.
0012 %    <options>      optional list of property-value pairs (see table below)
0013 %
0014 %    =========================================================================
0015 %     Properties    Values
0016 %    -------------------------------------------------------------------------
0017 %     'threshold'   minimum z-scored envelope tail (default = 2.5)
0018 %     'peak'        minimum z-scored envelope peak (default = 5)
0019 %     'durations'   min and max spindle durations in ms (default = [400 3500])
0020 %    =========================================================================
0021 %
0022 %  OUTPUT
0023 %
0024 %    spindles       for each spindle, start, peak and stop times, and z-scored
0025 %                   peak amplitude, in a Nx4 matrix [start peak_t end peak_z]
0026 %
0027 
0028 % Copyright (C) 2012-2016 Michaël Zugaro, 2012-2015 Nicolas Maingret, 2016 Ralitsa Todorova
0029 %
0030 % This program is free software; you can redistribute it and/or modify
0031 % it under the terms of the GNU General Public License as published by
0032 % the Free Software Foundation; either version 3 of the License, or
0033 % (at your option) any later version.
0034 
0035 % Defaults
0036 minDistance = 0.5;
0037 threshold = 2.5;
0038 minPeak = 5 ;
0039 minDuration = 0.4;
0040 maxDuration = 3.5;
0041 
0042 % Constants
0043 minBoutDuration = 0.050;
0044 window = 0.320;
0045 
0046 % Check number of parameters
0047 if nargin < 1 | mod(length(varargin),2) ~= 0,
0048   error('Incorrect number of parameters (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).');
0049 end
0050 
0051 % Check parameter sizes
0052 if ~isdmatrix(filtered,'@2'),
0053     error('Parameter ''filtered'' is not a Nx2 matrix (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).');
0054 end
0055 
0056 % Parse parameter list
0057 for i = 1:2:length(varargin),
0058     if ~ischar(varargin{i}),
0059         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).']);
0060     end
0061     switch(lower(varargin{i})),
0062         case 'threshold',
0063             threshold = varargin{i+1};
0064             if ~isdscalar(threshold,'>0'),
0065                 error('Incorrect value for property ''threshold'' (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).');
0066             end
0067         case 'peak',
0068             minPeak = varargin{i+1};
0069             if ~isdscalar(minPeak,'>0'),
0070                 error('Incorrect value for property ''peak'' (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).');
0071             end
0072         case 'durations',
0073             durations = varargin{i+1};
0074             if ~isdvector(durations,'#2','>0'),
0075                 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).');
0076             end
0077             minDuration = durations(1)/1000;
0078             maxDuration = durations(2)/1000;
0079         otherwise,
0080             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help FindSpindles">FindSpindles</a>'' for details).']);
0081     end
0082 end
0083 
0084 
0085 % Compute low-pass filtered squared envelope of z-scored signal
0086 time = filtered(:,1);
0087 frequency = 1/median(diff(time));
0088 windowLength = round(frequency*window);
0089 window = ones(windowLength,1)/windowLength;
0090 envelope = abs(hilbert(zscore(filtered(:,2)))).^2;
0091 envelope = filtfilt(window,1,envelope);
0092 
0093 % Find start/stop as indices of threshold crossings, and discard incomplete pairs
0094 crossings = envelope > threshold;
0095 crossings = diff(crossings);
0096 start = find(crossings>0);
0097 stop = find(crossings<0);
0098 if time(stop(1)) < time(start(1)), stop(1) = []; end
0099 if time(start(end)) > time(stop(end)), start(end) = []; end
0100 
0101 % Discard short events
0102 tooShort = time(stop) - time(start) < minBoutDuration;
0103 start(tooShort) = [];
0104 stop(tooShort) = [];
0105 
0106 % Merge events when they are too close
0107 while true,
0108     tooClose = find(time(start(2:end)) - time(stop(1:end-1)) < minDistance);
0109     if isempty(tooClose), break; end
0110     start(tooClose+1) = [];
0111     stop(tooClose) = [];
0112 end
0113 
0114 % Discard events when envelope peak is too small
0115 % For each timestamp of the LFP, get the id of the spindle
0116 [~,id] = InIntervals(time,time([start stop]));
0117 % Find peak amplitudes and indices
0118 % (add 1 so that Accumulate does not complain when id=0, i.e. time points that do not belong to any candidate interval)
0119 [peak_z,peak_i] = Accumulate(id+1,envelope,'mode','max');
0120 % (now get rid of id=0, i.e. time points that do not belong to any candidate interval)
0121 peak_z(1) = []; peak_i(1) = []; 
0122 % Keep only spindles with large enough amplitude
0123 tooSmall = peak_z < minPeak;
0124 start(tooSmall) = [];
0125 stop(tooSmall) = [];
0126 peak_z(tooSmall) = [];
0127 peak_i(tooSmall) = [];
0128 
0129 % Replace indices with timestamps
0130 start = time(start);
0131 stop = time(stop);
0132 peak_t = time(peak_i);
0133 
0134 % Spindle start may be inaccurate due to leading delta wave, we need to correct for this
0135 % 1) Update threshold to 1/3 peak (if this is higher than previous threshold)
0136 newThrehold = peak_z/3;
0137 % 2) Select spindles that need correction
0138 indices = find(newThrehold>threshold); 
0139 % 3) Find last threshold crossing before peak (vectorized code)
0140 % Candidate intervals for the new spindle start = between the previous spindle end and current spindle peak
0141 if indices(1) > 1, candidateIntervals = [stop(indices-1) peak_t(indices)];
0142 else candidateIntervals = [0 peak_t(indices(1)); stop(indices(2:end)-1) peak_t(indices(2:end))]; end
0143 % For each timestamp, get the id of the interval in which it falls
0144 [in,id] = InIntervals(time,candidateIntervals);
0145 belowThreshold = false(size(id));
0146 belowThreshold(in) = envelope(in) < newThrehold(indices(id(in)));
0147 dt = mode(diff(time));
0148 start(indices) = Accumulate(id(in&belowThreshold),time(in&belowThreshold),'mode','max')+dt;
0149 
0150 
0151 % Merge events when they are too close
0152 n = length(peak_z);
0153 while true,
0154     tooClose = find(start(2:end) - stop(1:end-1) < minDistance);
0155     if isempty(tooClose), break; end
0156     start(tooClose+1) = [];
0157     stop(tooClose) = [];
0158     [peak_z,i] = max([peak_z(tooClose) peak_z(tooClose+1)],[],2);
0159     peak_t = peak_t(tooClose+i-1);
0160 end
0161 
0162 % Discard events that are too long or too short
0163 duration = stop-start;
0164 bad = abs(duration)<minDuration | abs(duration)>maxDuration;
0165 start(bad) = [];
0166 stop(bad) = [];
0167 peak_z(bad) = [];
0168 peak_t(bad) = [];
0169 
0170 spindles = [start peak_t stop peak_z];
0171 spindles = unique(spindles,'rows');
FindSpindles

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
