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