smurf_util.py

#!/usr/bin/env python
#-----------------------------------------------------------------------------
# Title      : pysmurf util module - SmurfUtilMixin class
#-----------------------------------------------------------------------------
# File       : pysmurf/util/smurf_util.py
# Created    : 2018-08-29
#-----------------------------------------------------------------------------
# This file is part of the pysmurf software package. It is subject to
# the license terms in the LICENSE.txt file found in the top-level directory
# of this distribution and at:
#    https://confluence.slac.stanford.edu/display/ppareg/LICENSE.html.
# No part of the pysmurf software package, including this file, may be
# copied, modified, propagated, or distributed except according to the terms
# contained in the LICENSE.txt file.
#-----------------------------------------------------------------------------
from contextlib import contextmanager
import glob
import os
import threading
import time
import re

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

from pysmurf.client.base import SmurfBase
from pysmurf.client.command.sync_group import SyncGroup as SyncGroup
from pysmurf.client.util.SmurfFileReader import SmurfStreamReader
from pysmurf.client.util.pub import set_action

class SmurfUtilMixin(SmurfBase):

    @set_action()
    def take_debug_data(self, band, channel=None, nsamp=2**19, filename=None,
            IQstream=1, single_channel_readout=1, debug=False, rf_iq=False,
            write_log=True):
        """Takes raw debugging data.

        Args
        ----
        band : int
            The band to take data on.
        channel : int or None, optional, default None
            The channel to take debug data on in single_channel_mode.
        nsamp : int, optional, default 2**19
            The number of samples to take.
        filename : str or None, optional, default None
            The name of the file to save to.
        IQstream : int, optional, default 1
            Whether to take the raw IQ stream.
        single_channel_readout : int, optional, default 1
            Whether to look at one channel.
        debug : bool, optional, default False
            Whether to take data in debug mode.
        rf_iq : bool, optional, default False
            Return the RF IQ. Must provide channel.
        write_log : bool, optional, default True
            Whether to write low-level commands to the log file.

        Returns
        -------
        f : float array
            The frequency response.
        df : float array
            The frequency error.
        sync : float array
            The sync count.

        """
        # Set proper single channel readout
        if channel is not None:
            if rf_iq:
                IQstream = False
                single_channel_readout = 2
                self.set_rf_iq_stream_enable(band, 1)

            if single_channel_readout == 1:
                self.set_single_channel_readout(band, 1)
                self.set_single_channel_readout_opt2(band, 0)
            elif single_channel_readout == 2:
                self.set_single_channel_readout(band, 0)
                self.set_single_channel_readout_opt2(band, 1)
            else:
                self.log('single_channel_readout must be 1 or 2',
                    self.LOG_ERROR)
                raise ValueError('single_channel_readout must be 1 or 2')
            self.set_readout_channel_select(band, channel, write_log=write_log)
        else: # exit single channel otherwise
            self.set_single_channel_readout(band, 0, write_log=write_log)
            self.set_single_channel_readout_opt2(band, 0, write_log=write_log)

        # Set IQstream
        if IQstream == 1:
            self.set_iq_stream_enable(band, 1)
        else:
            self.set_iq_stream_enable(band, 0)

        # set filename
        if filename is not None:
            data_filename = os.path.join(self.output_dir, filename+'.dat')
            self.log(f'Writing to file : {data_filename}',
                self.LOG_USER)
        else:
            timestamp = self.get_timestamp()
            data_filename = os.path.join(self.output_dir, timestamp+'.dat')
            self.log(f'Writing to file : {data_filename}',
                self.LOG_USER)

        dtype = 'debug'
        dchannel = 0 # I don't really know what this means and I'm sorry -CY
        self.setup_daq_mux(dtype, dchannel, nsamp, band=band, debug=debug)
        self.log('Data acquisition in progress...', self.LOG_USER)
        char_array = [ord(c) for c in data_filename] # convert to ascii
        write_data = np.zeros(300, dtype=int)
        for j in np.arange(len(char_array)):
            write_data[j] = char_array[j]

        self.set_streamdatawriter_datafile(write_data) # write this

        #self.set_streamdatawriter_open('True') # str and not bool
        self.set_streamdatawriter_open(True)

        bay=self.band_to_bay(band)
        self.set_trigger_daq(bay, 1, write_log=True) # this seems to = TriggerDM

        time.sleep(.1) # maybe unnecessary

        done=False
        while not done:
            done=True
            for k in range(2):
                # see pysmurf issue 161.  This call is no longer used,
                # and causes take_debug_data to crash if
                # get_waveform_wr_addr is called before the
                # acquisition completes.
                #wr_addr = self.get_waveform_wr_addr(bay, engine=0)
                empty = self.get_waveform_empty(bay, engine=k)
                if not empty:
                    done=False
            time.sleep(1)

        time.sleep(.25) # do we need all of these?

        # Close the streamdatawriter
        self.set_streamdatawriter_close(True)

        self.log('Done taking data', self.LOG_USER)
        self.pub.register_file(data_filename, 'debug', format='dat')

        if rf_iq:
            self.set_rf_iq_stream_enable(band, 0)

        if single_channel_readout > 0:
            f, df, sync = self.decode_single_channel(data_filename)
        else:
            f, df, sync = self.decode_data(data_filename)

        return f, df, sync

    # the JesdWatchdog will check if an instance of the JesdWatchdog is already
    # running and kill itself if there is
    def start_jesd_watchdog(self):
        import pysmurf.client.watchdog.JesdWatchdog as JesdWatchdog
        import subprocess
        import sys
        subprocess.Popen([sys.executable,JesdWatchdog.__file__])

    # Shawn needs to make this better and add documentation.
    @set_action()
    def estimate_phase_delay(self, band, nsamp=2**19, make_plot=True,
            show_plot=True, save_plot=True, save_data=True, n_scan=5,
            timestamp=None, uc_att=None, dc_att=None, freq_min=-2.4E6, freq_max=2.4E6):
        """Estimates total system latency for requested band.

        Measures the analog and digital (=processing) phase delay (or
        latency) by sweeping the [`freq_min`, `freq_max`] interval of the
        requested 500 MHz band.  Estimates the best values for the
        `refPhaseDelay` and `refPhaseDelayFine` firmware registers to
        compensate for the measured delay.  Three steps:

        1. Measures the analog latency using the [`freq_min`, `freq_max`]
           sub-interval of data taken using the
           :func:`~pysmurf.client.tune.smurf_tune.SmurfTuneMixin.full_band_resp`
           routine.
        2. Measures the total latency using the
           :func:`~pysmurf.client.tune.smurf_tune.SmurfTuneMixin.find_freq`
           routine.
        3. Sets the `refPhaseDelay` and `refPhaseDelayFine` registers
           to the values computed to best compensate for the total
           measured latency and re-measures the total latency again
           using the
           :func:`~pysmurf.client.tune.smurf_tune.SmurfTuneMixin.find_freq`
           routine.

        On completion, the `refPhaseDelay` and `refPhaseDelayFine`
        registers are set to the estimated optimal values.

        Args
        ----
        band : int
           The band to estimate phase delay on.
        nsamp : int, optional, default 2**19
           The number of samples to take.
        make_plot : bool, optional, default True
           Whether or not to make plots.
        show_plot : bool, optional, default True
           Whether or not to show plots.
        save_plot : bool, optional, default True
           Whether or not to save plot to file.
        save_data : bool, optional, default True
           Whether or not to save data to file.
        n_scan : int, optional, default 5
           Number of scans to do to estimate analog phase delay,
           passed to
           :func:`~pysmurf.client.tune.smurf_tune.SmurfTuneMixin.full_band_resp`.
        timestamp : str or None, optional, default None
           ctime to timestamp the plot and data with (if saved to
           file).  If None, it gets the time stamp right before
           acquiring data.
        uc_att : int or None, optional, default None
           UC attenuator setting to use during measurements.  If None,
           uses currently programmed setting.
        dc_att : int or None, optional, default None
           DC attenuator setting to use during measurements.  If None,
           uses currently programmed setting.
        freq_min : float, optional, default -2.4E6
           Lower bound of the frequency interval used to estimate the
           phase delay, in Hz.  From the center of the 500 MHz band.
        freq_max : float, optional, default 2.4E6
           Upper bound of the frequency interval used to estimate the
           phase delay, in Hz.  From the center of the 500 MHz band.

        Returns
        -------
        processing_delay_us : float
           Estimated processing phase delay, in microseconds.
        dsp_corr_delay_us : float
           Residual total phase delay, for estimated `refPhaseDelay`
           and `refPhaseDelayFine`.

        """

        self.set_band_delay_us(band, 0)

        uc_att0 = self.get_att_uc(band)
        dc_att0 = self.get_att_dc(band)
        if uc_att is None:
            uc_att = self.get_att_uc(band)
        if dc_att is None:
            dc_att = self.get_att_dc(band)
        self.set_att_uc(band, uc_att, write_log=True)
        self.set_att_dc(band, dc_att, write_log=True)

        # only loop over dsp subbands in requested frequency range (to
        # save time)
        n_subbands = self.get_number_sub_bands(band)
        digitizer_frequency_mhz = self.get_digitizer_frequency_mhz(band)
        subband_half_width_mhz = digitizer_frequency_mhz/\
            n_subbands
        subbands,subband_centers=self.get_subband_centers(band)
        subband_freq_min=-subband_half_width_mhz/2.
        subband_freq_max=subband_half_width_mhz/2.
        dsp_subbands=[]
        for sb,sbc in zip(subbands,subband_centers):
            # ignore unprocessed sub-bands
            if sb not in subbands:
                continue
            lower_sb_freq=sbc+subband_freq_min
            upper_sb_freq=sbc+subband_freq_max
            if lower_sb_freq>=(freq_min/1.e6-subband_half_width_mhz) and \
                    upper_sb_freq<=(freq_max/1.e6+subband_half_width_mhz):
                dsp_subbands.append(sb)

        if timestamp is None:
            timestamp = self.get_timestamp()

        if make_plot:
            if show_plot:
                plt.ion()
            else:
                plt.ioff()

        bay=int(band/4)
        fw_abbrev_sha=self.get_fpga_git_hash_short()

        self.band_off(band)
        self.flux_ramp_off()

        self.log('Running full band resp')
        freq_cable, resp_cable = self.full_band_resp(
            band, nsamp=nsamp, make_plot=make_plot,
            save_data=save_data, n_scan=n_scan)

        idx_cable = np.where( (freq_cable > freq_min) & (freq_cable < freq_max) )

        cable_z = np.polyfit(freq_cable[idx_cable], np.unwrap(np.angle(resp_cable[idx_cable])), 1)
        cable_p = np.poly1d(cable_z)
        cable_delay_us=np.abs(1.e6*cable_z[0]/2/np.pi)

        freq_cable_subset=freq_cable[idx_cable]
        resp_cable_subset=resp_cable[idx_cable]
        #### done measuring cable delay

        #### start measuring dsp delay (cable+processing)
