diff --git a/.gitignore b/.gitignore
index 3dc927c..3ce2502 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,2 @@
 data/*
-./_pycache_/
\ No newline at end of file
+./__pycache__/
\ No newline at end of file
diff --git a/python/fig_pathway_stages.py b/python/fig_pathway_stages.py
index 80e55d4..73b0c68 100644
--- a/python/fig_pathway_stages.py
+++ b/python/fig_pathway_stages.py
@@ -6,10 +6,10 @@ from thunderhopper.modeltools import load_data
 from IPython import embed
 
 # GENERAL SETTINGS:
-data_paths = glob.glob('../data/*.npz')
+data_paths = glob.glob('../data/processed/*.npz')
+stages = ['filt', 'env', 'log', 'inv',
+          'conv', 'bi', 'feat']
 channel = 0
-stages_pre = ['raw', 'filt', 'env', 'log', 'inv']
-stages_feat = ['conv', 'bi', 'feat', 'norm']
 
 # PLOT SETTINGS:
 fig_kwargs = dict(
@@ -19,42 +19,115 @@ fig_kwargs = dict(
     sharey = 'row'
 )
 zoom_rel = np.array([0.4, 0.6])
-colors_pre = ['k' for _ in range(len(stages_pre))]
-colors_feat = ['k' for _ in range(len(stages_feat))]
+colors = dict(
+    filt='k',
+    env='k',
+    log='k',
+    inv='k',
+    conv='k',
+    bi='k',
+    feat='k'
+)
+linewidths = dict(
+    filt=0.25,
+    env=0.5,
+    log=0.5,
+    inv=0.5,
+    conv=0.5,
+    bi=1,
+    feat=1
+)
 
+# EXECUTION:
 for data_path in data_paths:
+    if 'Gomphocerippus' in data_path:
+        continue
+    print(f'Processing {data_path}')
+
     # Load overall data:
-    data, config = load_data(data_path)
-    t_full = np.arange(data['raw'].shape[0]) / config['rate']
+    data, config = load_data(data_path, stages)
+    t_full = np.arange(data['filt'].shape[0]) / config['rate']
 
     # Establish zoom frame:
     zoom_abs = zoom_rel * t_full[-1]
     zoom_mask = (t_full >= zoom_abs[0]) & (t_full <= zoom_abs[1])
     t_zoom = t_full[zoom_mask]
 
+
     # PART I: PREPROCESSING STAGE
-    fig, axes = plt.subplots(len(stages_pre), 2, **fig_kwargs)
+    fig, axes = plt.subplots(4, 2, **fig_kwargs)
+    fig.supylabel('amplitude', fontsize=plt.rcParams['axes.labelsize'])
+    fig.supxlabel('time [s]', fontsize=plt.rcParams['axes.labelsize'])
 
-    for i, stage in enumerate(stages_pre):
-        signal = data[stage][:, channel]
-        ax_full, ax_zoom = axes[i, :]
-        ax_full.plot(t_full, signal, c=colors_pre[i])
-        ax_full.set_ylim(signal.min(), signal.max())
+    # Bandpass-filtered signal:
+    signal = data['filt'][:, channel]
+    ax_full, ax_zoom = axes[0, :]
+    c, lw = colors['filt'], linewidths['filt']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask], c=c, lw=lw)
+    ax_full.set_ylim(signal.min(), signal.max())
 
-        ax_zoom.plot(t_zoom, signal[zoom_mask], c=colors_pre[i])
+    # Signal envelope:
+    signal = data['env'][:, channel]
+    ax_full, ax_zoom = axes[1, :]
+    c, lw = colors['env'], linewidths['env']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask], c=c, lw=lw)
+    ax_full.set_ylim(0, signal.max())
+
+    # Logarithmic envelope:
+    signal = data['log'][:, channel]
+    ax_full, ax_zoom = axes[2, :]
+    c, lw = colors['log'], linewidths['log']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask], c=c, lw=lw)
+    ax_full.set_ylim(signal.min(), 0)
+
+    # Adapted envelope:
+    signal = data['inv'][:, channel]
+    ax_full, ax_zoom = axes[3, :]
+    c, lw = colors['inv'], linewidths['inv']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask], c=c, lw=lw)
+    ax_full.set_ylim(signal.min(), signal.max())
+
+    # Posthoc adjustments:
     ax_full.set_xlim(t_full[0], t_full[-1])
+    ax_zoom.set_xlim(t_zoom[0], t_zoom[-1])
+
 
     # PART II: FEATURE EXTRACTION STAGE:
-    fig, axes = plt.subplots(len(stages_feat), 2, **fig_kwargs)
+    fig, axes = plt.subplots(3, 2, **fig_kwargs)
+    fig.supylabel('amplitude', fontsize=plt.rcParams['axes.labelsize'])
+    fig.supxlabel('time [s]', fontsize=plt.rcParams['axes.labelsize'])
 
-    for i, stage in enumerate(stages_feat):
-        signal = data[stage][:, ..., channel]
-        ax_full, ax_zoom = axes[i, :]
-        ax_full.plot(t_full, signal, c=colors_feat[i])
-        ax_full.set_ylim(signal.min(), signal.max())
+    # Convolutional filter responses:
+    signal = data['conv'][:, :, channel]
+    ax_full, ax_zoom = axes[0, :]
+    c, lw = colors['conv'], linewidths['conv']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask, :], c=c, lw=lw)
+    ax_full.set_ylim(signal.min(), signal.max())
 
-        ax_zoom.plot(t_zoom, signal[zoom_mask, ...], c=colors_feat[i])
+    # Binary responses:
+    signal = data['bi'][:, :, channel]
+    ax_full, ax_zoom = axes[1, :]
+    c, lw = colors['bi'], linewidths['bi']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask, :], c=c, lw=lw)
+    ax_full.set_ylim(signal.min(), signal.max())
+
+    # Finalized features:
+    signal = data['feat'][:, :, channel]
+    ax_full, ax_zoom = axes[2, :]
+    c, lw = colors['feat'], linewidths['feat']
+    ax_full.plot(t_full, signal, c=c, lw=lw)
+    ax_zoom.plot(t_zoom, signal[zoom_mask, :], c=c, lw=lw)
+    ax_full.set_ylim(0, 1)
+
+    # Posthoc adjustments:
     ax_full.set_xlim(t_full[0], t_full[-1])
+    ax_zoom.set_xlim(t_zoom[0], t_zoom[-1])
     plt.show()