Full Day Systematics Inspection

by Tyler Cox and Josh Dillon, last updated July 20, 2024

This notebooks is designed to inspect the full day delay spectra for a small number of redundantly-averaged baselines, by inpainting in frequency and time, and performing delay filtering, cross-talk filtering, and main-beam fringe rate filtering. It is based on a similar notebook written by Josh Dillon for H1C analysis.

More specifically, it looks at 1-unit (i.e. 14.6 m) EW baselines, 2-unit EW baselines, 4-unit EW baselines, and one inter-sector baseline. This is done for both polarizations (and eventually for Pseudo-Stokes I and Q). We look at delay-spectrum waterfalls vs. LST, delay-spectrum waterfalls vs. fringe-rate, and averaged delay-spectra for these baselines. Those spectra are plotted for two wide bands, one below FM and one above FM.

This notebook is divided into four sections:

• Section 1: Time + Frequency Inpainting of Redundantly Averaged Calibrated Visibilities

• Section 2: Delay-Filter and Cross-Talk Filter Redundantly-Averaged Visibilities

• Section 3: Form Psuedo-Stokes I + Q

• Section 4: Main-Beam Fringe-Rate Filter

import os
os.environ['HDF5_USE_FILE_LOCKING'] = 'FALSE'
import h5py
import hdf5plugin  # REQUIRED to have the compression plugins available
import glob
import uvtools
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from copy import deepcopy
from pyuvdata import uvbeam, uvdata
from hera_filters import dspec
from hera_cal import io, frf, vis_clean, redcal, utils
from hera_cal.datacontainer import DataContainer, RedDataContainer
from IPython.display import display, HTML

%matplotlib inline
display(HTML("<style>.container { width:100% !important; }</style>"))
_ = np.seterr(all='ignore')  # get rid of red warnings
%config InlineBackend.figure_format = 'retina'

# Get file names
SUM_FILE = os.environ.get("SUM_FILE", None)
#SUM_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/2459866/zen.2459869.25282.sum.uvh5'
RED_AVG_SUFFIX = os.environ.get("RED_AVG_SUFFIX", "sum.smooth_calibrated.red_avg.uvh5")

# Arguements for delay-filtering, fringe-rate filtering, cross-talk subtraction, and inpainting
XTALK_FR = float(os.environ.get("XTALK_FR", 0.025)) # Fringe rate half-width used for fringe rate filtering
FILT_MIN_DLY = float(os.environ.get("FILT_MIN_DLY", 150.)) # minimum delay for delay-filter in units of ns
INPAINT_MIN_DLY = float(os.environ.get("INPAINT_MIN_DLY", 500.)) # minimum delay for inpainting in units of ns
INPAINT_REGULARIZATION = float(os.environ.get("INPAINT_REGULARIZATION", 1e-5)) # reasonable values are between 1e-2 and 1e-5
STANDOFF = float(os.environ.get("STANDOFF", 50.)) # additional standoff added to min_dly for delay-filter in units of ns
INPAINT_FR = float(os.environ.get("INPAINT_FR", 2.5)) # Fringe-rate half-width in mHz for time inpainting
EIGENVAL_CUTOFF = float(os.environ.get("EIGENVAL_CUTOFF", 1e-12))

# Flag settings to get delay transform regions
FM_LOW_FREQ = float(os.environ.get("FM_LOW_FREQ", 87.5)) # in MHz
FM_HIGH_FREQ = float(os.environ.get("FM_HIGH_FREQ", 108.0)) # in MHz

# Maximum number of contiguous frequency channels in a time integration to attempt to inpaint over
# If number of contiguous flags is greater than this value, flagged channels will remain flagged
# and will attempt to be inpainted over along the time axis
MAX_CONTIGUOUS_FLAGS = int(os.environ.get("MAX_CONTIGUOUS_FLAGS", 20))

# Remove this number of channels on either side of the below/above FM bands to exclude from delay spectrum plots
SPECTRUM_CHAN_BUFFER = int(os.environ.get("SPECTRUM_CHAN_BUFFER", 25))

# Is the input LST-binned
IS_LST_BINNED = os.environ.get("IS_LST_BINNED", "FALSE").upper() == "TRUE"

for setting in ['XTALK_FR', 'FILT_MIN_DLY', 'INPAINT_MIN_DLY', 'STANDOFF', 'INPAINT_FR', 'EIGENVAL_CUTOFF', 
                'FM_LOW_FREQ', 'FM_HIGH_FREQ', 'MAX_CONTIGUOUS_FLAGS', 'SPECTRUM_CHAN_BUFFER', 'IS_LST_BINNED']:
        print(f'{setting} = {eval(setting)}')

XTALK_FR = 0.025
FILT_MIN_DLY = 150.0
INPAINT_MIN_DLY = 500.0
STANDOFF = 50.0
INPAINT_FR = 2.5
EIGENVAL_CUTOFF = 1e-12
FM_LOW_FREQ = 87.5
FM_HIGH_FREQ = 108.0
MAX_CONTIGUOUS_FLAGS = 20
SPECTRUM_CHAN_BUFFER = 25
IS_LST_BINNED = False

sum_glob = os.path.join(os.path.dirname(SUM_FILE), "*" + RED_AVG_SUFFIX)
red_avg_files = sorted(glob.glob(sum_glob))
print(f'Found {len(red_avg_files)} *.{RED_AVG_SUFFIX} files starting with {red_avg_files[0]}.')

Found 174 *.sum.smooth_calibrated.red_avg.uvh5 files starting with /mnt/sn1/data2/2460846/zen.2460846.21085.sum.smooth_calibrated.red_avg.uvh5.

