diff --git a/figures/fig_invariance_full.pdf b/figures/fig_invariance_full.pdf
index faa9b76..0fe2f55 100644
Binary files a/figures/fig_invariance_full.pdf and b/figures/fig_invariance_full.pdf differ
diff --git a/figures/fig_invariance_log_hp.pdf b/figures/fig_invariance_log_hp.pdf
index 72cab27..6ea4be6 100644
Binary files a/figures/fig_invariance_log_hp.pdf and b/figures/fig_invariance_log_hp.pdf differ
diff --git a/python/fig_invariance_full.py b/python/fig_invariance_full.py
index 944ba69..0e522f4 100644
--- a/python/fig_invariance_full.py
+++ b/python/fig_invariance_full.py
@@ -4,9 +4,11 @@ import numpy as np
 import matplotlib.pyplot as plt
 from itertools import product
 from thunderhopper.modeltools import load_data
+from misc_functions import get_saturation
 from color_functions import load_colors
 from plot_functions import hide_axis, ylimits, xlabel, ylabel, title_subplot,\
-                           plot_line, plot_barcode, strip_zeros, time_bar, super_xlabel
+                           plot_line, plot_barcode, strip_zeros, time_bar,\
+                           letter_subplot, letter_subplots
 from IPython import embed
 
 def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs):
@@ -20,52 +22,70 @@ def plot_bi_snippets(axes, time, snippets, **kwargs):
         plot_barcode(ax, time, snippets[:, ..., i], **kwargs)
     return None
 
+def plot_curves(ax, scales, measures, fill_kwargs={}, **kwargs):
+    if measures.ndim == 1:
+        ax.plot(scales, measures, **kwargs)[0]
+        return measures
+    median_measure = np.median(measures, axis=1)
+    spread_measure = [np.percentile(measures, 25, axis=1),
+                      np.percentile(measures, 75, axis=1)]
+    ax.plot(scales, median_measure, **kwargs)[0]
+    ax.fill_between(scales, *spread_measure, **fill_kwargs)
+    return median_measure
+
+def show_saturation(ax, scales, measures, high=0.95, **kwargs):
+    high_ind = get_saturation(measures, high=high)[1]
+    return ax.plot(scales[high_ind], 0, transform=ax.get_xaxis_transform(),
+                   marker='o', ms=10, zorder=6, clip_on=False, **kwargs)
+
 # GENERAL SETTINGS:
 target = 'Omocestus_rufipes'
 data_paths = glob.glob(f'../data/inv/full/{target}*.npz')
-stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat']
+ref_path = '../data/inv/full/ref_measures.npz'
+save_path = '../figures/fig_invariance_full.pdf'
+stages = ['filt', 'env', 'log', 'inv', 'conv', 'feat']
 load_kwargs = dict(
     files=stages,
     keywords=['scales', 'snip', 'measure']
 )
-save_path = '../figures/fig_invariance_full.pdf'
 
 # GRAPH SETTINGS:
 fig_kwargs = dict(
-    figsize=(32/2.54, 16/2.54),
+    figsize=(32/2.54, 20/2.54),
 )
 super_grid_kwargs = dict(
-    nrows=1,
-    ncols=3,
+    nrows=2,
+    ncols=1,
     wspace=0,
     hspace=0,
     left=0,
     right=1,
     bottom=0,
-    top=1
+    top=1,
+    height_ratios=[3, 2]
 )
 subfig_specs = dict(
-    snip=(slice(None), slice(0, -1)),
-    big=(slice(None), -1),
+    snip=(0, 0),
+    big=(1, 0),
 )
 snip_grid_kwargs = dict(
     nrows=len(stages),
     ncols=None,
     wspace=0.1,
     hspace=0.4,
-    left=0.15,
+    left=0.08,
     right=0.95,
     bottom=0.08,
     top=0.95
 )
 big_grid_kwargs = dict(
     nrows=1,
-    ncols=1,
-    wspace=0,
+    ncols=3,
+    wspace=0.2,
     hspace=0,
-    left=0.2,
+    left=snip_grid_kwargs['left'],
     right=0.96,
-    bottom=0.08,
+    bottom=0.2,
     top=0.95
 )
 
@@ -79,9 +99,7 @@ fs = dict(
     bar=16,
 )
 colors = load_colors('../data/stage_colors.npz')
