Hi Paul,

thank you a lot for the elaborate explanation! See my replies below.

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

Correct. This is why I showed the power spectra (and not their log-log version) in my first post. The “problem” (the difference in power) has nothing to do with the log, but is already present in the standard power spectra.

Quote
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

I changed the detrending option to "detrend='linear'", because I read that 3dPeriodogram uses linear detrending in the AFNI page of 3dPeriodogram. (Just as a side note here: the detrending option as well as the window option do not really change the power and its distribution in the standard power spectra (or in the log-log power-law spectra and their associated PLE=slope value), but they rather "smooth" the graph a bit differently).

Quote
+ 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.

I experienced that 3dPeriodogram does not allow an nfft number that is equal to the length (TRs) of the time-series, while normally that should be fine. As far as I remember (from my testing), 3dPeriodogam requires an nfft number that has the length of: time-series (N)+1 or N+2. I assume that this is related to the baseline/constant value you are talking about? I now always use a power of 2 (e.g. 1024 if the time-series has 560 time points).

Quote
+ 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 think this is the solution to my problem! As I wrote before, the PLE computation (a fitting of a linear regression line to the log-log power spectra) on power spectra created by sps.periodogram() and sps.welch() always yield approximately twice the PLE value than on power spectra created by 3dPeriodogram and 1dFFT.

If I apply the square root on the power output of sps.periodogram() and sps.welch() (e.g. np.sqrt(Power) in Python), the resulting PLE value turns out to be approximately the same as the PLE value for 3dPeriodogram and 1dFFT power spectra.
Conversely, if I square the power of the 3dPeriodogram results in Python (say "Power = Power**2 in pythonic), the resulting PLE value is again approximately the same as when using sps.periodogram().

There is a slight difference in the PLE value. But this difference, I assume, simply stems from the fact that in option 1 (3dPeriodogram and 1dFFT) power spectra were computed for every voxel, then averaged across the voxels with 3dmaskave, and finally the power-law is computed in Python.
While in option 2, all voxels’ time-series were extracted via 3dmaskave, and “one” power spectrum was computed in Python using sps.periodogram() (or sps.welch()). This, of course, can yield slightly different power spectra and consequently PLE values that differ a bit.

Quote
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).

As far as I remember, I have never seen people reporting in EEG or fMRI papers what kind they used to compute the PLE. This is surprising, because, as we saw, the resulting PLE can increase/decrease up to a 100%.

Philipp



Edited 1 time(s). Last edit at 05/28/2022 04:10AM by Philipp.