Skip to content

Off-by-one indexing issue in convolution #10

Description

@michaelhallquist

The fmri.stimulus command seems to have an indexing problem when going from the microtime resolution (upsampled) back to the original temporal resolution. The problem arises when calculating the indices for the subset:

ind <- unique((trunc(min(scale * TR, 1) * scale):scans) %/% (scale^2 * TR)) * scale^2 * TR.

The first value tends to be 0, which is not a valid index in R, and the result is one greater than scans (e.g., if I ask for 100 scans, ind will be 101). The invalid index of 0 is silently dropped by R, such that the first index is actually the very end of the first microtime bin. As a result, the convolved regressor tends to be off by 1 (shifted too early).

I believe the correct syntax should be:

ind <- seq(from = scale/2 * TR, to = scans, by = scale^2 * TR)

Here is a small reproducible example comparing the convolution result in the package to the regressor produced by convolution in AFNI:

TR <- 2
onsets <- TR*(seq(2.25, 80, by=5) - 1) # in seconds
durations <- 1*TR
f2.pkg <- fmri::fmri.stimulus(scans=100, onsets=onsets, durations=durations, times = TRUE, TR=TR)

system(glue::glue("/Users/hallquist/abin/3dDeconvolve \\
-overwrite \\
-nodata 100 {TR} \\
-polort -1 \\
-num_stimts 1 \\
-stim_times 1 '1D: {paste(onsets, collapse=' ')}' 'SPMG1({durations[1L]})' \\
-stim_label 1 mystim \\
-x1D X.xmat.1D -x1D_stop"))

afni_res <- read.table("X.xmat.1D")$V1
plot(scale(afni_res), type="l")
lines(scale(f2), type="l", col="blue")

print(ccf(afni_res, f2.pkg))
cor(afni_res, f2.pkg) # .55 without the fix

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions