RippleStats - Compute descriptive stats for ripples (100~200Hz oscillations).
USAGE
[maps,data,stats] = RippleStats(filtered,ripples,<options>)
filtered ripple-band filtered <a href="matlab:help samples">samples</a> (one channel)
ripples ripple timing information (obtained using <a href="matlab:help FindRipples">FindRipples</a>)
<options> optional list of property-value pairs (see table below)
=========================================================================
Properties Values
-------------------------------------------------------------------------
'frequency' sampling rate (in Hz) (default = 1250Hz)
'durations' durations before and after ripple peak (in s)
(default = [-0.5 0.5])
=========================================================================
OUTPUT
maps.ripples instantaneous amplitude (one ripple per row)
maps.frequency instantaneous frequency
maps.phase instantaneous phase
maps.amplitude enveloppe amplitude
data.peakFrequency frequency at peak
data.peakAmplitude amplitude at peak
data.duration durations
stats.durationAmplitude correlation between duration and amplitude (rho, p)
stats.durationFrequency correlation between duration and frequency (rho, p)
stats.amplitudeFrequency correlation between amplitude and frequency (rho, p)
stats.acg.data autocorrelogram data
stats.acg.t autocorrelogram time bins
SEE
See also FindRipples, SaveRippleEvents, PlotRippleStats.0001 function [maps,data,stats] = RippleStats(filtered,ripples,varargin) 0002 0003 %RippleStats - Compute descriptive stats for ripples (100~200Hz oscillations). 0004 % 0005 % USAGE 0006 % 0007 % [maps,data,stats] = RippleStats(filtered,ripples,<options>) 0008 % 0009 % filtered ripple-band filtered <a href="matlab:help samples">samples</a> (one channel) 0010 % ripples ripple timing information (obtained using <a href="matlab:help FindRipples">FindRipples</a>) 0011 % <options> optional list of property-value pairs (see table below) 0012 % 0013 % ========================================================================= 0014 % Properties Values 0015 % ------------------------------------------------------------------------- 0016 % 'frequency' sampling rate (in Hz) (default = 1250Hz) 0017 % 'durations' durations before and after ripple peak (in s) 0018 % (default = [-0.5 0.5]) 0019 % ========================================================================= 0020 % 0021 % OUTPUT 0022 % 0023 % maps.ripples instantaneous amplitude (one ripple per row) 0024 % maps.frequency instantaneous frequency 0025 % maps.phase instantaneous phase 0026 % maps.amplitude enveloppe amplitude 0027 % data.peakFrequency frequency at peak 0028 % data.peakAmplitude amplitude at peak 0029 % data.duration durations 0030 % stats.durationAmplitude correlation between duration and amplitude (rho, p) 0031 % stats.durationFrequency correlation between duration and frequency (rho, p) 0032 % stats.amplitudeFrequency correlation between amplitude and frequency (rho, p) 0033 % stats.acg.data autocorrelogram data 0034 % stats.acg.t autocorrelogram time bins 0035 % 0036 % SEE 0037 % 0038 % See also FindRipples, SaveRippleEvents, PlotRippleStats. 0039 0040 % Copyright (C) 2004-2011 by Michaël Zugaro 0041 % 0042 % This program is free software; you can redistribute it and/or modify 0043 % it under the terms of the GNU General Public License as published by 0044 % the Free Software Foundation; either version 3 of the License, or 0045 % (at your option) any later version. 0046 0047 % Default values 0048 samplingRate = 1250; 0049 durations = [-50 50]/1000; 0050 nCorrBins = 51; 0051 % corrBinSize = 400; 0052 corrDuration = 20; 0053 corrBinSize = 0.1; 0054 0055 % Check number of parameters 0056 if nargin < 2, 0057 error('Incorrect number of parameters (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).'); 0058 end 0059 0060 % Check parameter sizes 0061 if size(filtered,2) ~= 2, 0062 error('Parameter ''filtered'' is not a Nx2 matrix (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).'); 0063 end 0064 if size(ripples,2) ~= 4, 0065 error('Parameter ''ripples'' is not a Nx4 matrix (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).'); 0066 end 0067 0068 % Parse parameter list 0069 for i = 1:2:length(varargin), 0070 if ~ischar(varargin{i}), 0071 error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).']); 0072 end 0073 switch(lower(varargin{i})), 0074 case 'frequency', 0075 samplingRate = varargin{i+1}; 0076 if ~isdscalar(samplingRate,'>0'), 0077 error('Incorrect value for property ''frequency'' (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).'); 0078 end 0079 case 'durations', 0080 durations = varargin{i+1}; 0081 if ~isdvector(durations,'#2','<'), 0082 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).'); 0083 end 0084 otherwise, 0085 error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help RippleStats">RippleStats</a>'' for details).']); 0086 end 0087 end 0088 0089 nBins = floor(samplingRate*diff(durations)/2)*2+1; % must be odd 0090 nHalfCenterBins = 3; 0091 centerBin = ceil(nBins/2); 0092 centerBins = centerBin-nHalfCenterBins:centerBin+nHalfCenterBins; 0093 0094 % Compute instantaneous phase and amplitude 0095 h = hilbert(filtered(:,2)); 0096 phase(:,1) = filtered(:,1); 0097 phase(:,2) = angle(h); 0098 amplitude(:,1) = filtered(:,1); 0099 amplitude(:,2) = abs(h); 0100 unwrapped(:,1) = filtered(:,1); 0101 unwrapped(:,2) = unwrap(phase(:,2)); 0102 % Compute instantaneous frequency 0103 frequency = Diff(unwrapped,'smooth',0); 0104 frequency(:,2) = frequency(:,2)/(2*pi); 0105 0106 % Compute ripple map 0107 [r,i] = Sync(filtered,ripples(:,2),'durations',durations); 0108 maps.ripples = SyncMap(r,i,'durations',durations,'nbins',nBins,'smooth',0); 0109 0110 % Compute frequency Map 0111 [f,i] = Sync(frequency,ripples(:,2),'durations',durations); 0112 maps.frequency = SyncMap(f,i,'durations',durations,'nbins',nBins,'smooth',0); 0113 0114 % Compute phase map 0115 [p,i] = Sync(phase,ripples(:,2),'durations',durations); 0116 maps.phase = SyncMap(p,i,'durations',durations,'nbins',nBins,'smooth',0); 0117 0118 % Compute amplitude map 0119 [a,i] = Sync(amplitude,ripples(:,2),'durations',durations); 0120 maps.amplitude = SyncMap(a,i,'durations',durations,'nbins',nBins,'smooth',0); 0121 0122 % Ripple frequency and amplitude at peak 0123 data.peakFrequency = maps.frequency(:,centerBin); 0124 data.peakAmplitude = maps.amplitude(:,centerBin); 0125 0126 % Ripple durations 0127 data.duration = ripples(:,3)-ripples(:,1); 0128 0129 % Autocorrelogram and correlations 0130 if nargin >= 3, 0131 [stats.acg.data,stats.acg.t] = CCG(ripples(:,2),1,'binSize',corrBinSize,'nBins',nCorrBins); 0132 [stats.amplitudeFrequency.rho,stats.amplitudeFrequency.p] = corrcoef(data.peakAmplitude,data.peakFrequency); 0133 [stats.durationFrequency.rho,stats.durationFrequency.p] = corrcoef(data.duration,data.peakFrequency); 0134 [stats.durationAmplitude.rho,stats.durationAmplitude.p] = corrcoef(data.duration,data.peakAmplitude); 0135 end