## FIXME -- should be able to scan with "0" delay, not working
        self.set_band_delay_us(band, 1)

        self.log('Running find_freq')
        #freq_dsp,resp_dsp=self.find_freq(band, start_freq=freq_min, stop_freq=freq_max)
        freq_dsp,resp_dsp=self.find_freq(band,subband=dsp_subbands)

        # only preserve data in the subband half width
        freq_dsp_subset=[]
        resp_dsp_subset=[]
        est_delay=[]
        for sb,sbc in zip(subbands,subband_centers):
            freq_subband=freq_dsp[sb]-sbc
            idx = np.where( ( freq_subband > subband_freq_min ) &
                (freq_subband < subband_freq_max) )
            if len(idx[0]) > 0:
                dsp_z = np.polyfit(freq_dsp[sb][idx]*1e6, np.unwrap(np.angle(resp_dsp[sb][idx])), 1)
                dsp_p = np.poly1d(dsp_z)
                dsp_delay_us=np.abs(1.e6*dsp_z[0]/2/np.pi)
                dsp_delay_us=dsp_delay_us + self.get_band_delay_us(band)
                est_delay.append(dsp_delay_us)
            freq_dsp_subset.extend(freq_dsp[sb][idx])
            resp_dsp_subset.extend(resp_dsp[sb][idx])

        freq_dsp_subset=np.array(freq_dsp_subset)
        resp_dsp_subset=np.array(resp_dsp_subset)

        idx_dsp = np.where( (freq_dsp_subset > freq_min) &
            (freq_dsp_subset < freq_max) )

        # restrict to requested frequencies only
        freq_dsp_subset=freq_dsp_subset[idx_dsp]
        resp_dsp_subset=resp_dsp_subset[idx_dsp]

        # to Hz
        freq_dsp_subset=(freq_dsp_subset)*1.0E6

        # fit
        dsp_z = np.polyfit(freq_dsp_subset, np.unwrap(np.angle(resp_dsp_subset)), 1)
        dsp_p = np.poly1d(dsp_z)
        dsp_delay_us=np.abs(1.e6*dsp_z[0]/2/np.pi)
        dsp_delay_us=dsp_delay_us + self.get_band_delay_us(band)
        dsp_delay_us=np.mean(est_delay)

        processing_delay_us=dsp_delay_us-cable_delay_us

        print('-------------------------------------------------------')
        print(f'Estimated cable_delay_us={cable_delay_us}')
        print(f'Estimated dsp_delay_us={dsp_delay_us}')
        print(f'Estimated processing_delay_us={processing_delay_us}')
        print('-------------------------------------------------------')

        #### done measuring dsp delay (cable+processing)

        #### start measuring total (DSP + cable) delay with estimated correction applied
        self.set_band_delay_us(band, dsp_delay_us)

        self.log('Running find_freq')
        freq_dsp_corr,resp_dsp_corr=self.find_freq(band,subband=dsp_subbands)

        freq_dsp_corr_subset=[]
        resp_dsp_corr_subset=[]
        first = True
        for sb,sbc in zip(subbands,subband_centers):
            freq_subband=freq_dsp_corr[sb]-sbc
            idx = np.where( ( freq_subband > subband_freq_min ) & (freq_subband < subband_freq_max) )
            if len(idx[0]) > 0:
                if not first:
                    last_phase = np.angle(resp_dsp_corr_subset[-1])
                    new_phase  = np.angle(resp_dsp_corr[sb][idx[0]])
                    resp_dsp_corr[sb][idx] = resp_dsp_corr[sb][idx] * np.exp(1j*(last_phase - new_phase))

                if first:
                    first = False
            freq_dsp_corr_subset.extend(freq_dsp_corr[sb][idx])
            resp_dsp_corr_subset.extend(resp_dsp_corr[sb][idx])

        freq_dsp_corr_subset=np.array(freq_dsp_corr_subset)
        resp_dsp_corr_subset=np.array(resp_dsp_corr_subset)

        # restrict to requested frequency subset
        idx_dsp_corr = np.where( (freq_dsp_corr_subset > freq_min) & (freq_dsp_corr_subset < freq_max) )

        # restrict to requested frequencies only
        freq_dsp_corr_subset=freq_dsp_corr_subset[idx_dsp_corr]
        resp_dsp_corr_subset=resp_dsp_corr_subset[idx_dsp_corr]

        # to Hz
        freq_dsp_corr_subset=(freq_dsp_corr_subset)*1.0E6

        # fit
        dsp_corr_z = np.polyfit(freq_dsp_corr_subset, np.unwrap(np.angle(resp_dsp_corr_subset)), 1)
        dsp_corr_delay_us=np.abs(1.e6*dsp_corr_z[0]/2/np.pi)
        #### done measuring total (DSP) delay with estimated correction applied

        # plot unwraped phase in top panel, subtracted in bottom

        fig, ax = plt.subplots(3, figsize=(6,7.5), sharex=True)

        f_cable_plot = (freq_cable_subset) / 1.0E6
        cable_phase = np.unwrap(np.angle(resp_cable_subset))

        f_dsp_plot = (freq_dsp_subset) / 1.0E6
        dsp_phase = np.unwrap(np.angle(resp_dsp_subset))

        f_dsp_corr_plot = (freq_dsp_corr_subset) / 1.0E6
        dsp_corr_phase = np.unwrap(np.angle(resp_dsp_corr_subset))

        ax[0].set_title(f'AMC in Bay {bay}, Band {band} Cable Delay')
        ax[0].plot(f_cable_plot,cable_phase,label='Cable (full_band_resp)',
            c='g', lw=3)
        ax[0].plot(f_cable_plot,cable_p(f_cable_plot*1.0E6),'m--',
            label='Cable delay fit',lw=3)

        ax[1].set_title(f'AMC in Bay {bay}, Band {band} DSP Delay')
        ax[1].plot(f_dsp_plot,dsp_phase,label='DSP (find_freq)',c='c',lw=3)
        ax[1].plot(f_dsp_plot,dsp_p(f_dsp_plot*1.0E6), c='orange', ls='--',
                   label='DSP delay fit', lw=3)

        ax[0].set_ylabel("Phase [rad]")
        ax[0].set_xlabel('Frequency offset from band center [MHz]')

        ax[1].set_ylabel("Phase [rad]")
        ax[1].set_xlabel('Frequency offset from band center [MHz]')

        ax[0].legend(loc='lower left',fontsize=8)
        ax[1].legend(loc='lower left',fontsize=8)

        bbox = dict(boxstyle="round", ec='w', fc='w', alpha=.65)
        ax[0].text(.97, .90, f'cable delay={cable_delay_us:.5f} us',
                   transform=ax[0].transAxes, fontsize=10,
                   bbox=bbox,horizontalalignment='right')

        ax[1].text(.97, .90, f'dsp delay={dsp_delay_us:.5f} us',
                   transform=ax[1].transAxes, fontsize=10,
                   bbox=bbox,horizontalalignment='right')

        cable_residuals=cable_phase-(cable_p(f_cable_plot*1.0E6))
        ax[2].plot(f_cable_plot,cable_residuals-np.median(cable_residuals),
            label='Cable (full_band_resp)',c='g')
        ax[2].plot(f_dsp_corr_plot,dsp_corr_phase-np.median(dsp_corr_phase),
            label='DSP (find_freq)', c='c')
        ax[2].set_title(f'AMC in Bay {bay}, Band {band} Residuals'.format(bay,band))
        ax[2].set_ylabel("Residual [rad]")
        ax[2].set_xlabel('Frequency offset from band center [MHz]')
        ax[2].set_ylim([-5,5])

        ax[2].text(.97, .76,
                   f'processing delay={processing_delay_us:.5f} us (fw={fw_abbrev_sha})',
                   transform=ax[2].transAxes, fontsize=8,
                   bbox=bbox,horizontalalignment='right')
        ax[2].text(.97, .68, f'delay post-correction={dsp_corr_delay_us*1000.:.3f} ns',
                   transform=ax[2].transAxes, fontsize=8,
                bbox=bbox,horizontalalignment='right')

        ax[2].legend(loc='upper left',fontsize=8)

        plt.tight_layout()

        if save_plot:
            save_name = f'{timestamp}_b{band}_delay.png'

            path = os.path.join(self.plot_dir, save_name)
            plt.savefig(path,bbox_inches='tight')
            self.pub.register_file(path, 'delay', plot=True)

            if not show_plot:
                plt.close()


        self.log('Setting attenuator values back to original values')
        self.log(f'UC Att: {uc_att0}')
        self.log(f'DC Att: {dc_att0}')

        self.set_att_uc(band, uc_att0, write_log=True)
        self.set_att_dc(band, dc_att0, write_log=True)

        if show_plot:
            plt.show()

        return dsp_delay_us, dsp_corr_delay_us

    def process_data(self, filename, dtype=np.uint32):
        """ Reads a file taken with take_debug_data and processes it into data
        and header.

        Args
        ----
        filename : str
            Path to file
        dtype : numpy.dtype, optional, default numpy.uint32
            datatype to cast to.

        Returns
        -------
        header : numpy.ndarray
            The header information.
        data : numpy.ndarray
            The resonator data.
        """
        n_chan = 2 # number of stream channels
        #header_size = 4 # 8 bytes in 16-bit word

        rawdata = np.fromfile(filename, dtype='<u4').astype(dtype)

        # -1 is equiv to [] in Matlab
        rawdata = np.transpose(np.reshape(rawdata, (n_chan, -1)))

        if dtype==np.uint32:
            header = rawdata[:2, :]
            data = np.delete(rawdata, (0,1), 0).astype(dtype)
        elif dtype==np.int32:
            header = np.zeros((2,2))
            header[:,0] = rawdata[:2,0].astype(np.uint32)
            header[:,1] = rawdata[:2,1].astype(np.uint32)
            data = np.double(np.delete(rawdata, (0,1), 0))
        elif dtype==np.int16:
            header1 = np.zeros((4,2))
            header1[:,0] = rawdata[:4,0].astype(np.uint16)
            header1[:,1] = rawdata[:4,1].astype(np.uint16)
            header1 = np.double(header1)
            header = header1[::2] + header1[1::2] * (2**16) # what am I doing
        else:
            raise TypeError(f'Type {dtype} not yet supported!')

        if (header[1,1]>>24 == 2) or (header[1,1]>>24 == 0):
            header = np.fliplr(header)
            data = np.fliplr(data)

        return header, data

    @set_action()
    def decode_data(self, filename, swapFdF=False, recast=True, truncate=True):
        """ Take a dataset from take_debug_data and spit out results.

        Args
        ----
        filename : str
            Path to file.
        swapFdF : bool, optional, default False
            Whether the F and dF (or I/Q) streams are flipped.
        recast : bool, optional, default True
            Whether to recast from size n_channels_processed to
            n_channels.
        truncate : bool, optional, default True
            Truncates the data if the number of elements returned is
            not an integer multiple of the sample rate.

        Returns
        -------
        f : numpy.ndarray
            If iqStreamEnable = 0. f is the tracking frequency.
            Otherwise if iqStreamEnable = 1. f is the demodulated
            in-phase tracking component.
        df : numpy.ndarray
            If iqStreamEnable = 0. df is the tracking frequency error.
            Otherwise if iqStreamEnable = 1. f is the demodulated
            quadrature tracking component.
        flux_ramp_strobe : numpy.ndarray
            The synchronizing pulse.
        """
        n_proc = self.get_number_processed_channels()
        n_chan = self.get_number_channels()

        n_subbands = self.get_number_sub_bands()
        digitizer_frequency_mhz = self.get_digitizer_frequency_mhz()
        subband_half_width_mhz = (digitizer_frequency_mhz / n_subbands)

        header, rawdata = self.process_data(filename)

        # decode strobes
        strobes = np.floor(rawdata / (2**30))
        data = 0x00FFFFFF & rawdata
        ch0_strobe = np.remainder(strobes, 2)
        flux_ramp_strobe = np.floor((strobes - ch0_strobe) / 2)

        # decode frequencies
        ch0_idx = np.where(ch0_strobe[:,0] == 1)[0]
        f_first = ch0_idx[0]
        f_last = ch0_idx[-1]

        freqs = data[f_first:f_last, 0]
        neg = np.where(freqs >= 2**23)[0]
        f = np.double(freqs)
        if len(neg) > 0:
            f[neg] = f[neg] - 2**24

        if np.remainder(len(f), n_proc)!=0:
            if truncate:
                self.log(f'Number of points in f not a multiple of {n_proc}.' +
                    f' Truncating f to the nearest multiple of {n_proc}.',
                    self.LOG_USER)
                f=f[:(len(f)-np.remainder(len(f),n_proc))]
            else:
                self.log(f'Number of points in f not a multiple of {n_proc}.'+
                    ' Cannot decode', self.LOG_ERROR)
        f = np.reshape(f, (-1, n_proc)) * subband_half_width_mhz / 2**23

        # frequency errors
        ch0_idx_df = np.where(ch0_strobe[:,1] == 1)[0]
        if len(ch0_idx_df) > 0:
            d_first = ch0_idx_df[0]
            d_last = ch0_idx_df[-1]
            dfreq = data[d_first:d_last, 1]
            neg = np.where(dfreq >= 2**23)[0]
            df = np.double(dfreq)
            if len(neg) > 0:
                df[neg] = df[neg] - 2**24

            if np.remainder(len(df), n_proc)!=0:
                if truncate:
                    self.log('Number of points in df not a multiple of '+
                        f'{n_proc}. Truncating df to the nearest multiple ' +
                        f' of {n_proc}.', self.LOG_USER)
                    df=df[:(len(df)-np.remainder(len(df),n_proc))]
                else:
                    self.log(f'Number of points in df not a multiple of {n_proc}.' +
                        'Cannot decode', self.LOG_ERROR)
            df = np.reshape(df, (-1, n_proc)) * subband_half_width_mhz / 2**23

        else:
            df = []

        if recast:
            nsamp, nprocessed = np.shape(f)
            nsamp_df, _ = np.shape(df)
            if nsamp != nsamp_df:
                self.log('f and df are different sizes. Choosing the smaller'
                    ' value. Not sure why this is happening.')
                nsamp = np.min([nsamp, nsamp_df])

            ftmp = np.zeros((nsamp, n_chan))
            dftmp = np.zeros_like(ftmp)

            processed_ind = self.get_processed_channels()
            ftmp[:, processed_ind] = f[:nsamp]
            dftmp[:, processed_ind] = df[:nsamp]

            f = ftmp
            df = dftmp

        return f, df, flux_ramp_strobe

    @set_action()
    def decode_single_channel(self, filename, swapFdF=False):
        """
        decode take_debug_data file if in singlechannel mode

        Args
        ----
        filename : str
            Path to file to decode.
        swapFdF : bool, optional, default False
            Whether to swap f and df streams.

        Returns
        -------
        list
            [f, df, sync] if iq_stream_enable = False
            [I, Q, sync] if iq_stream_enable = True
        """

        n_subbands = self.get_number_sub_bands()
        digitizer_frequency_mhz = self.get_digitizer_frequency_mhz()
        subband_half_width_mhz = (digitizer_frequency_mhz / n_subbands)

        if swapFdF:
            nF = 1
            nDF = 0
        else:
            nF = 0
            nDF = 1

        header, rawdata = self.process_data(filename)

        # decode strobes
        strobes = np.floor(rawdata / (2**30))
        data = rawdata - (2**30)*strobes
        ch0_strobe = np.remainder(strobes, 2)
        flux_ramp_strobe = np.floor((strobes - ch0_strobe) / 2)

        # decode frequencies
        freqs = data[:,nF]
        neg = np.where(freqs >= 2**23)[0]
        f = np.double(freqs)
        if len(neg) > 0:
            f[neg] = f[neg] - 2**24

        f = np.transpose(f) * subband_half_width_mhz / 2**23

        dfreqs = data[:,nDF]
        neg = np.where(dfreqs >= 2**23)[0]
        df = np.double(dfreqs)
        if len(neg) > 0:
            df[neg] = df[neg] - 2**24

        df = np.transpose(df) * subband_half_width_mhz / 2**23

        return f, df, flux_ramp_strobe

    @set_action(action=None)
    def take_stream_data(self, meas_time, downsample_factor=None,
                         write_log=True, update_payload_size=True,
                         reset_unwrapper=True, reset_filter=True,
                         return_data=False, make_freq_mask=True,
                         register_file=False):
        """
        Takes streaming data for a given amount of time

        To do: move downsample_factor to config table

        Args
        ----
        meas_time : float
            The amount of time to observe for in seconds.
        downsample_factor : int or None, optional, default None
            The number of fast sample (the flux ramp reset rate -
            typically 4kHz) to skip between reporting. If None, does
            not update.
        write_log : bool, optional, default True
            Whether to write to the log file.
        update_payload_size : bool, optional, default True
            Whether to update the payload size (the number of channels
            written to disk). If the number of channels on is greater
            than the payload size, then only the first N channels are
            written. This bool will update the payload size to be the
            same as the number of channels on across all bands)
        reset_unwrapper : bool, optional, default True
            Whether to reset the unwrapper before taking data.
        reset_filter : bool, optional, default True
            Whether to reset the filter before taking data.
        return_data : bool, optional, default False
            Whether to return the data. If False, returns the full
            path to the data.
        make_freq_mask : bool, optional, default True
            Whether to write a text file with resonator frequencies.
        register_file : bool, optional, default False
            Whether to register the data file with the pysmurf
            publisher.


        Returns
        -------
        data_filename : str
            The fullpath to where the data is stored.
        """
        if write_log:
            self.log('Starting to take data.', self.LOG_USER)
        data_filename = self.stream_data_on(downsample_factor=downsample_factor,
            update_payload_size=update_payload_size, write_log=write_log,
            reset_unwrapper=reset_unwrapper, reset_filter=reset_filter,
            make_freq_mask=make_freq_mask)

        # Sleep for the full measurement time
        time.sleep(meas_time)

        # Stop acq
        self.stream_data_off(write_log=write_log, register_file=register_file)

        if write_log:
            self.log('Done taking data.', self.LOG_USER)

        if return_data:
            t, d, m = self.read_stream_data(data_filename)
            return t, d, m
        else:
            return data_filename

    @contextmanager
    def stream_data_cm(self, write_log=True, register_file=False,
                       **stream_on_kwargs):
        """
        Context manager for data streaming. If you intend to turn streaming
        on, do something, and then turn streaming off this is a safe way to make
        sure streaming is in fact stopped properly even if an error is raised.

        Args
        ----
        write_config : bool, optional, default False
            Whether to dump the entire config. Warning this can be
            slow.
        data_filename : str or None, optional, default None
            The full path to store the data. If None, it uses the
            timestamp.
        downsample_factor : int or None, optional, default None
            The number of fast samples to skip between sending.
        write_log : bool, optional, default True
            Whether to write to the log file.
        update_payload_size : bool, optional, default True
            Whether to update the payload size (the number of channels
            written to disk). If the number of channels on is greater
            than the payload size, then only the first N channels are
            written. This bool will update the payload size to be the
            same as the number of channels on across all bands)
        reset_filter : bool, optional, default True
            Whether to reset the filter before taking data.
        reset_unwrapper : bool, optional, default True
            Whether to reset the unwrapper before taking data.
        make_freq_mask : bool, optional, default True
            Whether to write a text file with resonator frequencies.
        register_file : bool, optional, default False
            If true, the stream data file will be registered through
            the publisher.

        Yields
        -------
        data_filename : str
            The fullpath to where the data is stored.
        """
        data_filename = self.stream_data_on(write_log=write_log, **stream_on_kwargs)
        try:
            yield data_filename
        finally:
            self.stream_data_off(write_log=write_log,
                                 register_file=register_file)

    @set_action()
    def stream_data_on(self, write_config=False, data_filename=None,
                       downsample_factor=None, write_log=True,
                       update_payload_size=True, reset_filter=True,
                       reset_unwrapper=True, make_freq_mask=True,
                       channel_mask=None, make_datafile=True,
                       filter_wait_time=0.1):
        """
        Turns on streaming data.

        Args
        ----
        write_config : bool, optional, default False
            Whether to dump the entire config. Warning this can be
            slow.
        data_filename : str or None, optional, default None
            The full path to store the data. If None, it uses the
            timestamp.
        downsample_factor : int or None, optional, default None
            The number of fast samples to skip between sending.
        write_log : bool, optional, default True
            Whether to write to the log file.
        update_payload_size : bool, optional, default True
            Whether to update the payload size (the number of channels
            written to disk). If this is True, will set the payload size to
            0, which tells rogue to automatically adjust it based on the
            channel count.
        reset_filter : bool, optional, default True
            Whether to reset the filter before taking data.
        reset_unwrapper : bool, optional, default True
            Whether to reset the unwrapper before taking data.
        make_freq_mask : bool, optional, default True
            Whether to write a text file with resonator frequencies.
        channel_mask : list or None, optional, default None
            Channel mask to set before streamig data. This should be an array
            of absolute smurf channels between 0 and
            ``nbands * chans_per_band``. If None will create the channel mask
            containing all channels with a non-zero tone amplitude.
        make_datafile : bool, optional, default True
            Whether to create a datafile.
        filter_wait_time : float, optional, default 0.1
            Time in seconds to wait after filter reset.

        Returns
        -------
        data_filename : str
            The fullpath to where the data is stored.
        """
        bands = self._bands

        if downsample_factor is not None:
            self.set_downsample_factor(downsample_factor)
        else:
            downsample_factor = self.get_downsample_factor()
            if write_log:
                self.log('Input downsample factor is None. Using '+
                     'value already in pyrogue:'+
                     f' {downsample_factor}')

        if update_payload_size:
            self.set_payload_size(0)

        # Check if flux ramp is non-zero
        ramp_max_cnt = self.get_ramp_max_cnt()
        if ramp_max_cnt == 0:
            self.log('Flux ramp frequency is zero. Cannot take data.',
                self.LOG_ERROR)
        else:
            # check which flux ramp relay state we're in
            # read_ac_dc_relay_status() should be 0 in DC mode, 3 in
            # AC mode.  this check is only possible if you're using
            # one of the newer C02 cryostat cards.
            flux_ramp_ac_dc_relay_status = self.C.read_ac_dc_relay_status()
            if flux_ramp_ac_dc_relay_status == 0:
                if write_log:
                    self.log("FLUX RAMP IS DC COUPLED.", self.LOG_USER)
            elif flux_ramp_ac_dc_relay_status == 3:
                if write_log:
                    self.log("Flux ramp is AC-coupled.", self.LOG_USER)
            else:
                self.log("flux_ramp_ac_dc_relay_status = " +
                         f"{flux_ramp_ac_dc_relay_status} " +
                         "- NOT A VALID STATE.", self.LOG_ERROR)

            if channel_mask is None:
                # Creates a channel mask with all channels that have enabled
                # tones
                smurf_chans = {}
                for b in bands:
                    smurf_chans[b] = self.which_on(b)

                channel_mask = self.make_channel_mask(bands, smurf_chans)
                self.set_channel_mask(channel_mask)
            else:
                channel_mask = np.atleast_1d(channel_mask)
                self.set_channel_mask(channel_mask)

            time.sleep(0.5)

            # start streaming before opening file
            # to avoid transient filter step
            self.set_stream_enable(1, write_log=False, wait_done=True)

            if reset_unwrapper:
                self.set_unwrapper_reset(write_log=write_log)
            if reset_filter:
                self.set_filter_reset(write_log=write_log)
            if reset_unwrapper or reset_filter:
                time.sleep(filter_wait_time)

            # Make the data file
            timestamp = self.get_timestamp()
            if data_filename is None:
                data_filename = os.path.join(self.output_dir, timestamp+'.dat')

            if make_datafile:
                self.set_data_file_name(data_filename)

            # Optionally write PyRogue configuration
            if write_config:
                config_filename=os.path.join(self.output_dir, timestamp+'.yml')
                if write_log:
                    self.log('Writing PyRogue configuration to file : ' +
                        f'{config_filename}', self.LOG_USER)
                self.write_config(config_filename)

                # short wait
                time.sleep(5.)
            if write_log:
                self.log(f'Writing to file : {data_filename}',
                         self.LOG_USER)


            # Save mask file as text file. Eventually this will be in the
            # raw data output
            mask_fname = os.path.join(data_filename.replace('.dat',
                '_mask.txt'))
            np.savetxt(mask_fname, channel_mask, fmt='%i')
            self.pub.register_file(mask_fname, 'mask')
            self.log(mask_fname)

            if make_freq_mask:
                if write_log:
                    self.log("Writing frequency mask.")
                freq_mask = self.make_freq_mask(channel_mask)
                np.savetxt(os.path.join(data_filename.replace('.dat',
                    '_freq.txt')), freq_mask, fmt='%4.4f')
                self.pub.register_file(
                    os.path.join(data_filename.replace('.dat', '_freq.txt')),
                    'mask', format='txt')

            if make_datafile:
                self.open_data_file(write_log=write_log)

            return data_filename

    @set_action()
    def stream_data_off(self, write_log=True, register_file=False):
        """
        Turns off streaming data.

        Args
        ----
        write_log : bool, optional, default True
            Whether to log the CA commands or not.
        register_file : bool, optional, default False
            If true, the stream data file will be registered through
            the publisher.
        """
        self.close_data_file(write_log=write_log)

        if register_file:
            datafile = self.get_data_file_name().tostring().decode()
            if datafile:
                self.log(f"Registering File {datafile}")
                self.pub.register_file(datafile, 'data', format='dat')

        self.set_stream_enable(0, write_log=write_log, wait_after=.15)

    @set_action()
    def read_stream_data(self, datafile, channel=None,
                         nsamp=None, array_size=None,
                         return_header=False,
                         return_tes_bias=False, write_log=True,
                         n_max=2048, make_freq_mask=False,
                         gcp_mode=False):
        """
        Loads data taken with the function stream_data_on.
        Gives back the resonator data in units of phase. Also
        can optionally return the header (which has things
        like the TES bias).

        Args
        ----
        datafile : str
            The full path to the data to read.
        channel : int or int array or None, optional, default None
            Channels to load.
        nsamp : int or None, optional, default None
            The number of samples to read.
        array_size : int or None, optional, default None
            The size of the output arrays. If 0, then the size will be
            the number of channels in the data file.
        return_header : bool, optional, default False
            Whether to also read in the header and return the header
            data. Returning the full header is slow for large
            files. This overrides return_tes_bias.
        return_tes_bias : bool, optional, default False
            Whether to return the TES bias.
        write_log : bool, optional, default True
            Whether to write outputs to the log file.
        n_max : int, optional, default 2048
            The number of elements to read in before appending the
            datafile. This is just for speed.
        make_freq_mask : bool, optional, default False
            Whether to write a text file with resonator frequencies.
        gcp_mode (bool) : Indicates that the data was written in GCP mode. This
            is the legacy data mode which was depracatetd in Rogue 4.

        Ret:
        ----
        t (float array): The timestamp data
        d (float array): The resonator data in units of radians
        m (int array): The maskfile that maps smurf num to gcp num
        h (dict) : A dictionary with the header information.
        """
        if gcp_mode:
            self.log('Data is in GCP mode.')
            return self.read_stream_data_gcp_save(datafile, channel=channel,
                unwrap=True, downsample=1, nsamp=nsamp)

        # Why were we globbing here?
        #try:
        #    datafile = glob.glob(datafile+'*')[-1]
        #except BaseException:
        #    self.log(f'datafile={datafile}')
        if not os.path.isfile(datafile):
            raise FileNotFoundError(f'datafile={datafile}')

        if write_log:
            self.log(f'Reading {datafile}')

        if channel is not None:
            self.log(f'Only reading channel {channel}')

        # Flag to indicate we are about the read the fist frame from the disk
        # The number of channel will be extracted from the first frame and the
        # data structures will be build based on that
        first_read = True
        with SmurfStreamReader(datafile,
                isRogue=True, metaEnable=True) as file:
            for header, data in file.records():
                if first_read:
                    # Update flag, so that we don't do this code again
                    first_read = False

                    # Read in all used channels by default
                    if channel is None:
                        channel = np.arange(header.number_of_channels)

                    channel = np.ravel(np.asarray(channel))
                    n_chan = len(channel)

                    # Indexes for input channels
                    channel_mask = np.zeros(n_chan, dtype=int)
                    for i, c in enumerate(channel):
                        channel_mask[i] = c

                    #initialize data structure
                    phase=list()
                    for _,_ in enumerate(channel):
                        phase.append(list())
                    for i,_ in enumerate(channel):
                        phase[i].append(data[i])
                    t = [header.timestamp]
                    if return_header or return_tes_bias:
                        tmp_tes_bias = np.array(header.tesBias)
                        tes_bias = np.zeros((0,16))

                    # Get header values if requested
                    if return_header or return_tes_bias:
                        tmp_header_dict = {}
                        header_dict = {}
                        for i, h in enumerate(header._fields):
                            tmp_header_dict[h] = np.array(header[i])
                            header_dict[h] = np.array([],
                                                      dtype=type(header[i]))
                        tmp_header_dict['tes_bias'] = np.array([header.tesBias])


                    # Already loaded 1 element
                    counter = 1
                else:
                    for i in range(n_chan):
                        phase[i].append(data[i])
                    t.append(header.timestamp)

                    if return_header or return_tes_bias:
                        for i, h in enumerate(header._fields):
                            tmp_header_dict[h] = np.append(tmp_header_dict[h],
                                                       header[i])
                        tmp_tes_bias = np.vstack((tmp_tes_bias, header.tesBias))

                    if counter % n_max == n_max - 1:
                        if write_log:
                            self.log(f'{counter+1} elements loaded')

                        if return_header:
                            for k in header_dict.keys():
                                header_dict[k] = np.append(header_dict[k],
                                                           tmp_header_dict[k])
                                tmp_header_dict[k] = \
                                    np.array([],
                                             dtype=type(header_dict[k][0]))
                            print(np.shape(tes_bias), np.shape(tmp_tes_bias))
                            tes_bias = np.vstack((tes_bias, tmp_tes_bias))
                            tmp_tes_bias = np.zeros((0, 16))

                        elif return_tes_bias:
                            tes_bias = np.vstack((tes_bias, tmp_tes_bias))
                            tmp_tes_bias = np.zeros((0, 16))

                    counter += 1

        phase=np.array(phase)
        t=np.array(t)

        if return_header:
            for k in header_dict.keys():
                header_dict[k] = np.append(header_dict[k],
                    tmp_header_dict[k])
            tes_bias = np.vstack((tes_bias, tmp_tes_bias))
            tes_bias = np.transpose(tes_bias)

        elif return_tes_bias:
            tes_bias = np.vstack((tes_bias, tmp_tes_bias))
            tes_bias = np.transpose(tes_bias)

        # rotate and transform to phase
        phase = phase.astype(float) / 2**15 * np.pi

        if np.size(phase) == 0:
            self.log("Only 1 element in datafile. This is often an indication" +
                "that the data was taken in GCP mode. Try running this"+
                " function again with gcp_mode=True")

        # make a mask from mask file
        #  regexp pattern to match any filename which ends in a
        #  . followed by a number, as occurs when MaxFileSize is
        #  nonzero and rogue rolls over files by appending an
        #  increasing number at the end after a .
        extpattern=re.compile('(.+?).dat.([0-9]|[1-9][0-9]+)$')
        extmatch=extpattern.match(datafile)
        if ".dat.part" in datafile:
            mask = self.make_mask_lookup(datafile.split(".dat.part")[0] +
                "_mask.txt")
        elif extmatch is not None:
            mask = self.make_mask_lookup(extmatch[1]+"_mask.txt")
        else:
            mask = self.make_mask_lookup(datafile.replace('.dat', '_mask.txt'),
                                         make_freq_mask=make_freq_mask)

        # If an array_size was defined, resize the phase array
        if array_size is not None:
            phase.resize(array_size, phase.shape[1])

        if return_header:
            header_dict['tes_bias'] = tes_bias
            return t, phase, mask, header_dict
        elif return_tes_bias:
            return t, phase, mask, tes_bias
        else:
            return t, phase, mask

    @set_action()
    def read_stream_data_gcp_save(self, datafile, channel=None,
            unwrap=True, downsample=1, nsamp=None):
        """
        Reads the special data that is designed to be a copy of the GCP data.
        This was the most common data writing mode until the Rogue 4 update.
        Maintining this function for backwards compatibility.

        Args
        ----
        datafile : str
            The full path to the data made by stream_data_on.

        channel : int or list of int or None, optional, default None
            Channels to load.
        unwrap : bool, optional, default True
            Whether to unwrap units of 2pi.
        downsample : int, optional, default 1
            The amount to downsample.
        nsamp : int or None, optional, default None
            The number of samples to read.

        Returns
        -------
        t : numpy.ndarray
            The timestamp data.
        d : numpy.ndarray
            The resonator data in units of phi0.
        m : numpy.ndarray
            The maskfile that maps smurf num to gcp num.
        """
        import struct
        try:
            datafile = glob.glob(datafile+'*')[-1]
        except ValueError:
            print(f'datafile={datafile}')

        self.log(f'Reading {datafile}')

        if channel is not None:
            self.log(f'Only reading channel {channel}')


        keys = ['protocol_version','crate_id','slot_number','number_of_channels',
                'rtm_dac_config0', 'rtm_dac_config1', 'rtm_dac_config2',
                'rtm_dac_config3', 'rtm_dac_config4', 'rtm_dac_config5',
                'flux_ramp_increment','flux_ramp_start', 'rate_since_1Hz',
                'rate_since_TM', 'nanoseconds', 'seconds', 'fixed_rate_marker',
                'sequence_counter', 'tes_relay_config', 'mce_word',
                'user_word0', 'user_word1', 'user_word2'
                ]

        data_keys = [f'data{i}' for i in range(528)]

        keys.extend(data_keys)
        keys_dict = dict(zip(keys, range(len(keys))))

        # Read in all channels by default
        if channel is None:
            n_channels = self.get_number_channels()
            channel = np.arange(n_channels)

        channel = np.ravel(np.asarray(channel))
        n_chan = len(channel)

        # Indices for input channels
        channel_mask = np.zeros(n_chan, dtype=int)
        for i, c in enumerate(channel):
            channel_mask[i] = keys_dict[f'data{c}']

        eval_n_samp = False
        if nsamp is not None:
            eval_n_samp = True

        # Make holder arrays for phase and timestamp
        phase = np.zeros((n_chan,0))
        timestamp2 = np.array([])
        counter = 0
        n = 20000  # Number of elements to load at a time
        tmp_phase = np.zeros((n_chan, n))
        tmp_timestamp2 = np.zeros(n)
        with open(datafile, mode='rb') as file:
            while True:
                chunk = file.read(2240)  # Frame size is 2240
                if not chunk:
                    # If frame is incomplete - meaning end of file
                    phase = np.hstack((phase, tmp_phase[:,:counter%n]))
                    timestamp2 = np.append(timestamp2, tmp_timestamp2[:counter%n])
                    break
                elif eval_n_samp:
                    if counter >= nsamp:
                        phase = np.hstack((phase, tmp_phase[:,:counter%n]))
                        timestamp2 = np.append(timestamp2,
                                               tmp_timestamp2[:counter%n])
                        break
                frame = struct.Struct('3BxI6Q8I5Q528i').unpack(chunk)

                # Extract detector data
                for i, c in enumerate(channel_mask):
                    tmp_phase[i,counter%n] = frame[c]

                # Timestamp data
                tmp_timestamp2[counter%n] = frame[keys_dict['rtm_dac_config5']]

                # Store the data in a useful array and reset tmp arrays
                if counter % n == n - 1 :
                    self.log(f'{counter+1} elements loaded')
                    phase = np.hstack((phase, tmp_phase))
                    timestamp2 = np.append(timestamp2, tmp_timestamp2)
                    tmp_phase = np.zeros((n_chan, n))
                    tmp_timestamp2 = np.zeros(n)
                counter = counter + 1

        phase = np.squeeze(phase)
        phase = phase.astype(float) / 2**15 * np.pi # where is decimal?  Is it in rad?

        rootpath = os.path.dirname(datafile)
        filename = os.path.basename(datafile)
        timestamp = filename.split('.')[0]

        mask = self.make_mask_lookup(os.path.join(rootpath,
                                                  f'{timestamp}_mask.txt'))

        return timestamp2, phase, mask

    @set_action()
    def header_to_tes_bias(self, header, as_volt=True,
                           n_tes_bias=15):
        """
        Takes the SmurfHeader returned from read_stream_data
        and turns it to a TES bias. The header is a 20 field,
        and each DAC is 18 bits signed. So the output of the
        data in the header is (dac_b - dac_a)/2. This function
        also takes care of the factor of 2 in the denominator.

        Args
        ----
        header : dict
            The header dictionary from read_stream_data.  This
            includes all the tes_byte data.

        as_volt : bool, optional, default True
            Whether to return the data as voltage. If False, returns
            as DAC units.
        n_tes_bias : int, optional, default 15
            The number of TES bias pairs.

        Returns
        -------
        bias : numpy.ndarray
            The tes bias data. (dac_b - dac_a) in voltage or DAC units
            depending on the as_volt opt arg.
        """
        # Numbr of total elements
        n_els = len(header['tes_byte_0'])

        # Pre-allocate array
        bias = np.zeros((n_tes_bias, n_els))

        # Iterate over bias groups
        for bias_group in np.arange(n_tes_bias):
            base_byte = int((bias_group*20) / 8)
            base_bit = int((bias_group*20) % 8)
            for i in np.arange(n_els):
                val = 0
                for idx, byte in enumerate(range(base_byte, base_byte+3)):
                    # Cast as type int instead of numpy.int64
                    val += int(header[f'tes_byte_{byte}'][i]) << idx*8

                # https://github.com/slaclab/pysmurf/blob/master/README.SmurfPacket.md
                # Dealing with the 16x20 bit in 10x32 bit words.
                tmp = (val >> base_bit) & 0xFFFFF
                if tmp & 0x80000:
                    tmp |= 0xF00000

                # Cast data into int
                ba = tmp.to_bytes(3, byteorder='little', signed=False)
                bias[bias_group,i] = int.from_bytes(ba, byteorder='little',
                                                    signed=True)

        # Take care of factor of 2 thrown away in writing to the header
        bias *= 2

        # Cast as voltage.
        if as_volt:
            bias *= self._rtm_slow_dac_bit_to_volt

        return bias

    @set_action()
    def make_mask_lookup(self, mask_file, make_freq_mask=False):
        """ Makes an n_band x n_channel array where the elements correspond
        to the smurf_to_mce mask number. In other words, mask[band, channel]
        returns the GCP index in the mask that corresonds to band, channel.

        Args
        ----
        mask_file : str
            The full path the a mask file.
        make_freq_mask : bool, optional, default False
            Whether to write a text file with resonator frequencies.

        Returns
        -------
        mask_lookup : int array
            An array with the GCP numbers.
        """
        # Look for .dat file and replace with mask file
        if ".dat" in mask_file:
            self.log("make_mask_lookup received a .dat file. " +
                "Replacing with mask path.")
            if ".dat.part" in mask_file:
                mask_file = mask_file.split(".dat.part")[0] + "_mask.txt"
            else:
                mask_file = mask_file.replace(".dat", "_mask.txt")

        n_channels = self.get_number_channels()
        mask = np.atleast_1d(np.loadtxt(mask_file))
        bands = np.unique(mask // n_channels).astype(int)
        ret = np.ones((np.max(bands)+1, n_channels), dtype=int) * -1
        if make_freq_mask:
            freq_mask_file = mask_file.replace("_mask.txt", "_freq.txt")
            freq_mask_ret = np.zeros_like(ret).astype(float)
            try:
                freq_mask = np.atleast_1d(np.loadtxt(freq_mask_file))
            except OSError:
                self.log(f'{freq_mask_file} does not exist.')
                make_freq_mask = False

        for gcp_chan, smurf_chan in enumerate(mask):
            b = int(smurf_chan//n_channels)
            ch = int((smurf_chan)%n_channels)
            ret[b,ch] = gcp_chan

            # fill corresponding elements with frequency
            if make_freq_mask:
                freq_mask_ret[b, ch] = freq_mask[gcp_chan]

        # Append freq data if requested
        if make_freq_mask:
            ret = (ret, freq_mask_ret)

        return ret

    @set_action()
    def read_stream_data_daq(self, data_length, bay=0, hw_trigger=False,
            write_log=False):
        """
        Reads the stream data from the DAQ.

        Args
        ----
        data_length : int
            The number of samples to process.
        bay : int, optional, default 0
            The AMC bay number.
        hw_trigger : bool, optional, default False
            Whether to trigger the start of the acquistion with a
            hardware trigger.
        write_log : bool, optional, default False
            Whether to write outputs to log.
        """
        # Ask mitch why this is what it is...
        if bay == 0:
            stream0 = self.epics_root + ":AMCc:Stream0"
            stream1 = self.epics_root + ":AMCc:Stream1"
        else:
            stream0 = self.epics_root + ":AMCc:Stream2"
            stream1 = self.epics_root + ":AMCc:Stream3"

        pvs = [stream0, stream1]
        sg  = SyncGroup(pvs, skip_first=True)

        # trigger PV
        if not hw_trigger:
            self.set_trigger_daq(bay, 1, write_log=write_log)
        else:
            self.set_arm_hw_trigger(bay, 1, write_log=write_log)

        time.sleep(.1)
        sg.wait()

        vals = sg.get_values()

        r0 = vals[pvs[0]]
        r1 = vals[pvs[1]]

        return r0, r1

    @set_action()
    def check_adc_saturation(self, band):
        """
        Reads data directly off the ADC.  Checks for input saturation.

        Args
        -----
        band : int
            Which band.  Assumes adc number is band%4.

        Returns
        -------
        saturated : bool
           True if ADC is saturated, otherwise False.
        """
        adc = self.read_adc_data(band, data_length=2**12, make_plot=False,
                  save_data=False, show_plot=False, save_plot=False)
        adc_max   = int(np.max((adc.real.max(), adc.imag.max())))
        adc_min   = int(np.min((adc.real.min(), adc.imag.min())))
        saturated = ((adc_max > 31000) | (adc_min < -31000))
        self.log(f'ADC{band} max count: {adc_max}')
        self.log(f'ADC{band} min count: {adc_min}')
        if saturated:
            self.log(f'\033[91mADC{band} saturated\033[00m') # color red
        else:
            self.log(f'\033[92mADC{band} not saturated\033[00m') # color green
        return saturated

    @set_action()
    def check_dac_saturation(self, band):
        """
        Reads data directly off the DAC.  Checks for input saturation.

        Args
        ----
        band : int
            Which band.  Assumes dac number is band%4.

        Returns
        -------
        saturated : bool
            Flag if DAC is saturated.
        """
        dac = self.read_dac_data(band, data_length=2**12, make_plot=False,
                  save_data=False, show_plot=False, save_plot=False)
        dac_max   = int(np.max((dac.real.max(), dac.imag.max())))
        dac_min   = int(np.min((dac.real.min(), dac.imag.min())))
        saturated = ((dac_max > 31000) | (dac_min < -31000))
        self.log(f'DAC{band} max count: {dac_max}')
        self.log(f'DAC{band} min count: {dac_min}')
        if saturated:
            self.log(f'\033[91mDAC{band} saturated\033[00m') # color red
        else:
            self.log(f'\033[92mDAC{band} not saturated\033[00m') # color green
        return saturated

    @set_action()
    def read_adc_data(self, band, data_length=2**19,
                      hw_trigger=False, make_plot=False, save_data=True,
                      timestamp=None, show_plot=True, save_plot=True,
                      plot_ylimits=[None,None]):
        """
        Reads data directly off the ADC.

        Args
        ----
        band : int
            Which band.  Assumes adc number is band%4.
        data_length : int, optional, default 2**19
            The number of samples.
        hw_trigger : bool, optional, default False
            Whether to use the hardware trigger. If False, uses an
            internal trigger.
        make_plot : bool, optional, default False
            Whether or not to plot.
        save_data : bool, optional, default True
            Whether or not to save the data in a time stamped file.
        timestamp : int or None, optional, default None
            ctime to timestamp the plot and data with (if saved to
            file).  If None, it gets the time stamp right before
            acquiring data.
        show_plot : bool, optional, default True
            Whether or not to show plots.
        save_plot : bool, optional, default True
            Whether or not to save plot to file.
        plot_ylimits : [float or None, float or None], optional, default [None,None]
            y-axis limit (amplitude) to restrict plotting over.

        Returns
        -------
        dat : int array
            The raw ADC data.
        """
        if timestamp is None:
            timestamp = self.get_timestamp()

        bay=self.band_to_bay(band)
        adc_number=band%4

        self.setup_daq_mux('adc', adc_number, data_length,band=band)

        res = self.read_stream_data_daq(data_length, bay=bay,
            hw_trigger=hw_trigger)
        dat = res[1] + 1.j * res[0]

        if make_plot:
            if show_plot:
                plt.ion()
            else:
                plt.ioff()

            import scipy.signal as signal
            digitizer_frequency_mhz = self.get_digitizer_frequency_mhz()
            f, p_adc = signal.welch(dat, fs=digitizer_frequency_mhz,
                nperseg=data_length/2, return_onesided=False, detrend=False)
            f_plot = f

            idx = np.argsort(f)
            f_plot = f_plot[idx]
            p_adc = p_adc[idx]

            plt.figure(figsize=(9,9))
            ax1 = plt.subplot(211)
            ax1.plot(np.real(dat), label='I')
            ax1.plot(np.imag(dat), label='Q')
            ax1.set_xlabel('Sample number')
            ax1.set_ylabel('Raw counts')
            ax1.set_title(f'{timestamp} Timeseries')
            ax1.legend()
            ax1.set_ylim((-2**15, 2**15))
            ax2 = plt.subplot(212)
            ax2.plot(f_plot, 10*np.log10(p_adc))
            ax2.set_ylabel(f'ADC{band}')
            ax2.set_xlabel('Frequency [MHz]')
            ax2.set_title(f'{timestamp} Spectrum')
            plt.grid(which='both')
            if plot_ylimits[0] is not None:
                plt.ylim(plot_ylimits[0],plt.ylim()[1])
            if plot_ylimits[1] is not None:
                plt.ylim(plt.ylim()[0],plot_ylimits[1])



            if save_plot:
                plot_fn = f'{self.plot_dir}/{timestamp}_adc{band}.png'
                plt.savefig(plot_fn)
                self.pub.register_file(plot_fn, 'adc', plot=True)
                self.log(f'ADC plot saved to {plot_fn}')

        if save_data:
            outfn=os.path.join(self.output_dir,
                               f'{timestamp}_adc{band}')
            self.log(f'Saving raw adc data to {outfn}',
                     self.LOG_USER)

            np.save(outfn, res)
            self.pub.register_file(outfn, 'adc', format='npy')

        return dat

    @set_action()
    def read_dac_data(self, band, data_length=2**19,
                      hw_trigger=False, make_plot=False, save_data=True,
                      timestamp=None, show_plot=True, save_plot=True,
                      plot_ylimits=[None,None]):
        """
        Read the data directly off the DAC.

        Args
        ----
        band : int
            Which band.  Assumes dac number is band%4.
        data_length : int, optional, default 2**19
            The number of samples.
        hw_trigger : bool, optional, default False
            Whether to use the hardware trigger. If False, uses an
            internal trigger.
        make_plot : bool, optional, default False
            Whether or not to plot.
        save_data : bool, optional, default True
            Whether or not to save the data in a time stamped file.
        timestamp : int or None, optional, default None
            ctime to timestamp the plot and data with (if saved to
            file).  If None, in which case it gets the time stamp
            right before acquiring data.
        show_plot : bool, optional, default True
            If make_plot is True, whether or not to show the plot.
        save_plot : bool, optional, default True
            Whether or not to save plot to file.
        plot_ylimits : list of float or list of None, optional, default [None,None]
            2-element list of y-axis limits (amplitude) to restrict
            plotting over.

        Returns
        -------
        dat : int array
            The raw DAC data.
        """
        if timestamp is None:
            timestamp = self.get_timestamp()

        bay=self.band_to_bay(band)
        dac_number=band%4

        self.setup_daq_mux('dac', dac_number, data_length, band=band)

        res = self.read_stream_data_daq(data_length, bay=bay, hw_trigger=hw_trigger)
        dat = res[1] + 1.j * res[0]

        if make_plot:
            if show_plot:
                plt.ion()
            else:
                plt.ioff()

            import scipy.signal as signal
            digitizer_frequency_mhz = self.get_digitizer_frequency_mhz()
            f, p_dac = signal.welch(dat, fs=digitizer_frequency_mhz,
                nperseg=data_length/2, return_onesided=False, detrend=False)
            f_plot = f

            idx = np.argsort(f)
            f_plot = f_plot[idx]
            p_dac = p_dac[idx]

            plt.figure(figsize=(9,9))
            ax1 = plt.subplot(211)
            ax1.plot(np.real(dat), label='I')
            ax1.plot(np.imag(dat), label='Q')
            ax1.set_xlabel('Sample number')
            ax1.set_ylabel('Raw counts')
            ax1.set_title(f'{timestamp} Timeseries')
            ax1.legend()
            ax1.set_ylim((-2**15, 2**15))
            ax2 = plt.subplot(212)
            ax2.plot(f_plot, 10*np.log10(p_dac))
            ax2.set_ylabel(f'DAC{band}')
            ax2.set_xlabel('Frequency [MHz]')
            ax2.set_title(f'{timestamp} Spectrum')
            plt.grid(which='both')
            if plot_ylimits[0] is not None:
                plt.ylim(plot_ylimits[0],plt.ylim()[1])
            if plot_ylimits[1] is not None:
                plt.ylim(plt.ylim()[0],plot_ylimits[1])


            if save_plot:
                plot_fn = f'{self.plot_dir}/{timestamp}_dac{band}.png'
                plt.savefig(plot_fn)
                self.pub.register_file(plot_fn, 'dac', plot=True)
                self.log(f'DAC plot saved to {plot_fn}')

        if save_data:
            outfn = os.path.join(self.output_dir,f'{timestamp}_dac{band}')
            self.log(f'Saving raw dac data to {outfn}', self.LOG_USER)

            np.save(outfn, res)
            self.pub.register_file(outfn, 'dac', format='npy')

        return dat

    @set_action()
    def setup_daq_mux(self, converter, converter_number, data_length,
                      band=0, debug=False, write_log=False):
        """
        Sets up for either ADC or DAC data taking.

        Args
        ----
        converter : str
            Whether it is the ADC or DAC. choices are 'adc', 'dac', or
            'debug'. The last one takes data on a single band.
        converter_number : int
            The ADC or DAC number to take data on.
        data_length : int
            The amount of data to take.
        band : int, optional, default 0
            which band to get data on.
        """

        bay=self.band_to_bay(band)

        if converter.lower() == 'adc':
            daq_mux_channel0 = (converter_number + 1)*2
            daq_mux_channel1 = daq_mux_channel0 + 1
        elif converter.lower() == 'dac':
            daq_mux_channel0 = (converter_number + 1)*2 + 10
            daq_mux_channel1 = daq_mux_channel0 + 1
        else:
            # In dspv3, daq_mux_channel0 and daq_mux_channel1 are now
            # the same for all eight bands.
            daq_mux_channel0 = 22
            daq_mux_channel1 = 23

        # setup buffer size
        self.set_buffer_size(bay, data_length, debug)

        # input mux select
        self.set_input_mux_sel(bay, 0, daq_mux_channel0,
                               write_log=write_log)
        self.set_input_mux_sel(bay, 1, daq_mux_channel1,
                               write_log=write_log)

        # which f,df stream to route to MUX, maybe?
        self.set_debug_select(bay, band%4, write_log=True)

    @set_action()
    def set_buffer_size(self, bay, size, debug=False,
                        write_log=False):
        """
        Sets the buffer size for reading and writing DAQs

        Args
        ----
        size : int
            The buffer size in number of points.
        """
        # Change DAQ data buffer size

        # Change waveform engine buffer size
        self.set_data_buffer_size(bay, size, write_log=True)
        for daq_num in np.arange(2):
            s = self.get_waveform_start_addr(bay, daq_num, convert=True,
                write_log=debug)
            e = s + 4*size
            self.set_waveform_end_addr(bay, daq_num, e, convert=True,
                write_log=debug)
            if debug:
                self.log(f'DAQ number {daq_num}: start {s} - end {e}')

    @set_action()
    def config_cryo_channel(self, band, channel, frequencyMHz, amplitude,
            feedback_enable, eta_phase, eta_mag):
        """
        Set parameters on a single cryo channel

        Args
        ----
        band : int
            The band for the channel.
        channel : int
            Which channel to configure.
        frequencyMHz : float
            The frequency offset from the subband center in MHz.
        amplitude : int
            Amplitude scale to set for the channel (0..15).
        feedback_enable : bool
            Whether to enable feedback for the channel.
        eta_phase : float
            Feedback eta phase, in degrees (-180..180).
        eta_mag : float
            Feedback eta magnitude.
        """

        n_subbands = self.get_number_sub_bands(band)
        digitizer_frequency_mhz = self.get_digitizer_frequency_mhz(band)
        subband_width = digitizer_frequency_mhz / (n_subbands / 2)

        # some checks to make sure we put in values within the correct ranges

        if frequencyMHz > subband_width / 2:
            self.log("frequencyMHz exceeds subband width! setting to top of subband")
            freq = subband_width / 2
        elif frequencyMHz < - subband_width / 2:
            self.log("frequencyMHz below subband width! setting to bottom of subband")
            freq = -subband_width / 2
        else:
            freq = frequencyMHz

        if amplitude > 15:
            self.log("amplitude too high! setting to 15")
            ampl = 15
        elif amplitude < 0:
            self.log("amplitude too low! setting to 0")
            ampl = 0
        else:
            ampl = amplitude

        # get phase within -180..180
        phase = eta_phase
        while phase > 180:
            phase = phase - 360
        while phase < -180:
            phase = phase + 360

        # now set all the PV's
        self.set_center_frequency_mhz_channel(band, channel, freq)
        self.set_amplitude_scale_channel(band, channel, ampl)
        self.set_eta_phase_degree_channel(band, channel, phase)
        self.set_eta_mag_scaled_channel(band, channel, eta_mag)

    @set_action()
    def which_on(self, band):
        """
        Finds all detectors that are on.

        Args
        ----
        band : int
            The band to search.

        Returns
        --------
        int array
            The channels that are on.
        """
        amps = self.get_amplitude_scale_array(band)
        return np.ravel(np.where(amps != 0))

    @set_action()
    def toggle_feedback(self, band, **kwargs):
        """
        Toggles feedbackEnable (->0->1) and lmsEnables1-3 (->0->1) for
        this band.  Only toggles back to 1 if it was 1 when asked to
        toggle, otherwise leaves it zero.

        Args
        ----
        band : int
           The band whose feedback to toggle.
        """

        # current vals?
        old_feedback_enable=self.get_feedback_enable(band)
        old_lms_enable1=self.get_lms_enable1(band)
        old_lms_enable2=self.get_lms_enable2(band)
        old_lms_enable3=self.get_lms_enable3(band)

        self.log(f'Before toggling feedback on band {band}, ' +
            f'feedbackEnable={old_feedback_enable}, ' +
            f'lmsEnable1={old_lms_enable1}, lmsEnable2={old_lms_enable2}, ' +
            f'and lmsEnable3={old_lms_enable3}.', self.LOG_USER)

        # -> 0
        self.log('Setting feedbackEnable=lmsEnable1=lmsEnable2=lmsEnable3=0'+
            ' (in that order).', self.LOG_USER)

        self.set_feedback_enable(band,0)
        self.set_lms_enable1(band,0)
        self.set_lms_enable2(band,0)
        self.set_lms_enable3(band,0)

        # -> 1
        logstr='Set '
        if old_feedback_enable:
            self.set_feedback_enable(band,1)
            logstr+='feedbackEnable='
        if old_lms_enable1:
            self.set_lms_enable1(band,1)
            logstr+='lmsEnable1='
        if old_lms_enable2:
            self.set_lms_enable2(band,1)
            logstr+='lmsEnable2='
        if old_lms_enable3:
            self.set_lms_enable3(band,1)
            logstr+='lmsEnable3='

        logstr+='1 (in that order).'
        self.log(logstr,
                 self.LOG_USER)

    @set_action()
    def band_off(self, band, **kwargs):
        """
        Turns off all tones in a band

        Args
        ----
        band : int
            The band that is to be turned off.
        """
        # Warning ; you might think using the
        # set_amplitude_scale_array function would be fast than this
        # but it is apparently not!
        self.set_amplitude_scales(band, 0, **kwargs)
        n_channels = self.get_number_channels(band)
        self.set_feedback_enable_array(
            band, np.zeros(n_channels, dtype=int), **kwargs)
        self.set_cfg_reg_ena_bit(0, wait_after=.2, **kwargs)

    def channel_off(self, band, channel, **kwargs):
        """
        Turns off the tone for a single channel by setting the amplitude to
        zero and disabling feedback.

        Args
        ----
        band : int
            The band that is to be turned off.
        channel : int
            The channel to turn off.
        """
        self.log(f'Turning off band {band} channel {channel}',
                 self.LOG_USER)
        self.set_amplitude_scale_channel(band, channel, 0, **kwargs)
        self.set_feedback_enable_channel(band, channel, 0, **kwargs)


    def set_feedback_limit_khz(self, band, feedback_limit_khz, **kwargs):
        """
        Sets the feedback limit

        Args
        ----
        band : int
            The band that is to be turned off.
        feedback_limit_khz : float
            The feedback rate in units of kHz.
        """
        digitizer_freq_mhz = self.get_digitizer_frequency_mhz(band)
        n_subband = self.get_number_sub_bands(band)

        subband_bandwidth = 2 * digitizer_freq_mhz / n_subband
        desired_feedback_limit_mhz = feedback_limit_khz/1000.

        if desired_feedback_limit_mhz > subband_bandwidth/2:
            desired_feedback_limit_mhz = subband_bandwidth/2

        desired_feedback_limit_dec = np.floor(desired_feedback_limit_mhz/
            (subband_bandwidth/2**16.))

        self.set_feedback_limit(band, desired_feedback_limit_dec, **kwargs)

    # if no guidance given, tries to reset both
    def recover_jesd(self, bay, recover_jesd_rx=True, recover_jesd_tx=True):
        """
        Attempts to recover the JESDs

        Args
        ----
        bay : int
            The AMC bay to recover
        recover_jesd_rx : bool, optional, default True
            Whether to attempt to recover the JESD RX
        recover_jesd_tx : bool, optional, default True
            Whether to attempt to recover the JESD TX
        """
        if recover_jesd_rx:
            #1. Toggle JesdRx:Enable 0x3F3 -> 0x0 -> 0x3F3
            self.set_jesd_rx_enable(bay, 0x0)
            self.set_jesd_rx_enable(bay, 0x3F3)

        if recover_jesd_tx:
            #1. Toggle JesdTx:Enable 0x3CF -> 0x0 -> 0x3CF
            self.set_jesd_tx_enable(bay, 0x0)
            self.set_jesd_tx_enable(bay, 0x3CF)

            #2. Toggle AMCcc:FpgaTopLevel:AppTop:AppCore:MicrowaveMuxCore[0]:
            # DAC[0]:JesdRstN 0x1 -> 0x0 -> 0x1
            self.set_jesd_reset_n(bay, 0, 0x0)
            self.set_jesd_reset_n(bay, 0, 0x1)

            #3. Toggle AMCcc:FpgaTopLevel:AppTop:AppCore:MicrowaveMuxCore[0]:
            # DAC[1]:JesdRstN 0x1 -> 0x0 -> 0x1
            self.set_jesd_reset_n(bay, 1, 0x0)
            self.set_jesd_reset_n(bay, 1, 0x1)

        # probably overkill...shouldn't call this function if you're not going
        # to do anything
        if (recover_jesd_rx or recover_jesd_tx):
            # powers up the SYSREF which is required to sync fpga and
            # adc/dac jesd
            self.run_pwr_up_sys_ref(bay)

        # check if Jesds recovered - enable printout
        (jesd_tx_ok, jesd_rx_ok, jesd_status) = self.check_jesd(
            bay, silent_if_valid=False)

        # raise exception if failed to recover
        if (jesd_rx_ok and jesd_tx_ok):
            self.log('Recovered Jesd.', self.LOG_USER)
        else:
            which_jesd_down='Jesd Rx and Tx are both down'
            if (jesd_rx_ok or jesd_tx_ok):
                which_jesd_down = ('Jesd Rx is down' if jesd_tx_ok else 'Jesd Tx is down')
            self.log('Failed to recover Jesds ...', self.LOG_ERROR)
            raise ValueError(which_jesd_down)


    def jesd_decorator(decorated):
        def jesd_decorator_function(self):
            # check JESDs
            (jesd_tx_ok0, jesd_rx_ok0, jesd_status) = self.check_jesd(silent_if_valid=True)

            # if either JESD is down, try to fix
            if not (jesd_rx_ok0 and jesd_tx_ok0):
                which_jesd_down0='Jesd Rx and Tx are both down'
                if (jesd_rx_ok0 or jesd_tx_ok0):
                    which_jesd_down0 = ('Jesd Rx is down' if
                        jesd_tx_ok0 else 'Jesd Tx is down')

                self.log(f'{which_jesd_down0} ... will attempt to recover.',
                         self.LOG_ERROR)

                # attempt to recover ; if it fails it will assert
                self.recover_jesd(recover_jesd_rx=(not jesd_rx_ok0),
                    recover_jesd_tx=(not jesd_tx_ok0))

                # rely on recover to assert if it failed
                self.log('Successfully recovered Jesd but may need to redo' +
                    ' some setup ... rerun command at your own risk.',
                    self.LOG_USER)

            # don't continue running the desired command by default.
            # just because Jesds are back doesn't mean we're in a sane
            # state.  User may need to relock/etc.
            if (jesd_rx_ok0 and jesd_tx_ok0):
                decorated()

        return jesd_decorator_function

    def check_jesd(self, bay, silent_if_valid=False,
                   max_timeout_sec=60, get_timeout_sec=5):
        """Checks JESD status for requested bay.

        Queries the Jesd tx and rx and compares the data_valid and
        enable bits.  If newer Rogue `AppTop.JesdHealth` method is
        available, returns its status.

        Args
        ----
        bay : int
            Which bay (0 or 1).
        silent_if_valid : bool, optional, default False
            If True, does not print anything if things are working.
        max_timeout_sec : float, optional, default 60.0
            Seconds to wait for JESD health check to complete before
            giving up.  Passed to
            :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.set_check_jesd`
        caget_timeout_sec : float, optional, default 5.0
            Seconds to wait for each poll of the JESD health check
            status register (see
            :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.get_jesd_status`).
            Passed to
            :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.set_check_jesd`

        Returns
        -------
        jesd_tx : bool
            True if DataValid = Enable for this bay's JesdTx,
            otherwise False.
        jesd_rx : bool
            True if DataValid = Enable for this bay's JesdRx,
            otherwise False.
        jesd_status : str or None
            System (all bays) JESD health status (see
            :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.get_jesd_status`
            for a description of the possible statuses).

        See Also
        --------
        :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.set_check_jesd` :
              Gets the status of the Rogue `AppTop.JesdHealth`
              method.

        """
        # JESD Tx
        jesd_tx_enable = self.get_jesd_tx_enable(bay)
        jesd_tx_valid = self.get_jesd_tx_data_valid(bay)
        jesd_tx_ok = (jesd_tx_enable==jesd_tx_valid)
        if not jesd_tx_ok:
            self.log("JESD Tx DOWN", self.LOG_ERROR)
        else:
            if not silent_if_valid:
                self.log("JESD Tx Okay", self.LOG_USER)

        # JESD Rx
        jesd_rx_enable = self.get_jesd_rx_enable(bay)
        jesd_rx_valid = self.get_jesd_rx_data_valid(bay)
        jesd_rx_ok = (jesd_rx_enable==jesd_rx_valid)
        if not jesd_rx_ok:
            self.log("JESD Rx DOWN", self.LOG_ERROR)
        else:
            if not silent_if_valid:
                self.log("JESD Rx Okay", self.LOG_USER)

        # New checks introduced
        status = self.set_check_jesd()

        return (jesd_tx_ok, jesd_rx_ok, status)

    def get_fpga_status(self):
        """
        Loads FPGA status checks if JESD is ok.

        Returns
        -------
        ret : dict
            A dictionary containing uptime, fpga_version, git_hash,
            build_stamp, jesd_tx_enable, and jesd_tx_valid
        """
        uptime = self.get_fpga_uptime()
        fpga_version = self.get_fpga_version()
        git_hash = self.get_fpga_git_hash()
        build_stamp = self.get_fpga_build_stamp()

        git_hash = ''.join([chr(y) for y in git_hash]) # convert from int to ascii
        build_stamp = ''.join([chr(y) for y in build_stamp])

        self.log("Build stamp: " + str(build_stamp), self.LOG_USER)
        self.log("FPGA version: Ox" + str(fpga_version), self.LOG_USER)
        self.log("FPGA uptime: " + str(uptime), self.LOG_USER)

        jesd_tx_enable = self.get_jesd_tx_enable()
        jesd_tx_valid = self.get_jesd_tx_data_valid()
        if jesd_tx_enable != jesd_tx_valid:
            self.log("JESD Tx DOWN", self.LOG_USER)
        else:
            self.log("JESD Tx Okay", self.LOG_USER)

        jesd_rx_enable = self.get_jesd_rx_enable()
        jesd_rx_valid = self.get_jesd_rx_data_valid()
        if jesd_rx_enable != jesd_rx_valid:
            self.log("JESD Rx DOWN", self.LOG_USER)
        else:
            self.log("JESD Rx Okay", self.LOG_USER)

        # dict containing all values
        ret = {
            'uptime' : uptime,
            'fpga_version' : fpga_version,
            'git_hash' : git_hash,
            'build_stamp' : build_stamp,
            'jesd_tx_enable' : jesd_tx_enable,
            'jesd_tx_valid' : jesd_tx_valid,
            'jesd_rx_enable': jesd_rx_enable,
            'jesd_rx_valid' : jesd_rx_valid,
        }

        return ret

    def which_bays(self):
        r"""Which carrier AMC bays are enabled.

        Returns which AMC bays were enabled on pysmurf server startup.
        Each SMuRF carrier has two AMC bays, indexed by an integer,
        either 0 or 1.  If looking at an installed carrier from the
        front of a crate, bay 0 is on the right and bay 1 is on the
        left.

        A bay is enabled if the `--disable-bay#` argument is not
        provided as a startup argument to the pysmurf server where #
        is the bay number, either 0 or 1.  The pysmurf server startup
        arguments are returned by the
        :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.get_smurf_startup_args`
        routine.

        Returns
        -------
        bays : list of int
            Which bays were enabled on pysmurf server startup.
        """
        if hasattr(self, '_cached_enabled_bays'):
            return self._cached_enabled_bays

        # New method of getting enabled bays:
        enabled_bays = self.get_enabled_bays()

        if enabled_bays is None:  # Then new rogue var doesn't exist
            # Old method of getting enabled bays
            smurf_startup_args = self.get_smurf_startup_args()

            # Bays are enabled unless --disable-bay{bay} is provided to
            # the pysmurf server on startup.
            enabled_bays = []
            for bay in [0, 1]:
                if f'--disable-bay{bay}' not in smurf_startup_args:
                    enabled_bays.append(bay)

        self._cached_enabled_bays = enabled_bays
        return enabled_bays

    def which_bands(self):
        """Which bands the carrier firmware was built for.

        Returns
        -------
        bands : list of int
            Which bands the carrier firmware was built for.
        """
        build_dsp_g=self.get_build_dsp_g()
        bands=[b for b,x in enumerate(bin(build_dsp_g)[2:]) if x=='1']
        return bands

    def freq_to_subband(self, band, freq):
        """
        Look up subband number of a channel frequency, and its subband
        frequency offset.

        Args
        ----
        band : float
            The band to place the resonator.
        freq : float
            Frequency in MHz.

        Returns
        -------
        subband_no : int
            Subband (0..n_subbands-1) of the frequency within the band.
        offset : float
            Offset from subband center.

        """
        subbands, subband_centers = self.get_subband_centers(band,
            as_offset=False)

        df = np.abs(freq - subband_centers)
        idx = np.ravel(np.where(df == np.min(df)))[0]

        subband_no = subbands[idx]
        offset = freq - subband_centers[idx]

        return subband_no, offset

    def channel_to_freq(self, band, channel=None, yml=None):
        """
        Gives the frequency of the channel.

        Args
        ----
        band : int
            The band the channel is in.
        channel : int, None, or array, optional, default None
            If None, will return the channel freqs of all enabled channels

        Returns
        -------
        freq : float
            The channel frequency in MHz or an array of values if
            channel is None. In the array format, the freq list is
            aligned with self.which_on(band).
        """
        if band is None and channel is None:
            return None

        band_center_mhz = self.get_band_center_mhz(band)
        subband_offset = self.get_tone_frequency_offset_mhz(band)
        channel_offset = self.get_center_frequency_array(band)
        channel_freqs = band_center_mhz + subband_offset + channel_offset

        if channel is None:
            return channel_freqs[self.which_on(band)]
        else:
            return channel_freqs[channel]

    def get_channel_order(self, band=None, channel_orderfile=None):
        """ produces order of channels from a user-supplied input file

        Args
        ----
        band : int or None, optional, default None
            Which band.  If None, assumes all bands have the same
            number of channels, and pulls the number of channels from
            the first band in the list of bands specified in the
            experiment.cfg.
        channelorderfile : str or None, optional, default None
            Path to a file that contains one channel per line.

        Returns
        -------
        channel_order : int array
            An array of channel orders.
        """

        if band is None:
            # assume all bands have the same channel order, and pull
            # the channel frequency ordering from the first band in
            # the list of bands specified in experiment.cfg.
            bands = self._bands
            band = bands[0]

        tone_freq_offset = self.get_tone_frequency_offset_mhz(band)
        freqs = np.sort(np.unique(tone_freq_offset))

        n_subbands = self.get_number_sub_bands(band)
        n_channels = self.get_number_channels(band)

        n_chanpersubband = int(n_channels / n_subbands)

        channel_order = np.zeros(len(tone_freq_offset), dtype=int)
        for i, f in enumerate(freqs):
            channel_order[n_chanpersubband*i:n_chanpersubband*(i+1)] = \
                np.ravel(np.where(tone_freq_offset == f))

        return channel_order

    def get_processed_channels(self, channel_orderfile=None):
        """
        take_debug_data, which is called by many functions including
        tracking_setup only returns data for the processed
        channels. Therefore every channel is not returned.

        Args
        ----
        channelorderfile : str or None, optional, default None
            Path to a file that contains one channel per line.
        """
        n_proc = self.get_number_processed_channels()
        n_chan = self.get_number_channels()
        n_cut = (n_chan - n_proc)//2
        return np.sort(self.get_channel_order(
            channel_orderfile=channel_orderfile)[n_cut:-n_cut])

    def get_subband_from_channel(self, band, channel, channelorderfile=None,
            yml=None):
        """Returns subband number given a channel number

        Args
        ----
        band : int
            Which band we're working in.
        channel : int
            Ranges 0..(n_channels-1), cryo channel number.
        channelorderfile : str or None, optional, default None
            Path to file containing order of channels.

        Returns
        -------
        subband : int
            The subband the channel lives in.
        """

        n_subbands = self.get_number_sub_bands(band, yml=yml)
        n_channels = self.get_number_channels(band, yml=yml)

        n_chanpersubband = n_channels / n_subbands

        if channel > n_channels:
            raise ValueError('channel number exceeds number of channels')

        if channel < 0:
            raise ValueError('channel number is less than zero!')

        chanOrder = self.get_channel_order(band,channelorderfile)
        idx = np.where(chanOrder == channel)[0]

        subband = idx // n_chanpersubband

        return int(subband)

    def get_subband_centers(self, band, as_offset=True, hardcode=False,
            yml=None):
        """ returns frequency in MHz of subband centers

        Args
        ----
        band : int
            Which band.
        as_offset : bool, optional, default True
            Whether to return as offset from band center.
        """

        if hardcode:
            #bandCenterMHz = 3.75 + 0.5*(band + 1)
            digitizer_frequency_mhz = 614.4
            n_subbands = 128
        else:
            digitizer_frequency_mhz = self.get_digitizer_frequency_mhz(band,
                yml=yml)
            n_subbands = self.get_number_sub_bands(band, yml=yml)

        subband_width_MHz = 2 * digitizer_frequency_mhz / n_subbands

        subbands = list(range(n_subbands))
        subband_centers = (np.arange(1, n_subbands + 1) - n_subbands/2) * \
            subband_width_MHz/2

        if not as_offset:
            subband_centers += self.get_band_center_mhz(band, yml=yml)

        return subbands, subband_centers

    def get_channels_in_subband(self, band, subband, channelorderfile=None):
        """
        Returns channels in subband

        Args
        ----
        band : int
            Which band.
        subband : int
            Subband number, ranges from 0..127.
        channelorderfile : str or None, optional, default None
            Path to file specifying channel order.

        Returns
        -------
        subband_chans : int array
            The channels in the subband.
        """

        n_subbands = self.get_number_sub_bands(band)
        n_channels = self.get_number_channels(band)
        n_chanpersubband = int(n_channels / n_subbands)

        if subband > n_subbands:
            raise ValueError("subband requested exceeds number of subbands")

        if subband < 0:
            raise ValueError("requested subband less than zero")

        chanOrder = self.get_channel_order(band,channelorderfile)
        subband_chans = chanOrder[subband * n_chanpersubband : subband *
            n_chanpersubband + n_chanpersubband]

        return subband_chans

    def iq_to_phase(self, i, q):
        """
        Changes IQ to phase
        """
        return np.unwrap(np.arctan2(q, i))


    def hex_string_to_int(self, s):
        """
        Converts hex string, which is an array of characters, into an int.

        Args
        ----
        s : character array
            An array of chars to be turned into a single int.

        Returns
        -------
        i : numpy.int64
           The 64 bit int.
        """
        return np.int64(int(''.join([chr(x) for x in s]),0))


    def int_to_hex_string(self, i):
        """
        Converts an int into a string of characters.

        Args
        ----
        i : int
            A 64 bit int to convert into hex.

        Returns
        -------
        s : char array
            A character array representing the int.
        """
        # Must be array length 300
        s = np.zeros(300, dtype=int)
        i_hex = hex(i)
        for j in np.arange(len(i_hex)):
            s[j] = ord(i_hex[j])

        return s


    def set_tes_bias_bipolar(self, bias_group, volt, do_enable=False,
            flip_polarity=False, **kwargs):
        """
        Set an individual TES bias group to the specified voltage, in
        volts.  Asserts if the requested bias group is not defined in
        the pysmurf configuration file.  The positive DAC in the bias
        group is set to +volt/2, while the negative DAC in the bias
        group is set to -volt/2.

        Args
        ----
        bias_group : int
           The bias group.
        volt : float
           The TES bias to command in volts.
        do_enable : bool, optional, default True
           Sets the enable bit. Only must be done once.
        flip_polarity : bool, optional, default False
           Sets the voltage to volt*-1.
        """

        # Make sure the requested bias group is in the list of defined
        # bias groups.
        bias_groups = self.bias_group_to_pair[:,0]
        assert (bias_group in bias_groups),\
            f'Bias group {bias_group} is not defined (available bias '+\
            f' groups are {bias_groups}).  Doing nothing!'

        bias_order = bias_groups
        dac_positives = self.bias_group_to_pair[:,1]
        dac_negatives = self.bias_group_to_pair[:,2]

        dac_idx = np.ravel(np.where(bias_order == bias_group))

        dac_positive = dac_positives[dac_idx][0]
        dac_negative = dac_negatives[dac_idx][0]

        volts_pos = volt / 2
        volts_neg = - volt / 2

        if flip_polarity:
            volts_pos *= -1
            volts_neg *= -1

        if do_enable:
            self.set_rtm_slow_dac_enable(dac_positive, 2, **kwargs)
            self.set_rtm_slow_dac_enable(dac_negative, 2, **kwargs)

        self.set_rtm_slow_dac_volt(dac_positive, volts_pos, **kwargs)
        self.set_rtm_slow_dac_volt(dac_negative, volts_neg, **kwargs)

    def set_tes_bias_bipolar_array(self, bias_group_volt_array, do_enable=False,
                                   **kwargs):
        """
        Set TES bipolar values for all DACs at once.  Set using a
        pyrogue array write, so should be much more efficient than
        setting each TES bias one at a time (a single register
        transaction vs. many).  Only DACs assigned to TES bias groups
        are touched by this function.  The enable status and output
        voltage of all DACs not assigned to a TES bias group are
        maintained.

        Args
        ----
        bias_group_volt_array : float array
            The TES bias to command in voltage for each bipolar TES
            bias group. Should be (n_bias_groups,).
        do_enable : bool, optional, default True
            Set the enable bit for both DACs for every TES bias group.
        """

        n_bias_groups = self._n_bias_groups

        # in this function we're only touching the DACs defined in TES
        # bias groups.  Need to make sure we carry along the setting
        # and enable of any DACs that are being used for something
        # else.
        dac_enable_array = self.get_rtm_slow_dac_enable_array()
        dac_volt_array = self.get_rtm_slow_dac_volt_array()

        if len(bias_group_volt_array) != n_bias_groups:
            self.log("Received the wrong number of biases. Expected " +
                     f"an array of n_bias_groups={n_bias_groups} voltages",
                     self.LOG_ERROR)
        else:
            for bg in np.arange(n_bias_groups):
                bias_order = self.bias_group_to_pair[:,0]
                dac_positives = self.bias_group_to_pair[:,1]
                dac_negatives = self.bias_group_to_pair[:,2]

                bias_group_idx = np.ravel(np.where(bias_order == bg))

                dac_positive = dac_positives[bias_group_idx][0] - 1 # freakin Mitch
                dac_negative = dac_negatives[bias_group_idx][0] - 1 # 1 vs 0 indexing

                volts_pos = bias_group_volt_array[bg] / 2
                volts_neg = - bias_group_volt_array[bg] / 2

                if do_enable:
                    dac_enable_array[dac_positive] = 2
                    dac_enable_array[dac_negative] = 2

                dac_volt_array[dac_positive] = volts_pos
                dac_volt_array[dac_negative] = volts_neg

            if do_enable:
                self.set_rtm_slow_dac_enable_array(dac_enable_array, **kwargs)

            self.set_rtm_slow_dac_volt_array(dac_volt_array, **kwargs)


    def set_tes_bias_off(self, **kwargs):
        """
        Turns off all of the DACs assigned to a TES bias group in the
        pysmurf configuration file.
        """
        self.set_tes_bias_bipolar_array(np.zeros(self._n_bias_groups), **kwargs)

    def get_tes_bias_bipolar(self, bias_group, return_raw=False, **kwargs):
        """
        Returns the bias voltage in units of Volts for the requested
        TES bias group.

        Args
        ----
        bias_group : int
            The number of the bias group.  Asserts if bias_group
            requested is not defined in the pysmurf configuration
            file.
        return_raw : bool, optional, default False
            If True, returns pos and neg terminal values.

        Returns
        -------
        val : float
            The bipolar output TES bias voltage for the requested bias
            group.  If return_raw=True, then returns a two element
            float array containing the output voltages of the two DACs
            assigned to the requested TES bias group.
        """
        # Make sure the requested bias group is in the list of defined
        # bias groups.
        bias_groups = self.bias_group_to_pair[:,0]
        assert (bias_group in bias_groups),\
            f'Bias group {bias_group} is not defined (available bias groups are {bias_groups}).  Doing nothing!'

        bias_order = bias_groups
        dac_positives = self.bias_group_to_pair[:,1]
        dac_negatives = self.bias_group_to_pair[:,2]

        dac_idx = np.ravel(np.where(bias_order == bias_group))
        dac_positive = dac_positives[dac_idx][0]-1
        dac_negative = dac_negatives[dac_idx][0]-1

        volt_array = self.get_rtm_slow_dac_volt_array(**kwargs)
        volts_pos = volt_array[dac_positive]
        volts_neg = volt_array[dac_negative]

        if return_raw:
            return volts_pos, volts_neg
        else:
            return volts_pos - volts_neg


    def get_tes_bias_bipolar_array(self, return_raw=False, **kwargs):
        """
        Returns array of bias voltages per bias group in units of volts.

        Args
        ----
        return_raw : bool, optional, default False
            If True, returns +/- terminal vals as separate arrays
            (pos, then negative)

        Returns
        -------
        bias_vals : float, array
            If return_raw is False, returns an array of size n_bias_group
            with the TES bias in units of volts. If return_raw is True,
            returns two arrays of size n_bias_group. The first is the positive
            bias DAC and the second is the negative bias DAC.
        """

        bias_order = self.bias_group_to_pair[:,0]
        dac_positives = self.bias_group_to_pair[:,1]
        dac_negatives = self.bias_group_to_pair[:,2]

        n_bias_groups = self._n_bias_groups

        bias_vals_pos = np.zeros((n_bias_groups,))
        bias_vals_neg = np.zeros((n_bias_groups,))

        volts_array = self.get_rtm_slow_dac_volt_array(**kwargs)

        for idx in np.arange(n_bias_groups):
            dac_idx = np.ravel(np.where(bias_order == idx))
            dac_positive = dac_positives[dac_idx][0] - 1
            dac_negative = dac_negatives[dac_idx][0] - 1

            bias_vals_pos[idx] = volts_array[dac_positive]
            bias_vals_neg[idx] = volts_array[dac_negative]

        if return_raw:
            return bias_vals_pos, bias_vals_neg
        else:
            return bias_vals_pos - bias_vals_neg

    def overbias_tes(self, bias_group, overbias_voltage=19.9, overbias_wait=1.,
                     tes_bias=19.9, cool_wait=20., high_current_mode=False,
                     flip_polarity=False, actually_overbias=True):
        """
        Overbiases requested bias group at overbias_voltage in high current mode
        for overbias_wait seconds. If high_current_mode=False,
        returns to low current mode, after which it biases the TESs at
        tes_bias.  Then waits cool_wait seconds before returning
        control.

        Args
        ----
        bias_group : int
            The bias group to overbias.  Asserts if not a valid bias
            group.
        overbias_voltage : float, optional, default 19.9
            The value of the TES bias in the high current mode.
        overbias_wait : float, optional, default 1.0
            The time to stay in high current mode in seconds.
        tes_bias : float, optional, default 19.9
            The value of the TES bias when put back in low current
            mode.
        cool_wait : float, optional, default 20.0
            The time to wait after setting the TES bias for transients
            to die off.
        high_current_mode : bool, optional, default False
            Whether to keep the TES bias in high current mode after
            the kick.
        flip_polarity : bool, optional, default False
            Whether to flip the TES bias bipolar DAC polarity.
        actually_overbias : bool, optional, default True
            Whether to actaully do the overbias.
        """
        bias_groups = self.bias_group_to_pair[:,0]
        assert (bias_group in bias_groups),\
            f'Bias group {bias_group} is not defined (available bias groups are {bias_groups}).  Doing nothing!'

        if actually_overbias:
            # drive high current through the TES to attempt to drive normal
            self.set_tes_bias_bipolar(bias_group, overbias_voltage,
                flip_polarity=flip_polarity)
            time.sleep(.1)

            self.set_tes_bias_high_current(bias_group)
            self.log('Driving high current through TES. ' +
                     f'Waiting {overbias_wait}', self.LOG_USER)

            time.sleep(overbias_wait)

        if not high_current_mode:
            self.set_tes_bias_low_current(bias_group)
            time.sleep(.1)

        self.set_tes_bias_bipolar(bias_group, tes_bias,
            flip_polarity=flip_polarity)
        self.log(f'Waiting {cool_wait:1.1f} seconds to cool', self.LOG_USER)
        time.sleep(cool_wait)

        self.log('Done waiting.', self.LOG_USER)


    def overbias_tes_all(self, bias_groups=None, overbias_voltage=19.9,
            overbias_wait=1.0, tes_bias=19.9, cool_wait=20.,
            high_current_mode=False, actually_overbias=True):
        """
        Overbiases all requested bias groups (specified by the
        bias_groups array) at overbias_voltage in high current mode
        for overbias_wait seconds.  If high_current_mode=False,
        returns to low current mode, after which it biases the TESs at
        tes_bias.  Then waits cool_wait seconds before returning
        control.

        Args
        ----
        bias_groups : array or None, optional, default None
            Which bias groups to overbias. defaults to all_groups.
            Asserts if any of the bias groups listed is not a defined
            bias group.
        overbias_voltage : float, optional, default 19.9
            The value of the TES bias in the high current mode.
        overbias_wait : float, optional, default 1.0
            The time to stay in high current mode in seconds.
        tes_bias : float, optional, default 19.9
            The value of the TES bias when put back in low current
            mode.
        cool_wait : float, optional, default 20.0
            The time to wait after setting the TES bias for transients
            to die off.
        high_current_mode : bool, optional, default False
            Whether to keep the TES bias in high current mode after
            the kick.
        actually_overbias : bool, optional, default True
            Whether to actaully do the overbias.
        """
        # drive high current through the TES to attempt to drive normal
        if bias_groups is None:
            bias_groups = self._all_groups
        else:
            # assert requires array
            bias_groups = np.atleast_1d(bias_groups)

        valid_bias_groups = self.bias_group_to_pair[:,0]

        assert (all(bg in valid_bias_groups for bg in bias_groups)),\
            'Some of the bias groups requested are not valid '+\
            f'(available bias groups are {valid_bias_groups}).  Doing nothing!'

        # Set the overbias voltage
        if actually_overbias:
            voltage_overbias_array = self.get_tes_bias_bipolar_array()
            voltage_overbias_array[bias_groups] = overbias_voltage
            self.set_tes_bias_bipolar_array(voltage_overbias_array)

            self.log(f'Driving {overbias_voltage} V in high current mode '+
                f'through bias groups {bias_groups}. ' +
                f'Waiting {overbias_wait}', self.LOG_USER)

            # Set high current mode
            self.set_tes_bias_high_current(bias_groups)

            time.sleep(overbias_wait)

        # Set to low current mode
        if not high_current_mode:
            self.log('setting to low current')
            self.set_tes_bias_low_current(bias_groups)

        # Set TES bias
        voltage_bias_array = self.get_tes_bias_bipolar_array()
        voltage_bias_array[bias_groups] = tes_bias
        self.set_tes_bias_bipolar_array(voltage_bias_array)

        # Cool wait
        self.log(f'Waiting {cool_wait:3.2f} seconds to cool',
                 self.LOG_USER)
        time.sleep(cool_wait)
        self.log('Done waiting.', self.LOG_USER)


    def set_tes_bias_high_current(self, bias_group, write_log=False):
        """
        Sets all bias groups to high current mode. Note that the bias group
        number is not the same as the relay number. It also does not matter,
        because Joe's code secretly flips all the relays when you flip one.

        Args
        ----
        bias_group : int
            The bias group(s) to set to high current mode.
        """
        old_relay = self.get_cryo_card_relays()
        old_relay = self.get_cryo_card_relays()  # querey twice to ensure update
        new_relay = np.copy(old_relay)
        if write_log:
            self.log(f'Old relay {bin(old_relay)}')

        n_bias_groups = self._n_bias_groups
        bias_group = np.ravel(np.array(bias_group))
        for bg in bias_group:
            if bg < n_bias_groups:
                r = np.ravel(self._pic_to_bias_group[
                    np.where(self._pic_to_bias_group[:,1]==bg)])[0]
            else:
                r = bg
            new_relay = (1 << r) | new_relay
        if write_log:
            self.log(f'New relay {bin(new_relay)}')
        self.set_cryo_card_relays(new_relay, write_log=write_log)
        self.get_cryo_card_relays()

    def get_tes_bias_high_current(self, bias_group):
        """
        Returns 1 if requested bias_group is in high_current_mode and 0
        otherwise.

        Args
        ----
        bias_group : int
            The bias group to query
        """
        relay = self.get_cryo_card_relays()
        relay = self.get_cryo_card_relays()  # querey twice to ensure update

        if bias_group >= self._n_bias_groups:
            raise ValueError("Biasgroup must be between 0 and {self._n_bias_groups}")

        r = np.ravel(self._pic_to_bias_group[np.where(
            self._pic_to_bias_group[:,1]==bias_group)])[0]
        return (relay>>r) & 1

    def set_tes_bias_low_current(self, bias_group, write_log=False):
        """
        Sets all bias groups to low current mode. Note that the bias group
        number is not the same as the relay number. It also does not matter,
        because Joe's code secretly flips all the relays when you flip one

        Args
        ----
        bias_group : int
            The bias group to set to low current mode.
        """
        old_relay = self.get_cryo_card_relays()
        old_relay = self.get_cryo_card_relays()  # querey twice to ensure update
        new_relay = np.copy(old_relay)

        n_bias_groups = self._n_bias_groups
        bias_group = np.ravel(np.array(bias_group))
        if write_log:
            self.log(f'Old relay {bin(old_relay)}')
        for bg in bias_group:
            if bg < n_bias_groups:
                r = np.ravel(self._pic_to_bias_group[np.where(
                    self._pic_to_bias_group[:,1]==bg)])[0]
            else:
                r = bg
            if old_relay & 1 << r != 0:
                new_relay = new_relay & ~(1 << r)
        if write_log:
            self.log(f'New relay {bin(new_relay)}')
        self.set_cryo_card_relays(new_relay, write_log=write_log)
        self.get_cryo_card_relays()


    def set_mode_dc(self, write_log=False):
        """
        Sets flux ramp to DC coupling

        Args
        ----
        write_log : bool, optional, default False
            Whether to write outputs to log.
        """
        # The 16th bit (0 indexed) is the AC/DC coupling
        # self.set_tes_bias_high_current(16)
        r = 16

        old_relay = self.get_cryo_card_relays()
        # query twice to ensure update
        old_relay = self.get_cryo_card_relays()
        self.log(f'Old relay {bin(old_relay)}')

        new_relay = np.copy(old_relay)
        new_relay = (1 << r) | new_relay
        self.log(f'New relay {bin(new_relay)}')
        self.set_cryo_card_relays(new_relay, write_log=write_log)
        self.get_cryo_card_relays()


    def set_mode_ac(self, write_log=False):
        """
        Sets flux ramp to AC coupling

        Args
        ----
        write_log : bool, optional, default False
            Whether to write outputs to log.
        """
        # The 16th bit (0 indexed) is the AC/DC coupling
        # self.set_tes_bias_low_current(16)
        old_relay = self.get_cryo_card_relays()
        old_relay = self.get_cryo_card_relays()  # querey twice to ensure update
        new_relay = np.copy(old_relay)

        r = 16
        if old_relay & 1 << r != 0:
            new_relay = new_relay & ~(1 << r)

        self.log(f'New relay {bin(new_relay)}')
        self.set_cryo_card_relays(new_relay)
        self.get_cryo_card_relays()

    def att_to_band(self, att):
        """500 MHz band associated with this attenuator number.

        Args
        ----
        att : int
            Attenuator number.

        Returns
        -------
        band : int
            The 500 MHz band associated with the attenuator.
        """
        return int(
            self._attenuator['band'][np.ravel(
                np.where(self._attenuator['att']==att))[0]])

    def band_to_att(self, band):
        """Attenuator number associated with this 500 MHz band.

        Args
        ----
        band : int
            500 MHz band number.

        Returns
        -------
        att : int
            The attenuator number associated with the band.
        """
        # for now, mod 4 ; assumes the band <-> att correspondence is
        # the same for the LB and HB AMCs.
        band=band%4
        return int(
            self._attenuator['att'][np.ravel(
                np.where(self._attenuator['band']==band))[0]])

    def flux_ramp_rate_to_PV(self, val):
        """
        Convert between the desired flux ramp reset rate and the PV number
        for the timing triggers.

        Hardcoded somewhere that we can't access; this is just a lookup table
        Allowed reset rates (kHz): 1, 2, 3, 4, 5, 6, 8, 10, 12, 15

        Returns:
        rate_sel (int): the rate sel PV for the timing trigger
        """

        rates_kHz = np.array([15, 12, 10, 8, 6, 5, 4, 3, 2, 1])

        try:
            idx = np.where(rates_kHz == val)[0][0] # weird numpy thing sorry
            return idx
        except IndexError:
            self.log("Reset rate not allowed! Look up help for allowed values")
            return


    def flux_ramp_PV_to_rate(self, val):
        """
        Convert between PV number in timing triggers and output flux ramp reset rate

        Returns:
        reset_rate (int): the flux ramp reset rate, in kHz
        """

        rates_kHz = [15, 12, 10, 8, 6, 5, 4, 3, 2, 1]
        return rates_kHz[val]


    def why(self):
        """
        Why not?
        """
        util_dir = os.path.dirname(__file__)
        aphorisms = open(os.path.join(util_dir, 'aphorism.txt'),'r').readlines()
        aph = np.random.choice(aphorisms).rstrip()
        self.log(aph)
        return (aph)

    def make_channel_mask(self, band=None, smurf_chans=None):
        """
        Makes the channel mask. Only the channels in the
        mask will be streamed or written to disk.

        If no optional arguments are given, mask will contain all channels
        that are on. If both band and smurf_chans are supplied, a mask
        in the input order is created.

        Args
        ----
        band : int array or None, optional, default None
            An array of band numbers. Must be the same length as
            smurf_chans
        smurf_chans : int_array or None, optional, default None
            An array of SMuRF channel numbers.  Must be the same
            length as band.

        Returns
        -------
        output_chans : int array
            The output channels.
        """
        output_chans = np.array([], dtype=int)

        # If no input, build one by querying pyrogue
        if smurf_chans is None and band is not None:
            band = np.ravel(np.array(band))
            n_chan = self.get_number_channels(band)
            output_chans = np.arange(n_chan) + n_chan*band
        # Take user inputs and make the channel map
        elif smurf_chans is not None:
            keys = smurf_chans.keys()  # the band numbers
            for k in keys:
                n_chan = self.get_number_channels(k)
                for ch in smurf_chans[k]:
                    output_chans = np.append(output_chans,
                                             ch + n_chan*k)

        return output_chans


    def make_freq_mask(self, mask):
        """
        Makes the frequency mask. These are the frequencies
        associated with the channels in the channel mask.

        Args
        ----
        mask : int array
            The channel mask file.

        Returns
        -------
        freqs : float array
            An array with frequencies associated with the mask file.
        """
        freqs = np.zeros(len(mask), dtype=float)
        bands = mask // 512
        chans = mask % 512

        for b in range(8):
            m = bands == b
            if not m.any():
                continue
            freqs[m] = self.channel_to_freq(b, chans[m])

        return freqs


    def set_downsample_filter(self, filter_order, cutoff_freq, write_log=False):
        """
        Sets the downsample filter. This is anti-alias filter
        that filters data at the flux_ramp reset rate, which is
        before the downsampler.

        Args
        ----
        filter_order : int
            The number of poles in the filter.
        cutoff_freq : float
            The filter cutoff frequency.
        """
        # Get flux ramp frequency
        flux_ramp_freq = self.get_flux_ramp_freq()*1.0E3

        # Get filter parameters
        b, a = signal.butter(filter_order,
                             2*cutoff_freq/flux_ramp_freq)

        # Set filter parameters
        self.set_filter_order(filter_order, write_log=write_log)
        self.set_filter_a(a, write_log=write_log)
        self.set_filter_b(b, write_log=write_log, wait_done=True)

        self.set_filter_reset(wait_after=.1, write_log=write_log)


    def get_filter_params(self):
        """
        Get the downsample filter parameters: filter order, filter
        gain, num averages, and the actual filter parameters.

        If filter order is -1, the downsampler is using a rectangula
        integrator. This will set filter_a, filter_b to None.

        Returns
        -------
        filter_params : dict
            A dictionary with the filter parameters.
        """
        # Get filter order, gain, and averages
        filter_order = self.get_filter_order()
        filter_gain = self.get_filter_gain()
        num_averages = self.get_downsample_factor()

        # Get filter order, gain, and averages

        if filter_order < 0:
            a = None
            b = None
        else:
            # Get filter parameters - (filter_order+1) elements
            a = self.get_filter_a()[:filter_order+1]
            b = self.get_filter_b()[:filter_order+1]

        # Cast into dictionary
        ret = {
            'filter_order' : filter_order,
            'filter_gain': filter_gain,
            'num_averages' : num_averages,
            'filter_a' : a,
            'filter_b' : b
        }

        return ret

    @set_action()
    def make_gcp_mask(self, band=None, smurf_chans=None,
                      gcp_chans=None, read_gcp_mask=True):
        """
        THIS FUNCTION WAS USED FOR BKUMUX DATA ACQUISITION.  IT'S
        COMPLETELY BROKEN NOW, POST ROGUE4 MIGRATION.

        Makes the gcp mask. Only the channels in this mask will be stored
        by GCP.

        If no optional arguments are given, mask will contain all channels
        that are on. If both band and smurf_chans are supplied, a mask
        in the input order is created.

        Args
        ----
        band : int array or None, optional, default None
            An array of band numbers. Must be the same length as
            smurf_chans
        smurf_chans : int_array or None, optional, default None
            An array of SMuRF channel numbers.  Must be the same
            length as band.
        gcp_chans : int_array or None, optional, default None
            A list of smurf numbers to be passed on as GCP channels.
        read_gcp_mask : bool, optional, default True
            Whether to read in the new GCP mask file.  If not read in,
            it will take no effect.

        """
        gcp_chans = np.array([], dtype=int)
        if smurf_chans is None and band is not None:
            band = np.ravel(np.array(band))
            n_chan = self.get_number_channels(band)
            gcp_chans = np.arange(n_chan) + n_chan*band
        elif smurf_chans is not None:
            keys = smurf_chans.keys()
            for k in keys:
                self.log(f'Band {k}')
                n_chan = self.get_number_channels(k)
                for ch in smurf_chans[k]:
                    gcp_chans = np.append(gcp_chans, ch + n_chan*k)

        n_channels = self.get_number_channels(band)
        if len(gcp_chans) > n_channels:
            self.log('WARNING: too many gcp channels!')
            return

        self.log(f'Generating gcp mask file. {len(gcp_chans)} ' +
                 'channels added')

        np.savetxt(self.smurf_to_mce_mask_file, gcp_chans, fmt='%i')

        if read_gcp_mask:
            self.read_smurf_to_gcp_config()
        else:
            self.log('Warning: new mask has not been read in yet.')

    @set_action()
    def bias_bump(self, bias_group, wait_time=.5, step_size=0.001,
                  duration=5.0, start_bias=None, make_plot=False,
                  skip_samp_start=10, high_current_mode=True,
                  skip_samp_end=10, plot_channels=None,
                  gcp_mode=False, gcp_wait=0.5, gcp_between=1.0,
                  dat_file=None, offset_percentile=2):
        """
        Toggles the TES bias high and back to its original state. From this, it
        calculates the electrical responsivity (sib), the optical responsivity (siq),
        and resistance.

        This is optimized for high_current_mode. For low current mode, you will need
        to step much slower. Try wait_time=1, step_size=.015, duration=10,
        skip_samp_start=50, skip_samp_end=50.

        Note that only the resistance is well defined now because the phase response
        has an un-set factor of -1. We will need to calibrate this out.

        Args
        ----
        bias_group : int of int array
            The bias groups to toggle. The response will return every
            detector that is on.
        wait_time : float, optional, default 0.5
            The time to wait between steps
        step_size : float, optional, default 0.001
            The voltage to step up and down in volts (for low current
            mode).
        duration : float, optional, default 5.0
            The total time of observation.
        start_bias : float or None, optional, default None
            The TES bias to start at. If None, uses the current TES
            bias.
        make_plot : bool, optional, default False
            Whether to make plots. Must set some channels in
            plot_channels.
        skip_samp_start : int, optional, default 10
            The number of samples to skip before calculating a DC
            level.
        high_current_mode : bool, optional, default True
            Whether to observe in high or low current mode.
        skip_samp_end : int, optional, default 10
            The number of samples to skip after calculating a DC
            level.
        plot_channels : int array or None, optional, default None
           The channels to plot.
        dat_file : str or None, optional, default None
            Filename to read bias-bump data from; if provided, data is
            read from file instead of being measured live.
        offset_percentile : float, optional, default 2.0
            Number between 0 and 100. Determines the percentile used
            to calculate the DC level of the TES data.

        Returns
        -------
        bands : int array
           The bands.
        channels : int array
           The channels.
        resistance : float array
            The inferred resistance of the TESs in Ohms.
        sib : float array
            The electrical responsivity. This may be incorrect until
            we define a phase convention. This is dimensionless.
        siq : float array
            The power responsivity. This may be incorrect until we
            define a phase convention. This is in uA/pW

        """
        if duration < 10* wait_time:
            self.log('Duration must bee 10x longer than wait_time for high enough' +
                     ' signal to noise.')
            return

        # Calculate sampling frequency
        # flux_ramp_freq = self.get_flux_ramp_freq() * 1.0E3
        # fs = flux_ramp_freq * self.get_downsample_factor()

        # Cast the bias group as an array
        bias_group = np.ravel(np.array(bias_group))

        # Fill in bias array if not provided
        if start_bias is not None:
            all_bias = self.get_tes_bias_bipolar_array()
            for bg in bias_group:
                all_bias[bg] += start_bias
            start_bias = all_bias
        else:
            start_bias = self.get_tes_bias_bipolar_array()

        step_array = np.zeros_like(start_bias)
        for bg in bias_group:
            step_array[bg] = step_size

        n_step = int(np.floor(duration / wait_time / 2))
        if high_current_mode:
            self.set_tes_bias_high_current(bias_group)

        if dat_file is None:
            filename = self.stream_data_on(make_freq_mask=False)

            if gcp_mode:
                self.log('Doing GCP mode bias bump')
                for j, bg in enumerate(bias_group):
                    self.set_tes_bias_bipolar(bg, start_bias[j] + step_size,
                                           wait_done=False)
                time.sleep(gcp_wait)
                for j, bg in enumerate(bias_group):
                    self.set_tes_bias_bipolar(bg, start_bias[j],
                                          wait_done=False)
                time.sleep(gcp_between)
                for j, bg in enumerate(bias_group):
                    self.set_tes_bias_bipolar(bg, start_bias[j] + step_size,
                                           wait_done=False)
                time.sleep(gcp_wait)
                for j, bg in enumerate(bias_group):
                    self.set_tes_bias_bipolar(bg, start_bias[j],
                                          wait_done=False)

            else:
                # Sets TES bias high then low
                for _ in np.arange(n_step):
                    self.set_tes_bias_bipolar_array(start_bias + step_array,
                                                    wait_done=False)
                    time.sleep(wait_time)

                    self.set_tes_bias_bipolar_array(start_bias,
                                                    wait_done=False)
                    time.sleep(wait_time)

            self.stream_data_off(register_file=True)  # record data
        else:
            filename = dat_file

        if gcp_mode:
            return

        t, d, m, v_bias = self.read_stream_data(filename,
                                                return_tes_bias=True)

        # flag region after step
        flag = np.ediff1d(v_bias[bias_group[0]],
                          to_end=0).astype(bool)
        flag = self.pad_flags(flag, after_pad=20,
                              before_pad=2)

        # flag first full step
        s, e = self.find_flag_blocks(flag)
        flag[0:s[1]] = np.nan

        v_bias *= -2 * self._rtm_slow_dac_bit_to_volt  # FBU to V
        d *= self._pA_per_phi0/(2*np.pi*1.0E6) # Convert to microamp
        i_amp = step_size / self._bias_line_resistance * 1.0E6 # also uA
        i_bias = v_bias[bias_group[0]] / self._bias_line_resistance * 1.0E6


        # Scale the currents for high/low current
        if high_current_mode:
            i_amp *= self._high_low_current_ratio
            i_bias *= self._high_low_current_ratio

        # Demodulation timeline
        demod = (v_bias[bias_group[0]] - np.min(v_bias[bias_group[0]]))
        _amp = (np.max(np.abs(v_bias[bias_group[0]])) -
                np.min(np.abs(v_bias[bias_group[0]])))
        demod /= (_amp/2)
        demod -= 1
        demod[flag] = np.nan

        bands, channels = np.where(m!=-1)
        resp = np.zeros(len(bands))
        sib = np.zeros(len(bands))*np.nan

        timestamp = filename.split('/')[-1].split('.')[0]

        # Needs to be an array for the check later
        if plot_channels is None:
            plot_channels = np.array([])

        for i, (b, c) in enumerate(zip(bands, channels)):
            mm = m[b, c]
            offset = (np.percentile(d[mm], 100-offset_percentile) +
                      np.percentile(d[mm], offset_percentile))/2
            d[mm] -= offset

            # Calculate response amplitude and S_IB
            resp[i] = np.nanmedian(2*d[mm]*demod/i_amp)
            sib[i] = resp[i] / i_amp

            if c in plot_channels:
                fig, ax = plt.subplots(2, figsize=(4.5, 3),
                                       sharex=True)
                ax[0].plot(d[mm], label='resp')
                ax[0].plot(i_bias-np.min(i_bias), label='bias')

                ax[1].plot(2*d[mm]*demod/i_amp, color='k')
                ax[0].legend(loc='upper right')
                ax[0].set_ylabel('Current [uA]')
                ax[1].set_ylabel('Resp [A/A]')

                ax[1].set_xlabel('Samples')
                ax[0].set_title(f'Bias bump - b{b}ch{c:03}')
                plt.tight_layout()

                # Make plot name path
                plot_fn = os.path.join(self.plot_dir,
                                       f'{timestamp}_bias_bump_b{b}ch{c:03}.png')
                plt.savefig(plot_fn)
                plt.close(fig)
                self.pub.register_file(plot_fn, 'bias_bump', plot=True)

        resistance = np.abs(self._R_sh * (1-1/sib))
        siq = (2*sib-1)/(self._R_sh*i_amp) * 1.0E6/1.0E12  # convert to uA/pW

        ret = {}
        for b in np.unique(bands):
            ret[b] = {}
            idx = np.where(bands == b)[0]
            for i in idx:
                c = channels[i]
                ret[b][c] = {}
                ret[b][c]['resp'] = resp[i]
                ret[b][c]['R'] = resistance[i]
                ret[b][c]['Sib'] = sib[i]
                ret[b][c]['Siq'] = siq[i]

        return ret


    def all_off(self):
        """
        Turns off everything. Does band off, flux ramp off, then TES bias off.
        """
        self.log('Turning off tones')
        bands = self._bands
        for b in bands:
            self.band_off(b)

        self.log('Turning off flux ramp')
        self.flux_ramp_off()

        self.log('Turning off all TES biases')
        n_bias_groups = self._n_bias_groups
        for bg in np.arange(n_bias_groups):
            self.set_tes_bias_bipolar(bg, 0)


    def mask_num_to_gcp_num(self, mask_num):
        """
        Goes from the smurf2mce mask file to a gcp number.
        Inverse of gcp_num_to_mask_num.

        Args
        ----
        mask_num : int
            The index in the mask file.

        Returns
        -------
        gcp_num : int
            The index of the channel in GCP.
        """
        return (mask_num*33)%528+mask_num//16


    def gcp_num_to_mask_num(self, gcp_num):
        """
        Goes from a GCP number to the smurf2mce index.
        Inverse of mask_num_to_gcp_num

        Args
        ----
        gcp_num : int
            The gcp index.

        Returns
        -------
        mask_num : int
            The index in the mask.
        """
        return (gcp_num*16)%528 + gcp_num//33


    def smurf_channel_to_gcp_num(self, band, channel):
        """
        Converts from smurf channel (band and channel) to a gcp number

        Args
        ----
        band : int
            The smurf band number.
        channel : int
            The smurf channel number.

        Returns
        -------
        gcp_num : int
            The GCP number.
        """
        mask = self.get_channel_mask()

        if mask[band, channel] == -1:
            self.log(f'Band {band} Ch {channel} not in mask')
            return None

        return self.mask_num_to_gcp_num(mask[band, channel])


    def gcp_num_to_smurf_channel(self, gcp_num):
        """
        Converts from gcp number to smurf channel (band and channel).

        Args
        ----
        gcp_num : int
            The GCP number.

        Returns
        -------
        band : int
            The smurf band number.
        channel : int
            The smurf channel number.
        """
        mask = self.get_channel_mask()
        n_channels = self.get_number_channels()

        mask_num = self.gcp_num_to_mask_num(gcp_num)
        return int(mask[mask_num]//n_channels), int(mask[mask_num]%n_channels)


    def play_sine_tes(self, bias_group, tone_amp, tone_freq, dc_amp=None):
        """
        Play a sine wave on the bias group pair.

        Tone file is in bias bit units. The bias is int20. The
        inputs of this function are in units of bias dac output
        voltage. The conversion from requested volts to bits
        is calculated in this function.

        Args
        ----
        bias_group : int
            The bias group to play a sine wave on.
        tone_amp : float
            The amplitude of the sine wave in units of out TES bias in
            volts.
        tone_freq : float
            The frequency of the tone in Hz.
        dc_amp : float or None, optional, default None
            The amplitude of the DC term of the sine wave.  If None,
            reads the current DC value and uses that.
        """
        if dc_amp is None:
            dc_amp = self.get_tes_bias_bipolar(bias_group)
            self.log(f"No dc_amp provided. Using current value: {dc_amp} V")

        # The waveform is played on 2 DACs, so amp/2. Then convert
        # to bits
        dc_amp /= (2*self._rtm_slow_dac_bit_to_volt)
        tone_amp /= (2*self._rtm_slow_dac_bit_to_volt)

        # Handles issue where it won't play faster than ~7 Hz
        freq_split = 5
        scale = 1
        if tone_freq > freq_split:
            scale = np.ceil(tone_freq / freq_split)

        # Make tone file. 2048 elements
        n_tes_samp = 2048
        sig = tone_amp * \
            np.cos(2*np.pi*scale*np.arange(n_tes_samp)/n_tes_samp) + dc_amp

        # Calculate frequency - 6.4ns * TimerSize between samples
        ts = int((tone_freq * n_tes_samp * 6.4E-9)**-1)
        ts *= scale
        self.set_rtm_arb_waveform_timer_size(ts, wait_done=True)

        self.play_tes_bipolar_waveform(bias_group, sig)


    def play_tone_file(self, band, tone_file=None, load_tone_file=True):
        """
        Plays the specified tone file on this band.  If no path provided
        for tone file, assumes the path to the correct tone file has
        already been loaded.

        Args
        ----
        band : int
            Which band to play tone file on.
        tone_file : str or None, optional, default None
            Path (including csv file name) to tone file.  If None,
            uses whatever's already been loaded.
        load_tone_file : bool, optional, default True
            Whether or not to load the tone file.  The tone file is
            loaded per DAC, so if you already loaded the tone file for
            this DAC you don't have to do it again.
        """
        # the bay corresponding to this band.
        bay = self.band_to_bay(band)

        # load the tone file
        if load_tone_file:
            self.load_tone_file(bay,tone_file)

        # play it!
        self.log(f'Playing tone file {tone_file} on band {band}',
                 self.LOG_USER)
        self.set_waveform_select(band, 1)


    def stop_tone_file(self, band):
        """
        Stops playing tone file on the specified band and reverts
        to DSP.

        Args
        ----
        band : int
            Which band to play tone file on.
        """

        self.set_waveform_select(band,0)

        # may need to do this, not sure.  Try without
        # for now.
        #self.set_dsp_enable(band,1)


    def get_gradient_descent_params(self, band):
        """
        Convenience function for getting all the serial
        gradient descent parameters

        Args
        ----
        band : int
            The band to query.

        Returns
        -------
        params : dict
            A dictionary with all the gradient descent parameters
        """
        ret = {}
        ret['averages'] = self.get_gradient_descent_averages(band)
        ret['beta'] = self.get_gradient_descent_beta(band)
        ret['converge_hz'] = self.get_gradient_descent_converge_hz(band)
        ret['gain'] = self.get_gradient_descent_gain(band)
        ret['max_iters'] = self.get_gradient_descent_max_iters(band)
        ret['momentum'] = self.get_gradient_descent_momentum(band)
        ret['step_hz'] = self.get_gradient_descent_step_hz(band)

        return ret


    def set_fixed_tone(self, freq_mhz, tone_power, write_log=False):
        """
        Places a fixed tone at the requested frequency.  Asserts
        without doing anything if the requested resonator frequency
        falls outside of the usable 500 MHz bands, or if there are no
        unassigned channels available in the subband the requested
        frequency falls into (where a channel is deemed "assigned" if
        it has non-zero amplitude).

        Args
        ----
        freq_mhz : float
            The frequency in MHz at which to place a fixed tone.
        tone_power : int
            The amplitude for the fixed tone (0-15 in recent fw
            revisions).
        write_log : bool, optional, default False
            Whether to write low-level commands to the log file.

        Returns
        -------
        band : int
            The band number in which a tone was turned on.
        channel : int
            The band channel number that was used to turn on the tone.
        """

        # Find which band the requested frequency falls into.
        bands=self.which_bands()
        band_centers_mhz=[self.get_band_center_mhz(b) for b in bands]

        band_idx=min(range(len(band_centers_mhz)), key=lambda i: abs(band_centers_mhz[i]-freq_mhz))
        band = bands[band_idx]
        band_center_mhz=band_centers_mhz[band_idx]

        # Confirm that the requested frequency falls into a 500 MHz
        # band that's usable in this fw.  If not, assert.
        assert (np.abs(freq_mhz-band_center_mhz)<250),\
            f'! Requested frequency (={freq_mhz:0.1f} MHz) outside of the ' + \
            '500 MHz band with the closest band center ' + \
            f'(={band_center_mhz:0.0f} MHz). Doing nothing!'

        # Find subband this frequency falls in, and its channels.
        subband, foff = self.freq_to_subband(band, freq_mhz)
        subband_channels = self.get_channels_in_subband(band, subband)

        # Which channels in the subband are unallocated?
        allocated_channels=self.which_on(band)
        unallocated_channels=[chan for chan in subband_channels if chan not in allocated_channels]
        # If no unallocated channels available in the subband, assert.
        assert (len(unallocated_channels)), \
            f'! No unallocated channels available in subband (={subband:d}).' + \
            ' Doing nothing!'

        # Take lowest channel number in the list of unallocated
        # channels for this subband.
        channel = sorted(unallocated_channels)[0]

        # Put a fixed tone at the requested frequency
        self.set_center_frequency_mhz_channel(band, channel, foff)
        self.set_amplitude_scale_channel(band, channel, tone_power)
        self.set_feedback_enable_channel(band, channel, 0)

        # Unless asked to be quiet, print where we're putting a fixed
        # tone.
        if write_log:
            self.log(f'Setting a fixed tone at {freq_mhz:.2f} MHz' +
                     f' and amplitude {tone_power}', self.LOG_USER)

        return band, channel

    # SHOULD MAKE A GET FIXED TONE CHANNELS FUNCTION - WOULD MAKE IT
    # EASIER TO CHANGE THINGS FAST USING THE ARRAY GET/SETS

    def turn_off_fixed_tones(self,band):
        """
        Turns off every channel which has nonzero amplitude but
        feedback set to zero.

        Args
        ----
        band : int
            The band to query.
        """
        amplitude_scale_array=self.get_amplitude_scale_array(band)
        feedback_enable_array=self.get_feedback_enable_array(band)

    # want to turn off all channels for which the amplitude is
    # nonzero, but feedback is not enabled.
        fixed_tone_channels=np.where((amplitude_scale_array*(1-feedback_enable_array))!=0)
        new_amplitude_scale_array=amplitude_scale_array.copy()
        new_amplitude_scale_array[fixed_tone_channels]=0

    # set by array, not by channel
        self.set_amplitude_scale_array(band,new_amplitude_scale_array)

    __hardware_logging_pause_event=None

    def pause_hardware_logging(self):
        """Pauses hardware logging thread.

        See Also
        --------
        resume_hardware_logging : Resumes hardware logging thread.
        start_hardware_logging : Starts hardware logging thread.
        stop_hardware_logging : Stops hardware logging thread.
        """
        self.__hardware_logging_pause_event.set()

    def resume_hardware_logging(self):
        """Resumes hardware logging thread.

        See Also
        --------
        pause_hardware_logging : Pauses hardware logging thread.
        start_hardware_logging : Starts hardware logging thread.
        stop_hardware_logging : Stops hardware logging thread.
        """
        self.__hardware_logging_pause_event.clear()

    __hardware_log_file=None

    def get_hardware_log_file(self):
        """Returns path to current hardware log file.

        Returns
        -------
        str or None
           Path on disk to current hardware log file.  If not
           currently hardware logging, returns None.

        See Also
        --------
        start_hardware_logging : Starts hardware logging thread.
        """
        return self.__hardware_log_file

    _hardware_logging_thread=None
    __hardware_logging_stop_event=None

    def start_hardware_logging(self,filename=None,wait_btw_sec=5.0):
        """Starts hardware logging in external thread.

        Args
        ----
        filename : str or None, optional, default None
           Name of file on disk to write hardware logging to
           (including path).  If None, file name is automatically
           generated as *CTIME_hwlog.dat* with CTIME the current unix
           epoch timestamp returned by
           :func:`~pysmurf.client.base.smurf_control.SmurfControl.get_timestamp`,
           and saved in the directory specified by the
           :class:`~pysmurf.client.base.smurf_control.SmurfControl`
           class attribute
           :attr:`~pysmurf.client.base.smurf_control.SmurfControl.output_dir`.
        wait_btw_sec : float, optional, default 5.0 Time to wait, in
           seconds, between each poll of the hardware registers being
           logged.

        See Also
        --------
        get_hardware_log_entry : Generates each row of hardware logging data written to file.
        pause_hardware_logging : Pauses hardware logging thread.
        resume_hardware_logging : Resumes hardware logging thread.
        stop_hardware_logging : Stops hardware logging thread.
        """
        # Just in case somewhere the enable got set to false,
        # explicitly enable here
        if filename is None:
            filename=os.path.join(
                self.output_dir,
                str(self.get_timestamp())+'_hwlog.dat')
        self.__hardware_log_file = os.path.join(filename)
        self.log('Starting hardware logging to file : ' +
                 f'{self.__hardware_log_file}',
                 self.LOG_USER)
        self.__hardware_logging_stop_event=threading.Event()
        self.__hardware_logging_pause_event=threading.Event()
        self._hardware_logging_thread = threading.Thread(target=self._hardware_logger,
            args=(
                self.__hardware_logging_pause_event,
                self.__hardware_logging_stop_event,wait_btw_sec))
        self._hardware_logging_thread.daemon = True
        self._hardware_logging_thread.start()

    def stop_hardware_logging(self):
        """Stops and cleans up hardware logging thread.

        See Also
        --------
        start_hardware_logging : Starts hardware logging thread.
        """
        self.__hardware_logging_stop_event.set()
        self._hardware_logging_thread.join()
        self._hardware_logging_thread=None
        self.__hardware_log_file=None

    def _hardware_logger(self,pause_event,stop_event,wait_btw_sec=5):
        """Hardware logging thread function.

        Args
        ----
        pause_event : :py:class:`threading.Event`
           :py:class:`threading.Event` object for pausing the hardware
           logging thread.
        stop_event : :py:class:`threading.Event`
           :py:class:`threading.Event` object for stopping the
           hardware logging thread.
        wait_btw_sec : float, optional, default 5.0 Time to wait, in
           seconds, between each poll of the hardware registers being
           logged.

        See Also
        --------
        pause_hardware_logging : Pauses hardware logging thread.
        resume_hardware_logging : Resumes hardware logging thread.
        start_hardware_logging : Starts hardware logging thread.
        stop_hardware_logging : Stops hardware logging thread.
        """
        filename=self.get_hardware_log_file()
        import fcntl
        #counter=0
        while not stop_event.wait(wait_btw_sec):
            if not pause_event.isSet():
                hdr,entry=self.get_hardware_log_entry()
                # only write header once, if file doesn't exist yet if
                # file *does* already exist, check to make sure header
                # will be the same, otherwise the resulting data won't
                # make sense if multiple carriers are logging to the same
                # file.
                if not os.path.exists(filename):
                    with open(filename,'a') as logf:
                        logf.write(hdr)
                else:
                    with open(filename) as logf:
                        hdr2=logf.readline()
                        if not hdr.rstrip().split() == hdr2.rstrip().split():
                            self.log('Attempting to hardware log to an ' +
                                'incompatible file.  Giving up without ' +
                                'logging any data!', self.LOG_ERROR)
                            return

                with open(filename,'a') as logf:
                    # file locking so multiple hardware loggers running in
                    # multiple pysmurf sessions can write to the same
                    # requested file if desired
                    fcntl.flock(logf, fcntl.LOCK_EX)
                    logf.write(entry)
                    fcntl.flock(logf, fcntl.LOCK_UN)
                #counter+=1

    def get_hardware_log_entry(self):
        """Returns hardware log file header and data.

        Returns
        -------
        hdr : str
           Header for hardware log file.
        row : str
           One row of data for hardware log file.  Measured values are
           polled once each time function is called.

        See Also
        --------
        start_hardware_logging : Starts hardware logging thread.
        """
        d={}
        d['epics_root']=lambda:self.epics_root
        d['ctime']=self.get_timestamp
        d['fpga_temp']=self.get_fpga_temp
        d['fpgca_vccint']=self.get_fpga_vccint
        d['fpgca_vccaux']=self.get_fpga_vccaux
        d['fpgca_vccbram']=self.get_fpga_vccbram
        d['cc_temp']=self.get_cryo_card_temp

        # probably should check for which AMCs are in in a smarter way
        bays=[]
        bands=self.which_bands()
        if 0 in bands:
            bays.append(0)
        if 4 in bands:
            bays.append(1)

        for bay in bays:
            for dac in [0,1]:
                d[f'bay{bay}_dac{dac}_temp']=(
                    lambda:self.get_dac_temp(bay,dac))

        # atca monitor
        d['atca_temp_fpga']=self.get_board_temp_fpga
        d['atca_temp_rtm']=self.get_board_temp_rtm
        d['atca_temp_amc0']=(
            lambda:self.get_board_temp_amc(bay=0))
        d['atca_temp_amc2']=(
            lambda:self.get_board_temp_amc(bay=1))
        d['atca_jct_temp_fpga']=self.get_junction_temp_fpga

        # regulator
        d['regulator_iout']=self.get_regulator_iout
        d['regulator_temp1']=self.get_regulator_temp1
        d['regulator_temp2']=self.get_regulator_temp2

        columns=[]
        names=[]
        fmt=''
        counter=0
        for key, value in d.items():
            columns.append(str(value()))
            names.append(key)
            fmt+=('{0['+f'{counter}'+']:<20}')
            counter+=1
        fmt+='\n'

        hdr=fmt.format(names)
        row=fmt.format(columns)
        return hdr,row

    def play_tes_bipolar_waveform(self, bias_group, waveform, do_enable=False,
            continuous=True, **kwargs):
        """ Play a bipolar waveform on the bias group.

        Args
        ----
        bias_group : int
            The bias group
        waveform : float array
            The waveform the play on the bias group.
        do_enable : bool, optional, default True
            Whether to enable the DACs (similar to what is required
            for TES bias).
        continuous : bool, optional, default True
            Whether to play the TES waveform continuously.
        """
        bias_order = self.bias_group_to_pair[:,0]
        dac_positives = self.bias_group_to_pair[:,1]
        dac_negatives = self.bias_group_to_pair[:,2]

        dac_idx = np.ravel(np.where(bias_order == bias_group))

        dac_positive = dac_positives[dac_idx][0]
        dac_negative = dac_negatives[dac_idx][0]

        # https://confluence.slac.stanford.edu/display/SMuRF/SMuRF+firmware#SMuRFfirmware-RTMDACarbitrarywaveforms
        # Target the two bipolar DACs assigned to this bias group:
        self.set_dac_axil_addr(0, dac_positive)
        self.set_dac_axil_addr(1, dac_negative)

        # Enable waveform generation (3=on both DACs)
        self.set_rtm_arb_waveform_enable(3)

        # Must enable the DACs (if not enabled already)
        if do_enable:
            self.set_rtm_slow_dac_enable(dac_positive, 2, **kwargs)
            self.set_rtm_slow_dac_enable(dac_negative, 2, **kwargs)

        # Load waveform into each DAC's LUT table.  Opposite sign so
        # they combine coherently
        self.set_rtm_arb_waveform_lut_table(0, waveform)
        self.set_rtm_arb_waveform_lut_table(1, -waveform)

        # Continous mode to play the waveform continuously
        if continuous:
            self.set_rtm_arb_waveform_continuous(1)
        else:
            self.set_rtm_arb_waveform_continuous(0)

    # Readback on which DACs are selected is broken right now,
    # so has to be specified.
    def stop_tes_bipolar_waveform(self, bias_group, **kwargs):
        """
        Stop the bipolar waveform being played on a bias group.

        Args
        ----
        bias_group : int
            The bias group.
        """
        # https://confluence.slac.stanford.edu/display/SMuRF/SMuRF+firmware#SMuRFfirmware-RTMDACarbitrarywaveforms
        # Target the two bipolar DACs assigned to this bias group:
        self.set_dac_axil_addr(0,0) # Disabled
        self.set_dac_axil_addr(1,0) # Disabled

        # Enable waveform generation (3=on both DACs)
        self.set_rtm_arb_waveform_enable(0)

        # Zero TES biases on this bias group
        self.set_tes_bias_bipolar(bias_group, 0)

    @set_action()
    def get_sample_frequency(self):
        """ Gives the data rate.

        Returns
        -------
        sample_frequency : float
            The data sample rate in Hz.
        """
        flux_ramp_freq = self.get_flux_ramp_freq() * 1.0E3
        downsample_factor = self.get_downsample_factor()

        return flux_ramp_freq / downsample_factor

    def identify_bias_groups(self, bias_groups=None,
            probe_freq=2.5, probe_time=3, probe_amp=.1, make_plot=False,
            show_plot=False, save_plot=True, cutoff_frac=.05,
            update_channel_assignment=True, high_current_mode=True):
        """ Identify bias groups of all the channels that are on. Plays
        a sine wave on a bias group and looks for a response. Does
        this with the TESs superconducting so it can look for an
        response is exactly the same amplitude as the input.

        Args
        ----
        bias_groups : int array or None, optional, default None
            The bias groups to search. If None, does the first 8 bias
            groups.
        probe_freq : float, optional, default 2.5
            The frequency of the probe tone.
        probe_time : float, optional, default 3
            The length of time to probe each bias group in seconds.
        probe_amp : float, optional, default 0.1
            Amplitude of the probe signal in volts.
        make_plot : bool, optional, default False
            Whether to make the plot.
        show_plot : bool, optional, default False
            Whether to show the plot.
        save_plot : bool, optional, default True
            Whether to save the plot.
        cutoff_frac : float, optional, default 0.05
            The fraction difference the response can be away from the
            expected amplitude.
        update_channel_assignment : bool, optional, default True
            Whether to update the master channels assignment to
            contain the new bias group information.
        high_current_mode : bool, optional, default True
            Whether to use high or low current mode.

        Returns
        -------
        channels_dict : dict of {int : dict of {str : numpy.ndarray} }
            A dictionary where the first key is the bias group that is
            being probed. In each is the band, channnel pairs, and
            frequency of the channels.
        """
        # Check if probe frequency is too high
        flux_ramp_freq = self.get_flux_ramp_freq() * 1.0E3
        fs = flux_ramp_freq * self.get_downsample_factor()

        # Calculate downsample filter transfer function
        filter_params = self.get_filter_params()
        w, h = signal.freqz(filter_params['filter_b'],
                            filter_params['filter_a'],
                            fs=flux_ramp_freq)
        df = np.abs(w - probe_freq)
        df_idx = np.ravel(np.where(df == np.min(df)))[0]
        if probe_freq > fs:
            self.log('Probe frequency is higher than sample rate. Exiting',
                     self.LOG_ERROR)
            return
        elif h[df_idx] < 1 - cutoff_frac/3:
            self.log('Downsample filter cutting into the signal too much.' +
                     ' Exiting.', self.LOG_ERROR)
            return

        # There should be something smarter than this
        if bias_groups is None:
            bias_groups = np.arange(self._n_bias_groups)

        channels_dict = {}

        # Get the cryocard settings before starting this script
        cryo_card_bits = self.get_cryo_card_relays()

        timestamp = self.get_timestamp()

        for bias_group in bias_groups:
            self.log(f"Working on bias group {bias_group}")

            # Work in high current mode to bypass filter
            if high_current_mode:
                self.set_tes_bias_high_current(bias_group)
            else:
                self.set_tes_bias_low_current(bias_group)

            # Play sine wave and take data
            self.play_sine_tes(bias_group, probe_amp, probe_freq, dc_amp=0)
            datafile = self.take_stream_data(probe_time, write_log=False)

            self.stop_tes_bipolar_waveform(bias_group)

            # Read back data
            t, d, mm = self.read_stream_data(datafile, make_freq_mask=True)
            m = mm[0]  # extract mask
            m_freq = mm[1]  #frequency mask
            freq_arr = m_freq[np.where(m!=-1)]
            d *= (self._pA_per_phi0/2/np.pi)  # convert to pA
            d = np.transpose(d.T - np.mean(d.T, axis=0))

            n_det, nsamp = np.shape(d)

            # currents on lines
            if high_current_mode:
                r_inline = self._bias_line_resistance / self._high_low_current_ratio
            else:
                r_inline = self._bias_line_resistance

            i_bias = probe_amp / r_inline * 1.0E12  # Bias current in pA

            # sine/cosine decomp templates
            s = np.sin(2*np.pi*np.arange(nsamp) / nsamp*probe_freq*probe_time)
            c = np.cos(2*np.pi*np.arange(nsamp) / nsamp*probe_freq*probe_time)
            s /= np.sum(s**2)
            c /= np.sum(c**2)

            # cosine/sine decomposition
            sa = np.zeros(n_det)
            ca = np.zeros(n_det)
            for ch in np.arange(n_det):
                sa[ch] = np.dot(d[ch], s)
                ca[ch]= np.dot(d[ch], c)
            amp = np.sqrt(sa**2 + ca**2)  # amplitude calculation

            # In superconducting, amplitude of response should be 1
            norm_amp = amp/i_bias
            idx = np.where(np.logical_and(norm_amp < 1+cutoff_frac,
                           norm_amp > 1-cutoff_frac))[0]

            bands = np.zeros(len(idx), dtype=int)
            channels = np.zeros(len(idx), dtype=int)
            freqs = np.zeros(len(idx))
            for i, ii in enumerate(idx):
                bands[i], channels[i] = np.ravel(np.where(m == ii))
                freqs[i] = m_freq[bands[i], channels[i]]

            channels_dict[bias_group] = {}
            channels_dict[bias_group]['band'] = bands
            channels_dict[bias_group]['channel'] = channels
            channels_dict[bias_group]['freq'] = freqs

            if make_plot:
                # Turn off interactive plot
                if show_plot:
                    plt.ion()
                else:
                    plt.ioff()

                fig, ax = plt.subplots(1, 2, figsize=(8.5, 3),
                                       sharey=True)
                # Plot timestreams
                ax[0].plot((t-t[0])*1.0E-9, d.T,
                         color='k', alpha=.1)
                ax[0].axhline(-i_bias, linestyle='--', color='r')
                ax[0].axhline(i_bias, linestyle='--', color='r')
                ax[1].axhline(-i_bias, linestyle='--', color='r')
                ax[1].axhline(i_bias, linestyle='--', color='r')
                ax[0].set_xlabel('Time [s]')
                ax[0].set_ylabel('Amp [pA]')

                current_mode_label = 'high current'
                if not high_current_mode:
                    current_mode_label = 'low current'
                ax[0].text(.02, .98, current_mode_label,
                    transform=ax[0].transAxes, va='top', ha='left')

                ax[1].plot(freq_arr, sa, 'x', color='b',
                           label='sine', alpha=.5)
                ax[1].plot(freq_arr, ca, '+', color='y',
                           label='cos', alpha=.5)
                ax[1].plot(freq_arr, amp, 'o', color='k',
                           label='amp')
                ax[1].legend(loc='lower right')
                ax[1].set_ylim((-1.5*i_bias, 1.5*i_bias))
                ax[1].set_xlabel('Res Freq [MHz]')
                plt.tight_layout()

                fig.suptitle(f'Bias Group {bias_group}')

                if save_plot:
                    savename = f'{timestamp}_identify_bg{bias_group:02}.png'
                    plt.savefig(os.path.join(self.plot_dir, savename),
                                bbox_inches='tight')
                    self.pub.register_file(
                        os.path.join(self.plot_dir, savename),
                        'identify_bg', plot=True)
                if not show_plot:
                    plt.close(fig)

            # Set relays back to original state
            self.set_cryo_card_relays(cryo_card_bits)

        # To do - add a check for band, channels that are on two different
        # bias groups.

        for bias_group in bias_groups:
            self.log(f'Bias Group {bias_group} : ')
            self.log(f"   Bands : {np.unique(channels_dict[bias_group]['band'])}")
            n_chan = len(channels_dict[bias_group]['channel'])
            self.log("   Number of channels : "+
                     f"{n_chan}")
            if n_chan > 0:
                ff = channels_dict[bias_group]['freq']
                self.log(f"   Between freq : {np.min(ff)} and {np.max(ff)}")


        if update_channel_assignment:
            self.log('Updating channel assignment')
            self.write_group_assignment(channels_dict)

        return channels_dict

    def find_probe_tone_gap(self, band, n_tone=2, n_scan=5,
                          make_plot=True,
                          save_plot=True, show_plot=False):
        """
        This finds the larges n_tone gaps in the band. These
        are used for finding gaps to put probe tones for checking
        JESD skips.

        Args
        ----
        band : int
            The 500 MHz frequency band to find gaps.
        n_tone : int, optional
            The number of probe tones. Or the number of gaps
            to find.
        n_scan : int, optional
            The number of times to measure the band using
            full_band_resp.
        make_plot : bool, optional, default True
            Whether to make a plot.
        save_plot : bool, optional, default True
            Whether to save the plot.
        show_plot : bool, optional, default False
            Whether to show the plot

        Returns
        -------
        gaps : float array
            An array of size [n_tone, 2], indicating the start and
            end of the gap. recomment using the middle value between
            the two edges. You can calulate this quickly
            with np.mean(gaps, axis=1).
        """
        timestamp = self.get_timestamp()

        # Run full band response to measure transfer function
        f, resp = self.full_band_resp(band, n_scan=n_scan,
                                      make_plot=make_plot,
                                      save_plot=save_plot,
                                      show_plot=False,
                                      timestamp=timestamp)

        f *= 1.0E-6  # Mitch wants units of MHz

        # Identify the resonator frequencies
        res_freq = self.find_peak(f, resp,
                                  make_plot=make_plot,
                                  save_plot=save_plot,
                                  show_plot=False)

        # Throw out frequencies not in the band
        res_freq = res_freq[np.logical_and(res_freq > -250,
                                           res_freq < 250)]

        df = np.diff(res_freq)
        df_vals = np.sort(df)[-n_tone:] # find the largest gaps

        # loop over gaps
        gap_freq = np.zeros((n_tone, 2))
        for i, dfv in enumerate(df_vals):
            idx = np.ravel(np.where(df == dfv))[0]
            gap_freq[i, 0] = res_freq[idx]
            gap_freq[i, 1] = res_freq[idx+1]


        if make_plot:
            plt.figure(figsize=(8,4))
            plt.plot(f, np.abs(resp))
            for i in range(n_tone):
                plt.axvspan(gap_freq[i,0], gap_freq[i,1], color='k', alpha=.1)

            plt.title(f'{timestamp} - find probe tone gaps')
            plt.ylabel(f'Band {band} resp')
            plt.xlabel('Freq [MHz]')
            plt.tight_layout()

            if save_plot:
                plt.savefig(os.path.join(self.plot_dir,
                                         f'{timestamp}_probe_tone_gap.png'),
                            bbox_inches='tight')

            if show_plot:
                plt.show()
            else:
                plt.close()

        return gap_freq

    def check_full_band_resp(self,
                             n_scan_per_band=5,
                             amc_auto_detect=True,
                             make_plot=True,
                             save_plot=True,
                             show_plot=False,
                             save_results=True,
                             **kwargs):
        r"""
        Measures RF transfer function for all configured AMCs.  This
        function is based on scratch/shawn/full_band_response.py.
        Usually we run this function with the AMCs connected in
        loopback with a short RF cable to verify that they are
        working.

        Args
        ----
        n_scan : int, optional, default 5
            The number of times to measure each band using
            full_band_resp.
        amc_auto_detect : bool, optional, default True
            If True, determines AMC RF bands from serial number.
            Otherwise, uses `get_band_center_mhz`.
        make_plot : bool, optional, default True
            Whether to make a plot.
        save_plot : bool, optional, default True
            Whether to save the plot.
        show_plot : bool, optional, default False
            Whether to show the plot
        save_results : bool, optional, default True
            Whether or not to save the results.
        \**kwargs
            Arbitrary keyword arguments.  Passed directly to the
            `full_band_resp` call.

        Returns
        -------
        results : dict
            Dictionary with results.

        """
        timestamp = self.get_timestamp()
        # Take which bands are configured from pysmurf configuration
        # file.
        bands=self._bands

        carrier_sn=self.get_carrier_sn(use_shell=True)

        results_dict={}
        results_dict['n_scan_per_band']=n_scan_per_band
        results_dict['amc_auto_detect']=amc_auto_detect
        results_dict['kwargs']=kwargs
        for bay in self.bays:
            results_dict[bay]={}
            bands=range(4*bay,4*(bay+1))
            for band in bands:
                print(' ')
                print(' ')
                print(f'Band {band}')
                print(' ')
                print(' ')
                results_dict[bay][band%4]={}
                amc_sn=self.get_amc_sn(bay=bay,use_shell=True)
                results_dict[bay]['amc_sn']=amc_sn
                amc_type=amc_sn.split('-')[1]

                # Sometimes the rogue zip files don't properly set the
                # correct band center frequency if e.g. a LB is
                # plugged into bay 1.
                if amc_auto_detect and amc_type in ['A01','A02']:
                    if amc_type=='A01': #LB
                        results_dict[bay][band%4]['fc']=4250.+500.*(band%4)
                    if amc_type=='A02': #HB
                        results_dict[bay][band%4]['fc']=6250.+500.*(band%4)

                else:
                    results_dict[bay][band%4]['fc']=self.get_band_center_mhz(band)

                f,resp=self.full_band_resp(band=band, make_plot=False, show_plot=False, n_scan=n_scan_per_band, timestamp=timestamp, save_data=False, **kwargs)
                results_dict[bay][band%4]['f']=f
                results_dict[bay][band%4]['resp']=resp

                # Also record the UC and DC attenuator settings for this
                # band
                results_dict[bay][band%4]['uc_att']=self.get_att_uc(band)
                results_dict[bay][band%4]['dc_att']=self.get_att_dc(band)

            # Plot results (if desired)
            fig, ax = plt.subplots(2, figsize=(6,7.5), sharex=True)
            plt.suptitle(f'slot={self.slot_number} / Carrier={carrier_sn} / AMC{bay*2}={amc_sn}')

            save_name = f'{timestamp}_full_band_resp_all_amc{bay*2}.png'
            ax[0].set_title(save_name)

            # Compute unwrapped phase for each band
            last_angle=None
            for band in bands:
                f_plot=results_dict[bay][band%4]['f']/1e6
                resp_plot=results_dict[bay][band%4]['resp']
                plot_idx = np.where(np.logical_and(f_plot>-250, f_plot<250))
                ax[0].plot(f_plot[plot_idx]+results_dict[bay][band%4]['fc'], np.log10(np.abs(resp_plot[plot_idx])),label=f'b{band%4}')
                angle = np.unwrap(np.angle(resp_plot))
                if last_angle is not None:
                    angle-=(angle[0]-last_angle)
                ax[1].plot(f_plot[plot_idx]+results_dict[bay][band%4]['fc'], angle[plot_idx],label=f'b{band%4}')
                results_dict[bay][band%4]['phase']=angle[plot_idx]
                last_angle=angle[plot_idx][-1]

            ax[0].legend(loc='lower left',fontsize=8)
            ax[0].set_ylabel("log10(abs(Response))")
            ax[0].set_xlabel('Frequency [MHz]')
            ax[0].set_ylim(-0.5,1.25)

            ax[1].legend(loc='lower left',fontsize=8)
            ax[1].set_ylabel("Phase [rad]")
            ax[1].set_xlabel('Frequency [MHz]')

            plt.tight_layout(rect=[0, 0.03, 1, 0.95])

            if save_plot:
                save_path = os.path.join(self.plot_dir, save_name)
                print(f'Saving plot to {save_path}.')
                plt.savefig(save_path,
                            bbox_inches='tight')
                self.pub.register_file(save_path, f'fbr{bay*2}_plot', plot=True)

                if not show_plot:
                    plt.close()

        if save_results:
            save_path = os.path.join(self.output_dir,f'{timestamp}_full_band_resp_all.npy')
            print(f'Saving data to {save_path}.')
            np.save(save_path,results_dict)
            self.pub.register_file(save_path, 'fbr_results', format='npy')

        if show_plot:
            plt.show()

        return results_dict

    def get_timing_mode(self):
        r"""Determines timing mode configuration.

        Returns the current timing mode configuration, or None if the
        system is not in one of these three known configurations:

        * "ext_ref" : locked to an external reference, or free running
          if an external reference is absent.
        * "backplane" : locked to external timing signals distributed
          over the crate backplane by a carrier in slot 2 of this
          carrier's crate receiving timing on its RTM's timing input.
          Only carriers in slot 2 of typical crates can distribute
          timing to other carriers through the crate backplane, so
          carriers configured in "backplane" timing mode must be in
          slots 3 or higher.
        * "fiber" : locked to external timing input from the carrier's
          RTM timing input and if the carrier is in slot 2 of a crate
          with a dual-start backplane, distributing the timing signals
          to all other carriers in the crate's backplane.

        The timing mode configuration is determined by polling the
        configuration of the timing crossbar, LMKs,
        triggers, and RTM.

        For systems configured in "fiber" or "backplane" modes, a
        warning is printed if no external timing data is being
        received.

        Returns
        -------
        mode : str or None
           Current timing mode configuration.  Returns None if system
           not in one of the three recognized configurations :
           "ext_ref", "backplane", or "fiber".

        See Also
        --------
        :func:`set_timing_mode` : Can be used to set the timing mode.
        :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.get_timing_link_up` : Is external timing data being received?
        """
        ## Poll all registers needed to determine which timing mode we're in.

        # Crossbar
        cbar = [self.get_crossbar_output_config(i) for i in range(4)]

        # RTM
        rsm = self.get_ramp_start_mode()

        # Timing triggers (only used with external timing system)
        ecre = self.get_evr_channel_reg_enable(0)
        etdt = self.get_evr_trigger_dest_type(0)
        te = self.get_trigger_enable(0)

        # LMKs
        lmks={}
        for bay in self.bays:
            lmks[bay]={}
            for reg in [0x146,0x147]:
                lmks[bay][reg]=self.get_lmk_reg(bay,reg)

        ## Check polled register values against known timing
        ## configurations

        # External reference timing mode configuration
        if ( cbar == [0x0, 0x0, 0x1, 0x1] and
             rsm == 0 and
             all([lmks[bay][0x146]==0x10 for bay in self.bays]) and
             all([lmks[bay][0x147]==0x1a for bay in self.bays]) ):
            return 'ext_ref'

        # Fiber or backplane timing mode configurations
        if ( rsm == 1 and
             ( ecre == 1 and etdt == 0 and te == 1 ) and
             all([lmks[bay][0x146]==0x8 for bay in self.bays]) and
             all([lmks[bay][0x147]==0xa for bay in self.bays]) ):

            # Fiber timing mode configuration
            if ( cbar == [0x0, 0x0, 0x0, 0x0] ):
                mode = 'fiber'

            # Backplane timing mode configuration
            if ( cbar == [0x0, 0x2, 0x1, 0x1] ):
                mode = 'backplane'

            # Check if receiving timing.  Checks if FPGA recovered
            # clock is receiving timing data.  Warn if we're
            # configured for backplane or fiber timing but not
            # receiving any external timing data.
            status = self.get_timing_link_up()
            if status != 1:
                self.log(f'\033[91mConfigured for {mode} timing but not receiving external timing data.\033[0m',
                         self.LOG_ERROR) # color red

            return mode

        # Timing configuration not recognized, return None
        return None

    def set_timing_mode(self, mode, write_log=False):
        r"""Sets timing mode.

        Use this function to configure the system's timing mode.  The
        currently supported configurations are:

        * "ext_ref" : lock system to an external reference, or free
          run if an external reference is absent.
        * "backplane" : lock system to external timing signals
          distributed over the crate backplane by a carrier receiving
          timing on its RTM's timing input in slot 2 running in
          "fiber" timing mode.  Only carriers in slots 3 or higher can
          lock to backplane distributed timing.
        * "fiber" : lock to external timing input from the carrier's
          RTM timing input and if the carrier is in slot 2 of a crate
          with a dual-start backplane, distribute the timing signals
          to all other carriers in the crate's backplane.

        The timing mode configuration adjusts the configuration of the
        timing crossbar, LMKs, triggers, and RTM.

        For systems configured in "fiber" or "backplane" modes, after
        configuration a warning is printed if no external timing data
        is being received after a one second wait after configuring
        the timing crossbar.

        Before changing the timing configuration, the current
        configuration is polled via :func:`get_timing_mode` and if the
        system is already in the requested mode, the function prints a
        message and does nothing.

        .. warning::
           If "backplane" mode is requested for a slot 2 system,
           `set_timing_mode` instead configures the system for "fiber"
           timing.  This may or may not be what was desired.

        Args
        ----
        mode : str
            The timing mode to configure the system with.  Currently
            valid options are "ext_ref", "backplane", or "fiber".
        write_log : bool, optional, default False
            Whether or not to print all EPICS calls, for logging or
            debugging purposes.  Passed through all register set calls
            to the underlying `epics.caput` calls.


        See Also
        --------
        :func:`get_timing_mode` : Can be used to set the timing mode.
        :func:`~pysmurf.client.command.smurf_command.SmurfCommandMixin.get_timing_link_up` : Is external timing data being received?
        """
        valid_timing_modes = ['ext_ref','backplane','fiber']
        assert (mode in valid_timing_modes), f'\033[91mRequested timing mode={mode} unknown.\033[0m'

        current_timing_mode = self.get_timing_mode()
        if mode == current_timing_mode:
            self.log(f'System already configured for {mode} timing, doing nothing.',self.LOG_USER)
            return

        self.log(f'System configured for {current_timing_mode} timing.  Reconfiguring slot {self.slot_number} for {mode} timing.',self.LOG_USER)

        # Make sure we can write to the LMK regs
        for bay in self.bays:
            self.set_lmk_enable(bay, 1, write_log=write_log)

        if mode == 'backplane' and self.slot_number == 2:
            # Complain if backplane timing mode is requested for a
            # carrier in slot 2 and configure for ext_ref timing
            # instead.  In backplane timing mode, the system must
            # receive timing through the backplane from a system in
            # slot 2 configured in fiber mode.
            self.log('\033[91mSystem is in slot 2, which cannot be configured for backplane timing.  Will configure for fiber timing instead.\033[0m',
                     self.LOG_ERROR) # color red
            mode = 'fiber'

        if mode == 'ext_ref':
            # Configure crossbar for ext_ref timing
            self.set_crossbar_output_config(0, 0x0, write_log=write_log)
            self.set_crossbar_output_config(1, 0x0, write_log=write_log)
            self.set_crossbar_output_config(2, 0x1, write_log=write_log)
            self.set_crossbar_output_config(3, 0x1, write_log=write_log)

            # Configure LMK
            for bay in self.bays:
                self.log(f'Select external reference for bay {bay}' +
                         ' or free running if there is no reference.')
                for bay in self.bays:
                    self.set_lmk_reg(bay, 0x146, 0x10, write_log=write_log)
                    self.set_lmk_reg(bay, 0x147, 0x1a, write_log=write_log)

            # Select internal flux ramp trigger source
            self.set_ramp_start_mode(0, write_log=write_log)

        if mode in ['backplane','fiber']:

            if mode == 'backplane':
                # Configure crossbar for backplane timing
                self.set_crossbar_output_config(0, 0x0, write_log=write_log)
                # OutputConfig[1] = 0x2 configures the SMuRF carrier's
                # FPGA to take the timing signals from the backplane
                # (TO_FPGA = FROM_BACKPLANE)
                self.set_crossbar_output_config(1, 0x2, write_log=write_log)
                self.set_crossbar_output_config(2, 0x1, write_log=write_log)
                self.set_crossbar_output_config(3, 0x1, write_log=write_log)

            if mode == 'fiber':
                # Configure crossbar for fiber timing
                self.set_crossbar_output_config(0, 0x0, write_log=write_log)
                self.set_crossbar_output_config(1, 0x0, write_log=write_log)
                self.set_crossbar_output_config(2, 0x0, write_log=write_log)
                self.set_crossbar_output_config(3, 0x0, write_log=write_log)

            # Configure triggering
            #
            # From Matt Weaver : The triggering firmware is broken
            # into two parts: (1) the event selection logic "Channel",
            # and (2) the trigger pulse generation "Trigger".  The
            # event selection logic consists of two parts: (1)
            # choosing either a rate marker or sequence bit, and (2)
            # optionally including (logical AND of) a selection on the
            # presence of beam (irrelevant for you).  The DestType is
            # this optional selection on presence of beam (to a
            # destination), and you are selecting "All" which is
            # better described as "DontCare". [The other options
            # should be "Inclusive" and "Exclusive"] The trigger pulse
            # generation has configuration of which "Channel" to
            # listen to, the delay, width, and polarity of the trigger
            # to generate.  The "Enable" registers just turn on each
            # of these two components.

            #  EvrV2CoreTriggers EvrV2ChannelReg[0] EnableReg True
            self.set_evr_channel_reg_enable(0, True, write_log=write_log)
            #  EvrV2CoreTriggers EvrV2ChannelReg[0] DestType All
            self.set_evr_trigger_dest_type(0, 0, write_log=write_log)
            #  EvrV2CoreTriggers EVrV2TriggerReg[0] Enable Trig True
            self.set_trigger_enable(0, True, write_log=write_log)

            # Select external flux ramp trigger source
            self.set_ramp_start_mode(1, write_log=write_log)

            # Configure LMK
            for bay in self.bays:
                self.set_lmk_reg(bay, 0x146, 0x08, write_log=write_log)
                self.set_lmk_reg(bay, 0x147, 0x0A, write_log=write_log)

            # Check if receiving timing.  Checks if FPGA recovered
            # clock is receiving timing data.  Warn if we've just
            # configured the system for backplane or fiber timing but
            # not receiving any external timing data.  Have to wait a
            # second for timing to come in after crossbar
            # configuration change.
            time.sleep(1)
            status = self.get_timing_link_up()
            if status != 1:
                self.log(f'\033[91mConfigured for {mode} timing but not receiving external timing data after waiting 1 second.\033[0m',
                         self.LOG_ERROR) # color red