Another update.

With massive amount of help from a friend of mine, I managed to finally calculate the PLE (for each voxel). I will show you my script here for two reasons:

1. Maybe you have comments on the script.
2. Maybe the script, even though it is a very basic and simple script at this point (and there is tons to improve), could help other beginners who search for how to calculate the PLE in the AFNI board.

Here it is:

FFT & PLE

FFT = Fast Fourier Transform (transformation of a time domain into its frequency domain)
PLE = Power Law Exponent (the PLE quantifies the shape/slope of the frequency domain in one value/number)

The frequency domain has different slopes everywhere, i.e., it is not linear. Consequently, three steps are required to calculate the PLE.

1. The time-series/domain needs to be transformed into its frequency domain via the FFT.
2. The non-linear frequency domain needs adjustments to make it linear. This is done via the logarithm of both frequency and power (log(P) and log(F)).
3. The applied logarithm results into a line which is still not exactly linear. To create a linear line/slope, a linear regression line between log(F) and log(P) is applied.
The slope of this perfect linear line is the power law exponent (PLE).



# Fast Fourier Transform (FFT) - Transformation of the time-series into its frequency domain
directory=/users/philipp/desktop/fmri/dataset/subjects
directory_PLE=/users/philipp/desktop/fmri/dataset/Info
for subject in Subject1
do
mkdir $directory/$subject/FFT_PLE_RestingState
for fMRIruns (errts.Subject1_Rest.anaticor+tlrc)
do
cd $directory/$subject/Preprocessing_RestingState
echo "Processing $subject ..."

3dPeriodogram \
-nfft 192 \
-prefix $directory/$subject/FFT_PLE_RestingState/FFT.$fMRIruns $fMRIruns

1deval \
-a $directory_PLE/FFT_ideal.1D'{4..95}' \
-expr 'log(a)' \
> $directory/$subject/FFT_PLE_RestingState/FFT.1D

done
done


This is the first part of the script. What does it do? First, 3dPeriodogram adjusts the preprocessed and regressed functional dataset (e.g. "errts.Subject1_Rest.anaticor+tlrc" in this example, a resting-state run) to match its index number to the index number of our FFT_ideal.1D input file. It is the latter file that we have to provide for AFNI. Example: Lets say our that our provided "FFT_ideal.1D" file contains 96 values. Then, for the -nfft command, we will multiply this number by 2 (96*2=192). We will consequently use the value of 192: -nfft 192. The result is a transformation of the time-series into its frequency domain.

Second, 1deval calculates the logarithm for the FFT_ideal.1D input file. The output is a .1D file which contains all the numbers/values of the FFT.


# Logarithm of Amplitude/Power log(P) y-axis
directory=/users/philipp/desktop/fmri/dataset/subjects
for subject in Subject1
do
for fMRIruns (FFT.errts.Subject1_Rest.anaticor+tlrc)
do
cd $directory/$subject/FFT_PLE_RestingState
echo "Processing $subject ..."

3dTcat \
-prefix BP.$fMRIruns \
$fMRIruns'[4..95]'

3dTsmooth \
-hamming 7 \
-prefix Smooth.BP.$fMRIruns \
BP.$fMRIruns

3dcalc \
-prefix Log_Y.Smooth.BP.$fMRIruns \
-a Smooth.BP.$fMRIruns \
-expr '-log(a)'

done
done

In this third part above we basically do the same as for the second part, but for the y-axis. Note: 3dTsmooth smoothes the power. 3dcalc takes the logarithm of the power. We are now ready to calculate the PLE in the final step below. By opening the output file of this final script, one can inspect the PLE of every voxel via the Graph window in AFNI’s Gui.

# Linear Regression Line between log(F) and log(P)
directory=/users/philipp/desktop/fmri/dataset/subjects
for subject in Subject1
do
for fMRIruns (Log_Y.Smooth.BP.FFT.errts.Subject1_Rest.anaticor+tlrc)
do
cd $directory/$subject/FFT_PLE_RestingState
echo "Processing $subject ..."
3dfim+ \
-input $fMRIruns \
-ideal_file FFT.1D \
-out 'Fit Coef' \
-bucket PLE.$fMRIruns
done
done

Philipp



Edited 2 time(s). Last edit at 09/28/2021 09:45AM by Philipp.