Model-Free Timescale Estimation
Performed via the function acw
in IntrinsicTimescales.jl. The acw
function calculates ACF or PSD depending on the acwtypes you specify. If there is no missing data (indicated by NaN
or missing
), acw
calculates ACF as the inverse fourier transform of the power spectrum, using comp_ac_fft
internally. Otherwise it calculates ACF as correlations between a time-series and its lag-shifted variants, using comp_ac_time_missing
. For PSD, it uses periodogram method (comp_psd
) in the case of no missing data and Lomb-Scargle method (comp_psd_lombscargle
) in the case of missing data.
acwresults = acw(data, fs; acwtypes=[:acw0, :acw50, :acweuler, :auc, :tau, :knee],
n_lags=nothing, freqlims=nothing, dims=ndims(data),
return_acf=true, return_psd=true,
average_over_trials=false, trial_dims=setdiff([1, 2], dims)[1],
skip_zero_lag=false,
max_peaks=1, oscillation_peak::Bool=true,
allow_variable_exponent::Bool=false,
constrained::Bool=false,
parallel::Bool=false)
Simple usage:
results = acw(data, fs)
acw_results = results.acw_results
acw_0 = acw_results[1]
acw_50 = acw_result[2]
# And so on...
Arguments
Mandatory arguments:
data
: Your time-series data as a vector or n-dimensional array.
If it is n-dimensional, by default, the dimension of time is assumed to be the last dimension. For example, if you have a 2D array where rows are subjects and columns are time points, acw
function will correctly assume that the last (2nd) dimension is time. If the dimension of time is any other dimension than the last, you can set it via dims
argument. For example:
data = randn(100, 200, 5) # 100 trials, 200 time points, 5 channels
acw(data, fs; dims=2)
fs
: Sampling rate of your data. A floating point number.
Optional arguments:
acwtypes
: A symbol or a vector of symbols denoting which ACW types to calculate.
In Julia, a symbol is written as :symbol
. Example:
acw(data, fs; acwtypes=:acw0)
acw(data, fs; acwtypes=[:acw0, :acw50])
Supported ACW types:
:acw0
: The lag where autocorrelation function crosses 0.
:acw50
: The lag where autocorrelation function crosses 0.5.
:acweuler
: The lag where autocorrelation function crosses $1/e$. Corresponds to the inverse decay rate of an exponential decay function.
:tau
: Fit an exponential decay function $e^{\frac{t}{\tau}}$ to the autocorrelation function and extract $\tau$, which is the inverse decay rate. The parameter skip_zero_lag
is used to specify whether to skip the zero lag for fitting an exponential decay function. When zero-lag is skipped, the function fits an exponential decay of the form $A (exp(-lags / tau) + B)$ where A is the amplitude and B is the offset. See below for details and references.
:auc
: Calculate the area under the autocorrelation function from lag 0 to the lag where autocorrelation function crosses 0.
:knee
: Fit a lorentzian function $\frac{A}{1 + (f/a)^2}$ to the power spectrum using an iterative FOOOF-style approach. By Wiener-Khinchine theorem, this is the power spectrum of a time-series with an autocorrelation function of exponential decay form. The parameter $a$ corresponds to the knee frequency. $\tau$ and $a$ has the relationship $\tau = \frac{1}{2 \pi a}$. The :knee
method uses this relationship to estimate $\tau$ from the knee frequency. In practice, first, an initial lorentzian fit is performed. Then, any oscillatory peaks are identified and fitted with gaussian functions. These gaussians are subtracted from the original power spectrum to ensure the remaining PSD is closer to a Lorentzian, and a final Lorentzian is fit to this "cleaned" spectrum. freqlims
is used to specify the frequency limits to fit the lorentzian function. You can set the maximum number of oscillatory peaks to fit with the max_peaks
argument. The argument oscillation_peak
is used to specify whether to fit the oscillatory peaks or not. If set to false
, just fit a Lorentzian and return the timescale estimated from the knee frequency. The parameter allow_variable_exponent
is used to specify whether to allow a variable exponent in the lorentzian fit (i.e. $b$ in $\frac{A}{1 + (f/a)^b}$). This might be useful for cases where the power spectrum does not conform to a simple Lorentzian. The parameter constrained
is used to specify whether to use constrained optimization when fitting a lorentzian to the PSD for knee frequency estimation (see below).
n_lags
: An integer. Only used when:tau
is inacwtypes
. The number of lags to be used for fitting an exponential decay function.
By default, this is set to 1.1*acw0
. The reason for cutting the autocorrelation function is due to increased measurement noise for autocorrelation function for later lags. Intuitively, when we perform correlation of a time-series with a shifted version of itself, we have less and less number of time points to calculate the correlation. This is why if you plot the autocorrelation function, the portion of it after ACW-0 looks more noisy. For more details, see [Practice 1].
freqlims
: Only used when:knee
is inacwtypes
. The frequency limits to fit the lorentzian function.
By default, the lowest and highest frequencies that can be estimated from your data. A tuple of two floating point numbers, for example, (freq_low, freq_high)
or (1.0, 50.0)
.
dims
: The dimension of time in your data. See above, thedata
argument for the explanation. An integer.return_acf
: Whether or not to return the autocorrelation function (ACF) in the results object. A boolean.return_psd
: Whether or not to return the power spectrum (psd) in the results object. A boolean.average_over_trials
: Whether or not to average the ACF or PSD across trials, as in [Honey et al., 2012].
Assuming that your data is stationary, averaging over trials can greatly reduce noise in your ACF/PSD estimations. By default, the dimension of trials is assumed to be the first dimension of your data. For example, if your data is two dimensional with rows as trials and columns as time points, the function will correctly infer the dimension of trials. If this is not the case, set the dimension of trials with the argument trial_dims
. Below is an example of a three dimensional data with time as second and trials as third argument:
data = randn(10, 1000, 20) # 10 subjects, 1000 time points, 20 trials
result = acw(data, fs; dims=2, average_over_trials=true, trial_dims=3)
trial_dims
: Dimension of trials to average over. See above (average_over_trials
) for explanation. An integer.skip_zero_lag
: Whether or not to skip the zero lag for fitting an exponential decay function. Default isfalse
. If true, the function will fit an exponential decay of the form $A (exp(-lags / tau) + B)$ where A is the amplitude and B is the offset. This can be useful for cases with very low sampling rate (e.g. fMRI). The technique is used in Ito et al., 2020 and Murray et al., 2014.max_peaks
: Maximum number of oscillatory peaks to fit when cleaning the PSD for knee frequency estimation. Default is 1.oscillation_peak
: Whether or not to fit the oscillatory peaks when cleaning the PSD for knee frequency estimation. Default istrue
.allow_variable_exponent
: Whether or not to allow variable exponent (PLE) when fitting a lorentzian to the PSD for knee frequency estimation. Default isfalse
. If true, the function will admit Lorentzian's of form $\frac{A}{1 + (f/a)^b}$ where $b$ is not confined to -2.constrained
: Whether or not to use constrained optimization when fitting a lorentzian to the PSD for knee frequency estimation. Default isfalse
. If true, the function will use constrained optimization via Optim.jl (using Optimization.jl as a frontend). The lower constraints for amplitude, knee frequency and exponent are 0, freqs[1], 0 respectively. The upper constraints are Inf, freqs[end], 5.0. For optimization, LBFGS method is used.parallel
: Whether or not to use parallel computation. Default isfalse
. If true, the function will use theOhMyThreads
library to parallelize the computation.
Returns
acwresults
: AnACWResults
object. It has the fieldsfs
,acw_results
,acwtypes
,n_lags
,freqlims
,acf
,psd
,freqs
,lags
,x_dim
. You can access these fields asacwresults.field
. The fieldacw_results
contains the ACW results indicated by the input argumentacwtypes
in the same order you specify. Each element ofacw_results
is an array of the same size of your data minus the dimension of time, which will be dropped. See below for details.
The reason to not return the results directly but return the ACWResults
object is 1) give access to ACF and PSDs where the calculations are performed as well as n_lags
and freqlims
if the user is using defaults, 2) make plotting easy. You can simply type acwplot(acwresults)
to plot ACF and PSDs.
Your primary interest should be the field acwresults.acw_results
. This is a vector of arrays. Easiest way to explain this is via an example:
data = randn(2, 1000, 10) # assume 2 subjects, 1000 time points and 10 trials
fs = 1.0
acwresults = acw(data, fs; acwtypes=[:acw0, :tau], dims=2)
acw_results = acwresults.acw_results
acw_results
is a two element vector containing the results with the same order of acwtypes
as you specify. Since we wrote :acw0
as the first element and :tau
as the second element, we can extract the results as
acw_0 = acw_results[1]
tau = acw_results[2]
Let's check the dimensionality of these results. Remember that we specified 2 subjects, 1000 time points and 10 trials. The result collapses the dimension of time and gives the result as an 2x10 matrix.
size(acw_0) # should be (2, 10)
size(tau) # should be (2, 10)
Other fields:
fs
: Sampling rate. Floating point number.acwtypes
: The ACW types you specified.n_lags
: The number of lags used to fit exponential decay function to the autocorrelation function. See above in the input arguments for details. An integer.freqlims
: Frequency limits to fit a lorentzian to the power spectrum. See above in the input arguments for details. A tuple of floating point numbers.acf
: Autocorrelation function(s). Has the same size of your data with the time dimension replaced by lag dimension withn_lags
elements.psd
: Power spectrum/spectra. Has the same size of your data with the time dimension replaced by frequency dimension withfreqlims
as lowest and highest frequencies.freqs
: Frequencies corresponding to PSD.lags
: Lags corresponding to ACF.x_dim
: The dimension corresponding to lags and frequencies. Used internally in plotting.
Plotting
The function acwplot
can plot power spectra and autocorrelation functions. Currently it supports only two dimensions (for example, subjects x time or trials x time).
p = acwplot(acwresults; only_acf=false, only_psd=false, show=true)
Mandatory Arguments
acwresults
:ACWResults
type obtained by running the functionacw
.
Optional Arguments
only_acf
/only_psd
: Plot only the ACF or only the PSD. Boolean.show
: Whether to show the plot or only return the variable that contains the plot.
Returns
p
: The plot for further modification using the Plots library.