Home > FMAToolbox > Analyses > RippleStats.m

RippleStats

PURPOSE ^

RippleStats - Compute descriptive stats for ripples (100~200Hz oscillations).

SYNOPSIS ^

function [maps,data,stats] = RippleStats(filtered,ripples,varargin)

DESCRIPTION ^

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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 16-Mar-2018 13:00:20 by m2html © 2005