diff --git a/neurodsp/rhythm/__init__.py b/neurodsp/rhythm/__init__.py
index 5d8872d4..342c4ada 100644
--- a/neurodsp/rhythm/__init__.py
+++ b/neurodsp/rhythm/__init__.py
@@ -2,3 +2,4 @@
 
 from .lc import compute_lagged_coherence
 from .swm import sliding_window_matching
+from .phase import pairwise_phase_consistency
diff --git a/neurodsp/rhythm/phase.py b/neurodsp/rhythm/phase.py
new file mode 100644
index 00000000..5d05bce5
--- /dev/null
+++ b/neurodsp/rhythm/phase.py
@@ -0,0 +1,123 @@
+"""Phase consistency measures."""
+
+import warnings
+from importlib import import_module
+from itertools import combinations, combinations_with_replacement
+
+import numpy as np
+
+###################################################################################################
+###################################################################################################
+
+def pairwise_phase_consistency(pha0, pha1=None, return_pairs=True, progress=None):
+    """Compute pairwise phase consistency.
+
+    Parameters
+    ----------
+    pha0 : 1d array
+        First phases from a wavelet analysis, in radians, from -pi to pi.
+    pha1 : 1d array, optional, default: None
+        Second phases from a wavelet analysis, in radians, from -pi to pi.
+    return_pairs : True
+        Returns distance pairs as a 1d array if True.
+    progress : {None, 'tqdm', 'tqdm.notebook}
+        Displays tqdm progress bar.
+
+    Returns
+    -------
+    distance_avg : float
+        Average pairwise circular distance index.
+    distances : 2d array, optional
+        Pairwise circular distance indices. Only returned if ``return_pairs` is True.
+
+    Notes
+    -----
+
+    - distance == -1: inverse phases
+    - distance ==  0: pi / 2 phase difference
+    - distance ==  1: equal phases
+
+
+    Reference
+    ---------
+    Vinck, M., van Wingerden, M., Womelsdorf, T., Fries, P., & Pennartz, C. M. A. (2010).
+    The pairwise phase consistency: A bias-free measure of rhythmic neuronal synchronization.
+    NeuroImage, 51(1), 112–122. https://doi.org/10.1016/j.neuroimage.2010.01.073
+    """
+
+    if pha0.ndim != 1:
+        raise ValueError("Phase array must be 1-dimensional.")
+
+    if pha1 is not None and pha0.shape != pha1.shape:
+        raise ValueError("Phase arrays must be the same length.")
+
+    # Pairwise indices generator
+    if pha1 is None:
+
+        # Number of pairwise combinations
+        n_combs = int((len(pha0) * (len(pha0) - 1)) / 2)
+
+        # Exclude self-combinations (i.e. ignore (0, 0), (1, 1)...)
+        iterable = enumerate(combinations(np.arange(len(pha0)), 2))
+
+    else:
+
+        # Include all combinations
+        n_combs = int(len(pha0) ** 2)
+
+        iterable = enumerate((row, col) for row in range(len(pha0)) for col in range(len(pha1)))
+
+    # Initialize variables
+    if return_pairs:
+        cumulative = None
+        distances = np.ones((len(pha0), len(pha0)))
+    else:
+        cumulative = 0
+        distances = None
+
+    # Optional progress bar
+    if progress is not None:
+        try:
+            tqdm = import_module(progress)
+            iterable = tqdm.tqdm(iterable, total=n_combs, dynamic_ncols=True,
+                                 desc='Computing Pairwise Distances')
+        except ImportError:
+            pass
+
+    # Compute distance indices
+    for idx, pair in iterable:
+
+        phi0 = pha0[pair[0]]
+
+        if pha1 is None:
+            phi1 = pha0[pair[1]]
+        else:
+            phi1 = pha1[pair[1]]
+
+        # Convert range from (-pi, pi) to (0, 2pi)
+        phi0 = phi0 + (2*np.pi) if phi0 < 0 else phi0
+        phi1 = phi1 + (2*np.pi) if phi1 < 0 else phi1
+
+        # Absolute angular distance
+        abs_dist = np.abs(phi0 - phi1)
+
+        # Take smaller angle (range 0 to pi)
+        abs_dist = (2*np.pi) - abs_dist if abs_dist > np.pi else abs_dist
+
+        # Pairwise circular distance index (PCDI)
+        distance = (np.pi - 2 * abs_dist) / np.pi
+
+        if isinstance(distances, np.ndarray) and pha1 is None:
+            distances[pair[0], pair[1]] = distance
+            distances[pair[1], pair[0]] = distance
+        elif isinstance(distances, np.ndarray) and pha1 is not None:
+            distances[pair[0], pair[1]] = distance
+        else:
+            cumulative += distance
+
+    distance_avg = cumulative / n_combs if distances is None else np.mean(distances)
+
+    if return_pairs:
+        return distance_avg, distances
+    else:
+        return distance_avg
diff --git a/neurodsp/tests/rhythm/test_phase.py b/neurodsp/tests/rhythm/test_phase.py
new file mode 100644
index 00000000..5fe649c4
--- /dev/null
+++ b/neurodsp/tests/rhythm/test_phase.py
@@ -0,0 +1,66 @@
+"""Tests for neurodsp.rhythm.phase."""
+
+from pytest import mark
+
+import numpy as np
+
+from neurodsp.tests.settings import FS, FREQ_SINE
+from neurodsp.timefrequency import phase_by_time
+from neurodsp.rhythm.phase import *
+
+###################################################################################################
+###################################################################################################
+
+@mark.parametrize('return_pairs', [True, False])
+@mark.parametrize('phase_shift', [0, .25, .5])
+def test_pairwise_phase_consistency(tsig_sine, return_pairs, phase_shift):
+
+    peaks = np.where(tsig_sine == 1.0)[0]
+
+    pha0 = phase_by_time(tsig_sine, FS)
+
+    # Phase shift
+    sig_shift = np.roll(tsig_sine, int((FS / FREQ_SINE) * phase_shift))
+    pha1 = phase_by_time(sig_shift, FS)
+
+    # Compute consistency
+    dist_avg = pairwise_phase_consistency(pha0[peaks], pha1[peaks], return_pairs, 'tqdm')
+
+    # Unpack results if needed
+    if return_pairs:
+
+        dist_avg, dists = dist_avg[0], dist_avg[1]
+
+        assert isinstance(dists, np.ndarray)
+        assert len(dists[0]) * len(dists[1]) == len(peaks) ** 2
+        assert np.mean(dists) == dist_avg
+
+    # Expected consistency
+    if phase_shift == 0:
+        expected = 1
+    elif phase_shift == .25:
+        expected = 0
+    elif phase_shift == .5:
+        expected = -1
+
+    assert isinstance(dist_avg, float)
+    assert round(dist_avg) == expected
+
+    # Test self-consistency
+    dist_avg, dists = pairwise_phase_consistency(pha0[peaks], return_pairs=True)
+
+    assert dist_avg == 1
+    assert len(dists[0]) == len(dists[1]) == len(peaks)
+
+    # Cases where arrays are invalid sizes
+    try:
+        pairwise_phase_consistency(pha0[peaks], pha1[peaks][:2], return_pairs, 'tqdm')
+        assert False
+    except ValueError:
+        pass
+
+    try:
+        pairwise_phase_consistency(np.zeros((2, 2)), return_pairs, 'tqdm')
+        assert False
+    except ValueError:
+        pass