Description of FindRipples

0001 function [ripples,sd,bad] = FindRipples(filtered,varargin)
0002 
0003 %FindRipples - Find hippocampal ripples (100~200Hz oscillations).
0004 %
0005 %  USAGE
0006 %
0007 %    [ripples,stdev,noise] = FindRipples(filtered,<options>)
0008 %
0009 %    Ripples are detected using the normalized squared signal (NSS) by
0010 %    thresholding the baseline, merging neighboring events, thresholding
0011 %    the peaks, and discarding events with excessive duration.
0012 %    Thresholds are computed as multiples of the standard deviation of
0013 %    the NSS. Alternatively, one can use explicit values, typically obtained
0014 %    from a previous call.
0015 %
0016 %    filtered       ripple-band filtered LFP <a href="matlab:help samples">samples</a> (one channel).
0017 %    <options>      optional list of property-value pairs (see table below)
0018 %
0019 %    =========================================================================
0020 %     Properties    Values
0021 %    -------------------------------------------------------------------------
0022 %     'thresholds'  thresholds for ripple beginning/end and peak, in multiples
0023 %                   of the stdev (default = [2 5])
0024 %     'durations'   minimum inter-ripple interval, and minimum and maximum
0025 %                   ripple durations, in ms (default = [30 20 100])
0026 %     'baseline'    interval used to compute normalization (default = all)
0027 %     'restrict'    same as 'baseline' (for backwards compatibility)
0028 %     'frequency'   sampling rate (in Hz) (default = 1250Hz)
0029 %     'stdev'       reuse previously computed stdev
0030 %     'show'        plot results (default = 'off')
0031 %     'noise'       noisy ripple-band filtered channel used to exclude ripple-
0032 %                   like noise (events also present on this channel are
0033 %                   discarded)
0034 %    =========================================================================
0035 %
0036 %  OUTPUT
0037 %
0038 %    ripples        for each ripple, [start_t peak_t end_t peakNormalizedPower]
0039 %    stdev          standard deviation of the NSS (can be reused subsequently)
0040 %    noise          ripple-like activity recorded simultaneously on the noise
0041 %                   channel (for debugging info)
0042 %
0043 %  SEE
0044 %
0045 %    See also FilterLFP, RippleStats, SaveRippleEvents, PlotRippleStats.
0046 
0047 % Copyright (C) 2004-2011 by Michaël Zugaro, 2016 Ralitsa Todorova (vectorized secondPass),
0048 % initial algorithm by Hajime Hirase
0049 %
0050 % This program is free software; you can redistribute it and/or modify
0051 % it under the terms of the GNU General Public License as published by
0052 % the Free Software Foundation; either version 3 of the License, or
0053 % (at your option) any later version.
0054 
0055 % Default values
0056 frequency = 1250;
0057 show = 'off';
0058 restrict = [];
0059 sd = [];
0060 lowThresholdFactor = 2; % Ripple envoloppe must exceed lowThresholdFactor*stdev
0061 highThresholdFactor = 5; % Ripple peak must exceed highThresholdFactor*stdev
0062 minInterRippleInterval = 30; % in ms
0063 minRippleDuration = 20; % in ms
0064 maxRippleDuration = 100; % in ms
0065 noise = [];
0066 
0067 % Check number of parameters
0068 if nargin < 1 | mod(length(varargin),2) ~= 0,
0069   error('Incorrect number of parameters (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0070 end
0071 
0072 % Check parameter sizes
0073 if ~isdmatrix(filtered) | size(filtered,2) ~= 2,
0074     error('Parameter ''filtered'' is not a Nx2 matrix (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0075 end
0076 
0077 % Parse parameter list
0078 for i = 1:2:length(varargin),
0079     if ~ischar(varargin{i}),
0080         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).']);
0081     end
0082     switch(lower(varargin{i})),
0083         case 'thresholds',
0084             thresholds = varargin{i+1};
0085             if ~isdvector(thresholds,'#2','>0'),
0086                 error('Incorrect value for property ''thresholds'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0087             end
0088             lowThresholdFactor = thresholds(1);
0089             highThresholdFactor = thresholds(2);
0090         case 'durations',
0091             durations = varargin{i+1};
0092             if ~isdvector(durations,'#2','>0') && ~isdvector(durations,'#3','>0'),
0093                 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0094             end
0095             if length(durations) == 2,
0096                 minInterRippleInterval = durations(1);
0097                 maxRippleDuration = durations(2);
0098             else
0099                 minInterRippleInterval = durations(1);
0100                 minRippleDuration = durations(2);
0101                 maxRippleDuration = durations(3);
0102             end
0103         case 'frequency',
0104             frequency = varargin{i+1};
0105             if ~isdscalar(frequency,'>0'),
0106                 error('Incorrect value for property ''frequency'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0107             end
0108         case 'show',
0109             show = varargin{i+1};
0110             if ~isastring(show,'on','off'),
0111                 error('Incorrect value for property ''show'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0112             end
0113         case {'baseline','restrict'},
0114             restrict = varargin{i+1};
0115             if ~isempty(restrict) & ~isdvector(restrict,'#2','<'),
0116                 error('Incorrect value for property ''restrict'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0117             end
0118         case 'stdev',
0119             sd = varargin{i+1};
0120             if ~isdscalar(sd,'>0'),
0121                 error('Incorrect value for property ''stdev'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0122             end
0123         case 'noise',
0124             noise = varargin{i+1};
0125             if ~isdmatrix(noise) | size(noise,1) ~= size(filtered,1) | size(noise,2) ~= 2,
0126                 error('Incorrect value for property ''noise'' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).');
0127             end
0128         otherwise,
0129             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help FindRipples">FindRipples</a>'' for details).']);
0130     end
0131 end
0132 
0133 % Parameters
0134 windowLength = round(frequency/1250*11);
0135 
0136 % Square and normalize signal
0137 time = filtered(:,1);
0138 signal = filtered(:,2);
0139 squaredSignal = signal.^2;
0140 window = ones(windowLength,1)/windowLength;
0141 keep = [];
0142 if ~isempty(restrict),
0143     keep = filtered(:,1)>=restrict(1)&filtered(:,1)<=restrict(2);
0144 end
0145 
0146 [normalizedSquaredSignal,sd] = unity(Filter0(window,sum(squaredSignal,2)),sd,keep);
0147 
0148 % Detect ripple periods by thresholding normalized squared signal
0149 thresholded = normalizedSquaredSignal > lowThresholdFactor;
0150 start = find(diff(thresholded)>0);
0151 stop = find(diff(thresholded)<0);
0152 % Exclude last ripple if it is incomplete
0153 if length(stop) == length(start)-1,
0154     start = start(1:end-1);
0155 end
0156 % Exclude first ripple if it is incomplete
0157 if length(stop)-1 == length(start),
0158     stop = stop(2:end);
0159 end
0160 % Correct special case when both first and last ripples are incomplete
0161 if start(1) > stop(1),
0162     stop(1) = [];
0163     start(end) = [];
0164 end
0165 firstPass = [start,stop];
0166 if isempty(firstPass),
0167     disp('Detection by thresholding failed');
0168     return
0169 else
0170     disp(['After detection by thresholding: ' num2str(length(firstPass)) ' events.']);
0171 end
0172 
0173 % Merge ripples if inter-ripple period is too short (unless this would yield too long a ripple)
0174 secondPass = firstPass;
0175 iri = time(secondPass(2:end,1)) - time(secondPass(1:end-1,2));
0176 duration = time(secondPass(2:end,2)) - time(secondPass(1:end-1,1));
0177 toMerge = iri<minInterRippleInterval/1000 & duration<maxRippleDuration/1000;
0178 while any(toMerge),
0179     % Get indices of first ripples in pairs to be merged
0180     rippleStart = strfind([0 toMerge'],[0 1])';
0181     % Incorporate second ripple into first in all pairs
0182     rippleEnd = rippleStart+1;
0183     secondPass(rippleStart,2) = secondPass(rippleEnd,2);
0184     % Remove second ripples and loop
0185     secondPass(rippleEnd,:) = [];
0186     iri = time(secondPass(2:end,1)) - time(secondPass(1:end-1,2));
0187     duration = time(secondPass(2:end,2)) - time(secondPass(1:end-1,1));
0188     toMerge = iri<minInterRippleInterval/1000 & duration<maxRippleDuration/1000;
0189 end
0190 
0191 if isempty(secondPass),
0192     disp('Ripple merge failed');
0193     return
0194 else
0195     disp(['After ripple merge: ' num2str(length(secondPass)) ' events.']);
0196 end
0197 
0198 % Discard ripples with a peak power < highThresholdFactor
0199 thirdPass = [];
0200 peakNormalizedPower = [];
0201 for i = 1:size(secondPass,1)
0202     maxValue = max(normalizedSquaredSignal([secondPass(i,1):secondPass(i,2)]));
0203     if maxValue > highThresholdFactor,
0204         thirdPass = [thirdPass ; secondPass(i,:)];
0205         peakNormalizedPower = [peakNormalizedPower ; maxValue];
0206     end
0207 end
0208 if isempty(thirdPass),
0209     disp('Peak thresholding failed.');
0210     return
0211 else
0212     disp(['After peak thresholding: ' num2str(length(thirdPass)) ' events.']);
0213 end
0214 
0215 % Detect negative peak position for each ripple
0216 peakPosition = zeros(size(thirdPass,1),1);
0217 for i=1:size(thirdPass,1),
0218     [~,minIndex] = min(signal(thirdPass(i,1):thirdPass(i,2)));
0219     peakPosition(i) = minIndex + thirdPass(i,1) - 1;
0220 end
0221 
0222 % Discard ripples that are way too short
0223 ripples = [time(thirdPass(:,1)) time(peakPosition) time(thirdPass(:,2)) peakNormalizedPower];
0224 duration = ripples(:,3)-ripples(:,1);
0225 ripples(duration<minRippleDuration/1000,:) = [];
0226 disp(['After min duration test: ' num2str(size(ripples,1)) ' events.']);
0227 
0228 % Discard ripples that are way too long
0229 duration = ripples(:,3)-ripples(:,1);
0230 ripples(duration>maxRippleDuration/1000,:) = [];
0231 disp(['After max duration test: ' num2str(size(ripples,1)) ' events.']);
0232 
0233 % If a noisy channel was provided, find ripple-like events and exclude them
0234 bad = [];
0235 if ~isempty(noise),
0236     % Square and pseudo-normalize (divide by signal stdev) noise
0237     squaredNoise = noise(:,2).^2;
0238     window = ones(windowLength,1)/windowLength;
0239     normalizedSquaredNoise = unity(Filter0(window,sum(squaredNoise,2)),sd,[]);
0240     excluded = logical(zeros(size(ripples,1),1));
0241     % Exclude ripples when concomittent noise crosses high detection threshold
0242     previous = 1;
0243     for i = 1:size(ripples,1),
0244         j = FindInInterval(noise,[ripples(i,1),ripples(i,3)],previous);
0245         previous = j(2);
0246         if any(normalizedSquaredNoise(j(1):j(2))>highThresholdFactor),
0247             excluded(i) = 1;
0248         end
0249     end
0250     bad = ripples(excluded,:);
0251     ripples = ripples(~excluded,:);
0252     disp(['After noise removal: ' num2str(size(ripples,1)) ' events.']);
0253 end
0254 
0255 % Optionally, plot results
0256 if strcmp(show,'on'),
0257     figure;
0258     if ~isempty(noise),
0259         MultiPlotXY([time signal],[time squaredSignal],[time normalizedSquaredSignal],[time noise(:,2)],[time squaredNoise],[time normalizedSquaredNoise]);
0260         nPlots = 6;
0261         subplot(nPlots,1,3);
0262          ylim([0 highThresholdFactor*1.1]);
0263         subplot(nPlots,1,6);
0264           ylim([0 highThresholdFactor*1.1]);
0265     else
0266         MultiPlotXY([time signal],[time squaredSignal],[time normalizedSquaredSignal]);
0267 %          MultiPlotXY(time,signal,time,squaredSignal,time,normalizedSquaredSignal);
0268         nPlots = 3;
0269         subplot(nPlots,1,3);
0270           ylim([0 highThresholdFactor*1.1]);
0271     end
0272     for i = 1:nPlots,
0273         subplot(nPlots,1,i);
0274         hold on;
0275           yLim = ylim;
0276         for j=1:size(ripples,1),
0277             plot([ripples(j,1) ripples(j,1)],yLim,'g-');
0278             plot([ripples(j,2) ripples(j,2)],yLim,'k-');
0279             plot([ripples(j,3) ripples(j,3)],yLim,'r-');
0280             if i == 3,
0281                 plot([ripples(j,1) ripples(j,3)],[ripples(j,4) ripples(j,4)],'k-');
0282             end
0283         end
0284         for j=1:size(bad,1),
0285             plot([bad(j,1) bad(j,1)],yLim,'k-');
0286             plot([bad(j,2) bad(j,2)],yLim,'k-');
0287             plot([bad(j,3) bad(j,3)],yLim,'k-');
0288             if i == 3,
0289                 plot([bad(j,1) bad(j,3)],[bad(j,4) bad(j,4)],'k-');
0290             end
0291         end
0292         if mod(i,3) == 0,
0293             plot(xlim,[lowThresholdFactor lowThresholdFactor],'k','linestyle','--');
0294             plot(xlim,[highThresholdFactor highThresholdFactor],'k-');
0295         end
0296     end
0297 end
0298 
0299 function y = Filter0(b,x)
0300 
0301 if size(x,1) == 1,
0302     x = x(:);
0303 end
0304 
0305 if mod(length(b),2) ~= 1,
0306     error('filter order should be odd');
0307 end
0308 
0309 shift = (length(b)-1)/2;
0310 
0311 [y0 z] = filter(b,1,x);
0312 
0313 y = [y0(shift+1:end,:) ; z(1:shift,:)];
0314 
0315 function [U,stdA] = unity(A,sd,restrict)
0316 
0317 if ~isempty(restrict),
0318     meanA = mean(A(restrict));
0319     stdA = std(A(restrict));
0320 else
0321     meanA = mean(A);
0322     stdA = std(A);
0323 end
0324 if ~isempty(sd),
0325     stdA = sd;
0326 end
0327 
0328 U = (A - meanA)/stdA;
0329
FindRipples

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
