AFNI program: 3dLombScargle

Output of -help

++ Reading in options.

  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

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
 
____________________________________________________________________________
This page auto-generated on Thu Oct 31 09:41:59 PM EDT 2024