AFNI Message Board

Dear AFNI users-

We are very pleased to announce that the new AFNI Message Board framework is up! Please join us at:

https://discuss.afni.nimh.nih.gov

Existing user accounts have been migrated, so returning users can login by requesting a password reset. New users can create accounts, as well, through a standard account creation process. Please note that these setup emails might initially go to spam folders (esp. for NIH users!), so please check those locations in the beginning.

The current Message Board discussion threads have been migrated to the new framework. The current Message Board will remain visible, but read-only, for a little while.

Sincerely, AFNI HQ

History of AFNI updates  

|
May 27, 2022 08:52PM
Hi, Philipp-

This reply talks about AFNI's 3dPeriodogram, and scipy's scipy.signal.periodogram(); I refer to the latter here as sps.periodogram().

So, neither of the programs involve logs---that is not an issue here (you can take log after).

You can taper approx. similarly in each program---in my example, I will not taper, for simplicity. Therefore, 3dPeriodogram divides the power values by N, following the help note "(if ntaper = 0, then P = npts).", to get a power density. As noted just below, sps.periodogram() will also scale in some way---setting kwarg scaling='density' matches with AFNI.

Each program by default detrends the time series: 3dPeriodogram detrends to first order (so, constant+line); sps.periodogram() by default to 'constant'. So, this is a difference in the two programs. At present, 3dPeriodogram has no option about changing this. In sps.periodogram(), you can set the kwarg detrend='linear', to match 3dPeriodogram.

In terms of power vs magnitude: 3dPeriodogram estimates power spectrum (only); sps.periodogram() can do different 'scaling' types, as they call it. The kwarg scaling='density' (default), what is called power spectral density, matches with 3dPeriodogram.

So, the following two commands should produce (*almost*) identical output:
+ AFNI:
3dPeriodogram                                            \
    -taper 0                                                    \
    -prefix DSET_FREQ                     \
    DSET_TIME
+ sps.periodogram():
fvals, ARR_FREQ = sps.periodogram(ARR_TIME, scaling='density', window=None,
                               detrend='linear')
if ARR_TIME and DSET_TIME contain the same time series.
-> Now, the *almost*:
+ for a time series of, say, N=128 time points, 3dPeriodogram will output a frequency series of N=64 time points, and sps.periodogram() will output a time series of 65 points. The difference in length is due to sps.periodogram() outputting the baseline/constant value: this will be default by zero, but one could turn off detrending, so this wouldn't be trivial. I guess because 3dPeriodogram does always detrend, that value is always zero, and hence not output.
+ for a time series of odd N, the lower frequencies match, but the higher ones differ slightly. I am not sure why---I need to ponder this. The two methods are doing something different to deal with the odd number of time points, and I have to ponder.
+ The sps.periodogram() power spectrum is exactly 2x the value of the 3dPeriodogram value in each associated frequency point (for even N). I believe what is happening here is: for N input time points, the Fourier relations give N output frequency points. However, for real-valued time series, there is a symmetry of positive and negative frequencies, so that there are only N/2 independent ones. Often, people only report the positive half of the spectrum, for wave number k=0...N/2 values (or k=1...N/2 values, in 3dPeriodogram's case). Now, in terms of power-at-a-given-frequency, one would add the power values of the positive and negative frequencies---I suspect that sps.periodogram is doing this, to give total power at the frequency, while 3dPeriodogram is returning the power of that positive half-spectrum. You can choose what is a more appropriate value for your context of power-law spectrum.

I hope that makes sense. As with all things, there is a lot of subtlety when moving between different programs, even for something as fundamental as the FFT.

--pt
Subject Author Posted

3dPeriodogram and 1dFFT - How exactly is the power computed? Attachments

Philipp May 27, 2022 10:36AM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

ptaylor May 27, 2022 03:28PM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

Philipp May 27, 2022 03:58PM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

ptaylor May 27, 2022 08:52PM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

Philipp May 28, 2022 03:36AM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

ptaylor May 28, 2022 07:11AM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed? Attachments

Philipp May 28, 2022 11:03AM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

ptaylor May 31, 2022 05:41PM

Re: 3dPeriodogram and 1dFFT - How exactly is the power computed?

Philipp June 01, 2022 12:24AM