Description of PlotRippleStats

0001 function PlotRippleStats(r,m,d,s,varargin)
0002 
0003 %PlotRippleStats - Plot descriptive stats for ripples (100~200Hz oscillations).
0004 %
0005 %  USAGE
0006 %
0007 %    PlotRippleStats(ripples,maps,data,stats,<options>)
0008 %
0009 %    Use the ouputs of <a href="matlab:help RippleStats">RippleStats</a> as input parameters. Using cell arrays of
0010 %    repeated measures will yield pairwise statistics.
0011 %
0012 %    ripples        ripple timing info (provided by <a href="matlab:help FindRipples">FindRipples</a>)
0013 %    maps           ripple instantaneous amplitude, frequency, and phase maps
0014 %    data           frequency and amplitude at peak, durations
0015 %    stats          autocorrelogram and correlations
0016 %    <options>      optional list of property-value pairs (see table below)
0017 %
0018 %    =========================================================================
0019 %     Properties    Values
0020 %    -------------------------------------------------------------------------
0021 %     'frequency'   sampling rate (in Hz) (default = 1250Hz)
0022 %     'durations'   durations before and after ripple peak (in s)
0023 %                   (default = [-0.5 0.5])
0024 %    =========================================================================
0025 %
0026 %  SEE
0027 %
0028 %    See also FindRipples, RippleStats, SaveRippleEvents.
0029 
0030 % Copyright (C) 2004-2011 by Michaël Zugaro
0031 %
0032 % This program is free software; you can redistribute it and/or modify
0033 % it under the terms of the GNU General Public License as published by
0034 % the Free Software Foundation; either version 3 of the License, or
0035 % (at your option) any later version.
0036 
0037 % Default values
0038 samplingRate = 1250;
0039 durations = [-50 50]/1000;
0040 nCorrBins = 1000;
0041 
0042 % Check number of parameters
0043 if nargin < 3,
0044   error('Incorrect number of parameters (type ''help <a href="matlab:help PlotRippleStats">PlotRippleStats</a>'' for details).');
0045 end
0046 
0047 % Parse parameter list
0048 for i = 1:2:length(varargin),
0049     if ~ischar(varargin{i}),
0050         error(['Parameter ' num2str(i+2) ' is not a property (type ''help <a href="matlab:help PlotRippleStats">PlotRippleStats</a>'' for details).']);
0051     end
0052     switch(lower(varargin{i})),
0053         case 'frequency',
0054             samplingRate = varargin{i+1};
0055             if ~isdscalar(samplingRate,'>0'),
0056                 error('Incorrect value for property ''frequency'' (type ''help <a href="matlab:help PlotRippleStats">PlotRippleStats</a>'' for details).');
0057             end
0058         case 'durations',
0059             durations = varargin{i+1};
0060             if ~isdvector(durations,'#2','<'),
0061                 error('Incorrect value for property ''durations'' (type ''help <a href="matlab:help PlotRippleStats">PlotRippleStats</a>'' for details).');
0062             end
0063         otherwise,
0064             error(['Unknown property ''' num2str(varargin{i}) ''' (type ''help <a href="matlab:help PlotRippleStats">PlotRippleStats</a>'' for details).']);
0065     end
0066 end
0067 
0068 if isa(m,'cell'),
0069     ripples = r;
0070     maps = m;
0071     data = d;
0072     stats = s;
0073 else
0074     ripples{1} = r;
0075     maps{1} = m;
0076     data{1} = d;
0077     stats{1} = s;
0078 end
0079 nElectrodes = length(maps);
0080 
0081 %  nBins = floor(samplingRate*diff(durations)/2)*2+1; % must be odd
0082 nBins =     length(maps{1}.ripples);
0083 nHalfCenterBins = 3;
0084 %  centerBin = durations(1)/sum(durations)*(nBins-1)+1;
0085 centerBin = ceil(nBins/2);
0086 centerBins = centerBin-nHalfCenterBins:centerBin+nHalfCenterBins;
0087 
0088 % Stats for individual electrodes
0089 
0090 for electrode = 1:nElectrodes,
0091 
0092     % Unsorted color plots
0093     f = figure;
0094     set(f,'name',['Ripples (unsorted) - ' int2str(electrode)]);
0095     dx = diff(durations);
0096     x = durations(1):dx/size(maps{electrode}.ripples,2):durations(2);
0097     subplot(2,2,1);PlotColorMap(maps{electrode}.ripples,1,'bar','on','cutoffs',[-1000 1000],'x',x);xlabel('Ripples');
0098     subplot(2,2,2);PlotColorMap(maps{electrode}.phase,1,'bar','on','type','circular','x',x);xlabel('Ripple Phase');
0099     subplot(2,2,3);PlotColorMap(maps{electrode}.frequency,1,'bar','on','cutoffs',[100 250],'x',x);xlabel('Ripple Frequency');
0100     subplot(2,2,4);PlotColorMap(maps{electrode}.amplitude,1,'bar','on','x',x);xlabel('Ripple Amplitude'); %,'cutoffs',[0 3500]
0101 
0102     % Sorted color plots (sort by ripple frequency)
0103     %  f = figure;
0104     %  set(f,'name',['Ripples (sorted by frequency) - ' int2str(electrode)]);
0105     %  subplot(2,2,1);PlotColorMap(maps.ripples{electrode}(order,:),1,'bar','on');xlabel('Ripples');
0106     %  subplot(2,2,2);PlotColorMap(maps.phase{electrode}(order,:),1,'bar','on','type','circular');xlabel('Ripple Phase');
0107     %  subplot(2,2,3);PlotColorMap(maps.frequency{electrode}(order,:),1,'bar','on','cutoffs',[100 250]);xlabel('Ripple Frequency');
0108     %  subplot(2,2,4);PlotColorMap(maps.amplitude{electrode}(order,:),1,'bar','on','cutoffs',[0 3500]);xlabel('Ripple Amplitude');
0109 
0110     % Ripples stats: ripples, peak frequency vs amplitude, autocorrelogram, peak frequency vs duration, peak amplitude vs duration
0111     f = figure;set(f,'name',['Ripple Stats - ' int2str(electrode)]);
0112 
0113     subplot(2,2,1);a = gca;hold on;
0114     plot(((1:nBins)'-ceil(nBins/2))/nBins*diff(durations),maps{electrode}.ripples,'b');
0115 
0116     subplot(2,2,2);
0117     b = bar(stats{electrode}.acg.t,stats{electrode}.acg.data);set(b,'FaceColor',[0 0 0]);xlabel('Autocorrelogram');
0118 %      b = bar(((0:nCorrBins)-nCorrBins/2)/1000,stats{electrode}.acg.data);xlim([-nCorrBins nCorrBins]/2000);set(b,'FaceColor',[0 0 0]);xlabel('Autocorrelogram');
0119 
0120     subplot(2,3,4);a = gca;
0121     PlotDistribution2(data{electrode}.peakAmplitude,data{electrode}.peakFrequency,'nbins',1000); %,'smooth',5
0122     axes(a);xlabel(['r=' num2str(stats{electrode}.amplitudeFrequency.rho(1,2)) ' p=' num2str(stats{electrode}.amplitudeFrequency.p(1,2))]);ylabel('Frequency vs Amplitude');
0123 
0124     subplot(2,3,5);a = gca;
0125     PlotDistribution2(data{electrode}.duration,data{electrode}.peakFrequency,'nbins',1000); %,'smooth',5
0126     axes(a);xlabel(['r=' num2str(stats{electrode}.durationFrequency.rho(1,2)) ' p=' num2str(stats{electrode}.durationFrequency.p(1,2))]);ylabel('Frequency vs Duration');
0127 
0128     subplot(2,3,6);a = gca;
0129     PlotDistribution2(data{electrode}.duration,data{electrode}.peakAmplitude,'nbins',1000); %,'smooth',5
0130     axes(a);xlabel(['r=' num2str(stats{electrode}.durationAmplitude.rho(1,2)) ' p=' num2str(stats{electrode}.durationAmplitude.p(1,2))]);ylabel('Amplitude vs Duration');
0131 end
0132 
0133 % Pairwise stats
0134 
0135 if nElectrodes == 2,
0136 
0137         % Match ripple pairs
0138         ripplePairs = MatchPairs(ripples{1}(:,2),ripples{2}(:,2),'error',0.01); % within 10 ms
0139         % Ripples on LFP 1 only
0140         only1 = isnan(ripplePairs(:,2));
0141         % Ripples on LFP 2 only
0142         only2 = isnan(ripplePairs(:,1));
0143         % Ripples on both electrodes
0144         both = ~only1 & ~only2;
0145         % Peak time on LFP 1 when there are ripples on both electrodes
0146         ripplePairTimes = ripplePairs(both,1);
0147         % For each ripple on LFP 1, whether there are ripples on both electrodes
0148         ripplesOn1AlsoOn2 = ~isnan(ripplePairs(~only2,2));
0149         % For each ripple on LFP 2, whether there are ripples on both electrodes
0150         ripplesOn2AlsoOn1 = ~isnan(ripplePairs(~only1,1));
0151         % For each ripple on LFP 1, whether there are no ripples on LFP 2
0152         ripplesOn1NotOn2 = isnan(ripplePairs(~only2,2));
0153         % For each ripple on LFP 2, whether there are no ripples on LFP 1
0154         ripplesOn2NotOn1 = isnan(ripplePairs(~only1,1));
0155 
0156 %          % Compute phase shift map
0157 %          size(maps{1}.phase)
0158 %          [p1,i1] = Sync(maps{1}.phase,ripplePairTimes,'durations',durations);
0159 %          [p2,i2] = Sync(maps{2}.phase,ripplePairTimes,'durations',durations);
0160 %          shift(:,1) = p1(:,1);
0161 %          shift(:,2) = p1(:,2)-p2(:,2);
0162 %          phaseShiftMap = SyncMap(shift,i1,'durations',durations,'nbins',nBins,'smooth',0);
0163 %          phaseShiftMap(phaseShiftMap<-pi) = phaseShiftMap(phaseShiftMap<-pi) + 2*pi;
0164 %          phaseShiftMap(phaseShiftMap>pi) = phaseShiftMap(phaseShiftMap>pi) - 2*pi;
0165 %
0166 %          % Determine average phase shift using center bins
0167 %          meanCenterPhaseShift = mean(phaseShiftMap(:,centerBins),2);
0168 %          [~,order] = sortrows(meanCenterPhaseShift);
0169 
0170         % Plot phase difference, frequency and amplitude (1 vs 2)
0171         f = figure;set(f,'name','Pairwise Stats');
0172 
0173 %          subplot(2,2,1);PlotColorMap(phaseShiftMap,1,'bar','on','type','circular');xlabel('Phase Shift');
0174 %          hold on;plot([1 1]*size(phaseShiftMap,2)/2,ylim,'w');
0175 %
0176 %          subplot(2,2,2);PlotColorMap(phaseShiftMap(order,:),1,'bar','on','type','circular');xlabel('Phase Shift');
0177 %          hold on;plot([1 1]*size(phaseShiftMap,2)/2,ylim,'w');
0178 
0179         subplot(2,3,4);a = gca;
0180         [rho,p] = corrcoef(data{1}.peakFrequency(ripplesOn1AlsoOn2),data{2}.peakFrequency(ripplesOn2AlsoOn1));
0181         PlotDistribution2(data{1}.peakFrequency(ripplesOn1AlsoOn2),data{2}.peakFrequency(ripplesOn2AlsoOn1),'nbins',1000,'smooth',50);
0182         axes(a);xlabel(['r=' num2str(rho(1,2)) ' p=' num2str(p(1,2))]);ylabel('Frequency (1 vs 2)');
0183 
0184         subplot(2,3,5);a = gca;
0185         [rho,p] = corrcoef(data{1}.peakAmplitude(ripplesOn1AlsoOn2),data{2}.peakAmplitude(ripplesOn2AlsoOn1));
0186         PlotDistribution2(data{1}.peakAmplitude(ripplesOn1AlsoOn2),data{2}.peakAmplitude(ripplesOn2AlsoOn1),'nbins',1000,'smooth',50);
0187         axes(a);xlabel(['r=' num2str(rho(1,2)) ' p=' num2str(p(1,2))]);ylabel('Amplitude (1 vs 2)');
0188 
0189         subplot(2,3,6);a = gca;
0190         rippleDurations1 = data{1}.duration(ripplesOn1AlsoOn2);
0191         rippleDurations2 = data{2}.duration(ripplesOn2AlsoOn1);
0192         discard = rippleDurations1>0.05|rippleDurations2>0.05;
0193         rippleDurations1 = rippleDurations1(~discard);
0194         rippleDurations2 = rippleDurations2(~discard);
0195         [rho,p] = corrcoef(rippleDurations1,rippleDurations2);
0196         PlotDistribution2(rippleDurations1,rippleDurations2,'nbins',1000,'smooth',50);
0197         axes(a);xlabel(['r=' num2str(rho(1,2)) ' p=' num2str(p(1,2))]);ylabel('Duration (1 vs 2)');
0198         % Plot info for unilateral vs bilateral ripples
0199         f = figure;set(f,'name','Unilateral vs Bilateral Stats');
0200 
0201         subplot(2,2,1);
0202         S_boxplot([data{1}.peakFrequency(ripplesOn1NotOn2);data{1}.peakFrequency(ripplesOn1AlsoOn2);data{2}.peakFrequency(ripplesOn2NotOn1);data{2}.peakFrequency(ripplesOn2AlsoOn1)],[ones(sum(ripplesOn1NotOn2),1);2*ones(sum(ripplesOn1AlsoOn2),1);3*ones(sum(ripplesOn2NotOn1),1);4*ones(sum(ripplesOn2AlsoOn1),1)]);
0203         p1 = S_ranksum(data{1}.peakFrequency(ripplesOn1NotOn2),data{1}.peakFrequency(ripplesOn1AlsoOn2));
0204         p2 = S_ranksum(data{2}.peakFrequency(ripplesOn2NotOn1),data{2}.peakFrequency(ripplesOn2AlsoOn1));
0205         xlabel(['p=' num2str(p1) ' N=' num2str(sum(ripplesOn1NotOn2)) ',' num2str(sum(ripplesOn1AlsoOn2)) ' p=' num2str(p2) ' N=' num2str(sum(ripplesOn2NotOn1)) ',' num2str(sum(ripplesOn2AlsoOn1))]);
0206         ylabel('Frequency (1-2,1+2,2-1,2+1)');
0207 
0208         subplot(2,2,2);
0209         S_boxplot([data{1}.peakAmplitude(ripplesOn1NotOn2);data{1}.peakAmplitude(ripplesOn1AlsoOn2);data{2}.peakAmplitude(ripplesOn2NotOn1);data{2}.peakAmplitude(ripplesOn2AlsoOn1)],[ones(sum(ripplesOn1NotOn2),1);2*ones(sum(ripplesOn1AlsoOn2),1);3*ones(sum(ripplesOn2NotOn1),1);4*ones(sum(ripplesOn2AlsoOn1),1)]);
0210         p1 = S_ranksum(data{1}.peakAmplitude(ripplesOn1NotOn2),data{1}.peakAmplitude(ripplesOn1AlsoOn2));
0211         p2 = S_ranksum(data{2}.peakAmplitude(ripplesOn2NotOn1),data{2}.peakAmplitude(ripplesOn2AlsoOn1));
0212         xlabel(['p=' num2str(p1) ' N=' num2str(sum(ripplesOn1NotOn2)) ',' num2str(sum(ripplesOn1AlsoOn2)) ' p=' num2str(p2) ' N=' num2str(sum(ripplesOn2NotOn1)) ',' num2str(sum(ripplesOn2AlsoOn1))]);
0213         ylabel('Amplitude (1-2,1+2,2-1,2+1)');
0214 
0215         subplot(2,2,3);
0216         S_boxplot([data{1}.duration(ripplesOn1NotOn2);data{1}.duration(ripplesOn1AlsoOn2);data{2}.duration(ripplesOn2NotOn1);data{2}.duration(ripplesOn2AlsoOn1)],[ones(sum(ripplesOn1NotOn2),1);2*ones(sum(ripplesOn1AlsoOn2),1);3*ones(sum(ripplesOn2NotOn1),1);4*ones(sum(ripplesOn2AlsoOn1),1)]);
0217         p1 = S_ranksum(data{1}.duration(ripplesOn1NotOn2),data{1}.duration(ripplesOn1AlsoOn2));
0218         p2 = S_ranksum(data{2}.duration(ripplesOn2NotOn1),data{2}.duration(ripplesOn2AlsoOn1));
0219         xlabel(['p=' num2str(p1) ' N=' num2str(sum(ripplesOn1NotOn2)) ',' num2str(sum(ripplesOn1AlsoOn2)) ' p=' num2str(p2) ' N=' num2str(sum(ripplesOn2NotOn1)) ',' num2str(sum(ripplesOn2AlsoOn1))]);
0220         ylabel('Duration (1-2,1+2,2-1,2+1)');
0221 end
PlotRippleStats

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
