diff --git a/plotting.py b/plotting.py
index 37542df18adcc51e7bad5f7132dcbdc53dfe334c..9513a8006ffb33ac602c163abd29e1aa9ce1d97d 100644
--- a/plotting.py
+++ b/plotting.py
@@ -525,3 +525,88 @@ def plot_cuml_divisions_per_cell(t, buddings, nboots=30, col="b", label=None):
     )
     plt.plot(t, cuml, color=col, label=label)
     plt.fill_between(t, cuml - err_cuml, cuml + err_cuml, color=col, alpha=0.2)
+
+
+def bud_to_bud_plot(tpt, signal, dl, nbins=None, return_signal=False):
+    """
+    Plot the median and percentiles of a signal between consecutive buddings.
+
+    The first budding event is assigned to be at time point 0 and the second
+    at time point 1.
+
+    Parameters
+    ----------
+    tpt: float
+        The time point of interest in hours. Data from a bud-to-bud
+        time-series for each cell that contains this time point will
+        be analysed.
+    signal: str
+        The signal to plot.
+    dl: dataloader object
+        A dataloader object with the data to be plotted.
+    nbins: integer (optional)
+        The number of time bins to partition the interval between the
+        first and the second budding event.
+    return_signal: boolean (optional)
+        If True, return the signal for each cell interpolated to the time
+        bins.
+
+    Example
+    -------
+    >>> from wela.plotting import bud_to_bud_plot
+    >>> bud_to_bud_plot(8.4, "bud_growth_rate", dl)
+    """
+    t, signal_data = dl.get_time_series(signal)
+    t, buddings = dl.get_time_series("buddings")
+    if np.max(t) > 48:
+        # convert to hours
+        t = t / 60
+    tpt_i = np.argmin((t - tpt) ** 2)
+    # get data for bud-to-bud around tpt for each cell
+    local_signals, local_times = [], []
+    for i in range(signal_data.shape[0]):
+        future_buddings = np.nonzero(buddings[i, :][tpt_i:])[0]
+        if np.any(future_buddings):
+            end_tpt_i = tpt_i + future_buddings[0]
+            past_buddings = np.nonzero(buddings[i, :][:tpt_i])[0]
+            if np.any(past_buddings):
+                start_tpt_i = past_buddings[-1]
+                local_signals.append(
+                    signal_data[i, start_tpt_i : end_tpt_i + 1]
+                )
+                local_times.append(t[start_tpt_i : end_tpt_i + 1])
+    # find bins for normalised time, between 0 and 1
+    nbins = int(np.median([len(local_time) for local_time in local_times]))
+    ntbins = np.linspace(0, 1, nbins)
+    # interpolate each local signal to make a new signal
+    new_signal = np.nan * np.ones((len(local_signals), nbins))
+    for i in range(len(local_signals)):
+        s = local_signals[i]
+        # normalise time between 0 and 1
+        nt = local_times[i] - local_times[i][0]
+        nt /= nt[-1]
+        # interpolate into the bins
+        new_signal[i, :] = np.interp(
+            ntbins,
+            nt[~np.isnan(s)],
+            s[~np.isnan(s)],
+            left=np.nan,
+            right=np.nan,
+        )
+    # plot median and percentiles
+    plt.figure()
+    plt.plot(ntbins, np.nanmedian(new_signal, axis=0), "b.-")
+    for lower, upper in zip([45, 40, 35], [55, 60, 65]):
+        plt.fill_between(
+            ntbins,
+            np.nanpercentile(new_signal, lower, axis=0),
+            np.nanpercentile(new_signal, upper, axis=0),
+            alpha=0.06,
+            color="b",
+        )
+    plt.xlabel("position between budding events")
+    plt.ylabel(signal.replace("_", " "))
+    plt.title(f"t={tpt}")
+    plt.show(block=False)
+    if return_signal:
+        return new_signal