<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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