Make a periodogram or amplitude-spectrum of a time series that has a
non-constant sampling rate. The spectra output by this program are
'one-sided', so that they represent the half-amplitude or power
associated with a frequency, and they would require a factor of 2 to
account for both the the right- and left-traveling frequency solutions
of the Fourier transform (see below 'OUTPUT' and 'NOTE').
Of particular interest is the application of this functionality to
resting state time series that may have been censored. The theory behind
the mathematics and algorithms of this is due to separate groups, mainly
in the realm of astrophysical applications: Vaníček (1969, 1971),
Lomb (1976), Scargle (1982), and Press & Rybicki (1989). Shoutout to them.
This particular implementation is due to Press & Rybicki (1989), by
essentially translating their published Fortran implementation into C,
while using GSL for the FFT, instead of NR's realft(), and making
several adjustments based on that.
The Lomb-Scargle adaption was done with fairly minimal changes here by
PA Taylor (v1.4, June, 2016).
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ USAGE:
Input a 4D volumetric time series (BRIK/HEAD or NIFTI data set)
as well as an optional 1D file of 0s and 1s that defines which points
to censor out (i.e., each 0 represents a point/volume to censor out);
if no 1D file is input, the program will check for volumes that are
uniformly zero and consider those to be censored.
The output is a LS periodogram, describing spectral magnitudes
up to some 'maximum frequency'-- the default max here is what
the Nyquist frequency of the time series *would have been* without
any censoring. (Interestingly, this analysis can actually be
legitimately applied in cases to estimate frequency content >Nyquist.
Wow!)
The frequency spectrum will be in the range [df, f_N], where:
df = 1/T, and T is the total duration of the uncensored time series;
f_N = 1/dt, and dt is the sampling time (i.e., TR);
and the interval of frequencies is also df.
These ranges and step sizes should be *independent* of the censoring
which is a nice property of the Lomb-Scargle-iness.
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ OUTPUT:
1) PREFIX_time.1D :a 1D file of the sampled time points (in units of
seconds) of the analyzed (and possibly censored)
data set.
2) PREFIX_freq.1D :a 1D file of the frequency sample points (in units
of 1/seconds) of the output periodogram/spectrum
data set.
3) PREFIX_amp+orig :volumetric data set containing a LS-derived
or amplitude spectrum (by default, named 'amp') or a
PREFIX_pow+orig power spectrum (see '-out_pow_spec', named 'pow')
one per voxel.
Please note that the output amplitude and power
spectra are 'one-sided', to represent the
*half* amplitude or power of a given frequency
(see the following note).
+ A NOTE ABOUT Fourier+Parseval matters (please forgive the awkward
formatting):
In the formulation used here, for a time series x[n] of length N,
the periodogram value S[k] is related to the amplitude value |X[k]|:
(1) S[k] = (|X[k]|)**2,
for each k-th harmonic.
Parseval's theorem relates time fluctuations to spectral amplitudes,
stating that (for real time series with zero mean):
(2) sum_n{ x[n]**2 } = (1/N) * sum_k{ |X[k]|**2 },
= (1/N) * sum_k{ S[k] },
where n=0,1,..,N-1 and k=0,1,..,N-1 (NB: A[0]=0, for zero mean
series). The LHS is essentially the variance of the time series
(times N-1). The above is derived from Fourier transform maths, and
the Lomb-Scargle spectra are approximations to Fourier, so the above
can be expected to approximately hold, if all goes well.
Another Fourier-related result is that for real, discrete time series,
the spectral amplitudes/power values are symmetric and periodic in N.
Therefore, |X[k]| = |X[-k]| = |X[N-k-1]| (in zero-base array
counting);
the distinction between positive- and negative-indexed frequencies
can be thought of as signifying right- and left-traveling waves, which
both contribute to the total power of a specific frequency.
The upshot is that one could write the Parseval formula as:
(3) sum_n{ x[n]**2 } = (2/N) * sum_l{ |X[l]|**2 },
= (2/N) * sum_l{ S[l] },
where n=0,1,..,N-1 and l=0,1,..,(N/2)-1 (note the factor of 2 now
appearing on the RHS relations). These symmetries/considerations
are the reason why ~N/2 frequency values are output here (we assume
that only real-valued time series are input), without any loss of
information.
Additionally, with a view toward expressing the overall amplitude
or power of a given frequency, which many people might want to use to
estimate spectral 'functional connectivity' parameters such as ALFF,
fALFF, RSFA, etc. (using, for example, 3dAmptoRSFC), we therefore
note that the *total* amplitude or power of a given frequency would
be:
A[k] = 2*|X[k]|
P[k] = 2*S[k] = 2*|X[k]|**2 = 0.5*A[k]**2
instead of just that of the left/right traveling part. These types of
quantities (A and P) are also referred to as 'two-sided' spectra. The
resulting Parseval relation could then be written:
(4) sum_n{ x[n]**2 } = (1/(2N)) * sum_l{ A[l]**2 },
= (1/N) * sum_l{ P[l] },
where n=0,1,..,N-1 and l=0,1,..,(N/2)-1. Somehow, it just seems easier
to output the one-sided values, X and S, so that the Parsevalian
summation rules look more similar.
With all of that in mind, the 3dLombScargle results are output as
follows. For amplitudes, the following approx. Parsevellian relation
should hold between the 'holey' time series x[m] of M points and
the frequency series Y[l] of L~M/2 points (where {|Y[l]|} approaches
the Fourier amplitudes {|X[l]|} as the number of censored points
decreases and M->N):
(5) sum_m{ x[m]**2 } = (1/L) * sum_l{ Y[l]**2 },
where m=0,1,..,M-1 and l=0,1,..,L-1. For the power spectrum T[l]
of L~M/2 values, then:
(6) sum_m{ x[m]**2 } = (1/L) * sum_l{ T[l] }
for the same ranges of summations.
So, please consider that when using the outputs of here. 3dAmpToRSFC
is prepared for this when calculating spectral parameters (from
amplitudes).
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ COMMAND: 3dLombScargle -prefix PREFIX -inset FILE \
{-censor_1D C1D} {-censor_str CSTR} \
{-mask MASK} {-out_pow_spec} \
{-nyq_mult N2} {-nifti}
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ RUNNING:
-prefix PREFIX :output prefix name for data volume, time point 1D file
and frequency 1D file.
-inset FILE :time series of volumes, a 4D volumetric data set.
-censor_1D C1D :single row or column of 1s (keep) and 0s (censored)
describing which volumes of FILE are kept in the
sampling and which are censored out, respectively. The
length of the list of numbers must be of the
same length as the number of volumes in FILE.
If not entered, then the program will look for subbricks
of all-zeros and assume those are censored out.
-censor_str CSTR :AFNI-style selector string of volumes to *keep* in
the analysis. Such as:
'[0..4,7,10..$]'
Why we refer to it as a 'censor string' when it is
really the list of volumes to keep... well, it made
sense at the time. Future historians can duel with
ink about it.
-mask MASK :optional, mask of volume to analyze; additionally, any
voxel with uniformly zero values across time will
produce a zero-spectrum.
-out_pow_spec :switch to output the amplitude spectrum of the freqs
instead of the periodogram. In the formulation used
here, for a time series of length N, the power spectral
value S is related to the amplitude value X as:
S = (X)**2. (Without this opt, default output is
amplitude spectrum.)
-nyq_mult N2 :L-S periodograms can include frequencies above what
would typically be considered Nyquist (here defined
as:
f_N = 0.5*(number of samples)/(total time interval)
By default, the maximum frequency will be what
f_N *would* have been if no censoring of points had
occurred. (This makes it easier to compare L-S spectra
across a group with the same scan protocol, even if
there are slight differences in censoring, per subject.)
Acceptable values are >0. (For those reading the
algorithm papers, this sets the 'hifac' parameter.)
If you don't have a good reason for changing this,
dooon't change it!
-nifti :switch to output *.nii.gz volume file
(default format is BRIK/HEAD).
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+ EXAMPLE:
3dLombScargle -prefix LSout -inset TimeSeries.nii.gz \
-mask mask.nii.gz -censor_1D censor_list.txt
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
