README (last revised: 05 November 2012) ======================================= Matlab scripts for ATvDFA as published in Journal of Neuroscience Methods 209:178–188, 2012. If you use those scripts and publish the results, please cite the paper (and drop me a note!). Please note that this is experimental code and is provided with no warranty! It is far from being optimised and has very limited error checking within it. I welcome any feedback / comments / bug reports. Archive contents: ----------------- - dfa.m and dfaData2.m are two functions that compute detrended fluctuation analysis. There are historical reasons why both of these functions exist. In fact, they are doing the same thing but have different headers. Sorry about this (to be fixed). Note that the code is based on the code by McSharry (used to reside at: http://www.eng.ox.ac.uk/samp/software/cardiodynamics/dfa.m"). - rsquare.m is called by dfa and computes the linear regression of the fluctuations in the loglog space. - progressbar.m is only used if setting pflag=1 in odfat9ckf6.m (see below) and provides a visual indicator of progress (computations can be long!). - ekfukf_1_3.tar.gz contains a publicly available Kalman filtering toolbox (available under GNU GPL 2 from:http://becs.aalto.fi/en/research/bayes/ekfukf/). This was not developed by me and is only included in full for convenience. Only a few of its functions are actually used by ATvDFA. Note that the authors of this toolbox may release newer versions of it and that compatibility of my code with any newer version is not to be assumed. - odfat9ckf6.m is the ATvDFA method (see below for detailed notes regarding its inputs and outputs). How to use the method: ---------------------- The function has the following header: function [expArray,r2Array,ovslope,ovr2,seq,lns,xV,xVs,xVs2]=odfat9ckf6(data,fs,winlength,inc,nmeas,mbs,pflag,q) Inputs: - data: the time series to be analysed. Note this function does not take multivariate time series. - fs: the sampling frequency - winlength: the length (in seconds, not in number of samples) of the window over which DFA estimates are obtained. If this is set to the full length of the data, then this is identical to applying DFA to the time series (Note that this is done in the function anyway). - inc: the increment (in seconds) between the windows over which DFA estimates are obtained. Because the method computes DFA a number of times (see nmeas below) within this increment, this increment should be kept large enough. I have typically used inc=1sec, which means that I will get a new estimate of DFA exponent every second. There is no problem (other than computational) making it smaller but of course one should expect the exponent to be pretty constant within a very short period (in fact, it is a fundamental assumption of the construction of the filter). The increment can be made larger obviously but this means lowering the temporal resolution of the tracked estimates and, more importantly, the above assumption may be violated, e.g., the exponents may change a lot within 3 seconds for example. - nmeas: the number of DFA estimates collected within the increment period in order to build an estimate of the covariance of the fluctuations (assuming a constant exponent). The higher the number, the better the estimate but the higher the computational cost. I typically use 10 which means that if inc is 1sec, a DFA calculation is performed 10 times per second, each time shifting the time series by 100ms. - mbs: the min box size (in seconds) for the DFA calculation. This is an oft-ignored free parameter of the DFA method. One wants to set mbs as small as possible (to increase the range over which long-range temporal correlations are assessed) but not too small that it varies with the signal itself. For example, for an oscillatory signals, one would want a sufficiently high number of cycles and one should beware of the relationship to any filter used on the signal. For example a 8-12Hz bandpass filter could be designed so that the filter used has a length of at least 3 cycles of the slowest frequency. This would mean 3*1/8s = 375ms. It would therefore be inappropriate to set mbs to a value less than 0.5s. With higher frequencies, this is less of an issue and one could bring mbs down to 0.25s. Note that strictly speaking the max box size should also be a free parameter (this is discussed in the paper). Typically it is recommended to use 1/10 of the data. However, as ATvDFA applies DFA to short segments (of length winlength), using 1/10 can be too restrictive. In the code provided, the value is hard-coded to 4 (line 63 of the code). This means that for winlength=20, one would be looking at LRTCs over [mbs,5]s. - pflag: use 1 for graphically verbose output, 0 otherwise. - q is a 2-component array and contains the remaining free parameters of the method (see paper). A default value of [0.000001,0.000001] is provided. The choice of this value will typically depend on the choice of winlength and the constancy of the exponent within the record. At this point in time, the authors do not have a perfect recipe for setting this up. Outputs (I must apologise for the very non-intuitive naming of the output variables) - expArray: the sequence of DFA exponents obtained using moving DFA, i.e., *without* Kalman filtering. In the paper, we refer to it as mDFA (purple lines in plots). - r2Array: the r2 values associated with each of the above exponents (see paper) - ovslope: the DFA exponent over the entire time series (i.e., the standard application of DFA when one does not expect the exponent to change within the record) - ovr2: the r2 value associated with the above exponent. - seq: the starting index for each window over which DFA estimates were obtained (not including the within-inc measurements). This is mostly useful for plotting. - lns: the log box sizes. This is mostly useful for plotting. - xV: the Kalman filtered sequence of DFA exponents without smoothing. In the paper, this corresponds to the blue line in our plots. - xVs: the Kalman filtered sequence of DFA exponents with Rauch-Tung-Striebel smoother (not shown in paper). - xVs2: the Kalman filtered sequence of DFA exponents with a two-filter based smoother. In the paper, this is shown by the yellow line and is meant to be *the* method. Feel free to get in touch should you encounter any difficulty. Obviously we would be happy to hear of any results/findings using this method. Luc Berthouze 05 November 2012