0001 function [maps,data,stats] = RippleStats(filtered,ripples,varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048 samplingRate = 1250;
0049 durations = [-50 50]/1000;
0050 nCorrBins = 51;
0051
0052 corrDuration = 20;
0053 corrBinSize = 0.1;
0054
0055
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
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
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;
0090 nHalfCenterBins = 3;
0091 centerBin = ceil(nBins/2);
0092 centerBins = centerBin-nHalfCenterBins:centerBin+nHalfCenterBins;
0093
0094
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
0103 frequency = Diff(unwrapped,'smooth',0);
0104 frequency(:,2) = frequency(:,2)/(2*pi);
0105
0106
0107 [r,i] = Sync(filtered,ripples(:,2),'durations',durations);
0108 maps.ripples = SyncMap(r,i,'durations',durations,'nbins',nBins,'smooth',0);
0109
0110
0111 [f,i] = Sync(frequency,ripples(:,2),'durations',durations);
0112 maps.frequency = SyncMap(f,i,'durations',durations,'nbins',nBins,'smooth',0);
0113
0114
0115 [p,i] = Sync(phase,ripples(:,2),'durations',durations);
0116 maps.phase = SyncMap(p,i,'durations',durations,'nbins',nBins,'smooth',0);
0117
0118
0119 [a,i] = Sync(amplitude,ripples(:,2),'durations',durations);
0120 maps.amplitude = SyncMap(a,i,'durations',durations,'nbins',nBins,'smooth',0);
0121
0122
0123 data.peakFrequency = maps.frequency(:,centerBin);
0124 data.peakAmplitude = maps.amplitude(:,centerBin);
0125
0126
0127 data.duration = ripples(:,3)-ripples(:,1);
0128
0129
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