AFNI program: 3dPeriodogram

Output of -help

Usage: 3dPeriodogram [options] dataset
Computes the periodogram of each voxel time series.
(Squared FFT = a crude estimate of the power spectrum)

--------
Options:
--------
 -prefix p = use string 'p' for the prefix of the
               output dataset [DEFAULT = 'pgram']
 -taper    = fraction of data to taper [DEFAULT = 0.1]
 -nfft L   = set FFT length to 'L' points
               (longer than the data ==> zero padding)
               (shorter than the data ==> data pruning)
------
Notes:
------
* Output is in float format; number of sub-bricks will be
   half the FFT length; sub-brick #0 = FFT bin #1, etc.
* Grid spacing in the frequency (sub-brick) dimension will
   be 1/(nfft*TR) where nfft=FFT length, TR=dataset timestep.
* There is no '-mask' option.  The hyper-clever user could
   use something like
     '3dcalc( -a dset+orig -b mask+orig -expr a*b )'
   to apply a binary mask on the command line.
* Data is not scaled exactly as in the AFNI Power plugin.
* Each time series is linearly detrended prior to FFT-ization.
* FFT length defaults to be the next legal length >= input dataset.
* The program can only do FFT lengths that are positive even integers.
  ++ '-nfft' with an illegal value will cause the program to fail.
* If you want to do smaller FFTs, then average the periodograms
   (to reduce random fluctuations), you can use 3dPeriodogram in
   a script with "[...]" sub-brick selectors, then average
   the results with 3dMean.
* Or you could use the full-length FFT, then smooth that FFT
   in the frequency direction (e.g., with 3dTsmooth).
* This is a really quick hack for DH and PB and SfN.

* Author = RWCox -- who doesn't want any bribe at all for this!
                 -- http://ethics.od.nih.gov/topics/gifts.htm

---------------------------------------------------
More Details About What 3dPeriodogram Actually Does
---------------------------------------------------
* Tapering is done with the Hamming window (if taper > 0):
    Define npts   = number of time points analyzed (<= nfft)
                    (i.e., the length of the input dataset)
           ntaper = taper * npts / 2        (0 < taper <= 1)
                  = number of points to taper on each end
           ktop   = npts - ntaper
           phi    = PI / ntaper
    Then the k-th point (k=0..nfft-1) is tapered by
      w(k) = 0.54 - 0.46 * cos(k*phi)           0    <= k < ntaper
      w(k) = 0.54 + 0.46 * cos((k-ktop+1)*phi)  ktop <= k < npts
      w(k) = 1.0                                otherwise
    Also define P = sum{ w(k)*w(k) } from k=0..npts-1
    (if ntaper = 0, then P = npts).

* The result is the squared magnitude of the FFT of w(k)*data(k),
  divided by P.  This division makes the result be the 'power',
  which is to say the data's sum-of-squares ('energy') per unit
  time (in units of 1/TR, not 1/sec) ascribed to each FFT bin.

* Normalizing by P also means that the values output for different
  amounts of tapering or different lengths of data are comparable.

* To be as clear as I can: this program does NOT do any averaging
  across multiple windows of the data (such as Welch's method does)
  to estimate the power spectrum.  This program:
  ++ tapers the data,
  ++ zero-pads it to the FFT length,
  ++ FFTs it (in time),
  ++ squares it and divides by the P factor.

* The number of output sub-bricks is nfft/2:
     sub-brick #0 = FFT bin #1 = frequency 1/(nfft*dt)
               #1 = FFT bin #2 = frequency 2/(nfft*dt)
  et cetera, et cetera, et cetera.

* If you desire to implement Welch's method for spectrum estimation
  using 3dPeriodogram, you will have to run the program multiple
  times, using different subsets of the input data, then average
  the results with 3dMean.
  ++ https://en.wikipedia.org/wiki/Welch's_method

++ Compile date = Oct 31 2024 {AFNI_24.3.06:linux_ubuntu_24_64}
This page auto-generated on Thu Oct 31 09:42:23 PM EDT 2024