diff --git a/brainiak/isfc.py b/brainiak/isfc.py
new file mode 100644
index 000000000..dca2a8670
--- /dev/null
+++ b/brainiak/isfc.py
@@ -0,0 +1,167 @@
+#  Copyright 2017 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+"""Intersubject analyses (ISC/ISFC)
+
+Functions for computing intersubject correlation (ISC) and intersubject
+functional correlation (ISFC), and utility functions for loading subject data
+files and masks
+
+Paper references:
+ISC: Hasson, U., Nir, Y., Levy, I., Fuhrmann, G. & Malach, R. Intersubject
+synchronization of cortical activity during natural vision. Science 303,
+1634–1640 (2004).
+
+ISFC: Simony E, Honey CJ, Chen J, Lositsky O, Yeshurun Y, Wiesel A, Hasson U
+(2016) Dynamic reconfiguration of the default mode network during narrative
+comprehension. Nat Commun 7.
+"""
+
+# Authors: Christopher Baldassano and Mor Regev
+# Princeton University, 2017
+
+from brainiak.fcma.util import compute_correlation
+import nibabel as nib
+import numpy as np
+from scipy import stats
+
+
+def isc(D, collapse_subj=True):
+    """Intersubject correlation
+
+    For each voxel, computes the correlation of each subject's timecourse with
+    the mean of all other subjects' timecourses. By default the result is
+    averaged across subjects, unless collapse_subj is set to False.
+
+    Parameters
+    ----------
+    D : voxel by time by subject ndarray
+        fMRI data for which to compute ISC
+
+    collapse_subj : bool, default:True
+        Whether to average across subjects before returning result
+
+    Returns
+    -------
+    ISC : voxel ndarray (or voxel by subject ndarray, if collapse_subj=False)
+        pearson correlation for each voxel, across subjects
+    """
+
+    n_vox = D.shape[0]
+    n_subj = D.shape[2]
+    ISC = np.zeros((n_vox, n_subj))
+
+    # Loop across choice of leave-one-out subject
+    for loo_subj in range(n_subj):
+        group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
+        subj = D[:, :, loo_subj]
+        for v in range(n_vox):
+            ISC[v, loo_subj] = stats.pearsonr(group[v, :], subj[v, :])[0]
+
+    if collapse_subj:
+        ISC = np.mean(ISC, axis=1)
+    return ISC
+
+
+def isfc(D, collapse_subj=True):
+    """Intersubject functional correlation
+
+    Computes the correlation between the timecoure of each voxel in each
+    subject with the average of all other subjects' timecourses in *all*
+    voxels. By default the result is averaged across subjects, unless
+    collapse_subj is set to False.
+
+    Uses the high performance compute_correlation routine from fcma.util
+
+    Parameters
+    ----------
+    D : voxel by time by subject ndarray
+        fMRI data for which to compute ISFC
+
+    collapse_subj : bool, default:True
+        Whether to average across subjects before returning result
+
+    Returns
+    -------
+    ISFC : voxel by voxel ndarray
+        (or voxel by voxel by subject ndarray, if collapse_subj=False)
+        pearson correlation between all pairs of voxels, across subjects
+    """
+
+    n_vox = D.shape[0]
+    n_subj = D.shape[2]
+    ISFC = np.zeros((n_vox, n_vox, n_subj))
+
+    # Loop across choice of leave-one-out subject
+    for loo_subj in range(D.shape[2]):
+        group = np.mean(D[:, :, np.arange(n_subj) != loo_subj], axis=2)
+        subj = D[:, :, loo_subj]
+        ISFC[:, :, loo_subj] = compute_correlation(group, subj)
+
+        # Symmetrize matrix
+        ISFC[:, :, loo_subj] = (ISFC[:, :, loo_subj] +
+                                ISFC[:, :, loo_subj].T) / 2
+
+    if collapse_subj:
+        ISFC = np.mean(ISFC, axis=2)
+    return ISFC
+
+
+def load_subjects_nii(data_files, mask_file, mask_func=None):
+    """Loading masked nifti data into matrix
+
+    Given a list of subject data files and a mask, loads voxel timecourses
+    inside mask into a matrix, which is returned along with the voxel
+    coordinates of the mask. If mask_func is given, it specifies which mask
+    values should be included (e.g. for thresholding a continuous-valued mask)
+
+    Parameters
+    ----------
+    data_files : list of filenames of subject nii files
+        each nii file should be 4D (space+time), with all dimensions identical
+        for all subjects
+
+    mask_file : filename of mask nii
+        mask should be 3D (space), with dimensions identical to the first
+        three dimensions of the data files
+
+    mask_func : Callable[[ndarray], bool] : default x>0
+
+    Returns
+    -------
+    D : voxel by time by subject ndarray
+        all data within mask, from all subjects
+    coords : tuple of 3 ndarrays
+        x,y,z (as provided by nibabel) coordinates of mask voxel locations
+    """
+
+    mask_nii = nib.load(mask_file)
+    mask = mask_nii.get_data()
+    if mask_func is not None:
+        mask = mask_func(mask)
+    else:
+        mask = mask > 0
+    mask_shape = mask_nii.shape
+    coords = np.where(mask)
+
+    data_shape = nib.load(data_files[0]).shape
+    D = np.zeros((np.sum(mask), data_shape[3], len(data_files)))
+    for s in range(len(data_files)):
+        nii = nib.load(data_files[s])
+        if nii.shape != data_shape:
+            raise ValueError("Data has different shapes across subjects")
+        if nii.shape[:3] != mask_shape:
+            raise ValueError("Data and mask have different shapes")
+        D[:, :, s] = nii.get_data()[mask]
+
+    return (D, coords)
diff --git a/examples/isfc/download_data.sh b/examples/isfc/download_data.sh
new file mode 100755
index 000000000..19d0f05be
--- /dev/null
+++ b/examples/isfc/download_data.sh
@@ -0,0 +1,4 @@
+#!/bin/sh
+curl --location -o isfc.zip https://api.onedrive.com/v1.0/shares/s!Aobi2ryypFQCgqQOm2Zhs-TSLu9leQ/root/content
+unzip -qo isfc.zip
+rm -f isfc.zip
diff --git a/examples/isfc/isfc.py b/examples/isfc/isfc.py
new file mode 100644
index 000000000..9d0bfd2f1
--- /dev/null
+++ b/examples/isfc/isfc.py
@@ -0,0 +1,65 @@
+#  Copyright 2017 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+"""Example of intersubject analyses (ISC/ISFC)
+
+Computes ISC for all voxels within a brain mask, and
+computes ISFC for voxels with high ISC
+"""
+
+# Authors: Christopher Baldassano and Mor Regev
+# Princeton University, 2017
+
+import brainiak.isfc
+import nibabel as nib
+import numpy as np
+from matplotlib import pyplot as plt
+from scipy.cluster.hierarchy import fcluster, linkage
+import sys
+import os
+
+curr_dir = os.path.dirname(__file__)
+
+brain_fname = os.path.join(curr_dir,'avg152T1_gray_3mm.nii.gz')
+fnames = [os.path.join(curr_dir,
+          'sub-0' + format(subj, '02') + '-task-intact1.nii.gz') for
+          subj in np.arange(1, 5)]
+
+print('Loading data from ', len(fnames), ' subjects...')
+D, coords = brainiak.isfc.load_subjects_nii(fnames, brain_fname,
+                                            lambda x: x > 50)
+
+print('Calculating ISC on ', D.shape[0], ' voxels')
+ISC = brainiak.isfc.isc(D)
+ISC[np.isnan(ISC)] = 0
+
+print('Writing ISC map to file...')
+brain_nii = nib.load(brain_fname)
+ISC_vol = np.zeros(brain_nii.shape)
+ISC_vol[coords] = ISC
+ISC_nifti = nib.Nifti1Image(ISC_vol, brain_nii.affine, brain_nii.header)
+nib.save(ISC_nifti, 'ISC.nii.gz')
+
+ISC_mask = ISC > 0.2
+print('Calculating ISFC on ', np.sum(ISC_mask), ' voxels...')
+D_masked = D[ISC_mask, :, :]
+ISFC = brainiak.isfc.isfc(D_masked)
+
+print('Clustering ISFC...')
+Z = linkage(ISFC, 'ward')
+z = fcluster(Z, 2, criterion='maxclust')
+clust_inds = np.argsort(z)
+
+# Show the ISFC matrix, sorted to show the two main clusters
+plt.imshow(ISFC[np.ix_(clust_inds,clust_inds)])
+plt.show()
diff --git a/examples/isfc/requirements.txt b/examples/isfc/requirements.txt
new file mode 100644
index 000000000..f0b259309
--- /dev/null
+++ b/examples/isfc/requirements.txt
@@ -0,0 +1,2 @@
+matplotlib
+
diff --git a/tests/isfc/mask.nii.gz b/tests/isfc/mask.nii.gz
new file mode 100644
index 000000000..83b31ffbd
Binary files /dev/null and b/tests/isfc/mask.nii.gz differ
diff --git a/tests/isfc/subj1.nii.gz b/tests/isfc/subj1.nii.gz
new file mode 100644
index 000000000..9cf26af98
Binary files /dev/null and b/tests/isfc/subj1.nii.gz differ
diff --git a/tests/isfc/subj2.nii.gz b/tests/isfc/subj2.nii.gz
new file mode 100644
index 000000000..ec6bb9a44
Binary files /dev/null and b/tests/isfc/subj2.nii.gz differ
diff --git a/tests/isfc/test_isfc.py b/tests/isfc/test_isfc.py
new file mode 100644
index 000000000..c16d66b8d
--- /dev/null
+++ b/tests/isfc/test_isfc.py
@@ -0,0 +1,41 @@
+import brainiak.isfc
+import numpy as np
+import sys, os
+
+def test_ISC():
+    D = np.zeros((2,5,3))
+    D[:,:,0] = \
+        [[-0.36225433, -0.43482456,  0.26723158,  0.16461712, -0.37991465], \
+         [-0.62305959, -0.46660116, -0.50037994,  1.81083754,  0.23499509]]
+    D[:,:,1] = \
+        [[-0.31484153, -0.47486988,  0.18966625, -0.12568572, -0.40535156], \
+         [-0.78433298, -0.35875085, -0.6121344 ,  1.68267639,  0.28603493]]
+    D[:,:,2] = \
+        [[-0.26593192, -0.56914734,  0.28397317,  0.30276589, -0.42637472], \
+         [-0.67598379, -0.51549055, -0.64196342,  1.60686666,  0.04127293]]
+
+    ISC = brainiak.isfc.isc(D)
+    assert np.isclose(ISC, [0.9540602, 0.99585304]).all(), \
+        "Calculated ISC does not match ground truth"
+
+def test_loading_and_ISFC():
+    curr_dir = os.path.dirname(__file__)
+
+    mask_fname = os.path.join(curr_dir,'mask.nii.gz')
+    fnames = [os.path.join(curr_dir,'subj1.nii.gz'),
+              os.path.join(curr_dir,'subj2.nii.gz')]
+
+    D, coords = brainiak.isfc.load_subjects_nii(fnames, mask_fname)
+
+    assert D.shape == (4,5,2), "Loaded data has incorrect shape"
+
+    ISFC = brainiak.isfc.isfc(D)
+    
+    ground_truth = \
+        [[ 1,  1,  0, -1],
+         [ 1,  1,  0, -1],
+         [ 0,  0,  1,  0],
+         [-1, -1,  0,  1]]
+
+    assert np.isclose(ISFC, ground_truth).all(), \
+        "Calculated ISFC does not match ground truth"