code: def phase_randomize(data, voxelwise=False, random_state=None):
      """Randomize phase of time series across subjects
  
      For each subject, apply Fourier transform to voxel time series
      and then randomly shift the phase of each frequency before inverting
      back into the time domain. This yields time series with the same power
      spectrum (and thus the same autocorrelation) as the original time series
      but will remove any meaningful temporal relationships among time series
      across subjects. By default (voxelwise=False), the same phase shift is
      applied across all voxels; however if voxelwise=True, different random
      phase shifts are applied to each voxel. The typical input is a time by
      voxels by subjects ndarray. The first dimension is assumed to be the
      time dimension and will be phase randomized. If a 2-dimensional ndarray
      is provided, the last dimension is assumed to be subjects, and different
      phase randomizations will be applied to each subject.
  
      The implementation is based on the work in [Lerner2011]_ and
      [Simony2016]_.
  
      Parameters
      ----------
      data : ndarray (n_TRs x n_voxels x n_subjects)
          Data to be phase randomized (per subject)
  
      voxelwise : bool, default: False
          Apply same (False) or different (True) randomizations across voxels
  
      random_state : RandomState or an int seed (0 by default)
          A random number generator instance to define the state of the
          random permutations generator.
  
      Returns
      ----------
      shifted_data : ndarray (n_TRs x n_voxels x n_subjects)
          Phase-randomized time series
      """
  
      # Check if input is 2-dimensional
      data_ndim = data.ndim
  
      # Get basic shape of data
      data, n_TRs, n_voxels, n_subjects = _check_timeseries_input(data)
  
      # Random seed to be deterministically re-randomized at each iteration
      if isinstance(random_state, np.random.RandomState):
          prng = random_state
      else:
          prng = np.random.RandomState(random_state)
  
      # Get randomized phase shifts
      if n_TRs % 2 == 0:
          # Why are we indexing from 1 not zero here? n_TRs / -1 long?
          pos_freq = np.arange(1, data.shape[0] // 2)
          neg_freq = np.arange(
              data.shape[0] - 1, data.shape[0] // 2, -1
          )
      else:
          pos_freq = np.arange(1, (data.shape[0] - 1) // 2 + 1)
          neg_freq = np.arange(
              data.shape[0] - 1, (data.shape[0] - 1) // 2, -1
          )
  
      if not voxelwise:
          phase_shifts = (
              prng.rand(len(pos_freq), 1, n_subjects) * 2 * np.math.pi
          )
      else:
          phase_shifts = (
              prng.rand(len(pos_freq), n_voxels, n_subjects)
              * 2
              * np.math.pi
          )
  
      # Fast Fourier transform along time dimension of data
      fft_data = fft(data, axis=0)
  
      # Shift pos and neg frequencies symmetrically, to keep signal real
      fft_data[pos_freq, :, :] *= np.exp(1j * phase_shifts)
      fft_data[neg_freq, :, :] *= np.exp(-1j * phase_shifts)
  
      # Inverse FFT to put data back in time domain
      shifted_data = np.real(ifft(fft_data, axis=0))
  
      # Go back to 2-dimensions if input was 2-dimensional
      if data_ndim == 2:
          shifted_data = shifted_data[:, 0, :]
  
      return shifted_data
  
imports: import numpy
  import re
  import warnings
  import os.path
  import psutil
  from .fmrisim import generate_stimfunction, _double_gamma_hrf, convolve_hrf
  from scipy.fftpack import fft, ifft
  import logging
label_desc: Use image node_id 12 for calling a global function with example usage: phase_randomize(data, voxelwise, random_state) and returns: shifted_data