-colors['raw'] = "#000000"
 lw = dict(
-    raw=0.25,
     filt=0.25,
     env=0.25,
     log=0.25,
@@ -92,25 +110,17 @@ lw = dict(
     big=3
 )
 xlabels = dict(
-    snip='time [s]',
     big='scale $\\alpha$',
 )
 ylabels = dict(
-    raw='$x$',
     filt='$x_{\\text{filt}}$',
     env='$x_{\\text{env}}$',
-    log='$x_{\\text{log}}$',
+    log='$x_{\\text{db}}$',
     inv='$x_{\\text{inv}}$',
     conv='$c_i$',
     bi='$b_i$',
     feat='$f_i$',
-    big='norm. intensity measure'
-)
-xlab_snip_kwargs = dict(
-    y=0,
-    fontsize=fs['lab_norm'],
-    ha='center',
-    va='bottom',
+    big=['intensity', 'rel. intensity', 'norm. intensity']
 )
 xlab_big_kwargs = dict(
     y=0,
@@ -126,18 +136,17 @@ ylab_snip_kwargs = dict(
     va='center'
 )
 ylab_big_kwargs = dict(
-    x=0,
+    x=-0.12,
     fontsize=fs['lab_norm'],
     ha='center',
-    va='top',
+    va='bottom',
 )
 yloc = dict(
-    raw=500,
-    filt=500,
-    env=250,
-    log=25,
-    inv=10,
-    conv=1,
+    filt=3000,
+    env=1000,
+    log=50,
+    inv=20,
+    conv=2,
     feat=1,
 )
 title_kwargs = dict(
@@ -148,20 +157,18 @@ title_kwargs = dict(
     fontsize=fs['tit_norm'],
 )
 letter_snip_kwargs = dict(
-    x=0.02,
-    y=1,
+    x=0,
+    yref=0.5,
     ha='left',
-    va='top',
+    va='center',
     fontsize=fs['letter'],
-    fontweight='bold'
 )
 letter_big_kwargs = dict(
     x=0,
     y=1,
     ha='left',
-    va='top',
+    va='bottom',
     fontsize=fs['letter'],
-    fontweight='bold'
 )
 bar_time = 5
 bar_kwargs = dict(
@@ -181,6 +188,8 @@ bar_kwargs = dict(
     )
 )
 
+# PREPARATION:
+ref_data = dict(np.load(ref_path))
 
 # EXECUTION:
 for data_path in data_paths:
@@ -188,7 +197,7 @@ for data_path in data_paths:
 
     # Load invariance data:
     data, config = load_data(data_path, **load_kwargs)
-    t_full = np.arange(data['snip_raw'].shape[0]) / config['rate']
+    t_full = np.arange(data['snip_filt'].shape[0]) / config['rate']
 
     # Adjust grid parameters:
     snip_grid_kwargs['ncols'] = data['example_scales'].size
@@ -204,78 +213,91 @@ for data_path in data_paths:
     for i, j in product(range(snip_grid.nrows), range(snip_grid.ncols)):
         ax = snip_subfig.add_subplot(snip_grid[i, j])
         ax.set_xlim(t_full[0], t_full[-1])
+        ax.yaxis.set_major_locator(plt.MultipleLocator(yloc[stages[i]]))
         hide_axis(ax, 'bottom')
         if i == 0:
-            title = f'$\\alpha={strip_zeros(data["example_scales"][j])}$'
-            title_subplot(ax, title, ref=snip_subfig, **title_kwargs)
+            title = title_subplot(ax, f'$\\alpha={strip_zeros(data["example_scales"][j])}$',
+                                  ref=snip_subfig, **title_kwargs)
         if j == 0:
             ylabel(ax, ylabels[stages[i]], **ylab_snip_kwargs, transform=snip_subfig.transSubfigure)
         else:
             hide_axis(ax, 'left')
-        if stages[i] != 'bi':
-            ax.yaxis.set_major_locator(plt.MultipleLocator(yloc[stages[i]]))
         snip_axes[i, j] = ax
     time_bar(snip_axes[-1, -1], **bar_kwargs)
+    letter_subplot(snip_subfig, 'a', ref=title, **letter_snip_kwargs)
 
-    # Prepare single analysis axis:
+    # Prepare analysis axes:
     big_subfig = fig.add_subfigure(super_grid[subfig_specs['big']])
     big_grid = big_subfig.add_gridspec(**big_grid_kwargs)
-    big_ax = big_subfig.add_subplot(big_grid[0, 0])
-    big_ax.set_xlim(data['scales'].min(), data['scales'].max())
-    big_ax.set_xscale('symlog', linthresh=data['scales'][1], linscale=0.5)
-    big_ax.set_yscale('symlog', linthresh=0.01, linscale=0.1)
-    xlabel(big_ax, xlabels['big'], **xlab_big_kwargs, transform=big_subfig.transSubfigure)
-    ylabel(big_ax, ylabels['big'], **ylab_big_kwargs, transform=big_subfig.transSubfigure)
+    big_axes = np.zeros((big_grid.ncols,), dtype=object)
+    for i in range(big_grid.ncols):
+        ax = big_subfig.add_subplot(big_grid[0, i])
+        ax.set_xlim(data['scales'][0], data['scales'][-1])
+        ax.set_xscale('symlog', linthresh=data['scales'][1], linscale=0.5)
+        ax.set_yscale('symlog', linthresh=0.01, linscale=0.1)
+        xlabel(ax, xlabels['big'], transform=big_subfig, **xlab_big_kwargs)
+        ylabel(ax, ylabels['big'][i], **ylab_big_kwargs)
+        big_axes[i] = ax
+    letter_subplots(big_axes, 'bc', **letter_big_kwargs)
 
-    # Plot raw snippets:
-    plot_snippets(snip_axes[0, :], t_full, data['snip_raw'],
-                  c=colors['raw'], lw=lw['raw'])
+    # # Plot filtered snippets:
+    # plot_snippets(snip_axes[0, :], t_full, data['snip_filt'],
+    #               c=colors['filt'], lw=lw['filt'])
 
-    # Plot filtered snippets:
-    plot_snippets(snip_axes[1, :], t_full, data['snip_filt'],
-                  c=colors['filt'], lw=lw['filt'])
+    # # Plot envelope snippets:
+    # plot_snippets(snip_axes[1, :], t_full, data['snip_env'],
+    #               ymin=0, c=colors['env'], lw=lw['env'])
 
-    # Plot envelope snippets:
-    plot_snippets(snip_axes[2, :], t_full, data['snip_env'],
-                  ymin=0, c=colors['env'], lw=lw['env'])
+    # # Plot logarithmic snippets:
+    # plot_snippets(snip_axes[2, :], t_full, data['snip_log'],
+    #               c=colors['log'], lw=lw['log'])
 
-    # Plot logarithmic snippets:
-    plot_snippets(snip_axes[3, :], t_full, data['snip_log'],
-                  ymax=0, c=colors['log'], lw=lw['log'])
+    # # Plot invariant snippets:
+    # plot_snippets(snip_axes[3, :], t_full, data['snip_inv'],
+    #               c=colors['inv'], lw=lw['inv'])
 
-    # Plot invariant snippets:
-    plot_snippets(snip_axes[4, :], t_full, data['snip_inv'],
-                  c=colors['inv'], lw=lw['inv'])
+    # # Plot kernel response snippets:
+    # plot_snippets(snip_axes[4, :], t_full, data['snip_conv'],
+    #               c=colors['conv'], lw=lw['conv'])
 
-    # Plot kernel response snippets:
-    plot_snippets(snip_axes[5, :], t_full, data['snip_conv'],
-                  c=colors['conv'], lw=lw['conv'])
-
-    # Plot binary snippets:
-    plot_bi_snippets(snip_axes[6, :], t_full, data['snip_bi'],
-                     color=colors['bi'], lw=lw['bi'])
-
-    # Plot feature snippets:
-    plot_snippets(snip_axes[7, :], t_full, data['snip_feat'],
-                  ymin=0, ymax=1, c=colors['feat'], lw=lw['feat'])
+    # # Plot feature snippets:
+    # plot_snippets(snip_axes[5, :], t_full, data['snip_feat'],
+    #               ymin=0, ymax=1, c=colors['feat'], lw=lw['feat'])
 
     # Analysis results:
+    scales_rel = data['scales'] - data['scales'][0]
+    scales_rel /= scales_rel[-1]
     for stage in stages:
-        key = f'measure_{stage}'
-        if stage == 'bi':
-            continue
-        # Min-max normalization:
-        base_ind = np.argmin(data['scales'])
-        data[key] -= data[key][base_ind, ...]
-        data[key] /= data[key].max(axis=0)
+        measure = data[f'measure_{stage}']
 
-        # Condense measure:
-        if stage in ['conv', 'feat']:
-            data[key] = np.nanmedian(data[key], axis=1)
+        # Plot unmodified intensity measures:
+        curve = plot_curves(big_axes[0], data['scales'], measure, c=colors[stage], lw=lw['big'],
+                             fill_kwargs=dict(color=colors[stage], alpha=0.25))
+        if stage in ['log', 'inv', 'conv', 'feat']:
+            show_saturation(big_axes[0], data['scales'], curve, c=colors[stage])
 
-        # Plot measure over scales:
-        big_ax.plot(data['scales'], data[key],
-                    c=colors[stage], lw=lw['big'])
+        # # Relate to pure-noise reference:
+        # norm_measure = measure / ref_data[stage]
+
+        # # Plot noise-related intensity measures:
+        # big_axes[1].plot(data['scales'], norm_measure, c=colors[stage], lw=lw['big'])
+
+        # Normalize measure to [0, 1]:
+        min_measure = measure.min(axis=0)
+        max_measure = measure.max(axis=0)
+        norm_measure = (measure - min_measure) / (max_measure - min_measure)
+
+        # Plot normalized intensity measures:
+        curve = plot_curves(big_axes[1], data['scales'], norm_measure, c=colors[stage], lw=lw['big'],
+                             fill_kwargs=dict(color=colors[stage], alpha=0.25))
+        if stage in ['log', 'inv', 'conv', 'feat']:
+            show_saturation(big_axes[1], data['scales'], curve, c=colors[stage])
+
+        # # Plot over relative scales:
+        # plot_curves(big_axes[2], scales_rel, norm_measure, c=colors[stage], lw=lw['big'],
+        #             fill_kwargs=dict(color=colors[stage], alpha=0.25))
+        # scales_rel = curve - curve.min()
+        # scales_rel /= scales_rel.max()
 
     if save_path is not None:
         fig.savefig(save_path)
diff --git a/python/fig_invariance_log-hp.py b/python/fig_invariance_log-hp.py
index afe367f..3b0396c 100644
--- a/python/fig_invariance_log-hp.py
+++ b/python/fig_invariance_log-hp.py
@@ -4,10 +4,11 @@ import matplotlib.pyplot as plt
 from itertools import product
 from thunderhopper.filetools import search_files
 from thunderhopper.modeltools import load_data
+from misc_functions import shorten_species, get_kde, get_saturation
 from color_functions import load_colors
-from plot_functions import hide_axis, ylimits, xlabel, ylabel, hide_ticks,\
+from plot_functions import hide_axis, ylimits, super_xlabel, ylabel, hide_ticks,\
                            plot_line, strip_zeros, time_bar, zoom_inset,\
-                           letter_subplot, title_subplot
+                           letter_subplot, letter_subplots, title_subplot
 from IPython import embed
 
 def add_snip_axes(fig, grid_kwargs):
@@ -26,39 +27,86 @@ def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs):
         handles.extend(plot_line(ax, time, snippet, ymin=ymin, ymax=ymax, **kwargs))
     return handles
 
+def plot_dist_shifted(ax, data, axis, pdf=None, sigma=0.1, which='x',
+                      base=None, cap=None, add_pdf=False, shifted=False, **kwargs):
+    if pdf is None:
+        pdf, axis = get_kde(data, sigma, axis)
+    if base is None:
+        base = pdf.min()
+    if cap is None:
+        cap = pdf.max()
+    pdf = (pdf - pdf.min()) / (pdf.max() - pdf.min()) * (cap - base) + base
+
+    if which == 'x':
+        transform = ax.get_xaxis_transform()
+    elif which == 'y':
+        transform = ax.get_yaxis_transform()
+    else:
+        transform = ax.transData
+
+    rng = np.random.default_rng()
+    handles = []
+    for value in data:
+        ind = np.nonzero(axis == value)[0][0]
+        offset = base if not shifted else rng.uniform(base, pdf[ind])
+        variables = (offset, value) if which=='y' else (value, offset)
+        handles.extend(ax.plot(*variables, transform=transform, **kwargs))
+    if add_pdf:
+        variables = (pdf, axis) if which=='y' else (axis, pdf)
+        pdf_handle = ax.plot(*variables, transform=transform, c='k', lw=1)
+        return handles, pdf_handle
+    return handles
+
+def zalpha(handles, background='w', down=1):
+    twins = []
+    for handle in handles:
+        twin = handle.copy()
+        twin.set(color=background, alpha=1)
+        twin.set_zorder(handle.get_zorder() - down)
+        twins.append(twin)
+    return twins
 
 # GENERAL SETTINGS:
 target = 'Omocestus_rufipes'
 data_paths = search_files(target, excl='noise', dir='../data/inv/log_hp/')
-species_paths = search_files('*', incl='noise', dir='../data/inv/log_hp/')
+ref_path = '../data/inv/log_hp/ref_measures.npz'
+save_path = '../figures/fig_invariance_log_hp.pdf'
+target_species = [
+    'Omocestus_rufipes',
+    'Chorthippus_biguttulus',
+    'Chorthippus_mollis',
+    'Chrysochraon_dispar',
+    'Gomphocerippus_rufus',
+    'Pseudochorthippus_parallelus',
+    ]
 stages = ['env', 'log', 'inv']
 load_kwargs = dict(
     files=stages,
     keywords=['scales', 'snip', 'measure']
 )
-save_path = '../figures/fig_invariance_log_hp.pdf'
 compute_ratios = True
 show_diag = True
-show_noise = True
+show_plateaus = True
 
 # GRAPH SETTINGS:
 fig_kwargs = dict(
-    figsize=(32/2.54, 16/2.54),
+    figsize=(32/2.54, 32/2.54),
 )
 super_grid_kwargs = dict(
-    nrows=2,
-    ncols=3,
+    nrows=3,
+    ncols=1,
     wspace=0,
     hspace=0,
     left=0,
     right=1,
     bottom=0,
-    top=1
+    top=1,
+    height_ratios=[1, 1, 1]
 )
 subfig_specs = dict(
-    pure=(0, slice(0, -1)),
-    noise=(1, slice(0, -1)),
-    big=(slice(None), -1),
+    pure=(0, slice(None)),
+    noise=(1, slice(None)),
+    big=(2, slice(None)),
 )
 block_height = 0.8
 edge_padding = 0.08
@@ -67,7 +115,7 @@ pure_grid_kwargs = dict(
     ncols=None,
     wspace=0.1,
     hspace=0.15,
-    left=0.16,
+    left=0.11,
     right=0.95,
     bottom=1 - block_height - edge_padding,
     top=1 - edge_padding,
@@ -76,23 +124,23 @@ pure_grid_kwargs = dict(
 noise_grid_kwargs = dict(
     nrows=len(stages),
     ncols=None,
-    wspace=0.1,
-    hspace=0.15,
-    left=0.16,
-    right=0.95,
+    wspace=pure_grid_kwargs['wspace'],
+    hspace=pure_grid_kwargs['hspace'],
+    left=pure_grid_kwargs['left'],
+    right=pure_grid_kwargs['right'],
     bottom=edge_padding,
     top=edge_padding + block_height,
     height_ratios=[1, 2, 1]
 )
 big_grid_kwargs = dict(
-    nrows=2,
-    ncols=1,
-    wspace=0,
-    hspace=0.3,
-    left=0.19,
-    right=0.96,
-    bottom=0.09,
-    top=0.98
+    nrows=1,
+    ncols=3,
+    wspace=0.3,
+    hspace=0,
+    left=pure_grid_kwargs['left'],
+    right=pure_grid_kwargs['right'],
+    bottom=0.05,
+    top=1
 )
 anchor_kwargs = dict(
     aspect='equal',
@@ -110,8 +158,14 @@ fs = dict(
     bar=16,
 )
 colors = load_colors('../data/stage_colors.npz')
-lw_snippets = 1
-lw_big = 3
+species_colors = load_colors('../data/species_colors.npz')
+noise_colors = [(0.5, 0.5, 0.5), (0.7, 0.7, 0.7)]
+lw = dict(
+    snip=1,
+    big=4,
+    spec=2,
+    plateau=1,
+)
 xlabels = dict(
     big='scale $\\alpha$',
 )
@@ -135,7 +189,7 @@ ylab_snip_kwargs = dict(
     va='center',
 )
 ylab_big_kwargs = dict(
-    x=0.05,
+    x=0,
     fontsize=fs['lab_tex'],
     ha='center',
     va='top',
@@ -160,10 +214,10 @@ letter_snip_kwargs = dict(
     fontsize=fs['letter'],
 )
 letter_big_kwargs = dict(
-    x=0.05,
-    yref=letter_snip_kwargs['yref'],
+    x=0,
+    y=1,
     ha='left',
-    va='center',
+    va='bottom',
     fontsize=fs['letter'],
 )
 zoom_inset_bounds = [0.1, 0.2, 0.8, 0.6]
@@ -204,33 +258,77 @@ bar_kwargs = dict(
         va='center',
     )
 )
+leg_kwargs = dict(
+    ncols=2,
+    loc='upper right',
+    bbox_to_anchor=(0, 0.6, 1, 0.4),
+    frameon=False,
+    prop=dict(
+        size=12,
+        style='italic',
+    ),
+    borderpad=0,
+    borderaxespad=0,
+    handlelength=1,
+    columnspacing=1,
+)
 diag_kwargs = dict(
     c=(0.75, 0.75, 0.75),
     lw=2,
     ls='--',
     zorder=1.9,
 )
-noise_rel_thresh = 0.95
-noise_kwargs = dict(
-    fc=(0.9, 0.9, 0.9),
+plateau_settings = dict(
+    low=0.05,
+    high=0.95,
+    first=True,
+    last=True,
+    condense=None,
+)
+plateau_rect_kwargs = dict(
     ec='none',
     lw=0,
     zorder=1.5,
 )
+plateau_line_kwargs = dict(
+    lw=lw['plateau'],
+    ls='--',
+    zorder=1,
+)
+plateau_dot_kwargs = dict(
+    marker='o',
+    markersize=10,
+    markeredgecolor='k',
+    markeredgewidth=1,
+    # alpha=1,
+    zorder=6,
+    clip_on=False,
+    # base=0,
+    # cap=0.15,
+    # add_pdf=True,
+)
+kde_kwargs = dict(
+    sigma=0.1,
+)
 
 # PREPARATION:
 if compute_ratios:
-    ref_data = load_data('../data/processed/white_noise_sd-1.npz', files=stages)[0]
-    ref_measures = {k: v.std() for k, v in ref_data.items() if not k.endswith('rate')}
+    ref_measures = dict(np.load(ref_path))
 
-species_measures = []
-for species_path in species_paths:
-    species_data, _ = load_data(species_path, **load_kwargs)
-    species_measure = species_data['measure_inv']
+species_measures = {}
+thresh_inds = np.zeros((len(target_species),), dtype=int)
+thresh_scales = np.zeros((len(target_species),), dtype=float)
+for i, species in enumerate(target_species):
+    path = search_files(species, incl='noise', dir='../data/inv/log_hp/')[0]
+    species_data = load_data(path, **load_kwargs)[0]
+    measure = species_data['measure_inv']
+    scales = species_data['scales']
     if compute_ratios:
-        species_measure /= ref_measures['inv']
-    species_measures.append(species_measure)
-species_measures = np.array(species_measures).T
+        measure /= ref_measures['inv']
+    species_measures[species] = measure
+    thresh_inds[i] = get_saturation(measure, **plateau_settings)[1]
+    thresh_scales[i] = scales[thresh_inds[i]]
+thresh_pdf, pdf_axis = get_kde(thresh_scales, axis=scales, **kde_kwargs)
 
 # EXECUTION:
 for data_path in data_paths:
@@ -273,7 +371,7 @@ for data_path in data_paths:
                transform=noise_subfig.transSubfigure)
     for ax, scale in zip(noise_axes[0, :], noise_data['example_scales']):
         noise_title = title_subplot(ax, f'$\\alpha={strip_zeros(scale)}$', **title_kwargs)
-    letter_subplot(noise_subfig, 'c', ref=noise_title, **letter_snip_kwargs)
+    letter_subplot(noise_subfig, 'b', ref=noise_title, **letter_snip_kwargs)
     noise_inset = noise_axes[0, 0].inset_axes(zoom_inset_bounds)
     noise_inset.spines[:].set(visible=True, lw=zoom_kwargs['lw'])
     noise_inset.tick_params(**inset_tick_kwargs)
@@ -282,51 +380,49 @@ for data_path in data_paths:
     # Prepare analysis axes:
     big_subfig = fig.add_subfigure(super_grid[subfig_specs['big']])
     big_grid = big_subfig.add_gridspec(**big_grid_kwargs)
-    big_axes = np.zeros((big_grid.nrows,), dtype=object)
-    for i, scales in enumerate([pure_scales, noise_scales]):
-        ax = big_subfig.add_subplot(big_grid[i, 0])
+    big_axes = np.zeros((big_grid.ncols,), dtype=object)
+    for i, scales in enumerate([pure_scales, noise_scales, noise_scales]):
+        ax = big_subfig.add_subplot(big_grid[0, i])
         ax.set_xlim(scales[0], scales[-1])
         ax.set_ylim(scales[0], scales[-1])
         ax.set_xscale('symlog', linthresh=scales[1], linscale=0.5)
         ax.set_yscale('symlog', linthresh=scales[1], linscale=0.5)
         ax.set_aspect(**anchor_kwargs)
-        ylabel(ax, ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
-        if i == 0:
-            hide_ticks(ax, 'bottom')
-            letter_subplot(big_subfig, 'b', ref=pure_title, **letter_big_kwargs)
-        else:
-            xlabel(ax, xlabels['big'], transform=big_subfig.transSubfigure, **xlab_big_kwargs)
-            letter_subplot(big_subfig, 'd', ref=noise_title, **letter_big_kwargs)
+        if i > 0:
+            hide_ticks(ax, 'left')
         big_axes[i] = ax
+    ylabel(big_axes[0], ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
+    super_xlabel(xlabels['big'], big_subfig, big_axes[0], big_axes[-1], **xlab_big_kwargs)
+    letter_subplots(big_axes, 'cde', **letter_big_kwargs)
 
     # Plot pure-song envelope snippets:
     handle = plot_snippets(pure_axes[0, :], t_full, pure_data['snip_env'],
-                           ymin=0, c=colors['env'], lw=lw_snippets)[0]
+                           ymin=0, c=colors['env'], lw=lw['snip'])[0]
     zoom_inset(pure_axes[0, 0], pure_inset, handle, transform=pure_axes[0, 0].transAxes, **zoom_kwargs)
 
     # Plot pure-song logarithmic snippets:
     plot_snippets(pure_axes[1, :], t_full, pure_data['snip_log'],
-                  c=colors['log'], lw=lw_snippets)
+                  c=colors['log'], lw=lw['snip'])
 
     # Plot pure-song invariant snippets:
     plot_snippets(pure_axes[2, :], t_full, pure_data['snip_inv'],
-                  c=colors['inv'], lw=lw_snippets)
+                  c=colors['inv'], lw=lw['snip'])
 
     # Plot noise-song envelope snippets:
     ymin, ymax = pure_axes[0, 0].get_ylim()
     handle = plot_snippets(noise_axes[0, :], t_full, noise_data['snip_env'],
-                           ymin, ymax, c=colors['env'], lw=lw_snippets)[0]
+                           ymin, ymax, c=colors['env'], lw=lw['snip'])[0]
     zoom_inset(noise_axes[0, 0], noise_inset, handle, transform=noise_axes[0, 0].transAxes, **zoom_kwargs)
 
     # Plot noise-song logarithmic snippets:
     ymin, ymax = pure_axes[1, 0].get_ylim()
     plot_snippets(noise_axes[1, :], t_full, noise_data['snip_log'],
-                  ymin, ymax, c=colors['log'], lw=lw_snippets)
+                  ymin, ymax, c=colors['log'], lw=lw['snip'])
 
     # Plot noise-song invariant snippets:
     ymin, ymax = pure_axes[2, 0].get_ylim()
     plot_snippets(noise_axes[2, :], t_full, noise_data['snip_inv'],
-                  ymin, ymax, c=colors['inv'], lw=lw_snippets)
+                  ymin, ymax, c=colors['inv'], lw=lw['snip'])
 
     # Indicate time scale:
     time_bar(noise_axes[-1, -1], **bar_kwargs)
@@ -342,34 +438,46 @@ for data_path in data_paths:
         noise_data['measure_inv'] /= ref_measures['inv']
 
     # Plot pure-song measures (ideal):
-    big_axes[0].plot(pure_scales, pure_data['measure_env'], c=colors['env'], lw=lw_big)
-    big_axes[0].plot(pure_scales, pure_data['measure_log'], c=colors['log'], lw=lw_big)
-    big_axes[0].plot(pure_scales, pure_data['measure_inv'], c=colors['inv'], lw=lw_big)
+    big_axes[0].plot(pure_scales, pure_data['measure_env'], c=colors['env'], lw=lw['big'])
+    big_axes[0].plot(pure_scales, pure_data['measure_log'], c=colors['log'], lw=lw['big'])
+    big_axes[0].plot(pure_scales, pure_data['measure_inv'], c=colors['inv'], lw=lw['big'])
 
     # Plot noise-song measures (limited):
-    big_axes[1].plot(noise_scales, noise_data['measure_env'], c=colors['env'], lw=lw_big)
-    big_axes[1].plot(noise_scales, noise_data['measure_log'], c=colors['log'], lw=lw_big)
-    big_axes[1].plot(noise_scales, noise_data['measure_inv'], c=colors['inv'], lw=lw_big)
-
-    # Plot species measures:
-    big_axes[1].plot(noise_scales, species_measures, 'k', lw=lw_big, zorder=2.1)
+    big_axes[1].plot(noise_scales, noise_data['measure_env'], c=colors['env'], lw=lw['big'])
+    big_axes[1].plot(noise_scales, noise_data['measure_log'], c=colors['log'], lw=lw['big'])
+    big_axes[1].plot(noise_scales, noise_data['measure_inv'], c=colors['inv'], lw=lw['big'])
 
     if show_diag:
         # Indicate diagonal:
         big_axes[0].plot(pure_scales, pure_scales, **diag_kwargs)
         big_axes[1].plot(noise_scales, noise_scales, **diag_kwargs)
 
-    if show_noise:
-        # Indicate noise floor:
-        if compute_ratios:
-            span_measure = noise_data['measure_inv'][-1] - ref_measures['inv']
-            thresh_measure = ref_measures['inv'] + noise_rel_thresh * span_measure
-        else:     
-            span_measure = noise_data['measure_inv'][-1] - noise_data['measure_inv'][0]
-            thresh_measure = noise_data['measure_inv'][0] + noise_rel_thresh * span_measure
-        thresh_ind = np.nonzero(noise_data['measure_inv'] < thresh_measure)[0][-1]
-        thresh_scale = noise_scales[thresh_ind]
-        big_axes[1].axvspan(noise_scales[0], thresh_scale, **noise_kwargs)
+    if show_plateaus:
+        # Indicate low and high plateaus of noise invariance curve:
+        low_ind, high_ind = get_saturation(noise_data['measure_inv'], **plateau_settings)
+        big_axes[1].axvspan(noise_scales[0], noise_scales[low_ind],
+                            fc=noise_colors[0], **plateau_rect_kwargs)
+        big_axes[1].axvspan(noise_scales[low_ind], noise_scales[high_ind],
+                            fc=noise_colors[1], **plateau_rect_kwargs)
+
+    # Plot species-specific noise-song measures:
+    for i, (species, measure) in enumerate(species_measures.items()):
+        color = species_colors[species]
+        ind, scale = thresh_inds[i], thresh_scales[i]
+        big_axes[2].plot(noise_scales, measure, label=shorten_species(species),
+                         c=color, lw=lw['spec'])
+        big_axes[2].plot(scale, 0, c='w', alpha=1, zorder=5.5, **plateau_dot_kwargs,
+                         transform=big_axes[2].get_xaxis_transform()) 
+        handle = big_axes[2].plot(scale, 0, c=color, alpha=0.5, **plateau_dot_kwargs,
+                                  transform=big_axes[2].get_xaxis_transform())
+        big_axes[2].vlines(scale, big_axes[2].get_ylim()[0], measure[ind],
+                           color=color, **plateau_line_kwargs)
+    big_axes[2].legend(**leg_kwargs)
+
+    # handles = plot_dist_shifted(big_axes[2], species_threshs, axis=pdf_axis,
+    #                             pdf=thresh_pdf, **plateau_dot_kwargs)[0]
+    # [h.set_color(species_colors[s]) for h, s in zip(handles, target_species)]
+                      
 
     if save_path is not None:
         fig.savefig(save_path, bbox_inches='tight')
diff --git a/python/fig_invariance_log-hp_backup.py b/python/fig_invariance_log-hp_backup.py
deleted file mode 100644
index acda691..0000000
--- a/python/fig_invariance_log-hp_backup.py
+++ /dev/null
@@ -1,407 +0,0 @@
-import plotstyle_plt
-import numpy as np
-import matplotlib.pyplot as plt
-from itertools import product
-from thunderhopper.filetools import search_files
-from thunderhopper.modeltools import load_data
-from misc_functions import shorten_species, get_saturation
-from color_functions import load_colors
-from plot_functions import hide_axis, ylimits, super_xlabel, ylabel, hide_ticks,\
-                           plot_line, strip_zeros, time_bar, zoom_inset,\
-                           letter_subplot, letter_subplots, title_subplot
-from IPython import embed
-
-def add_snip_axes(fig, grid_kwargs):
-    grid = fig.add_gridspec(**grid_kwargs)
-    axes = np.zeros((grid.nrows, grid.ncols), dtype=object)
-    for i, j in product(range(grid.nrows), range(grid.ncols)):
-        axes[i, j] = fig.add_subplot(grid[i, j])
-    [hide_axis(ax, 'left') for ax in axes[:, 1:].flatten()]
-    [hide_axis(ax, 'bottom') for ax in axes.flatten()]
-    return axes
-
-def plot_snippets(axes, time, snippets, ymin=None, ymax=None, **kwargs):
-    ymin, ymax = ylimits(snippets, minval=ymin, maxval=ymax, pad=0.05)
-    handles = []
-    for ax, snippet in zip(axes, snippets.T):
-        handles.extend(plot_line(ax, time, snippet, ymin=ymin, ymax=ymax, **kwargs))
-    return handles
-
-
-# GENERAL SETTINGS:
-target = 'Omocestus_rufipes'
-data_paths = search_files(target, excl='noise', dir='../data/inv/log_hp/')
-target_species = [
-    'Omocestus_rufipes',
-    'Chorthippus_biguttulus',
-    'Chorthippus_mollis',
-    'Chrysochraon_dispar',
-    'Gomphocerippus_rufus',
-    'Pseudochorthippus_parallelus',
-    ]
-stages = ['env', 'log', 'inv']
-load_kwargs = dict(
-    files=stages,
-    keywords=['scales', 'snip', 'measure']
-)
-save_path = '../figures/fig_invariance_log_hp.pdf'
-compute_ratios = True
-show_diag = True
-show_plateaus = True
-
-# GRAPH SETTINGS:
-fig_kwargs = dict(
-    figsize=(32/2.54, 32/2.54),
-)
-# snip_rows = 1
-# big_rows = 1
-super_grid_kwargs = dict(
-    nrows=3,
-    ncols=1,
-    wspace=0,
-    hspace=0,
-    left=0,
-    right=1,
-    bottom=0,
-    top=1,
-    height_ratios=[1, 1, 1]
-)
-subfig_specs = dict(
-    pure=(0, slice(None)),
-    noise=(1, slice(None)),
-    big=(2, slice(None)),
-)
-block_height = 0.8
-edge_padding = 0.08
-pure_grid_kwargs = dict(
-    nrows=len(stages),
-    ncols=None,
-    wspace=0.1,
-    hspace=0.15,
-    left=0.11,
-    right=0.95,
-    bottom=1 - block_height - edge_padding,
-    top=1 - edge_padding,
-    height_ratios=[1, 2, 1]
-)
-noise_grid_kwargs = dict(
-    nrows=len(stages),
-    ncols=None,
-    wspace=pure_grid_kwargs['wspace'],
-    hspace=pure_grid_kwargs['hspace'],
-    left=pure_grid_kwargs['left'],
-    right=pure_grid_kwargs['right'],
-    bottom=edge_padding,
-    top=edge_padding + block_height,
-    height_ratios=[1, 2, 1]
-)
-big_grid_kwargs = dict(
-    nrows=1,
-    ncols=3,
-    wspace=0.3,
-    hspace=0,
-    left=pure_grid_kwargs['left'],
-    right=pure_grid_kwargs['right'],
-    bottom=0.05,
-    top=1
-)
-anchor_kwargs = dict(
-    aspect='equal',
-    adjustable='box',
-    anchor=(0.5, 0.5)
-)
-
-# PLOT SETTINGS:
-fs = dict(
-    lab_norm=16,
-    lab_tex=20,
-    letter=22,
-    tit_norm=16,
-    tit_tex=20,
-    bar=16,
-)
-colors = load_colors('../data/stage_colors.npz')
-species_colors = load_colors('../data/species_colors.npz')
-noise_colors = [(0.5, 0.5, 0.5), (0.7, 0.7, 0.7)]
-lw_snippets = 1
-lw_big = 3
-xlabels = dict(
-    big='scale $\\alpha$',
-)
-ylabels = dict(
-    env='$x_{\\text{env}}$',
-    log='$x_{\\text{dB}}$',
-    inv='$x_{\\text{adapt}}$',
-    big='$\\sigma_{\\alpha}\\,/\\,\\sigma_{\\eta}$',
-)
-xlab_big_kwargs = dict(
-    y=0,
-    fontsize=fs['lab_norm'],
-    ha='center',
-    va='bottom',
-)
-ylab_snip_kwargs = dict(
-    x=0,
-    fontsize=fs['lab_tex'],
-    rotation=0,
-    ha='left',
-    va='center',
-)
-ylab_big_kwargs = dict(
-    x=0,
-    fontsize=fs['lab_tex'],
-    ha='center',
-    va='top',
-)
-yloc = dict(
-    env=1000,
-    log=40,
-    inv=20
-)
-title_kwargs = dict(
-    x=0.5,
-    y=1,
-    ha='center',
-    va='bottom',
-    fontsize=fs['tit_norm'],
-)
-letter_snip_kwargs = dict(
-    x=0,
-    yref=0.5,
-    ha='left',
-    va='center',
-    fontsize=fs['letter'],
-)
-letter_big_kwargs = dict(
-    x=0,
-    y=1,
-    ha='left',
-    va='bottom',
-    fontsize=fs['letter'],
-)
-zoom_inset_bounds = [0.1, 0.2, 0.8, 0.6]
-zoom_kwargs = dict(
-    x0=0.45,
-    x1=0.55,
-    y0=0,
-    y1=0.0006,
-    low_left=True,
-    low_right=True,
-    ec='k',
-    lw=1,
-    alpha=1,
-)
-inset_tick_kwargs = dict(
-    axis='y',
-    length=3,
-    pad=1,
-    left=False,
-    labelleft=False,
-    right=True,
-    labelright=True,
-)
-bar_time = 5
-bar_kwargs = dict(
-    dur=bar_time,
-    y0=-0.25,
-    y1=-0.1,
-    xshift=1,
-    color='k',
-    lw=0,
-    clip_on=False,
-    text_pos=(-0.1, 0.5),
-    text_str=f'${bar_time}\\,\\text{{s}}$',
-    text_kwargs=dict(
-        fontsize=fs['bar'],
-        ha='right',
-        va='center',
-    )
-)
-leg_kwargs = dict(
-    ncols=2,
-    loc='upper right',
-    bbox_to_anchor=(0, 0.6, 1, 0.4),
-    frameon=False,
-    prop=dict(
-        size=12,
-        style='italic',
-    ),
-    borderpad=0,
-    borderaxespad=0,
-    handlelength=1,
-    columnspacing=1,
-)
-diag_kwargs = dict(
-    c=(0.75, 0.75, 0.75),
-    lw=2,
-    ls='--',
-    zorder=1.9,
-)
-plateau_settings = dict(
-    low=0.05,
-    high=0.95,
-    first=True,
-    last=True,
-    condense=None,
-)
-plateau_kwargs = dict(
-    ec='none',
-    lw=0,
-    zorder=1.5,
-)
-
-# PREPARATION:
-if compute_ratios:
-    ref_data = load_data('../data/processed/white_noise_sd-1.npz', files=stages)[0]
-    ref_measures = {k: v.std() for k, v in ref_data.items() if not k.endswith('rate')}
-
-species_measures = {}
-for species in target_species:
-    path = search_files(species, incl='noise', dir='../data/inv/log_hp/')[0]
-    measure = load_data(path, **load_kwargs)[0]['measure_inv']
-    if compute_ratios:
-        measure /= ref_measures['inv']
-    species_measures[species] = measure
-
-# EXECUTION:
-for data_path in data_paths:
-    print(f'Processing {data_path}')
-
-    # Load invariance data:
-    pure_data, config = load_data(data_path, **load_kwargs)
-    noise_data, _ = load_data(data_path.replace('.npz', '_noise.npz'), **load_kwargs)
-    pure_scales, noise_scales = pure_data['scales'], noise_data['scales']
-    t_full = np.arange(pure_data['snip_env'].shape[0]) / config['env_rate']
-
-    # Prepare overall graph:
-    fig = plt.figure(**fig_kwargs)
-    super_grid = fig.add_gridspec(**super_grid_kwargs)
-    fig.canvas.draw()
-
-    # Prepare pure-song snippet axes:
-    pure_grid_kwargs['ncols'] = pure_data['example_scales'].size
-    pure_subfig = fig.add_subfigure(super_grid[subfig_specs['pure']])
-    pure_axes = add_snip_axes(pure_subfig, pure_grid_kwargs)
-    for ax, stage in zip(pure_axes[:, 0], stages):
-        ax.yaxis.set_major_locator(plt.MultipleLocator(yloc[stage]))
-        ylabel(ax, ylabels[stage], **ylab_snip_kwargs,
-               transform=pure_subfig.transSubfigure)
-    for ax, scale in zip(pure_axes[0, :], pure_data['example_scales']):
-        pure_title = title_subplot(ax, f'$\\alpha={strip_zeros(scale)}$', **title_kwargs)
-    letter_subplot(pure_subfig, 'a', ref=pure_title, **letter_snip_kwargs)
-    pure_inset = pure_axes[0, 0].inset_axes(zoom_inset_bounds)
-    pure_inset.spines[:].set(visible=True, lw=zoom_kwargs['lw'])
-    pure_inset.tick_params(**inset_tick_kwargs)
-    hide_ticks(pure_inset, 'bottom', ticks=False)
-
-    # Prepare noise-song snippet axes:
-    noise_grid_kwargs['ncols'] = noise_data['example_scales'].size
-    noise_subfig = fig.add_subfigure(super_grid[subfig_specs['noise']])
-    noise_axes = add_snip_axes(noise_subfig, noise_grid_kwargs)
-    for ax, stage in zip(noise_axes[:, 0], stages):
-        ax.yaxis.set_major_locator(plt.MultipleLocator(yloc[stage]))
-        ylabel(ax, ylabels[stage], **ylab_snip_kwargs,
-               transform=noise_subfig.transSubfigure)
-    for ax, scale in zip(noise_axes[0, :], noise_data['example_scales']):
-        noise_title = title_subplot(ax, f'$\\alpha={strip_zeros(scale)}$', **title_kwargs)
-    letter_subplot(noise_subfig, 'b', ref=noise_title, **letter_snip_kwargs)
-    noise_inset = noise_axes[0, 0].inset_axes(zoom_inset_bounds)
-    noise_inset.spines[:].set(visible=True, lw=zoom_kwargs['lw'])
-    noise_inset.tick_params(**inset_tick_kwargs)
-    hide_ticks(noise_inset, 'bottom', ticks=False)
-
-    # Prepare analysis axes:
-    big_subfig = fig.add_subfigure(super_grid[subfig_specs['big']])
-    big_grid = big_subfig.add_gridspec(**big_grid_kwargs)
-    big_axes = np.zeros((big_grid.ncols,), dtype=object)
-    for i, scales in enumerate([pure_scales, noise_scales, noise_scales]):
-        ax = big_subfig.add_subplot(big_grid[0, i])
-        ax.set_xlim(scales[0], scales[-1])
-        ax.set_ylim(scales[0], scales[-1])
-        ax.set_xscale('symlog', linthresh=scales[1], linscale=0.5)
-        ax.set_yscale('symlog', linthresh=scales[1], linscale=0.5)
-        ax.set_aspect(**anchor_kwargs)
-        big_axes[i] = ax
-    ylabel(big_axes[0], ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
-    super_xlabel(xlabels['big'], big_subfig, big_axes[0], big_axes[-1], **xlab_big_kwargs)
-    letter_subplots(big_axes, 'cde', **letter_big_kwargs)
-
-    # Plot pure-song envelope snippets:
-    handle = plot_snippets(pure_axes[0, :], t_full, pure_data['snip_env'],
-                           ymin=0, c=colors['env'], lw=lw_snippets)[0]
-    zoom_inset(pure_axes[0, 0], pure_inset, handle, transform=pure_axes[0, 0].transAxes, **zoom_kwargs)
-
-    # Plot pure-song logarithmic snippets:
-    plot_snippets(pure_axes[1, :], t_full, pure_data['snip_log'],
-                  c=colors['log'], lw=lw_snippets)
-
-    # Plot pure-song invariant snippets:
-    plot_snippets(pure_axes[2, :], t_full, pure_data['snip_inv'],
-                  c=colors['inv'], lw=lw_snippets)
-
-    # Plot noise-song envelope snippets:
-    ymin, ymax = pure_axes[0, 0].get_ylim()
-    handle = plot_snippets(noise_axes[0, :], t_full, noise_data['snip_env'],
-                           ymin, ymax, c=colors['env'], lw=lw_snippets)[0]
-    zoom_inset(noise_axes[0, 0], noise_inset, handle, transform=noise_axes[0, 0].transAxes, **zoom_kwargs)
-
-    # Plot noise-song logarithmic snippets:
-    ymin, ymax = pure_axes[1, 0].get_ylim()
-    plot_snippets(noise_axes[1, :], t_full, noise_data['snip_log'],
-                  ymin, ymax, c=colors['log'], lw=lw_snippets)
-
-    # Plot noise-song invariant snippets:
-    ymin, ymax = pure_axes[2, 0].get_ylim()
-    plot_snippets(noise_axes[2, :], t_full, noise_data['snip_inv'],
-                  ymin, ymax, c=colors['inv'], lw=lw_snippets)
-
-    # Indicate time scale:
-    time_bar(noise_axes[-1, -1], **bar_kwargs)
-
-    if compute_ratios:
-        # Relate pure-song measures to zero scale:
-        pure_data['measure_env'] /= ref_measures['env']
-        pure_data['measure_log'] /= ref_measures['log']
-        pure_data['measure_inv'] /= ref_measures['inv']
-        # Relate noise-song measures to zero scale:
-        noise_data['measure_env'] /= ref_measures['env']
-        noise_data['measure_log'] /= ref_measures['log']
-        noise_data['measure_inv'] /= ref_measures['inv']
-
-    # Plot pure-song measures (ideal):
-    big_axes[0].plot(pure_scales, pure_data['measure_env'], c=colors['env'], lw=lw_big)
-    big_axes[0].plot(pure_scales, pure_data['measure_log'], c=colors['log'], lw=lw_big)
-    big_axes[0].plot(pure_scales, pure_data['measure_inv'], c=colors['inv'], lw=lw_big)
-
-    # Plot noise-song measures (limited):
-    big_axes[1].plot(noise_scales, noise_data['measure_env'], c=colors['env'], lw=lw_big)
-    big_axes[1].plot(noise_scales, noise_data['measure_log'], c=colors['log'], lw=lw_big)
-    big_axes[1].plot(noise_scales, noise_data['measure_inv'], c=colors['inv'], lw=lw_big)
-
-
-    if show_diag:
-        # Indicate diagonal:
-        big_axes[0].plot(pure_scales, pure_scales, **diag_kwargs)
-        big_axes[1].plot(noise_scales, noise_scales, **diag_kwargs)
-
-    if show_plateaus:
-        # Indicate low and high plateaus of noise invariance curve:
-        low_ind, high_ind = get_saturation(noise_data['measure_inv'], **plateau_settings)
-        big_axes[1].axvspan(noise_scales[0], noise_scales[low_ind],
-                             fc=noise_colors[0], **plateau_kwargs)
-        big_axes[1].axvspan(noise_scales[low_ind], noise_scales[high_ind],
-                             fc=noise_colors[1], **plateau_kwargs)
-
-    # Plot species-specific noise-song measures:
-    for species, measure in species_measures.items():
-        label = shorten_species(species)
-        big_axes[2].plot(noise_scales, measure, label=label,
-                         c=species_colors[species], lw=lw_big)
-    big_axes[2].legend(**leg_kwargs)
-
-    if save_path is not None:
-        fig.savefig(save_path, bbox_inches='tight')
-    plt.show()
-
-print('Done.')
-embed()
diff --git a/python/fig_invariance_thresh-lp_single.py b/python/fig_invariance_thresh-lp_single.py
index 274acf8..e794ab8 100644
--- a/python/fig_invariance_thresh-lp_single.py
+++ b/python/fig_invariance_thresh-lp_single.py
@@ -5,6 +5,7 @@ from itertools import product
 from thunderhopper.filetools import search_files
 from thunderhopper.modeltools import load_data
 from thunderhopper.filtertools import find_kern_specs
+from misc_functions import get_saturation
 from color_functions import load_colors, shade_colors
 from plot_functions import shift_subplot, hide_axis, ylimits, xlabel, ylabel,\
                            super_ylabel, plot_line, plot_barcode, strip_zeros,\
@@ -64,7 +65,7 @@ load_kwargs = dict(
     files=stages,
     keywords=['scales', 'snip', 'measure', 'thresh']
 )
-save_path = '../figures/fig_invariance_thresh_lp_single.pdf'
+save_path = None#'../figures/fig_invariance_thresh_lp_single.pdf'
 exclude_zero = True
 
 # GRAPH SETTINGS:
@@ -79,7 +80,7 @@ super_grid_kwargs = dict(
     left=0,
     right=1,
     bottom=0,
-    top=1
+    top=1,
 )
 input_rows = 1
 snip_rows = 2
@@ -111,10 +112,10 @@ input_grid_kwargs = dict(
     top=0.75,
 )
 big_grid_kwargs = dict(
-    nrows=1,
+    nrows=2,
     ncols=1,
     wspace=0,
-    hspace=0,
+    hspace=0.3,
     left=0.17,
     right=0.96,
     bottom=0.1,
@@ -141,7 +142,8 @@ lw = dict(
     big=4,
 )
 xlabels = dict(
-    big='scale $\\alpha$',
+    alpha='scale $\\alpha$',
+    sigma='$\\sigma_{\\text{adapt}}$',
 )
 ylabels = dict(
     inv='$x_{\\text{adapt}}$',
@@ -150,11 +152,17 @@ ylabels = dict(
     feat='$f_i$',
     big='$\\mu_f$',
 )
-xlab_big_kwargs = dict(
-    y=0,
+xlab_alpha_kwargs = dict(
+    y=-0.15,
     fontsize=fs['lab_norm'],
     ha='center',
-    va='bottom',
+    va='top',
+)
+xlab_sigma_kwargs = dict(
+    y=-0.12,
+    fontsize=fs['lab_tex'],
+    ha=xlab_alpha_kwargs['ha'],
+    va=xlab_alpha_kwargs['va'],
 )
 ylab_snip_kwargs = dict(
     x=0.08,
@@ -178,7 +186,7 @@ ylab_big_kwargs = dict(
 ypad = dict(
     inv=0.05,
     conv=0.05,
-    big=0.075
+    big=0.1
 )
 yloc = dict(
     inv=(2, 200),
@@ -242,6 +250,13 @@ leg_kwargs = dict(
     handlelength=1.5,
     columnspacing=1,
 )
+plateau_settings = dict(
+    low=0.05,
+    high=0.95,
+    first=True,
+    last=True,
+    condense=None,
+)
 kern_specs = np.array([
     [1, 0.008],
     [2, 0.004],
@@ -281,6 +296,7 @@ for data_path in data_paths:
         # Reduce to nonzero scales:
         nonzero_inds = scales > 0
         scales = scales[nonzero_inds]
+        noise_data['measure_inv'] = noise_data['measure_inv'][nonzero_inds]
         noise_data['measure_feat'] = noise_data['measure_feat'][nonzero_inds, :]
         pure_data['measure_feat'] = pure_data['measure_feat'][nonzero_inds, :]
 
@@ -293,7 +309,7 @@ for data_path in data_paths:
     )
 
     # Adjust grid parameters to loaded data:
-    super_grid_kwargs['nrows'] = snip_rows * thresh_rel.size + 1
+    super_grid_kwargs['nrows'] = snip_rows * thresh_rel.size + input_rows
     input_grid_kwargs['ncols'] = plot_scales.size
     snip_grid_kwargs['ncols'] = plot_scales.size
 
@@ -325,8 +341,6 @@ for data_path in data_paths:
             ax1.yaxis.set_major_locator(plt.MultipleLocator(yloc[stage][0]))
             ax2.yaxis.set_major_locator(plt.MultipleLocator(yloc[stage][1]))
             ylabel(ax1, ylabels[stage], transform=snip_subfig.transSubfigure, **ylab_snip_kwargs)
-        # for ax, scale in zip(axes[0, :], plot_scales):
-        #     title_subplot(ax, f'$\\alpha={strip_zeros(scale)}$', ref=snip_subfig, **title_kwargs)
         if i == thresh_rel.size - 1:
             axes[-1, -1].set_xlim(t_full[0], t_full[-1])
             time_bar(axes[-1, -1], **bar_kwargs)
@@ -334,17 +348,27 @@ for data_path in data_paths:
         snip_axes.append(axes)
     letter_subplots(snip_subfigs, 'bcd', **letter_snip_kwargs)
 
-    # Prepare analysis axis:
+    # Prepare analysis axes:
     big_subfig = fig.add_subfigure(super_grid[subfig_specs['big']])
     big_grid = big_subfig.add_gridspec(**big_grid_kwargs)
-    big_ax = big_subfig.add_subplot(big_grid[0, 0])
-    big_ax.set_xlim(scales[0], scales[-1])
-    big_ax.set_xscale('symlog', linthresh=scales[scales > 0][0], linscale=0.5)
-    ylimits(noise_data['measure_feat'], big_ax, minval=0, pad=ypad['big'])
-    big_ax.yaxis.set_major_locator(plt.MultipleLocator(yloc['big']))
-    xlabel(big_ax, xlabels['big'], transform=big_subfig.transSubfigure, **xlab_big_kwargs)
-    ylabel(big_ax, ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
-    letter_subplot(big_subfig, 'e', **letter_big_kwargs, ref=input_subfig)
+
+    alpha_ax = big_subfig.add_subplot(big_grid[0, 0])
+    alpha_ax.set_xlim(scales[0], scales[-1])
+    alpha_ax.set_xscale('symlog', linthresh=scales[scales > 0][0], linscale=0.5)
+    ylimits(pure_data['measure_feat'], alpha_ax, minval=0, pad=ypad['big'])
+    alpha_ax.yaxis.set_major_locator(plt.MultipleLocator(yloc['big']))
+    xlabel(alpha_ax, xlabels['alpha'], **xlab_alpha_kwargs)
+    ylabel(alpha_ax, ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
+
+    sigma_ax = big_subfig.add_subplot(big_grid[1, 0])
+    sigma_ax.set_xlim(noise_data['measure_inv'].min(), noise_data['measure_inv'].max())
+    # sigma_ax.set_xscale('log')
+    sigma_ax.set_xlim(scales[0], scales[-1])
+    sigma_ax.set_xscale('symlog', linthresh=scales[scales > 0][0], linscale=0.5)
+    ylimits(pure_data['measure_feat'], sigma_ax, minval=0, pad=ypad['big'])
+    sigma_ax.yaxis.set_major_locator(plt.MultipleLocator(yloc['big']))
+    xlabel(sigma_ax, xlabels['sigma'], **xlab_sigma_kwargs)
+    ylabel(sigma_ax, ylabels['big'], transform=big_subfig.transSubfigure, **ylab_big_kwargs)
 
     # Plot intensity-adapted snippets:
     plot_snippets(input_axes, t_full, noise_data['snip_inv'],
@@ -375,18 +399,25 @@ for data_path in data_paths:
                                 ymin=0, ymax=1, c=shaded['feat'][i], lw=lw['feat'])
         [set_clip_box(h[0], ax, bounds=[[0, -0.05], [1, 1.05]]) for h, ax in zip(handles, axes[2, :])]
 
-    # Plot pure-song analysis results:
-    handles = big_ax.plot(scales, pure_data['measure_feat'], lw=lw['big'], ls='dotted')
-    [h.set_color(c) for h, c in zip(handles, shaded['feat'])]
+    # Get threshold-specific saturation:
+    for i in range(thresh_rel.size):
+        ind = get_saturation(noise_data['measure_feat'][:, i], **plateau_settings)[1]
 
-    # Plot noise-song analysis results:
-    handles = big_ax.plot(scales, noise_data['measure_feat'], lw=lw['big'])
-    [h.set_color(c) for h, c in zip(handles, shaded['feat'])]
+    # Plot analysis results:
+    for ax, x in zip([alpha_ax, sigma_ax], [scales, noise_data['measure_inv']]):
+        # Plot pure-song analysis results:
+        handles = ax.plot(x, pure_data['measure_feat'], lw=lw['big'], ls='dotted')
+        [h.set_color(c) for h, c in zip(handles, shaded['feat'])]
 
-    # Add proxy legend:
-    h1 = big_ax.plot([], [], c='k', lw=lw['big'], label='$\\alpha\\cdot s(t) + \\eta(t)$')[0]
-    h2 = big_ax.plot([], [], c='k', lw=lw['big'], ls='dotted', label='$\\alpha\\cdot s(t)$')[0]
-    big_ax.legend(handles=[h1, h2], **leg_kwargs)
+        # Plot noise-song analysis results:
+        handles = ax.plot(x, noise_data['measure_feat'], lw=lw['big'])
+        [h.set_color(c) for h, c in zip(handles, shaded['feat'])]
+
+        # Add proxy legend:
+        if ax == alpha_ax:
+            h1 = ax.plot([], [], c='k', lw=lw['big'], label='$\\alpha\\cdot s(t) + \\eta(t)$')[0]
+            h2 = ax.plot([], [], c='k', lw=lw['big'], ls='dotted', label='$\\alpha\\cdot s(t)$')[0]
+            ax.legend(handles=[h1, h2], **leg_kwargs)
 
     if save_path is not None:
         fig.savefig(save_path)
diff --git a/python/fig_temp.py b/python/fig_temp.py
new file mode 100644
index 0000000..d93268a
--- /dev/null
+++ b/python/fig_temp.py
@@ -0,0 +1,131 @@
+import glob
+import numpy as np
+import matplotlib.pyplot as plt
+from thunderhopper.modeltools import load_data
+from thunderhopper.filtertools import find_kern_specs
+from misc_functions import unsort_unique
+from color_functions import sample_cmap
+from IPython import embed
+
+# Settings:
+target = 'Omocestus_rufipes'
+data_path = glob.glob(f'../data/inv/full/{target}*.npz')[0]
+stages = ['conv', 'feat']
+load_kwargs = dict(
+    files=stages,
+    keywords=['scales', 'measure']
+)
+
+# Subset settings:
+all_types = np.array([1, -1, 2, -2, 3, -3, 4, -4, 5, -5,
+                      6, -6, 7, -7, 8, -8, 9, -9, 10, -10]).astype(float)
+all_sigmas = np.array([0.001, 0.002, 0.004, 0.008, 0.016, 0.032])
+kerns = None
+types = np.array([1, -1, 2, -2, 3, -3, 4, -4, 5, -5,
+                  6, -6, 7, -7, 8, -8, 9, -9, 10, -10])
+sigmas = np.array([0.008, 0.016])
+
+# Plot settings:
+line_kwargs = dict(
+    linewidth=2,
+)
+median_kwargs = dict(
+    linewidth=4,
+    c='k',
+    linestyle='--'
+)
+
+# Load invariance data:
+data, config = load_data(data_path, **load_kwargs)
+scales = data['scales']
+
+# Reduce to kernel subset:
+if any(var is not None for var in [kerns, types, sigmas]):
+    subset_inds = find_kern_specs(config['k_specs'], kerns, types, sigmas)
+    data['measure_conv'] = data['measure_conv'][:, subset_inds]
+    data['measure_feat'] = data['measure_feat'][:, subset_inds]
+    config['kernels'] = config['kernels'][:, subset_inds]
+    config['k_specs'] = config['k_specs'][subset_inds, :]
+kern_types = unsort_unique(config['k_specs'][:, 0])
+kern_sigmas = unsort_unique(config['k_specs'][:, 1])
+
+# Prepare colors:
+type_colors = {t: c for t, c in zip(all_types, sample_cmap('turbo', all_types.size))}
+sigma_colors = {s: c for s, c in zip(all_sigmas, sample_cmap('turbo', all_sigmas.size))}
+
+# Prepare graph:
+fig, axes = plt.subplots(2, 4, figsize=(16, 16), layout='constrained', sharex=True)
+axes[0, 0].set_xlim(scales[0], scales[-1])
+axes[0, 0].set_xscale('log')
+axes[0, 0].set_ylabel('conv')
+axes[1, 0].set_ylabel('feat')
+
+# Condense across kernels:
+median_conv = np.median(data['measure_conv'], axis=1)
+median_feat = np.median(data['measure_feat'], axis=1)
+
+# Coded by type:
+leg_handles = []
+for kern_type in kern_types:
+    color = type_colors[kern_type]
+    inds = find_kern_specs(config['k_specs'], types=kern_type)
+    leg_handles.append(axes[0, 0].plot(scales, data['measure_conv'][:, inds],
+                       c=color, label=f'{kern_type}', **line_kwargs)[0])
+    axes[0, 0].plot(scales, median_conv, **median_kwargs)
+    axes[1, 0].plot(scales, data['measure_feat'][:, inds], c=color, **line_kwargs)
+    axes[1, 0].plot(scales, median_feat, **median_kwargs)
+axes[0, 0].legend(handles=leg_handles, loc='upper left')
+axes[0, 0].set_title('Coded by type')
+
+# Coded by sigma:
+leg_handles = []
+for kern_sigma in kern_sigmas:
+    color = sigma_colors[kern_sigma]
+    inds = find_kern_specs(config['k_specs'], sigmas=kern_sigma)
+    leg_handles.append(axes[0, 1].plot(scales, data['measure_conv'][:, inds],
+                       c=color, label=f'{kern_sigma}', **line_kwargs)[0])
+    axes[0, 1].plot(scales, median_conv, **median_kwargs)
+    axes[1, 1].plot(scales, data['measure_feat'][:, inds], c=color, **line_kwargs)
+    axes[1, 1].plot(scales, median_feat, **median_kwargs)
+axes[0, 1].legend(handles=leg_handles, loc='upper left')
+axes[0, 1].set_title('Coded by sigma')
+
+# Normalize measures:
+data['measure_conv'] -= data['measure_conv'].min(axis=0)
+data['measure_conv'] /= data['measure_conv'].max(axis=0)
+data['measure_feat'] -= data['measure_feat'].min(axis=0)
+data['measure_feat'] /= data['measure_feat'].max(axis=0)
+
+# Condense across kernels:
+median_conv = np.median(data['measure_conv'], axis=1)
+median_feat = np.median(data['measure_feat'], axis=1)
+
+# Coded by type:
+leg_handles = []
+for kern_type in kern_types:
+    color = type_colors[kern_type]
+    inds = find_kern_specs(config['k_specs'], types=kern_type)
+    leg_handles.append(axes[0, 2].plot(scales, data['measure_conv'][:, inds],
+                       c=color, label=f'{kern_type}', **line_kwargs)[0])
+    axes[0, 2].plot(scales, median_conv, **median_kwargs)
+    axes[1, 2].plot(scales, data['measure_feat'][:, inds], c=color, **line_kwargs)
+    axes[1, 2].plot(scales, median_feat, **median_kwargs)
+axes[0, 2].legend(handles=leg_handles, loc='upper left')
+axes[0, 2].set_title('Coded by type')
+
+
+# Coded by sigma:
+leg_handles = []
+for kern_sigma in kern_sigmas:
+    color = sigma_colors[kern_sigma]
+    inds = find_kern_specs(config['k_specs'], sigmas=kern_sigma)
+    leg_handles.append(axes[0, 3].plot(scales, data['measure_conv'][:, inds],
+                       c=color, label=f'{kern_sigma}', **line_kwargs)[0])
+    axes[0, 3].plot(scales, median_conv, **median_kwargs)
+    axes[1, 3].plot(scales, data['measure_feat'][:, inds], c=color, **line_kwargs)
+    axes[1, 3].plot(scales, median_feat, **median_kwargs)
+axes[0, 3].legend(handles=leg_handles, loc='upper left')
+axes[0, 3].set_title('Coded by sigma')
+
+plt.show()
+embed()
diff --git a/python/misc_functions.py b/python/misc_functions.py
index efd0ed1..ee9aedb 100644
--- a/python/misc_functions.py
+++ b/python/misc_functions.py
@@ -1,9 +1,20 @@
 import numpy as np
+from scipy.stats import gaussian_kde
 
 def shorten_species(name):
     genus, species = name.split('_')
     return genus[0] + '. ' + species
 
+def unsort_unique(array):
+    values, inds = np.unique(array, return_index=True)
+    return values[np.argsort(inds)]
+
+def get_kde(data, sigma, axis=None, n=1000, pad=10):
+    if axis is None:
+        axis = np.linspace(data.min() - pad * sigma, data.max() + pad * sigma, n)
+    pdf = gaussian_kde(data, sigma)(axis)
+    return pdf, axis
+
 def get_saturation(sigmoid, low=0.05, high=0.95, first=True, last=True,
                    condense=None):
     if condense == 'norm' and sigmoid.ndim == 2:
@@ -16,17 +27,17 @@ def get_saturation(sigmoid, low=0.05, high=0.95, first=True, last=True,
     low_value = min_value + low * span
     high_value = min_value + high * span
 
-    low_mask = sigmoid >= low_value
-    high_mask = sigmoid >= high_value
+    low_mask = sigmoid <= low_value
+    high_mask = sigmoid <= high_value
     if sigmoid.ndim == 1:
-        low_ind = np.nonzero(low_mask)[0][0]
-        high_ind = np.nonzero(high_mask)[0][0]
+        low_ind = np.nonzero(low_mask)[0][-1]
+        high_ind = np.nonzero(high_mask)[0][-1]
     elif condense == 'all':
-        low_ind = np.nonzero(low_mask.all(axis=1))[0][0]
-        high_ind = np.nonzero(high_mask.all(axis=1))[0][0]
+        low_ind = np.nonzero(low_mask.all(axis=1))[0][-1]
+        high_ind = np.nonzero(high_mask.all(axis=1))[0][-1]
     else:
         low_ind, high_ind = [], []
         for i in range(sigmoid.shape[1]):
-            low_ind.append(np.nonzero(low_mask[:, i])[0][0])
-            high_ind.append(np.nonzero(high_mask[:, i])[0][0])
-    return low_ind, high_ind
\ No newline at end of file
+            low_ind.append(np.nonzero(low_mask[:, i])[0][-1])
+            high_ind.append(np.nonzero(high_mask[:, i])[0][-1])
+    return low_ind, high_ind
diff --git a/python/save_inv_data_full.py b/python/save_inv_data_full.py
index 6efcc42..a1800e5 100644
--- a/python/save_inv_data_full.py
+++ b/python/save_inv_data_full.py
@@ -4,19 +4,23 @@ import matplotlib.pyplot as plt
 from thunderhopper.modeltools import load_data, save_data
 from thunderhopper.filetools import crop_paths
 from thunderhopper.filtertools import find_kern_specs
-from thunderhopper.model import process_signal
+from thunderhopper.model import process_signal, convolve_kernels
 from IPython import embed
 
 # GENERAL SETTINGS:
 target = 'Omocestus_rufipes'
 data_paths = glob.glob(f'../data/processed/{target}*.npz')
-stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat']
+noise_path = '../data/processed/white_noise_sd-1.npz'
+stages = ['filt', 'env', 'log', 'inv', 'conv', 'feat']
 save_path = '../data/inv/full/'
 
 # ANALYSIS SETTINGS:
-example_scales = np.array([0, 1, 10, 50])
-scales = np.geomspace(0.01, 100, 100)
+example_scales = np.array([0.1, 1, 10, 30, 100, 300])
+scales = np.geomspace(0.01, 10000, 100)
 scales = np.unique(np.concatenate((scales, example_scales)))
+thresh_rel = 3
+
+# SUBSET SETTINGS:
 kernels = np.array([
     [1, 0.002],
     [-1, 0.002],
@@ -29,20 +33,29 @@ kernels = None
 types = None#np.array([-1])
 sigmas = None#np.array([0.001, 0.002, 0.004, 0.008, 0.016, 0.032])
 
+# PREPARATION:
+noise_data = np.load(noise_path)
+pure_noise = noise_data['raw']
+
 # EXECUTION:
 for data_path, name in zip(data_paths, crop_paths(data_paths)):
     print(f'Processing {name}')
 
-    # Get song recording:
+    # Get song recording (prior to anything):
     data, config = load_data(data_path, files='raw')
     song, rate = data['raw'], config['rate']
 
+    if thresh_rel is not None:
+        # Get noise-bound kernel-specific thresholds:
+        config['feat_thresh'] = noise_data['conv'].std(axis=0) * thresh_rel
+
     # Reduce to kernel subset:
-    kern_inds = find_kern_specs(config['k_specs'], kernels, types, sigmas)
-    config['kernels'] = config['kernels'][:, kern_inds]
-    config['k_specs'] = config['k_specs'][kern_inds, :]
-    config['k_props'] = [config['k_props'][i] for i in kern_inds]
-    config['feat_thresh'] = config['feat_thresh'][kern_inds]
+    if any(var is not None for var in [kernels, types, sigmas]):
+        kern_inds = find_kern_specs(config['k_specs'], kernels, types, sigmas)
+        config['kernels'] = config['kernels'][:, kern_inds]
+        config['k_specs'] = config['k_specs'][kern_inds, :]
+        config['k_props'] = [config['k_props'][i] for i in kern_inds]
+        config['feat_thresh'] = config['feat_thresh'][kern_inds]
 
     # Get song segment to be analyzed:
     time = np.arange(song.shape[0]) / rate
@@ -52,22 +65,19 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
     # Normalize song component:
     song /= song[segment].std(axis=0)
 
-    # Get normalized noise:
-    rng = np.random.default_rng()
-    noise = rng.normal(size=song.shape[0])
+    # Get normalized noise component:
+    noise = pure_noise[:song.shape[0]]
     noise /= noise[segment].std()
 
     # Prepare snippet storage:
     shape_low = (song.shape[0], example_scales.size)
     shape_high = (song.shape[0], config['k_specs'].shape[0], example_scales.size)
     snippets = dict(
-        snip_raw=np.zeros(shape_low, dtype=float),
         snip_filt=np.zeros(shape_low, dtype=float),
         snip_env=np.zeros(shape_low, dtype=float),
         snip_log=np.zeros(shape_low, dtype=float),
         snip_inv=np.zeros(shape_low, dtype=float),
         snip_conv=np.zeros(shape_high, dtype=float),
-        snip_bi=np.zeros(shape_high, dtype=float),
         snip_feat=np.zeros(shape_high, dtype=float)
     )
 
@@ -75,7 +85,6 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
     shape_low = (scales.size,)
     shape_high = (scales.size, config['k_specs'].shape[0])
     measures = dict(
-        measure_raw=np.zeros(shape_low, dtype=float),
         measure_filt=np.zeros(shape_low, dtype=float),
         measure_env=np.zeros(shape_low, dtype=float),
         measure_log=np.zeros(shape_low, dtype=float),
@@ -96,36 +105,18 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
                                         signal=scaled, rate=rate)
         # Store results:
         for stage in stages:
-            key = f'measure_{stage}'
+            mkey, skey = f'measure_{stage}', f'snip_{stage}'
 
             # Log snippet data:
             if scale in example_scales:
                 scale_ind = np.nonzero(example_scales == scale)[0][0]
-                snippets[f'snip_{stage}'][:, ..., scale_ind] = signals[stage]
+                snippets[skey][:, ..., scale_ind] = signals[stage]
 
             # Log intensity measure per stage (excluding binary):
             if stage in ['raw', 'filt', 'env', 'log', 'inv', 'conv']:
-                measures[key][i] = signals[stage][segment, ...].std(axis=0)
+                measures[mkey][i] = signals[stage][segment, ...].std(axis=0)
             elif stage == 'feat':
-                measures[key][i] = signals[stage][segment, :].mean(axis=0)
-    
-    # thresh_y = np.percentile(measures['measure_feat'], 99, axis=0)
-    # kern_types = np.unique()
-    # thresh_x = np.zeros(thresh_y.shape, dtype=float)
-    # for i, thresh in enumerate(thresh_y):
-    #     if thresh < 0.1:
-    #         thresh_x[i] = scales[-1]
-    #         continue
-    #     mask = (measures['measure_feat'][:, i] < thresh)
-    #     thresh_x[i] = scales[np.nonzero(mask)[0][-1]]
-    # inds = np.argsort(thresh_x)
-    # print(config['k_specs'][inds, :])
-
-    # fig, axes = plt.subplots(1, 2)
-    # axes[0].plot(snippets['snip_feat'][:, inds, -1])
-    # axes[1].plot(scales, measures['measure_feat'][:, inds])
-    # plt.show()
-    # embed()
+                measures[mkey][i] = signals[stage][segment, :].mean(axis=0)
 
     # Save analysis results:
     if save_path is not None:
diff --git a/python/save_inv_data_log-hp.py b/python/save_inv_data_log-hp.py
index 78689d0..68e3f4c 100644
--- a/python/save_inv_data_log-hp.py
+++ b/python/save_inv_data_log-hp.py
@@ -7,15 +7,19 @@ from IPython import embed
 # GENERAL SETTINGS:
 target = ['Omocestus_rufipes', '*'][0]
 data_paths = search_files(target, excl='noise', dir='../data/processed/')
+noise_path = '../data/processed/white_noise_sd-1.npz'
 save_path = '../data/inv/log_hp/'
 
 # ANALYSIS SETTINGS:
-add_noise = False
+add_noise = target == '*' or False
 save_snippets = target == 'Omocestus_rufipes'
 example_scales = np.array([0.1, 1, 10, 30, 100, 300])
 scales = np.geomspace(0.1, 10000, 500)
 scales = np.unique(np.concatenate((scales, example_scales)))
 
+# PREPARATION:
+pure_noise = np.load(noise_path)['filt']
+
 # EXECUTION:
 for data_path, name in zip(data_paths, crop_paths(data_paths)):
     print(f'Processing {name}')
@@ -36,9 +40,8 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
     mix = song[:, None] * scales[None, :]
 
     if add_noise:
-        # Add normalized envelopenoise:
-        rng = np.random.default_rng()
-        noise = rng.normal(scale=1, size=song.shape)
+        # Add normalized noise component:
+        noise = pure_noise[:song.shape[0]]
         noise /= noise[segment].std()
         mix += noise[:, None]
 
diff --git a/python/save_inv_data_thresh-lp.py b/python/save_inv_data_thresh-lp.py
index 8d247bb..fbe9fe4 100644
--- a/python/save_inv_data_thresh-lp.py
+++ b/python/save_inv_data_thresh-lp.py
@@ -8,13 +8,13 @@ from thunderhopper.model import convolve_kernels
 from IPython import embed
 
 # GENERAL SETTINGS:
-target = ['Omocestus_rufipes', '*'][1]
+target = ['Omocestus_rufipes', '*'][0]
 data_paths = search_files(target, excl='noise', dir='../data/processed/')
 noise_path = '../data/processed/white_noise_sd-1.npz'
 save_path = '../data/inv/thresh_lp/'
 
 # ANALYSIS SETTINGS:
-add_noise = True
+add_noise = False
 save_snippets = add_noise and (target == 'Omocestus_rufipes')
 plot_results = False
 example_scales = np.array([0, 1, 10, 30, 100])
@@ -62,8 +62,8 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
     thresh_abs = ref_conv[segment, :].std(axis=0, keepdims=True) * thresh_rel[:, None]
 
     # Prepare measure storage:
-    shape = (scales.size, kern_specs.shape[0], thresh_rel.size)
-    measure_feat = np.zeros(shape, dtype=float)
+    measure_inv = np.zeros((scales.size,), dtype=float)
+    measure_feat = np.zeros((scales.size, kern_specs.shape[0], thresh_rel.size), dtype=float)
     if save_snippets:
         # Prepare snippet storage:
         snip_inv = np.zeros((song.size, example_scales.size), dtype=float)
@@ -81,6 +81,9 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
         if add_noise:
             # Add noise:
             scaled_song += noise
+        
+        # Log input intensity measure:
+        measure_inv[i] = scaled_song[segment].std()
 
         # Process mixture:
         scaled_conv = convolve_kernels(scaled_song, config['kernels'], config['k_specs'])
@@ -130,6 +133,7 @@ for data_path, name in zip(data_paths, crop_paths(data_paths)):
         data = dict(
             scales=scales,
             example_scales=example_scales,
+            measure_inv=measure_inv,
             measure_feat=measure_feat,
             thresh_rel=thresh_rel,
             thresh_abs=thresh_abs,
diff --git a/python/save_noise_data.py b/python/save_noise_data.py
index 9c4fd99..fc3ed92 100644
--- a/python/save_noise_data.py
+++ b/python/save_noise_data.py
@@ -7,7 +7,7 @@ from IPython import embed
 
 # General:
 save_path = '../data/processed/white_noise'
-stages = ['filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat']
+stages = ['raw', 'filt', 'env', 'log', 'inv', 'conv', 'bi', 'feat']
 sds = [1]
 dur = 60
 
@@ -23,9 +23,9 @@ types = [1, -1, 2, -2, 3, -3, 4, -4, 5, -5,
          6, -6, 7, -7, 8, -8, 9, -9, 10, -10]
 config = configuration(env_rate, feat_rate, types=types, sigmas=sigmas)
 config.update({
-    'bp_fcut': None,
     'rate_ratio': None,
     'env_fcut': 250,
+    'db_ref': 1,
     'inv_fcut': 5,
     'feat_thresh': np.load('../data/kernel_thresholds.npy') * 0.2,
     'feat_fcut': 0.5,
diff --git a/python/save_ref_measures.py b/python/save_ref_measures.py
new file mode 100644
index 0000000..8d61843
--- /dev/null
+++ b/python/save_ref_measures.py
@@ -0,0 +1,65 @@
+import numpy as np
+from thunderhopper.filters import decibel, sosfilter
+from thunderhopper.model import convolve_kernels, process_signal
+from thunderhopper.modeltools import load_data
+from IPython import embed
+
+## SETTINGS:
+
+# General:
+mode = ['log_hp', 'thresh_lp', 'full'][2]
+noise_path = '../data/processed/white_noise_sd-1.npz'
+save_path = '../data/inv/'
+pad = np.array([0.1, 0.9])
+
+stages = dict(
+    log_hp=['filt', 'env', 'log', 'inv'],
+    thresh_lp=['inv', 'conv', 'feat'],
+    full=['raw', 'filt', 'env', 'log', 'inv', 'conv', 'feat']
+)[mode]
+
+# PROCESSING:
+
+print(f'Fetching references for {mode} invariance...')
+
+# Load pure-noise starter representation:
+noise_data, config = load_data(noise_path, stages[0])
+starter = noise_data[stages[0]]
+
+# Prepare buffered measurement segment:
+pad = (pad * starter.shape[0]).astype(int)
+segment = np.arange(starter.shape[0])[pad[0]:pad[1]]
+
+# Normalize starter:
+starter /= starter[segment].std()
+
+# Run pipeline:
+if mode == 'log_hp':
+    data = {'filt': starter}
+    data['env'] = sosfilter(np.abs(data['filt']), config['rate'], config['env_fcut'],
+                            'lp', padtype='even', padlen=config['padlen'])
+    data['log'] = decibel(data['env'], ref=1)
+    data['inv'] = sosfilter(data['log'], config['env_rate'], config['inv_fcut'],
+                            'hp', padtype='constant', padlen=config['padlen'])
+elif mode == 'thresh_lp':
+    data = {'inv': starter}
+    data['conv'] = convolve_kernels(data['inv'], config['kernels'], config['k_specs'])
+    data['feat'] = sosfilter((data['conv'] > config['feat_thresh']).astype(float),
+                             config['env_rate'], config['feat_fcut'], 'lp',
+                             padtype='fixed', padlen=config['padlen'])
+elif mode == 'full':
+    data = process_signal(config, stages, signal=starter, rate=config['rate'])[0]
+
+# Get measures:
+measures = {}
+for stage in stages:
+    if stage == 'feat':
+        measures[stage] = data[stage][segment, :].mean(axis=0)
+    else:
+        measures[stage] = data[stage][segment, ...].std(axis=0)
+
+# Save results:
+np.savez(save_path + f'{mode}/ref_measures.npz', **measures)
+
+print('Done.')
+embed()
diff --git a/python/save_snippet_data.py b/python/save_snippet_data.py
index e52b759..0869123 100644
--- a/python/save_snippet_data.py
+++ b/python/save_snippet_data.py
@@ -16,7 +16,7 @@ if False:
 
 # Interactivity:
 reload_saved = False
-gui = False
+gui = True
 
 # Processing:
 env_rate = 44100.0
@@ -29,6 +29,7 @@ config.update({
     'channel': 0,
     'rate_ratio': None,
     'env_fcut': 250,
+    'db_ref': 1,
     'inv_fcut': 5,
     'feat_thresh': np.load('../data/kernel_thresholds.npy') * 0.2,
     'feat_fcut': 0.5,