# Select baselines for plotting
for i in range(len(red_avg_files) // 2):
    hd = io.HERADataFastReader(red_avg_files[(len(red_avg_files) // 2) + i])
    bls_to_plot = []

    # Find 1, 2, and 4 units EW auto-pol baselines 
    for bl in hd.bls:
        if utils.split_pol(bl[2])[0] == utils.split_pol(bl[2])[1]:
            bl_vec = (hd.antpos[bl[1]] - hd.antpos[bl[0]])
            if (np.abs(bl_vec[1]) < 1) and int(np.round(bl_vec[0] / 14.6)) in [1, 2, 4]:
                bls_to_plot.append(bl)

    # Get intersector baseline
    _, _, nsamples = hd.read()
    nsamples = RedDataContainer(nsamples, antpos=hd.antpos)
    # try a bunch of possible intersector baselines, hoping one will work
    for ant2 in list(range(226, 233)) + list(range(207, 215)) + list(range(187, 196)):
        if ((144, ant2, 'ee') in nsamples) and ((144, ant2, 'nn') in nsamples):
            if (np.median(nsamples[(144, ant2, 'nn')]) > 0) & (np.median(nsamples[(144, ant2, 'ee')]) > 1):
                bls_to_plot.append(nsamples.get_ubl_key((144, ant2, 'nn')))
                bls_to_plot.append(nsamples.get_ubl_key((144, ant2, 'ee')))
                break
    if len(bls_to_plot) >= 8:
        break

# Identify indices
FM_LOW_FREQ_IND = np.argmin(np.abs(hd.freqs - FM_LOW_FREQ * 1e6))
FM_HIGH_FREQ_IND = np.argmin(np.abs(hd.freqs - FM_HIGH_FREQ * 1e6))

hd = io.HERADataFastReader(red_avg_files)
d, f, n = hd.read(pols=['nn', 'ee'], bls=bls_to_plot)

# Create a truncated version of the data to deal with edge times that are flagged
data = deepcopy(d)
flags = deepcopy(f)
nsamples = deepcopy(n)

# Identify unflagged edge channels
unflagged_chans = np.where(~np.all(flags[bls_to_plot[0]], axis=1))[0]
inds = np.arange(np.min(unflagged_chans), np.max(unflagged_chans) + 1)

# Fix data
for bl in data:
    data[bl] = d[bl][inds]
    flags[bl] = f[bl][inds]
    nsamples[bl] = n[bl][inds]
    
# Fix metadata
data.times = d.times[inds]
data.lsts = d.lsts[inds]

Section 1: Time + Frequency Inpainting of Redundantly Averaged Calibrated Visibilities

unit = int(np.round(np.linalg.norm(data.antpos[bls_to_plot[0][1]] - data.antpos[bls_to_plot[0][0]]) / 14.6, 1))
display(HTML(f'<h2>Figure 1: Amplitude and Phase of Calibrated {unit}-unit EW Redundantly-Averaged Baseline</h2>'))

# Get the LSTs from the data
lsts_unwrapped = np.where(data.lsts > data.lsts[-1], data.lsts - 2 * np.pi, data.lsts) * 12 / np.pi  # in hours
extent = [data.freqs[0] / 1e6, data.freqs[-1] / 1e6, lsts_unwrapped[-1], lsts_unwrapped[0]]

if not IS_LST_BINNED:
    fig, axs = plt.subplots(2, 2, figsize=(18, 14), sharex=True, sharey=True, gridspec_kw={'hspace': .03, 'wspace': .03})
else:
    fig, axs = plt.subplots(2, 3, figsize=(18, 14), sharex=True, sharey=True, gridspec_kw={'hspace': .03, 'wspace': .03})
    
fig.subplots_adjust(hspace=0.1, wspace=0.1)

# Plot the amplitude of ee and nn
im1 = axs[0, 0].imshow(np.where(flags[bls_to_plot[0]], np.nan, np.abs(data[bls_to_plot[0]])), aspect='auto', interpolation='none',
                   norm=matplotlib.colors.LogNorm(), cmap='inferno', extent=extent)
axs[1, 0].imshow(np.where(flags[bls_to_plot[1]], np.nan, np.abs(data[bls_to_plot[1]])), aspect='auto', interpolation='none',
                   norm=matplotlib.colors.LogNorm(vmin=im1.get_clim()[0], vmax=im1.get_clim()[1]), cmap='inferno', extent=extent)

# Plot the phase of ee and nn
im2 = axs[0, 1].imshow(
    np.where(flags[bls_to_plot[0]], np.nan, np.angle(data[bls_to_plot[0]])), 
    aspect='auto', cmap='twilight', interpolation='None', extent=extent,
    vmin=-np.pi, vmax=np.pi
)
axs[1, 1].imshow(
    np.where(flags[bls_to_plot[1]], np.nan, np.angle(data[bls_to_plot[1]])), 
    aspect='auto', cmap='twilight', interpolation='None', extent=extent,
    vmin=-np.pi, vmax=np.pi
)

fig.colorbar(im1, ax=axs[:, 0], label='$|V_{ij}|$ (Jy)', location='top', pad=0.01)
fig.colorbar(im2, ax=axs[:, 1], label='Angle($V_{ij})$ (radians)', location='top', pad=0.01)

if IS_LST_BINNED:
    # Plot the amplitude of ee and nn
    im3 = axs[0, 2].imshow(np.where(flags[bls_to_plot[0]], np.nan, nsamples[bls_to_plot[0]]), aspect='auto', interpolation='none',
                       cmap='viridis', extent=extent)
    axs[1, 2].imshow(np.where(flags[bls_to_plot[1]], np.nan, nsamples[bls_to_plot[1]]), aspect='auto', interpolation='none',
                       vmin=im3.get_clim()[0], vmax=im3.get_clim()[1], cmap='viridis', extent=extent)
    fig.colorbar(im3, ax=axs[:, 2], label='Number of Samples', location='top', pad=0.01)
    
    bls_for_labeling = [bls_to_plot[0][2]] * 3 + [bls_to_plot[1][2]] * 3
else:
    bls_for_labeling = [bls_to_plot[0][2]] * 2 + [bls_to_plot[1][2]] * 2

for ax, pol in zip(axs.ravel(), bls_for_labeling):
    ax.text(0.03, 0.97, f"{unit}-unit EW ({pol})", 
             horizontalalignment='left', verticalalignment='top', 
             color='black', transform=ax.transAxes, fontsize=12, 
             bbox=dict(facecolor='w', alpha=.7, boxstyle='round'))
    
for i in range(axs.shape[1]):
    axs[1, i].set_xlabel('Frequency (MHz)')

for i in range(axs.shape[0]):
    axs[i, 0].set_ylabel('LST (hours)')

plt.tight_layout()

Figure 1: Amplitude and Phase of Calibrated 1-unit EW Redundantly-Averaged Baseline

This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.

def get_flagged_regions(flags, max_flagged_values=25):
    new_flags = np.zeros_like(flags)
    
    for i in range(flags.shape[0]):
        if np.all(flags[i]):
            new_flags[i] = True
        elif np.all(flags[i, low_band]) ^ np.all(flags[i, high_band]):
            new_flags[i, :] = True
            new_flags[i, high_band] = np.all(flags[i, high_band])
            new_flags[i, low_band] = np.all(flags[i, low_band])
        else:
            z = np.concatenate(([False], flags[i], [False]))
            starts = np.flatnonzero(~z[:-1] & z[1:])   
            ends = np.flatnonzero(z[:-1] & ~z[1:])
            for start, end in zip(starts, ends):
                if np.abs(end - start) > max_flagged_values:
                    new_flags[i, start:end] = True
                
    return new_flags

def sym_log_norm(to_plot):
    '''Shortcut for symmetric sym log norm.'''
    return matplotlib.colors.SymLogNorm(10, vmin=-np.nanmax(np.abs(to_plot)), vmax=np.nanmax(np.abs(to_plot)))

# Copy data
data_inpaint = deepcopy(data)
data_delay_filt = deepcopy(data)

# Peform delay filter separately above and below FM
FM_ind = np.argmin(np.abs(hd.freqs - 100e6))
unflagged_chans = np.argwhere(~np.all([np.all(flags[bl], axis=0) for bl in bls_to_plot], axis=0)).squeeze()
low_band = slice(np.min(unflagged_chans), np.max(unflagged_chans[unflagged_chans < FM_ind]) + 1)
high_band = slice(np.min(unflagged_chans[unflagged_chans > FM_ind]), np.max(unflagged_chans) + 1)

delay_bounds = {}
cache = {}

for inpaint, min_dly in zip([True, False], [INPAINT_MIN_DLY, FILT_MIN_DLY]):
    for bl in data:
        if np.all(data[bl] == 0.0) and np.all(flags[bl]):
            # don't bother delay filtering all 0s.
            continue

        wgts = np.where(flags[bl], 0, 1)

        # calculate filter properties
        bl_vec = (hd.antpos[bl[1]] - hd.antpos[bl[0]])
        bl_len = np.linalg.norm(bl_vec[:2]) / vis_clean.constants.c.value
        filter_centers, filter_half_widths = vis_clean.gen_filter_properties(ax='freq', horizon=1, standoff=STANDOFF, 
                                                                             min_dly=min_dly, bl_len=bl_len)
        delay_bounds[bl[:2]] = (filter_centers[0] - filter_half_widths[0], filter_centers[0] + filter_half_widths[0])

        # run delay filter for each band individually
        d_mdl = np.zeros_like(data[bl])
        for band in [low_band, high_band]:
            d_mdl[:, band], _, info = dspec.fourier_filter(hd.freqs[band], data[bl][:, band], wgts=wgts[:, band], filter_centers=filter_centers, 
                                                           filter_half_widths=filter_half_widths, mode='dpss_solve', 
                                                           ridge_alpha=(INPAINT_REGULARIZATION if inpaint else 0), fit_intercept=inpaint,
                                                           eigenval_cutoff=[EIGENVAL_CUTOFF], max_contiguous_edge_flags=len(hd.freqs[band]), cache=cache)
        if not inpaint:
            data_delay_filt[bl] = np.where(flags[bl], 0, data_inpaint[bl] - d_mdl)
        elif inpaint:
            data_inpaint[bl] = np.where(flags[bl], d_mdl, data_inpaint[bl])
            _flags = get_flagged_regions(flags[bl], max_flagged_values=MAX_CONTIGUOUS_FLAGS)
            _wgts = np.where(_flags, 0, 1)

            # now inpaint in time
            d_mdl = np.zeros_like(data[bl])
            for band in [low_band, high_band]:
                d_mdl[:, band], _, info = dspec.fourier_filter(data.times * 24 * 60 * 60, data_inpaint[bl][:, band], wgts=_wgts[:, band], filter_centers=filter_centers, 
                                                               filter_half_widths=[INPAINT_FR / 1000.], mode='dpss_solve', filter_dims=0, ridge_alpha=INPAINT_REGULARIZATION,
                                                               eigenval_cutoff=[EIGENVAL_CUTOFF], max_contiguous_edge_flags=len(data.times), cache=cache)


            data_inpaint[bl] = np.where(_flags, d_mdl, data_inpaint[bl])

# Indices to use for plotting
low_rng = np.arange(low_band.start + SPECTRUM_CHAN_BUFFER, low_band.stop - SPECTRUM_CHAN_BUFFER, dtype=int)
high_rng = np.arange(high_band.start + SPECTRUM_CHAN_BUFFER, high_band.stop - SPECTRUM_CHAN_BUFFER, dtype=int)

def plot_delay_vs_lst(data, rng, min_dly=FILT_MIN_DLY, vmin=None, vmax=None, title=""):
    
    # Add search-able title to figure
    display(HTML(f'<h2>{title}</h2>'))
    
    # plot in delay space
    fig, axs = plt.subplots(2, 4, figsize=(14,14), sharex=True, sharey=True)

    # Organize baselines
    plot_bls = list(set([bl[:2] for bl in data.keys() if bl[0] != bl[1]]))
    plot_bls = sorted(plot_bls, key=lambda k: np.linalg.norm(data.antpos[k[1]] - data.antpos[k[0]]))
    plot_pols = sorted(list(set([bl[-1] for bl in data.keys()])))

    # Fringe-rates
    lsts_unwrapped = np.where(data.lsts > data.lsts[-1], data.lsts - 2 * np.pi, data.lsts) * 12 / np.pi  # in hours

    for bi, blkey in enumerate(plot_bls):
        for pi, pol in enumerate(plot_pols):
            bl = blkey + (pol, )
            dfft = uvtools.utils.FFT(data[bl][:, rng], axis=1, taper='bh7')
            delays = uvtools.utils.fourier_freqs(hd.freqs[rng]) * 1e9
            to_plot = np.real(dfft)

            # Scale vmin and vmax 
            if bi == 0 and pi == 0: 
                _to_plot = np.copy(to_plot)

            im = axs[pi, bi].imshow(to_plot, interpolation='none', aspect='auto', cmap='bwr', norm=sym_log_norm(_to_plot),
                       extent=[delays[0], delays[-1], lsts_unwrapped[-1], lsts_unwrapped[0]])
            axs[pi, bi].set_xlim([-1499, 1499])

            # Get baseline-lengths
            blvec = data.antpos[bl[0]] - data.antpos[bl[1]]
            bllen = np.linalg.norm(blvec)
            if bi < 3:
                bltype = f"{int(np.round(bllen / 14.6))}-unit EW" 
            else:
                bltype = f"{np.abs(blvec[0]):.0f} m {'E' if blvec[0] > 0 else 'W'}, {np.abs(blvec[1]):.0f} m {'N' if blvec[1] > 0 else 'S'}"

            axs[pi, bi].text(0.05, 0.98, f"{bltype} ({pol})\n{np.round(data.freqs[rng[0]] / 1e6, 1)}—{np.round(data.freqs[rng[-1]] / 1e6, 1)} MHz", 
                             horizontalalignment='left', verticalalignment='top', 
                             color='black', transform=axs[pi, bi].transAxes, fontsize=12, 
                             bbox=dict(facecolor='w', alpha=.7, boxstyle='round'))

            # Label axes
            if pi == 1:
                axs[pi, bi].set_xlabel('Delay (ns)')
            if bi == 0:
                axs[pi, bi].set_ylabel('LST [hr]')
                
            bl_vec = (hd.antpos[blkey[1]] - hd.antpos[blkey[0]])
            bl_len = np.linalg.norm(bl_vec[:2]) / vis_clean.constants.c.value
            filter_centers, filter_half_widths = vis_clean.gen_filter_properties(
                ax='freq', horizon=1, standoff=STANDOFF, min_dly=min_dly, bl_len=bl_len
            )

            axs[pi, bi].axvline(-filter_half_widths[0] * 1e9, color='k', ls='-.')
            axs[pi, bi].axvline(filter_half_widths[0] * 1e9, color='k', ls='-.')
    
    plt.tight_layout()

    # Colorbar
    fig.subplots_adjust(right=0.9)
    cbar_ax = fig.add_axes([0.92, 0.04, 0.02, 0.95])
    fig.colorbar(im, cax=cbar_ax, label='$|\widetilde{V}_{ij}|$ (Jy)', extend='both')
    return {'vmin': vmin, 'vmax': vmax}

invalid escape sequence '\w'
invalid escape sequence '\w'
invalid escape sequence '\w'

def plot_delay_fringe_rate(data, rng, vmin=None, vmax=None, delay_bounds=None, fringe_bounds=None, title=""):
    
    # Add search-able title to figure
    display(HTML(f'<h2>{title}</h2>'))
    
    # plot in delay space
    fig, axs = plt.subplots(2, 4, figsize=(14, 8), sharex=True, sharey=True)

    # Organize baselines
    plot_bls = list(set([bl[:2] for bl in data.keys() if bl[0] != bl[1]]))
    plot_bls = sorted(plot_bls, key=lambda k: np.linalg.norm(data.antpos[k[1]] - data.antpos[k[0]]))
    plot_pols = sorted(list(set([bl[-1] for bl in data.keys()])))

    # Fringe-rates
    frates = uvtools.utils.fourier_freqs(data.times * 24 * 60 * 60) * 1000

    for bi, blkey in enumerate(plot_bls):
        for pi, pol in enumerate(plot_pols):
            bl = blkey + (pol, )
            dfft = uvtools.utils.FFT(data[bl][:, rng], axis=1, taper='bh7')
            delays = uvtools.utils.fourier_freqs(hd.freqs[rng]) * 1e9
            dfft2 = uvtools.utils.FFT(dfft, axis=0, taper='bh7')

            # Scale vmin and vmax 
            if bi == 0 and pi == 0: 
                vmin = np.median(np.abs(dfft2)) if vmin is None else vmin
                vmax = np.max(np.abs(dfft2)) if vmax is None else vmax

            im = axs[pi, bi].imshow(np.abs(dfft2), interpolation='none', aspect='auto', cmap='turbo', 
                                    norm=matplotlib.colors.LogNorm(vmin=vmin, vmax=vmax),
                                    extent=[delays[0], delays[-1], frates[-1], frates[0]])
            axs[pi, bi].set_xlim([-799, 799])
            axs[pi, bi].set_ylim([-5, 5])

            # Get baseline-lengths
            blvec = data.antpos[bl[0]] - data.antpos[bl[1]]
            bllen = np.linalg.norm(blvec)
            if bi < 3:
                bltype = f"{int(np.round(bllen / 14.6))}-unit EW" 
            else:
                bltype = f"{np.abs(blvec[0]):.0f} m {'E' if blvec[0] > 0 else 'W'}, {np.abs(blvec[1]):.0f} m {'N' if blvec[1] > 0 else 'S'}"

            axs[pi, bi].text(0.04, 0.96, f"{bltype} ({pol})\n{np.round(data.freqs[rng[0]] / 1e6, 1)}—{np.round(data.freqs[rng[-1]] / 1e6, 1)} MHz", 
                             horizontalalignment='left', verticalalignment='top', 
                             color='black', transform=axs[pi, bi].transAxes, fontsize=12, 
                             bbox=dict(facecolor='w', alpha=.7, boxstyle='round'))

            # Label axes
            if pi == 1:
                axs[pi, bi].set_xlabel('Delay (ns)')
            if bi == 0:
                axs[pi, bi].set_ylabel('Fringe Rate [mHz]')
                
            if delay_bounds:
                axs[pi, bi].axvline(delay_bounds[bl[:2]][0] * 1e9, color='white', ls='-')
                axs[pi, bi].axvline(delay_bounds[bl[:2]][1] * 1e9, color='white', ls='-')
            if fringe_bounds:
                axs[pi, bi].axhline(fringe_bounds[bl[:2]][0], color='white', ls='-.')
                axs[pi, bi].axhline(fringe_bounds[bl[:2]][1], color='white', ls='-.')

    plt.tight_layout()

    # Colorbar
    fig.subplots_adjust(right=0.9)
    cbar_ax = fig.add_axes([0.92, 0.04, 0.02, 0.95])
    fig.colorbar(im, cax=cbar_ax, label='$|\widetilde{V}_{ij}|$ (Jy)', extend='both')
    
    return {'vmin': vmin, 'vmax': vmax}

invalid escape sequence '\w'
invalid escape sequence '\w'
invalid escape sequence '\w'

plot_delay_vs_lst(
    data_inpaint, low_rng, min_dly=INPAINT_MIN_DLY,
    title=f'Figure 2: Real Part of Redundantly-Averaged, DPSS Time + Frequency Inpainted, \
    Delay Transformed, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
);

Figure 2: Real Part of Redundantly-Averaged, DPSS Time + Frequency Inpainted, Delay Transformed, 50.0—84.2 MHz

plot_delay_vs_lst(
    data_inpaint, high_rng, min_dly=INPAINT_MIN_DLY,
    title=f'Figure 3: Real Part of Redundantly-Averaged, DPSS Time + Frequency Inpainted, \
    Delay Transformed, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
);

Figure 3: Real Part of Redundantly-Averaged, DPSS Time + Frequency Inpainted, Delay Transformed, 111.4—230.1 MHz

low_band_meta = plot_delay_fringe_rate(
    data_inpaint, low_rng, title=f'Figure 4: Redundantly-Averaged, DPSS Time + Frequency Inpainted,\
    Delay + Fringe Transformed Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 4: Redundantly-Averaged, DPSS Time + Frequency Inpainted, Delay + Fringe Transformed Visibilities, 50.0—84.2 MHz

high_band_meta = plot_delay_fringe_rate(
    data_inpaint, high_rng, title=f'Figure 5: Redundantly-Averaged, DPSS Time + Frequency Inpainted, \
    Delay + Fringe Transformed Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 5: Redundantly-Averaged, DPSS Time + Frequency Inpainted, Delay + Fringe Transformed Visibilities, 111.4—230.1 MHz

Section 2: Delay-Filter and Cross-Talk Filter Redundantly-Averaged Visibilities

data_xtalk_filt = deepcopy(data_delay_filt)

cache = {}
for bl in data_delay_filt:
    # don't bother delay filtering all 0s.
    if np.all(data_delay_filt[bl] == 0.0) and np.all(flags[bl]):
        continue

    # run delay filter for each band individually
    d_res = np.zeros_like(data_delay_filt[bl])

    # Get weights
    wgts = (~flags[bl]).astype(float)

    for band in [low_band, high_band]:
        _, d_res[:, band], info = dspec.fourier_filter(data.times * 24 * 60 * 60, data_delay_filt[bl][:, band], 
                                                       wgts=wgts[:, band], filter_centers=filter_centers, 
                                                       filter_half_widths=[XTALK_FR / 1000], mode='dpss_solve', 
                                                       eigenval_cutoff=[EIGENVAL_CUTOFF], max_contiguous_edge_flags=len(data.times), cache=cache, filter_dims=0)

    data_xtalk_filt[bl] = d_res

_ = plot_delay_vs_lst(
    data_xtalk_filt, low_rng,
    title=f'Figure 6: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay Transformed Visibility, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 6: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay Transformed Visibility, 50.0—84.2 MHz

_ = plot_delay_vs_lst(
    data_xtalk_filt, high_rng,
    title=f'Figure 7: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay Transformed Visibility, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 7: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay Transformed Visibility, 111.4—230.1 MHz

low_band_meta = plot_delay_fringe_rate(
    data_xtalk_filt, low_rng, delay_bounds=delay_bounds, 
    title=f'Figure 8: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay + Fringe Transformed Visibility, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 8: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay + Fringe Transformed Visibility, 50.0—84.2 MHz

high_band_meta = plot_delay_fringe_rate(
    data_xtalk_filt, high_rng, delay_bounds=delay_bounds,
    title=f'Figure 9: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay + Fringe Transformed Visibility, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 9: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay + Fringe Transformed Visibility, 111.4—230.1 MHz

def calculate_avg_norm_spectra(dfft, flags, nsamples):
    '''Calculates incoherently-averaged spectrum, using ses flags and frequency-averaged nsamples as weights.'''   
    wgts = np.where(np.all(flags, axis=1, keepdims=True), np.nan, np.sum(nsamples, axis=1, keepdims=True))
    avg_norm_spectra = np.nanmean(np.where(np.all(flags, axis=1, keepdims=True), np.nan, np.abs(dfft) * wgts), axis=0) / np.nanmean(wgts)
    return avg_norm_spectra

def coherently_average_visibilities(dc_list, flags, nsamples, rng, nint=20):
    """
    """
    avg_vis = []
    delays = uvtools.utils.fourier_freqs(dc_list[0].freqs[rng]) * 1e9
    flags_for_avg = np.repeat(np.all(_flags, axis=1, keepdims=True), hd.freqs.shape[0], axis=1)
    
    for dc in dc_list:
        avg_delay_spec = {}
        for bl in dc.keys():
            bl_vec = (hd.antpos[bl[1]] - hd.antpos[bl[0]])

            # Coherently average
            ad, af, an, _, _ = frf.timeavg_waterfall(
                dc[bl], nint, nsamples=nsamples[bl], flags=flags_for_avg,
                rephase=True, lsts=dc.lsts, freqs=dc.freqs, bl_vec=bl_vec,
                extra_arrays=dict(times=dc.times), wgt_by_nsample=True, verbose=False
            )

            # Delay Transform
            dfft_avg = uvtools.utils.FFT(ad[:, rng], axis=1, taper='bh7')

            # Incoherently Average - Cohorently Averaged Visibilities
            avg_delay_spec[bl] = calculate_avg_norm_spectra(
                dfft_avg, af[:, rng], an[:, rng]
            )
            
        avg_vis.append(avg_delay_spec)
        
    return delays, avg_vis

def plot_coherently_averaged_1d(delays, dc_list, plot_labels, rng, min_dly=FILT_MIN_DLY, ymin=3e-1, ymax=3e4, title=""):
    # Add search-able title to figure
    display(HTML(f'<h2>{title}</h2>'))
    
    fig, axs = plt.subplots(2, 4, figsize=(14, 10), sharex=True, sharey=True, gridspec_kw={'wspace': 0, 'hspace': 0})

    plot_bls = list(set([bl[:2] for bl in dc_list[0].keys() if bl[0] != bl[1]]))
    plot_bls = sorted(plot_bls, key=lambda k: np.linalg.norm(data.antpos[k[1]] - data.antpos[k[0]]))
    plot_pols = sorted(list(set([bl[-1] for bl in dc_list[0].keys()])))
    
    # Lists for setting vmin and vmax
    ymins = []
    ymaxs = []
        
    for bi, blkey in enumerate(plot_bls):
        for pi, pol in enumerate(plot_pols):
            
            # Store plot output for legend
            plot_objs = []
            
            bl = blkey + (pol, )
            for fi, avg_filt in enumerate(dc_list):
                line = axs[pi, bi].semilogy(delays, avg_filt[bl], label=plot_labels[fi])
                plot_objs.append(line[0])
                
                # Get vmins and maxes
                ymins.append(np.median(avg_filt[bl]))
                ymaxs.append(np.max(avg_filt[bl]))

            axs[pi, bi].set_xlim([-2900, 2900])
            axs[pi, bi].grid()

            # Set widths
            bl_vec = (hd.antpos[bl[1]] - hd.antpos[bl[0]])
            bl_len = np.linalg.norm(bl_vec[:2]) / vis_clean.constants.c.value
            filter_centers, filter_half_widths = vis_clean.gen_filter_properties(ax='freq', horizon=1, standoff=STANDOFF, 
                                                                                 min_dly=min_dly, bl_len=bl_len)
            axs[pi, bi].axvline(-filter_half_widths[0] * 1e9, color='k', ls=':')
            fw_line = axs[pi, bi].axvline(filter_half_widths[0] * 1e9, color='k', ls=':', label='Filter Width')
            plot_objs.append(fw_line)

            axs[pi, bi].axvline(-bl_len * 1e9, color='k', ls='-', lw=0.75)
            horizon_line = axs[pi, bi].axvline(bl_len * 1e9, color='k', ls='-', label='Horizon', lw=0.75)
            plot_objs.append(horizon_line)

            # Get baseline-lengths
            blvec = data.antpos[bl[0]] - data.antpos[bl[1]]
            bllen = np.linalg.norm(blvec)
            if bi < 3:
                bltype = f"{int(np.round(bllen / 14.6))}-unit EW" 
            else:
                bltype = f"{np.abs(blvec[0]):.0f} m {'E' if blvec[0] > 0 else 'W'}, {np.abs(blvec[1]):.0f} m {'N' if blvec[1] > 0 else 'S'}"

            axs[pi, bi].text(0.05, 0.97, f"{bltype} ({pol})\n{np.round(data.freqs[rng[0]] / 1e6, 1)}—{np.round(data.freqs[rng[-1]] / 1e6, 1)} MHz", 
                             horizontalalignment='left', verticalalignment='top', 
                             color='black', transform=axs[pi, bi].transAxes, fontsize=12, 
                             bbox=dict(facecolor='w', alpha=.7, boxstyle='round'))

            if bi == 0:
                axs[pi, bi].set_ylabel(r'$|\widetilde{V}| \ [\rm Jy]$')
            if pi == 1:
                axs[pi, bi].set_xlabel(r'$\tau \ [\rm ns]$')
                
    
    for bi, blkey in enumerate(plot_bls):
        for pi, pol in enumerate(plot_pols):
            axs[pi, bi].set_ylim([np.min(ymins) * 0.1, np.max(ymaxs) * 10])

    fig.legend([pobj for pobj in plot_objs], [pobj._label for pobj in plot_objs], loc='upper center', ncol=len(plot_objs), 
               bbox_to_anchor=(0.5, 1.025))
    plt.tight_layout()
    plt.show()
    

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_inpaint, data_delay_filt, data_xtalk_filt], flags, nsamples, low_rng
)

# Plot coherenetly-averaged visibilities
plot_labels = ['DPSS Inpainted', 'DPSS Delay-Filtered', f'DPSS Delay-Filtered + DPSS Notch-Filtered: {XTALK_FR} mHz']
plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, low_rng,
    title=f'Figure 10: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities,\
    {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 10: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities, 50.0—84.2 MHz

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_inpaint, data_delay_filt, data_xtalk_filt], flags, nsamples, high_rng
)

# Plot coherenetly-averaged visibilities
plot_labels = ['DPSS Inpainted', 'DPSS Delay-Filtered', f'DPSS Delay-Filtered + DPSS Notch-Filtered: {XTALK_FR} mHz']
v = plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, high_rng,
    title=f'Figure 11: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities,\
    {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 11: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities, 111.4—230.1 MHz

Section 3: Form Psuedo-Stokes I + Q

# Create dictionaries for polarized data
data_p_xtalk_filt, flags_p_xtalk_filt, nsamples_p_xtalk_filt = {}, {}, {}
data_p_delay_filt, flags_p_delay_filt, nsamples_p_delay_filt = {}, {}, {}
data_p_inpaint, flags_p_inpaint, nsamples_p_inpaint = {}, {}, {}

for bl in data_inpaint.antpairs():
    # Stokes-I - data, flags, namples
    data_p_xtalk_filt[bl + ('pI',)] = (data_xtalk_filt[bl + ('nn',)] + data_xtalk_filt[bl + ('ee',)]) / 2
    data_p_delay_filt[bl + ('pI',)] = (data_delay_filt[bl + ('nn',)] + data_delay_filt[bl + ('ee',)]) / 2
    data_p_inpaint[bl + ('pI',)] = (data_inpaint[bl + ('nn',)] + data_inpaint[bl + ('ee',)]) / 2
    
    flags_p_xtalk_filt[bl + ('pI',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    flags_p_delay_filt[bl + ('pI',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    flags_p_inpaint[bl + ('pI',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    
    nsamples_p_xtalk_filt[bl + ('pI',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    nsamples_p_delay_filt[bl + ('pI',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    nsamples_p_inpaint[bl + ('pI',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    
    # Stokes Q - data, flags, namples
    data_p_xtalk_filt[bl + ('pQ',)] = (data_xtalk_filt[bl + ('ee',)] - data_xtalk_filt[bl + ('nn',)]) / 2
    data_p_delay_filt[bl + ('pQ',)] = (data_delay_filt[bl + ('ee',)] - data_delay_filt[bl + ('nn',)]) / 2
    data_p_inpaint[bl + ('pQ',)] = (data_inpaint[bl + ('ee',)] - data_inpaint[bl + ('nn',)]) / 2
    
    flags_p_xtalk_filt[bl + ('pQ',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    flags_p_delay_filt[bl + ('pQ',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    flags_p_inpaint[bl + ('pQ',)] = flags[bl + ('nn',)] | flags[bl + ('ee',)]
    
    nsamples_p_xtalk_filt[bl + ('pQ',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    nsamples_p_delay_filt[bl + ('pQ',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    nsamples_p_inpaint[bl + ('pQ',)] = (nsamples[bl + ('nn',)] + nsamples[bl + ('ee',)]) / 2
    
# Attach metadata
data_p_xtalk_filt = DataContainer(data_p_xtalk_filt)
data_p_xtalk_filt.times = data.times
data_p_xtalk_filt.freqs = data.freqs
data_p_xtalk_filt.antpos = data.antpos
data_p_xtalk_filt.lsts = data.lsts

data_p_delay_filt = DataContainer(data_p_delay_filt)
data_p_delay_filt.times = data.times
data_p_delay_filt.freqs = data.freqs
data_p_delay_filt.antpos = data.antpos
data_p_delay_filt.lsts = data.lsts

data_p_inpaint = DataContainer(data_p_inpaint)
data_p_inpaint.times = data.times
data_p_inpaint.freqs = data.freqs
data_p_inpaint.antpos = data.antpos
data_p_inpaint.lsts = data.lsts

_ = plot_delay_vs_lst(
    data_p_xtalk_filt, low_rng,
    title=f'Figure 12: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay-Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 12: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay-Transformed, Psuedo-Stokes I and Q Visibilities, 50.0—84.2 MHz

_ = plot_delay_vs_lst(
    data_p_xtalk_filt, high_rng,
    title=f'Figure 13: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay-Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 13: Real Part of Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay-Transformed, Psuedo-Stokes I and Q Visibilities, 111.4—230.1 MHz

# TODO: Add extend arrow and horizon lines 
_ = plot_delay_fringe_rate(
    data_p_xtalk_filt, low_rng, delay_bounds=delay_bounds, 
    title=f'Figure 14: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz',
    **low_band_meta
)

Figure 14: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, 50.0—84.2 MHz

_ = plot_delay_fringe_rate(
    data_p_xtalk_filt, high_rng, delay_bounds=delay_bounds,
    title=f'Figure 15: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, \
    Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz',
    **high_band_meta
)

Figure 15: Redundantly-Averaged, DPSS Delay + Cross-Talk Filtered, Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, 111.4—230.1 MHz

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_p_inpaint, data_p_delay_filt, data_p_xtalk_filt], flags_p_inpaint, nsamples_p_inpaint, low_rng
)

# Plot coherently-averaged visibilities
plot_labels = [
    'DPSS Inpainted', 'DPSS Delay-Filtered', f'DPSS Delay-Filtered + DPSS Notch-Filtered: {XTALK_FR} mHz'
]
plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, low_rng,
    title=f'Figure 16: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities,\
    Psuedo-Stokes I and Q, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 16: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities, Psuedo-Stokes I and Q, 50.0—84.2 MHz

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_p_inpaint, data_p_delay_filt, data_p_xtalk_filt], flags_p_inpaint, nsamples_p_inpaint, high_rng
)

# Plot coherenetly-averaged visibilities
plot_labels = [
    'DPSS Inpainted', 'DPSS Delay-Filtered', f'DPSS Delay-Filtered + DPSS Notch-Filtered: {XTALK_FR} mHz'
]
plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, high_rng,
    title=f'Figure 17: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities,\
    Psuedo-Stokes I and Q, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 17: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged Visibilities, Psuedo-Stokes I and Q, 111.4—230.1 MHz

Section 4: Main-Beam Fringe-Rate Filter

# Load one file as pyuvdata object to get metadata
uvd = uvdata.UVData()
uvd.read(red_avg_files[len(red_avg_files) // 2], axis='blt', bls=bls_to_plot)
uvd.use_future_array_shapes()

# Load uvbeam file
uvbeam_file_lustre = "/lustre/aoc/projects/hera/h6c-analysis/IDR2/beams/NF_HERA_Vivaldi_efield_beam_healpix.fits"
uvbeam_file_karoo = "/mnt/sn1/data1/beam_models/NF_HERA_Vivaldi_efield_beam_healpix.fits"
uvb = uvbeam.UVBeam()
try:
    uvb.read(uvbeam_file_lustre)
except:
    uvb.read(uvbeam_file_karoo)
uvb.use_future_array_shapes()

# Extract metadata
dt = np.diff(data.times)[0] * 60 * 60 * 24
nfr = data.times.shape[0]
dfr = np.diff(uvtools.utils.fourier_freqs(dt * np.arange(nfr)))[0] * 1000

# Get fringe rates centers and widths for each baseline
frate_centers, frate_half_widths = frf.get_fringe_rate_limits(uvd=uvd, uvb=uvb, nfr=nfr, dfr=dfr)
frate_bounds = {
    k[:2]: (frate_centers[k] - frate_half_widths[k], frate_centers[k] + frate_half_widths[k]) for k in frate_centers
}

data_beam_filt = deepcopy(data_p_xtalk_filt)
cache = {}
for bl in data_p_xtalk_filt:
    # don't bother delay filtering all 0s.
    if np.all(data_p_xtalk_filt[bl] == 0.0) and np.all(flags[bl]):
        continue

    d_mdl = np.zeros_like(data_p_xtalk_filt[bl])

    # Get weights
    wgts = (~np.repeat(np.all(flags_p_inpaint[bl], axis=1, keepdims=True), data.freqs.shape[0], axis=1)).astype(float)

    # Run main-beam filter for each band individually
    for band in [low_band, high_band]:
        d_mdl[:, band], _, info = dspec.fourier_filter(data.times * 24 * 60 * 60, data_p_xtalk_filt[bl][:, band], 
                                                       wgts=wgts[:, band], filter_centers=[frate_centers[bl[:2] + ('nn',)] / 1000], 
                                                       filter_half_widths=[frate_half_widths[bl[:2] + ('nn',)] / 1000], mode='dpss_solve',
                                                       eigenval_cutoff=[EIGENVAL_CUTOFF], max_contiguous_edge_flags=len(data.times), 
                                                       cache=cache, filter_dims=0)

    data_beam_filt[bl] = d_mdl

_ = plot_delay_vs_lst(
    data_beam_filt, low_rng,
    title=f'Figure 18: Real Part of Redundantly-Averaged, DPSS Delay + Main-Beam Filtered, \
    Delay-Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 18: Real Part of Redundantly-Averaged, DPSS Delay + Main-Beam Filtered, Delay-Transformed, Psuedo-Stokes I and Q Visibilities, 50.0—84.2 MHz

_ = plot_delay_vs_lst(
    data_beam_filt, high_rng,
    title=f'Figure 19: Real Part of Redundantly-Averaged, DPSS Delay + Main-Beam Filtered, \
    Delay-Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 19: Real Part of Redundantly-Averaged, DPSS Delay + Main-Beam Filtered, Delay-Transformed, Psuedo-Stokes I and Q Visibilities, 111.4—230.1 MHz

_ = plot_delay_fringe_rate(
    data_beam_filt, low_rng, delay_bounds=delay_bounds, fringe_bounds=frate_bounds,
    title=f'Figure 20: Redundantly-Averaged, DPSS Delay + Cross-Talk + Main-Beam Filtered, \
    Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz',
    **low_band_meta
)

Figure 20: Redundantly-Averaged, DPSS Delay + Cross-Talk + Main-Beam Filtered, Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, 50.0—84.2 MHz

_ = plot_delay_fringe_rate(
    data_beam_filt, high_rng, delay_bounds=delay_bounds, fringe_bounds=frate_bounds,
    title=f'Figure 21: Redundantly-Averaged, DPSS Delay + Cross-Talk + Main-Beam Filtered, \
    Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz',
    **high_band_meta
)

Figure 21: Redundantly-Averaged, DPSS Delay + Cross-Talk + Main-Beam Filtered, Delay + Fringe Transformed, Psuedo-Stokes I and Q Visibilities, 111.4—230.1 MHz

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_p_inpaint, data_p_delay_filt, data_p_xtalk_filt, data_beam_filt], flags_p_inpaint, nsamples_p_inpaint, low_rng
)

plot_labels = [
    'DPSS Inpainted', 'DPSS Delay-Filtered', 'DPSS Delay + Cross Talk-Filtered', 
    f'DPSS Delay-Filtered + Main-Beam Filter'
]
plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, low_rng,
    title=f'Figure 22: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged,\
    Psuedo-Stokes I and Q, Main-Beam Filtered Visibilities, {np.round(data.freqs[low_rng[0]] / 1e6, 1)}—{np.round(data.freqs[low_rng[-1]] / 1e6, 1)} MHz'
)

Figure 22: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged, Psuedo-Stokes I and Q, Main-Beam Filtered Visibilities, 50.0—84.2 MHz

# Coherently-average visibilities
delays, avg_vis = coherently_average_visibilities(
    [data_p_inpaint, data_p_delay_filt, data_p_xtalk_filt, data_beam_filt], 
    flags_p_inpaint, nsamples_p_inpaint, high_rng
)

plot_labels = [
    'DPSS Inpainted', 'DPSS Delay-Filtered', 'DPSS Delay + Cross Talk-Filtered', 
    f'DPSS Delay-Filtered + Main-Beam Filter'
]
plot_coherently_averaged_1d(
    delays, avg_vis, plot_labels, high_rng,
    title=f'Figure 23: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged,\
    Psuedo-Stokes I and Q, Main-Beam Filtered Visibilities, {np.round(data.freqs[high_rng[0]] / 1e6, 1)}—{np.round(data.freqs[high_rng[-1]] / 1e6, 1)} MHz'
)

Figure 23: Coherently-Averaged, Delay-Transformed, Incoherently-Averaged, Psuedo-Stokes I and Q, Main-Beam Filtered Visibilities, 111.4—230.1 MHz

Metadata

for repo in ['numpy', 'scipy', 'astropy', 'hera_cal', 'hera_qm', 'hera_filters', 'hera_notebook_templates', 'pyuvdata']:
    exec(f'from {repo} import __version__')
    print(f'{repo}: {__version__}')

numpy: 2.2.2
scipy: 1.15.1
astropy: 7.0.1
hera_cal: 3.7.1.dev45+g4a0c6f1
hera_qm: 2.2.1.dev2+ga535e9e
hera_filters: 0.1.6.dev9+gf165ec1

hera_notebook_templates: 0.1.dev989+gee0995d
pyuvdata: 3.1.3