tod_ops submodule

This submodule includes functions for processing time-ordered signals, including time and Fourier domain filtering, PCA, gap-filling, etc.

Fourier space filters

Applying Fourier space filters is facilitated by the fourier_filter function:

sotodlib.tod_ops.fourier_filter(tod, filt_function, detrend=None, resize='zero_pad', axis_name='samps', signal_name='signal', time_name='timestamps', **kwargs)

Return a filtered tod.signal_name along the axis axis_name.

Does not change the data in the axis manager.

Parameters:

• tod – axis manager

• filt_function – function( freqs, tod ) function that takes a set of frequencies and the axis manager and returns the filter in fourier space

• detrend – Method of detrending to be done before ffting. Can be ‘linear’, ‘mean’, or None. Note that detrending here can be slow for large arrays

• resize – How to resize the axis to increase fft speed. ‘zero_pad’ will increase to the next nice number (a product of small primes compatible with the FFT implementation). ‘trim’ will eliminate samples from the end so that axis has a nice length for FFTs. None will not change the axis length and might be quite slow. Trim will be kinda weird here, because signal will not be returned as the same size as it is input

• axis_name – name of axis you would like to fft along

• signal_name – name of the variable in tod to fft

• time_name – name for getting time of data (in seconds) from tod

Returns:

filtered tod.signal_name

Return type:

signal

Here’s an example, assuming you have time-ordered data in an AxisManager tod. Note how we combine two Fourier space filters using the multiplication operator:

# Filter with a smooth-edged band, passing 1 Hz to 10 Hz.
filt = tod_ops.filters.high_pass_sine(1.) * tod_ops.filters.low_pass_sine2(10.)
tod_ops.filters.fourier_filter(tod, filt)

Several filters are described in sotodlib.tod_ops.filters. See the source code for conventions on defining and wrapping new filters to be compatible with this system.

tod_ops.filters

The composable filters that can be passed into fourier_filter should be auto-documented here.

sotodlib.tod_ops.filters.low_pass_butter4(fc)

4th-order low-pass filter with f3db at fc (Hz).

sotodlib.tod_ops.filters.high_pass_butter4(fc)

4th-order high-pass filter with f3db at fc (Hz).

sotodlib.tod_ops.filters.timeconst_filter(timeconst=None, invert=False)

One-pole time constant filter for fourier_filter.

Builds filter for applying or removing time constants from signal data.

Parameters:

• timeconst – Array of time constant values (one per detector). Alternately, a string indicating what member of tod to use for the time constants array. Defaults to ‘timeconst’.

• invert (bool) – If true, returns the inverse transfer function, to deconvolve the time constants.

Example:

# Deconvolve time constants.
fourier_filter(tod, timeconst_filter(invert=True),
               detrend='linear', resize='zero_pad')

sotodlib.tod_ops.filters.timeconst_filter_single(timeconst, invert=False)

One-pole time constant filter for fourier_filter.

This version accepts a single time constant value, in seconds. To use different time constants for each detector, see timeconst_filter.

Example:

# Apply a 1ms time constant.
fourier_filter(tod, timeconst_filter_single(timeconst=0.001),
               detrend=None, resize='zero_pad')

sotodlib.tod_ops.filters.tau_filter(tau_name='timeconst', do_inverse=True)

tau_filter is deprecated; use timeconst_filter.

sotodlib.tod_ops.filters.gaussian_filter(fc=0.0, f_sigma=None, gain=1.0, t_sigma=None)

Gaussian bandpass filter

Parameters:

• fc (float0) – Central frequency of the filter (peak of passband), in Hz.

• f_sigma (float) – Standard deviation of the filter kernel, in Hz.

• gain (float) – Gain of the filter.

• t_sigma (float) – Instead of f_sigma, set t_sigma and f_sigma = 1/(2 pi t_sigma) will be used.

The filter kernel has the shape of a normal distribution, centered on fc with standard deviation f_sigma, and peak height gain.

sotodlib.tod_ops.filters.low_pass_sine2(cutoff, width=None)

Low-pass filter. Response falls from 1 to 0 between frequencies (cutoff - width/2, cutoff + width/2), with a sine-squared shape.

sotodlib.tod_ops.filters.high_pass_sine2(cutoff, width=None)

High-pass filter. Response rises from 0 to 1 between frequencies (cutoff - width/2, cutoff + width/2), with a sine-squared shape.

sotodlib.tod_ops.filters.iir_filter(b=None, a=None, fscale=1.0, iir_params=None, invert=False)

Infinite impulse response (IIR) filter. This sort of filter is used in digital applications as a low-pass filter prior to decimation. The Smurf and MCE readout filters can both be expressed in this form.

Parameters:

• b – numerator polynomial filter coefficients (z^0,z^1, …)

• a – denominator coefficients

• fscale – scalar used to compute z = exp(-2j*pi*freqs*fscale). In general this should be equal to 1/f_orig, where f_orig is the original sampling frequency (before downsampling).

• iir_params – Alternative way to specify b, a, and fscale (see notes).

• invert – If true, returns denom/num instead of num/denom.

Notes

The b and a coefficients are as implemented in scipy.signal.freqs, scipy.signal.butter, etc. The “angular frequencies”, w, are computed as 2*pi*freqs*fscale.

If the filter parameters (b, a, fscale) are not passed in explicitly, they will be extracted from an AxisManager based on the argument iir_params, which must be a dict or AxisManager with keys “b”, “a”, and “fscale”, but note that:

• If iir_params is a string, tod[iir_params] is used (and must be an AxisManager or dict).

• If iir_params is None, that’s the same as passing iir_params=’iir_params’.

Some types of low, high, and band pass filters can be derived by wrapper functions below.

sotodlib.tod_ops.filters.get_lpf(cfg)

Returns a low-pass filter based on the configuration.

Parameters:

cfg (dict) – A dictionary containing the low-pass filter configuration. It must have the following keys: - “type”: A string specifying the type of low-pass filter. Supported values are “identity”, “butter4” and “sine2”. - “cutoff”: A float specifying the cutoff frequency of the low-pass filter. - “trans_width”: A float specifying the transition width of the low-pass filter (only for “sine2” type).

Returns:

the low-pass filter.

Return type:

filters.fourier_filter

sotodlib.tod_ops.filters.get_hpf(cfg)

Returns a high-pass filter based on the configuration.

Parameters:

cfg (dict) – A dictionary containing the high-pass filter configuration. It must have the following keys: - “type”: A string specifying the type of high-pass filter. Supported values are “identity”, “butter4” and “sine2”. - “cutoff”: A float specifying the cutoff frequency of the high-pass filter. - “trans_width”: A float specifying the transition width of the high-pass filter (only for “sine2” type).

Returns:

the high-pass filter.

Return type:

filters.fourier_filter

sotodlib.tod_ops.filters.get_bpf(cfg)

Returns a band-pass filter based on the configuration.

Parameters:

cfg (dict) – A dictionary containing the band-pass filter configuration. It must have the following keys: - “type”: A string specifying the type of band-pass filter. Supported values are “identity”, “butter4” and “sine2”. - “center”: A float specifying the center frequency of the band-pass filter. - “width”: A float specifying the width of the band-pass filter. - “trans_width”: A float specifying the transition width of the band-pass filter (only for “sine2” type).

Returns:

the band-pass filter.

Return type:

filters.fourier_filter

tod_ops.gapfill

Tutorial

The gapfill submodule includes functions and classes for patching short segments of a TOD signal. For example, after identifying glitches in a timestream, one might want to eliminate those pathological data points by simply interpolating between good data segments at each end of the bad sample range.

The get_gap_fill function will fit low-order polynomial models to flagged data segments, and return vectors with modeled data that can be used to fill those gaps. For example, suppose we have an AxisManager tod pre-populated with signal and flags of some kind:

>>> tod.signal.shape
(619, 165016)
>>> tod.source_flags
RangesMatrix(619,165016)
>>> gaps = tod_ops.get_gap_fill(tod, flags=tod.source_flags)
>>> gaps
ExtractMatrix(619,165016@2.0%)
>>> gaps.swap(tod)

Note that by default the get_gap_fill function does not patch the signal data, and that is why we called gaps.swap(tod). (See the swap argument, however.) The object returned by get_gap_fill is an ExtractMatrix (see reference below). This contains only the gap-filled data – i.e. new values for each sample marked in tod.source_flags. The repr string tells us that the ExtractMatrix is compatible with a signal array of shape (619, 165016), and that it is only actually storing data for 2.0% of the full 2d matrix. The ExtractMatrix and Extract classes have several methods to help move data back and forth between compressed and expanded representation; see the full class reference below.

Superior gap-filling can sometimes be achieved with a mode-based model. For example, when there are large gaps in a single detector timestream, the signal common mode might be a better model for filling the gap than a polynomial interpolation based on only the single detector’s good amples. So if one has a model, such as the kind returned by the tod_ops.pca.get_pca_model function, the function get_gap_model can be used to lookup the values of the model at some set of flagged samples:

>>> model
AxisManager(weights[dets,eigen], modes[eigen,samps],
    dets:LabelAxis(619), eigen:IndexAxis(10),
    samps:OffsetAxis(165016))
>>> gaps = tod_ops.get_gap_model(tod, model, flags=tod.source_flags)
>>> gaps
ExtractMatrix(619,165016@2.0%)
>>> gaps.swap(tod)

Note that get_gap_model returns the same sort of object as get_gap_fill, and again in this example we have swapped the model into the tod.

Reference

Class and function references should be auto-generated here.

sotodlib.tod_ops.get_gap_fill(tod, nbuf=10, order=1, swap=False, signal=None, flags=None, _method='fast')

See get_gap_fill_single for meaning of arguments not described here.

Parameters:

• tod – AxisManager with (dets, samps) axes.

• signal – signal to pass to get_gap_fill_single as data argument; defaults to tod.signal

• flags – flags to pass to get_gap_fill_single; defaults to tod.flags.

Returns:

The ExtractMatrix object with per-detector Extracts from get_gap_fill_single.

sotodlib.tod_ops.get_gap_fill_single(data, flags, nbuf=10, order=1, swap=False)

Computes samples to fill the gaps in data identified by flags. Each flagged segment is modeled with a polynomial of the order specified, based on up to nbuf points on each side of the segment.

Parameters:

• data – 1d vector of samples.

• flags – Ranges object consistent with data size.

• nbuf – Maximum number of samples on each side of a flagged span to use for model fit. If the unflagged span between two flagged spans is shorter than nbuf, then it is still used to anchor one end of the fit.

• order – Maximum order of polynomial model used in interpolation. The actual order may be smaller if insufficient data are available to constrain the coefficients. Care should be taken when going higher than linear (order=1); this will tend to be unstable unless the gaps are much smaller than the anchors used to constrain the poly on each end.

• swap – If False, do not modify the input data vector. If True, patch the data with the model.

Returns:

An Extract object containing the modeled data. If swap is True, then the input data vector is patched with the model and the returned object, instead, contains the samples from data that were changed.

sotodlib.tod_ops.get_gap_model(tod, model, flags=None, weights=None, modes=None)

Calls get_gap_model_single on each detector and its corresponding model weights.

Parameters:

• tod – AxisManager with (dets, samps) axes.

• model – AxisManager with (dets, eigen, samps) axes.

• flags – flags to pass to get_gap_model_single; defaults to tod.flags.

• weights – array of mode couplings (dets, eigen); defaults to model.weights.

• modes – array with time-dependent modes (eigen, samps); defaults to model.modes.

Returns:

The ExtractMatrix object with per-detector models from get_gap_model_single.

sotodlib.tod_ops.get_gap_model_single(weights, modes, flags)

Computes samples to fill gaps in data identified by flags, based on weights and modes (such as would be contained in a PCA model).

Parameters:

• weights – 1-d array of weights (n_mode) to apply to each mode.

• modes – 2-d array of modes (n_mode, n_samps)

• flags – so3g.proj.Ranges object with count == n_samps.

Returns:

An Extract object containing the modeled data.

class sotodlib.tod_ops.gapfill.ExtractMatrix(items=None)

This class is a simple container for a list (of length n_dets) of Extract objects of equal length (n_samps), to provide abstract access to data with shape (n_dets, n_samps). Each child Extract can be accessed with [item] indexing. Most methods have (tod, signal=None) signature, and are simple loops that call the method of the same name on each child Extract object.

__init__(items=None)

The complete list of Extract objects should be passed in.

expand(fill_value=0)

Expands the extract into a full-size data array, filling missing bits with fill_value.

swap(tod, signal=None)

Swaps the current extract with the tracked samples of full data vector signal.

patch(tod, signal=None)

Copies the extract into the signal vector. Untracked samples are not modified.

accumulate(tod, signal, scale)

Adds the extract, scaled by factor scale, into the full length signal vector. Untracked samples are not modified.

class sotodlib.tod_ops.gapfill.Extract(ranges, init_data=True)

Container for storage of sparse sub-segments of a vector. This ties together a Ranges object and a data array with the same length as the total length of the Ranges segments.

This is useful in cases where a relatively small number of samples need to be excised from a signal vector, but stored and perhaps restored at a later time.

In docstrings below a “full data vector” is a 1-d array with self.ranges.count elements; the “tracked samples” are a subset consisting of self.ranges.mask.sum() elements.

ranges

The so3g.proj.Ranges object mapping tracked samples into full data vector.

n_ex

The number of tracked samples.

data

The data vector containing tracked samples.

__init__(ranges, init_data=True)

Constructor.

Parameters:

• ranges – a Ranges object; all positive ranges will be tracked.

• init_data – array or boolean. If False, self.data is initialized to None. If True, self.data is initialized to zeros. Otherwise, init_data is stored in self.data. Note this is a reference, not a copy.

offset_iter()

Returns an iterator for use in expanding / collapsing extracts. Each iteration returns a tuple (ex_lo, ex_hi, full_lo, full_hi) where ex_lo:ex_hi is a range in extracted data vector, and full_lo:full_hi is a range in the full data vector.

expand(fill_value=0)

Expands the extract into a full-length data array, filling missing bits with fill_value.

swap(signal)

Swaps the current extract with the tracked samples of full data vector signal.

patch(signal)

Copies the extract into the signal vector. Untracked samples are not modified.

accumulate(signal, scale)

Adds the extract, scaled by factor scale, into the full length signal vector. Untracked samples are not modified.

sotodlib.tod_ops.gapfill.get_contaminated_ranges(good_flags, bad_flags)

Determine what intervals in good_flags are contaminated (overlap with) intervals in bad_flags. Note this isn’t as simple as good_flags * bad_flags, because any contiguous region in good_flags is considered contaminated if even one sample of it is touched by bad_flags.

Parameters:

• good_flags (RangesMatrix) – The flags to check for contamination. Must have shape (dets, samps).

• bad_flags (RangesMatrix) – The flags marking bad data, which should be considered as a contaminant to good_flags. Same shape as good_flags.

Returns:

RangesMatrix with same shape as inputs, indicating the intervals from good_flags that overlap at all with some interval of bad_flags.

Example

Move contaminated intervals into bad_flags:

contam = get_contaminated_ranges(source_flags, glitch_flags)
source_flags *= ~contam
glitch_flags += contam

tod_ops.pca

Support for principal component analysis, useful in the signal-processing context, is provided in this submodule.

sotodlib.tod_ops.pca.get_pca_model(tod=None, pca=None, n_modes=None, signal=None, wrap=None)

Convert a PCA decomposition into the signal basis, i.e. into time-dependent modes that one might use for cleaning or calibrating.

The generalization of “common mode” computation is to convert the eigen-decomposition of the covariance matrix into a limited set of modes that explain a high fraction of the power in the input signals.

Here we select the strongest eigenmodes from the PCA and compute the time-dependent modes and the projection of each detector signal onto each mode. An approximation for the input signal may then be computed as

signal_model = weights . modes

where weights has shape (dets, eigen), modes has shape (eigen, samps), and . indicates matrix multipliaction. The size of the eigen axis can be as large as the number of dets, but is typically smaller and in this routine is set by n_modes.

Parameters:

• tod – AxisManager with dets and samps axes.

• pca – AxisManager holding a PCA decomposition, as returned by get_pca. If not specified, it will be computed from the data.

• n_modes – integer specifying the number of modes to compute; the strongest modes are taken and this sets the size of the “eigen” axis in the output. Defaults to len(dets), but beware that the resulting data object will be the same size as the input signal.

• signal – array of shape (dets, samps) that is used to construct the requested eigen modes. If pca is not passed in, this signal is also used to compute the covariance for PCA.

• wrap – string; if specified then the returned result is also stored in tod under that name.

Returns:

An AxisManager with (dets, eigen, samps) axes. The field ‘weights’ has shape (dets, eigen) and the field ‘modes’ has shape (eigen, samps).

sotodlib.tod_ops.pca.get_pca(tod=None, cov=None, signal=None, wrap=None)

Compute a PCA decomposition of the kind useful for signal analysis. A symmetric non-negative matrix cov of shape(n_dets, n_dets) can be decomposed into matrix R (same shape) and vector E (length n_dets) such that

cov = R . diag(E) . R^T

with . denoting matrix multiplication and T denoting matrix transposition.

Parameters:

• tod – AxisManager with dets and samps axes.

• cov – covariance matrix to decompose; if None then cov is computed from tod.signal (or signal).

• signal – array of shape (dets, samps). If cov is not provided, it will be computed from this matrix. Defaults to tod.signal.

• wrap – string; if set then the returned result is also stored in tod under this name.

Returns:

AxisManager with axes ‘dets’ and ‘eigen’ (of the same length), containing fields ‘R’ of shape (dets, eigen) and ‘E’ of shape (eigen). The eigenmodes are sorted from strongest to weakest.

sotodlib.tod_ops.pca.add_model(tod, model, scale=1.0, signal=None, modes=None, weights=None)

Adds modeled modes, multiplied by some scale factor, into signal.

Given a matrix of weights (dets, eigen) and modes (eigen, samps), this computes model signal:

model_signal = weights . modes

and adds the result to signal, scaled by scale:

signal += scale * (weights . modes)

The intended use is for PCA or common mode removal, for which user should pass scale = -1.

Parameters:

• tod – AxisManager with dets and samps axes.

• model – AxisManager with shape (dets, eigen, dets) containing ‘modes’ and ‘weights’ arrays. This is the type of object returned by get_pca_model, for example.

• scale – Factor by which to scale model_signal before accumulating into signal.

• signal – array into which to accumulate model_signal; defaults to tod.signal.

• modes – array of modes with shape (eigen, samps), defaults to model.modes.

• weights – array of weights with shape (dets, eigen), defaults to model.weights.

Returns:

The signal array into which modes were accumulated.

Notes

If you only want to operate on a subset of detectors, the most efficient solution is to pass in weights explicitly, forced to 0 for the dets you want to omit. There’s an optimization to skip the computation for a detector if all the mode coupling weights are zero.

sotodlib.tod_ops.pca.get_trends(tod, remove=False, size=1, signal=None)

Computes trends for each detector signal that remove the slope connecting first and last points, as well as the mean of the signal. The returned object can be treated like PCA model (e.g., it can be passed as the model input to add_model).

Parameters:

• tod – AxisManager with dets and samps axes.

• remove – boolean, if True then the computed trends (and means) are removed from the signal.

• size – the number of samples on each end of the signal to use for trend level computation. Defaults to 1.

• signal – array of shape (dets, samps) to compute trends on. Defaults to tod.signal.

Returns:

An AxisManager with (dets, eigen, samps) axes. The field ‘weights’ has shape (dets, eigen) and the field ‘modes’ has shape (eigen, samps). There are two modes, which always have the same form: index 0 is all ones, and index1 is a smooth line from -0.5 to +0.5.

tod_ops.detrend

Here are some detrending functions. Note detrending is applied automatically by sotodlib.tod_ops.fourier_filter(), by default, using this submodule.

sotodlib.tod_ops.detrend.detrend_tod(tod, method='linear', axis_name='samps', signal_name='signal', in_place=True, wrap_name=None, count=10)

Returns detrended data. Detrends data in place by default but pass in_place=False if you would like a copied array (such as if you’re just looking to use this in an FFT).

Using this with method =’mean’ and axis_name=’dets’ will remove a common mode from the detectors Using this with method =’median’ and axis_name=’dets’ will remove a common mode from the detectors with the median rather than the mean

Parameters:

• tod (axis manager)

• method (str) – method of detrending can be ‘linear’, ‘mean’, or median

• axis_name (str) – the axis along which to detrend. default is ‘samps’

• signal_name (str) – the name of the signal to detrend. defaults to ‘signal’. Can have any shape as long as axis_name can be resolved.

• in_place (bool.) – If False it makes a copy of signal before detrending and returns the copy.

• wrap_name (str or None.) – If not None, wrap the detrended data into tod with this name.

• count (int) – Number of samples to use, on each end, when measuring mean level for ‘linear’ detrend. Values larger than 1 suppress the influence of white noise.

Returns:

signal – Detrended signal. Done in place or on a copy depend on in_place argument.

Return type:

array of type tod[signal_name]

tod_ops.jumps

Functions to find and fix jumps in data. There are currently three jump finders included here each targeting a different use case:

• sotodlib.tod_ops.twopi_jumps(): For jumps that are a multiple of 2pi, these are a byproduct of the readout.

• sotodlib.tod_ops.slow_jumps(): For things that look like jumps when zoomed out, but happen on longer time scales (ie: a short unlock, timing issues, etc.).

• sotodlib.tod_ops.find_jumps(): For all other jumps, note here that we expect these to be “true” jumps and happen in a small number of samples.

When running multiple jump finders in a row it is advisible for performance reasons to cache the noise level used to compute the thresholds (ie: the output of sotodlib.tod_ops.std_est()) to avoid recomputing it. This can then be used to pass in variables like min_size and abs_thresh. Another way to increase performance is the set the NUM_FUTURES system variable to control how many futures are used to parallelize jump finding and fixing (currently only used by sotodlib.tod_ops.find_jumps() and sotodlib.tod_ops.jumpfix_subtract_heights()). Remember that at a certain point the overhead from increased parallelization will be greater than the speedup, so use with caution.

Historically the jump finder included options to do multiple passes with recursion and/or iteration. This was because in very jumpy data smaller jumps can be hard enough to find smaller jumps near large jumps. This was removed for the practical reason that the detectors that this helps with usually have low enough data quality that they should be cut anyways. If you would like to use iterative jump finding you can do so with a for loop or similar but remember to jumpfix and gapfill between each itteration.

sotodlib.tod_ops.jumps.std_est(x: ndarray[Any, dtype[floating]], ds: int = 1, axis: int = -1, method: str = 'median_unbiased') → ndarray[Any, dtype[floating]]

Estimate white noise standard deviation of data. More robust to jumps and 1/f then np.std()

Parameters:

• x – Data to compute standard deviation of.

• ds – Downsample factor to use, does a naive slicing.

• axis – The axis to compute along.

• method – The method to pass to np.quantile.

Returns:

The estimated white noise standard deviation of x.

Return type:

stdev

sotodlib.tod_ops.jumps.jumpfix_subtract_heights(x: ndarray[Any, dtype[floating]], jumps: RangesInt32 | RangesMatrix | ndarray[Any, dtype[bool_]], inplace: bool = False, heights: ndarray[Any, dtype[floating]] | csr_array | None = None, **kwargs) → ndarray[Any, dtype[floating]]

Naive jump fixing routine where we subtract known heights between jumps. Note that you should exepect a glitch at the jump locations. Works best if you buffer the jumps mask by a bit.

Parameters:

• x – Data to jumpfix on, expects 1D or 2D.

• jumps – Boolean mask or Ranges(Matrix) of jump locations. Should be the same shape at x.

• inplace – Whether of not x should be fixed inplace.

• heights – Array of jump heights, can be sparse. If None will be computed.

• **kwargs – Additional arguments to pass to estimate_heights if heights is None.

Returns:

x with jumps removed.

If inplace is True this is just a reference to x.

Return type:

x_fixed

sotodlib.tod_ops.jumps.estimate_heights(signal: ndarray[Any, dtype[floating]], jumps: ndarray[Any, dtype[bool_]], win_size: int = 20, twopi: bool = False, make_step: bool = False, diff_buffed: ndarray[Any, dtype[floating]] | None = None) → ndarray[Any, dtype[floating]]

Simple jump estimation routine.

Parameters:

• signal – The signal with jumps.

• jumps – Boolean mask of jump locations in signal.

• win_size – Number of samples to buffer when estimating heights.

• twopi – If True, heights will be rounded to the nearest 2*pi

• make_step – If True jump ranges will be turned into clean step functions.

• diff_buffed – Difference between signal and a signal shifted by win_size. If None will be computed.

Returns:

Array of jump heights.

Return type:

heights

sotodlib.tod_ops.jumps.twopi_jumps(aman, signal=None, win_size=20, nsigma=5.0, atol=None, fix: Literal[True] = True, inplace=False, merge=True, overwrite=False, name='jumps_2pi', **filter_pars) → Tuple[RangesMatrix, csr_array, ndarray[Any, dtype[floating]]]

sotodlib.tod_ops.jumps.twopi_jumps(aman, signal=None, win_size=20, nsigma=5.0, atol=None, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps_2pi', **filter_pars) → Tuple[RangesMatrix, csr_array]

Find and optionally fix jumps that are height ~N*2pi. TOD is expected to have detectors with high trends already cut. Data is assumed to have units of phase here.

Parameters:

• aman – The axis manager containing signal to find jumps on.

• signal – Signal to jumpfind on. If None than aman.signal is used.

• win_size – Size of window to use when looking for jumps. This should be set to something of order the width of the jumps.

• nsigma – How many multiples of the white noise level to set to use to compute atol. Only used if atol is None.

• atol – How close the jump height needs to be to N*2pi to count. If set to None, then nsigma times the WN level of the TOD is used. Note that in general this is faster than nsigma.

• fix – If True the jumps will be fixed by adding N*2*pi at the jump locations.

• inplace – If True jumps will be fixed inplace.

• merge – If True will wrap ranges matrix into aman.flags.<name>

• overwrite – If True will overwrite existing content of aman.flags.<name>

• name – String used to populate field in flagmanager if merge is True.

• **filter_pars – Parameters to pass to _filter

Returns:

RangesMatrix containing jumps in signal,

if signal is 1D Ranges in returned instead. Buffered to win_size.

heights: csr_array of jump heights.

fixed: signal with jump fixed. Only returned if fix is set.

Return type:

jumps

sotodlib.tod_ops.jumps.slow_jumps(aman, signal=None, win_size=800, thresh=20, abs_thresh=True, fix: Literal[True] = True, inplace=False, merge=True, overwrite=False, name='jumps_slow', **filter_pars) → Tuple[RangesMatrix, csr_array, ndarray[Any, dtype[floating]]]

sotodlib.tod_ops.jumps.slow_jumps(aman, signal=None, win_size=800, thresh=20, abs_thresh=True, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps_slow', **filter_pars) → Tuple[RangesMatrix, csr_array]

Find and optionally fix slow jumps. This is useful for catching things that aren’t really jumps but change the DC level in a jumplike way (ie: a short unlock period).

Parameters:

• aman – The axis manager containing signal to find jumps on.

• signal – Signal to jumpfind on. If None than aman.signal is used.

• win_size – Size of window to use when looking for jumps. This should be set to something of order the width of the jumps.

• thresh – ptp value at which to flag things as jumps. Default value is in phase units. You can also pass a percentile to use here if abs_thresh is False.

• abs_thresh – If True thresh is an absolute threshold. If False it is a percential in range [0, 1).

• fix – If True the jumps will be fixed by adding N*2*pi at the jump locations.

• inplace – If True jumps will be fixed inplace.

• merge – If True will wrap ranges matrix into aman.flags.<name>

• overwrite – If True will overwrite existing content of aman.flags.<name>

• name – String used to populate field in flagmanager if merge is True.

• **filter_pars – Parameters to pass to _filter

Returns:

RangesMatrix containing jumps in signal,

if signal is 1D Ranges in returned instead. Buffered to win_size.

heights: csr_array of jump heights.

fixed: signal with jump fixed. Only returned if fix is set.

Return type:

jumps

sotodlib.tod_ops.jumps.find_jumps(aman, signal=None, max_iters=1, min_sigma=None, min_size=None, win_size=20, nsigma=25, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps', **filter_pars) → Tuple[RangesMatrix, csr_array]

sotodlib.tod_ops.jumps.find_jumps(aman, signal=None, max_iters=1, min_sigma=None, min_size=None, win_size=20, nsigma=25, fix: Literal[True] = True, inplace=False, merge=True, overwrite=False, name='jumps', **filter_pars) → Tuple[RangesMatrix, csr_array, ndarray[Any, dtype[floating]]]

Find jumps in aman.signal_name with a matched filter for edge detection. Expects aman.signal_name to be 1D of 2D.

Parameters:

• aman – axis manager.

• signal – Signal to jumpfind on. If None than aman.signal is used.

• min_sigma – Number of standard deviations to count as a jump, note that the standard deviation here is computed by std_est and is the white noise standard deviation, so it doesn’t include contributions from jumps or 1/f. If min_size is provided it will be used instead of this.

• min_size – The smallest jump size counted as a jump. By default this is set to None and min_sigma is used instead, if set this will override min_sigma. If both min_sigma and min_size are None then the IQR is used as min_size.

• win_size – Size of window used by SG filter when peak finding. Also used for height estimation, should be of order jump width.

• nsigma – Number of sigma above the mean for something to be a peak.

• fix – Set to True to fix.

• inplace – Whether of not signal should be fixed inplace.

• merge – If True will wrap ranges matrix into aman.flags.<name>

• overwrite – If True will overwrite existing content of aman.flags.<name>

• name – String used to populate field in flagmanager if merge is True.

• **filter_pars – Parameters to pass to _filter

Returns:

RangesMatrix containing jumps in signal,

if signal is 1D Ranges in returned instead. There is some uncertainty on order of a few samples. Jumps within a few samples of each other may not be distinguished.

heights: csr_array of jump heights.

fixed: signal with jump fixed. Only returned if fix is set.

Return type:

jumps

sotodlib.tod_ops.jumps.jumps_aman(aman: AxisManager, jump_ranges: RangesMatrix, heights: csr_array) → AxisManager

Helper to wrap the jumpfinder outputs into a AxisManager for use with preprocess.

Parameters:

• aman – AxisManager to steam axis information from.

• jump_ranges – RangesMatrix containing the jump flag.

• heights – csr_array of jump heights.

Returns:

AxisManager containing the jump information wrapped in

’jump_flag’ and ‘jump_heights’.

Return type:

jumps_aman

tod_ops.azss

Function for binning signal by azimuth and fitting it with legendre polynomials in az-vs-signal space. This has been used in ABS and other experiments to get an Azimuth synchronous signal (AzSS) template largely due to polarized ground signal to remove from the data.

sotodlib.tod_ops.azss.get_azss(aman, signal=None, az=None, range=None, bins=100, flags=None, method='interpolate', max_mode=None, subtract_in_place=False, merge_stats=True, azss_stats_name='azss_stats', merge_model=True, azss_model_name='azss_model')

Derive azss (Azimuth Synchronous Signal) statistics and model from the given axismanager data. NOTE: This function does not modify the signal unless subtract_in_place = True.

Parameters:

• aman (TOD) – core.AxisManager

• signal (array-like, optional) – A numpy array representing the signal to be used for azss extraction. If not provided, the signal is taken from aman.signal.

• az (array-like, optional) – A 1D numpy array representing the azimuth angles. If not provided, the azimuth angles are taken from aman.boresight.az.

• range (list, optional) – A list specifying the range of azimuth angles to consider for binning. Defaults to [-np.pi, np.pi]. If None, [min(az), max(az)] will be used for binning.

• bins (integer) – The number of bins to divide the azimuth angle range into. Defaults to 360.

• flags (RangesMatrix, optinal) – Flag indicating whether to exclude flagged samples when binning the signal. Default is no mask applied.

• method (str) – The method to use for azss modeling. Options are ‘interpolate’ and ‘fit’. In ‘interpolate’, binned signal is used directly. In ‘fit’, fitting is applied to the binned signal. Defaults to ‘interpolate’.

• max_mode (integer, optinal) – The number of Legendre modes to use for azss when method is ‘fit’. Required when method is ‘fit’.

• subtract_in_place (bool) – If True, it subtract the modeled tod from original signal. The aman.signal will be modified.

• merge_stats (boolean, optional) – Boolean flag indicating whether to merge the azss statistics with aman. Defaults to True.

• azss_stats_name (string, optional) – The name to assign to the merged azss statistics. Defaults to ‘azss_stats’.

• merge_model (boolean, optional) – Boolean flag indicating whether to merge the azss model with the aman. Defaults to True.

• azss_model_name (string, optional) – The name to assign to the merged azss model. Defaults to ‘azss_model’.

Returns:

• azss_stats: core.AxisManager

  • azss statistics including: azumith bin centers, binned signal, std of each detector-az bin, std of each detector.

  • If method=fit then also includes: binned legendre model, legendre bin centers, fit coefficients, reduced chi2.

• model_sig_tod: numpy.array

  • azss model as a function of time either from fits or interpolation depending on method argument.

Return type:

Tuple

sotodlib.tod_ops.azss.subtract_azss(aman, signal=None, azss_template=None, subtract_name='azss_remove')

Subtract the scan synchronous signal (azss) template from the signal in the given axis manager.

Parameters:

• aman (AxisManager) – The axis manager containing the signal to which the azss template will be applied.

• signal (ndarray, optional) – The signal from which the azss template will be subtracted. If signal is None (default), the signal contained in the axis manager will be used.

• azss_template (ndarray, optional) – The azss template to be subtracted from the signal. If azss_template is None (default), the azss template stored in the axis manager under the key azss_extract will be used.

• subtract_name (str, optional) – The name of the output axis manager field that will contain the azss-subtracted signal. Defaults to azss_remove.

Return type:

None

tod_ops.binning

Function for binning signal along specified axis (i.e. azimuth, time, hwp angle).

sotodlib.tod_ops.binning.bin_signal(aman, bin_by, signal=None, range=None, bins=100, flags=None)

Bin time-ordered data by the bin_by and return the binned signal and its standard deviation.

Parameters:

• aman (TOD) – The Axismanager object to be binned.

• bin_by (array-like) – The array by which signal is binned. It should has the same samps as aman.

• signal (str, optional) – The name of the signal to be binned. Defaults to aman.signal if not specified.

• range (list or None) – A list specifying the bin range ([min, max]). Default is None, which means bin range is set to [min(bin_by), max(bin_by)].

• bins (int or sequence of scalars) – If bins is an int, it defines the number of equal-width bins in the given range (100, by default). If bins is a sequence, it defines the bin edges, including the rightmost edge, allowing for non-uniform bin widths. If bins is a sequence, bins overwrite range.

• flags (RangesMatrix, optional) – Flag indicating whether to exclude flagged samples when binning the signal. Default is no mask applied.

Returns:

• bin_edges (dict key): float array of bin edges length(bin_centers)+1.

• bin_centers (dict key): center of each bin.

• binned_signal (dict key): binned signal.

• binned_signal_sigma (dict key): estimated sigma of binned signal.

Return type:

Dictionary

tod_ops.sub_polyf

Function for remove low order polynominal component in each subscan.

sotodlib.tod_ops.sub_polyf.subscan_polyfilter(aman, degree, signal=None, exclude_turnarounds=False, mask=None, in_place=True)

Apply polynomial filtering to subscan segments in a data array. This function applies polynomial filtering to subscan segments within signal for each detector. Subscan segments are defined based on the presence of flags such as ‘left_scan’ and ‘right_scan’. Polynomial filtering is used to remove low-degree polynomial trends within each subscan segment.

Parameters:

• aman (AxisManager)

• degree (int) – The degree of the polynomial to be removed.

• signal (array-like, optional) – The TOD signal to use. If not provided, aman.signal will be used.

• exclude_turnarounds (bool) – Optional. If True, turnarounds are excluded from subscan identification. Default is False.

• mask (str or RangesMatrix) – Optional. A mask used to select specific data points for filtering. Default is None. If None, no mask is applied. If the mask is given in str, aman.flags['mask'] is used as mask. Arbitrary mask can be specified in the style of RangesMatrix.

• in_place (bool) – Optional. If True, aman.signal is overwritten with the processed signal.

Returns:

signal – The processed signal.

Return type:

array-like

tod_ops.flags

Module containing functions for generating flags for cuts and tod masking.

sotodlib.tod_ops.flags.get_det_bias_flags(aman, detcal=None, rfrac_range=(0.1, 0.7), psat_range=(0, 15), merge=True, overwrite=True, name='det_bias_flags', full_output=False)

Function for selecting detectors in appropriate bias range.

Parameters:

• aman (AxisManager) – Input axis manager.

• detcal (AxisManager) – AxisManager containing detector calibration information from bias steps and IVs. If None defaults to aman.det_cal.

• rfrac_range (Tuple) – Tuple (lower_bound, upper_bound) for rfrac det selection.

• psat_range (Tuple) – Tuple (lower_bound, upper_bound) for P_SAT from IV analysis. P_SAT in the IV analysis is the bias power at 90% Rn in pW.

• merge (bool) – If true, merges the generated flag into aman.

• overwrite (bool) – If true, write over flag. If false, don’t.

• name (str) – Name of flag to add to aman.flags if merge is True.

• full_output (bool) – If true, returns the full output with separated RangesMatrices

Returns:

msk_aman – AxisManager containing RangesMatrix shaped N_dets x N_samps that is True if the detector is flagged to be cut and false if it should be kept based on the rfrac, and psat ranges. To create a boolean mask from the RangesMatrix that can be used for aman.restrict() use keep = ~has_all_cut(mask) and then restrict with aman.restrict('dets', aman.dets.vals[keep]). If full_output is True, this will contain multiple RangesMatrices.

Return type:

AxisManager

sotodlib.tod_ops.flags.get_turnaround_flags(aman, az=None, method='scanspeed', name='turnarounds', merge=True, merge_lr=True, overwrite=True, t_buffer=2.0, kernel_size=400, peak_threshold=0.1, rel_distance_peaks=0.3, truncate=False, qlim=1)

Compute turnaround flags for a dataset.

Parameters:

• aman (AxisManager) – Input axis manager.

• az (Array) – (Optional). Azimuth data for turnaround flag computation. If not provided, it uses aman.boresight.az.

• method (str) – (Optional). The method for computing turnaround flags. Options are az or scanspeed.

• name (str) – (Optional). The name of the turnaround flag in aman.flags. Default is turnarounds

• merge (bool) – (Optional). Merge the computed turnaround flags into aman.flags if True.

• merge_lr (bool) – (Optional). Merge left and right scan flags as aman.flags.left_scan and aman.flags.right_scan if True.

• overwrite (bool) – (Optional). Overwrite an existing flag in aman.flags with the same name.

• t_buffer (float) – (Optional). Buffer time (in seconds) for flagging turnarounds in scanspeed method.

• kernel_size (int) – (Optional). Size of the step-wise matched filter kernel used in scanspeed method.

• peak_threshold (float) – (Optional). Peak threshold for identifying peaks in the matched filter response. It is a value used to determine the minimum peak height in the signal.

• rel_distance_peaks (float) – (Optional). Relative distance between peaks. It specifies the minimum distance between peaks as a fraction of the approximate number of samples in one scan period.

• truncate (bool) – (Optional). Truncate unstable scan segments if True in scanspeed method.

• qlim (float) – (Optional). Azimuth threshold percentile for az method turnaround detection.

Returns:

Ranges – The turnaround flags as a Ranges object.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_glitch_flags(aman, t_glitch=0.002, hp_fc=0.5, n_sig=10, buffer=200, detrend=None, signal_name=None, merge=True, overwrite=False, name='glitches', full_output=False, edge_guard=2000)

Find glitches with fourier filtering. Translation from moby2 as starting point

Parameters:

• aman (AxisManager) – The tod.

• t_glitch (float) – Gaussian filter width.

• hp_fc (float) – High pass filter cutoff.

• n_sig (int or float) – Significance of detection.

• buffer (int) – Amount to buffer flags around found location

• detrend (str) – Detrend method to pass to fourier_filter

• signal_name (str) – Field name in aman to detect glitches on if None, defaults to signal

• merge (bool)) – If true, add to aman.flags

• name (string) – Name of flag to add to aman.flags

• overwrite (bool) – If true, write over flag. If false, raise ValueError if name already exists in AxisManager

• full_output (bool) – If true, return sparse matrix with the significance of the detected glitches

• edge_guard (int) – Number of samples at the beginning and end of the tod to exclude from the returned glitch RangesMatrix. Defaults to 2000 samples (10 sec).

Returns:

flag – RangesMatrix object containing glitch mask.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_trending_flags(aman, max_trend=1.2, n_pieces=1, max_samples=500, signal=None, timestamps=None, merge=True, overwrite=True, name='trends', full_output=False)

Flag Detectors with trends larger than max_trend. This function can be used to find unlocked detectors. Note that this is a rough cut and unflagged detectors can still have poor tracking.

Parameters:

• aman (AxisManager) – The tod

• max_trend (float) – Slope at which detectors are unlocked. The default is for use with phase units.

• n_pieces (int) – Number of pieces to cut the timestream in to to look for trends.

• max_samples (int) – Maximum samples to compute the slope with.

• signal (array) – (Optional). Signal to use to generate flags, if None default is aman.signal.

• timestamps (array) – (Optional). Timestamps to use to generate flags, default is aman.timestamps.

• merge (bool) – If true, merges the generated flag into aman.

• overwrite (bool) – If true, write over flag. If false, don’t.

• name (str) – Name of flag to add to aman.flags if merge is True.

• full_output (bool) – If true, returns calculated slope sizes

Returns:

• cut (RangesMatrix) – RangesMatrix of trending regions

• trends (AxisManager) – If full_output is true, calculated slopes and the sample edges where they were calculated.

tod_ops.fft_ops

Module containing functions for FFTs and related operations.

FFTs and related operations

sotodlib.tod_ops.fft_ops.rfft(aman, detrend='linear', resize='zero_pad', window=<function hanning>, axis_name='samps', signal_name='signal', delta_t=None)

Return the real fft of aman.signal_name along the axis axis_name.

Does not change the data in the axis manager.

Parameters:

• aman – axis manager

• detrend – Method of detrending to be done before ffting. Can be ‘linear’, ‘mean’, or None. Note that detrending here can be slow for large arrays.

• resize – How to resize the axis to increase fft speed. ‘zero_pad’ will increase to the next 2**N. ‘trim’ will cut out so the factorization of N contains only low primes. None will not change the axis length and might be quite slow.

• window – a function that takes N are returns an fft window Can be None if no windowing

• axis_name – name of axis you would like to fft along

• signal_name – name of the variable in aman to fft

• delta_t – if none, it will look for ‘timestamps’ in the axis manager and will otherwise assume 1. if not None, it should be the sampling rate along the axis you’re ffting

Returns:

the fft’d data

freqs: the frequencies it is value at (since resizing is an option)

Return type:

fft

sotodlib.tod_ops.fft_ops.build_rfft_object(n_det, n, direction='FFTW_FORWARD', **kwargs)

Build PyFFTW object for fft-ing

Parameters:

• n_det – number of detectors (or just the arr.shape[0] for the array you are going to fft)

• n – number of samples in timestream

• direction – fft direction. Can be FFTW_FORWARD, FFTW_BACKWARD, or BOTH

• kwargs – additional arguments to pass to pyfftw.FFTW

Returns:

array for the real valued side of the fft

b: array for the the complex side of the fft

t_fun: function for performing FFT (two are returned if direction==’BOTH’)

Return type:

a

sotodlib.tod_ops.fft_ops.find_inferior_integer(target, primes=[2, 3, 5, 7, 11, 13])

Find the largest integer less than or equal to target whose prime factorization contains only the integers listed in primes.

sotodlib.tod_ops.fft_ops.find_superior_integer(target, primes=[2, 3, 5, 7, 11, 13])

Find the smallest integer less than or equal to target whose prime factorization contains only the integers listed in primes.

sotodlib.tod_ops.fft_ops.calc_psd(aman, signal=None, timestamps=None, max_samples=262144, prefer='center', freq_spacing=None, merge=False, overwrite=True, **kwargs)

Calculates the power spectrum density of an input signal using signal.welch(). Data defaults to aman.signal and times defaults to aman.timestamps. By default the nperseg will be set to power of 2 closest to the 1/50th of the samples used, this can be overridden by providing nperseg or freq_spacing.

Parameters:

• aman (AxisManager) – with (dets, samps) OR (channels, samps)axes.

• signal (float ndarray) – data signal to pass to scipy.signal.welch().

• timestamps (float ndarray) – timestamps associated with the data signal.

• max_samples (int) – maximum samples along sample axis to send to welch.

• prefer (str) – One of [‘left’, ‘right’, ‘center’], indicating what part of the array we would like to send to welch if cuts are required.

• freq_spacing (float) – The approximate desired frequency spacing of the PSD. If None the default nperseg of ~1/50th the signal length is used. If an nperseg is explicitly passed then that will be used.

• merge (bool) – if True merge results into axismanager.

• overwrite (bool) – if true will overwrite f, pxx axes.

• **kwargs – keyword args to be passed to signal.welch().

Returns:

array of frequencies corresponding to PSD calculated from welch. Pxx: array of PSD values.

Return type:

freqs

sotodlib.tod_ops.fft_ops.calc_wn(aman, pxx=None, freqs=None, low_f=5, high_f=10)

Function that calculates the white noise level as a median PSD value between two frequencies. Defaults to calculation of white noise between 5 and 10Hz. Defaults frequency information to a wrapped “freqs” field in aman.

Parameters:

• (AxisManager) (aman) – Uses aman.freq as frequency information associated with the PSD, pxx.

• array) (freqs (1d Float) – Psd information to calculate white noise. Defaults to aman.pxx

• array) – frequency information related to the psd. Defaults to aman.freqs

• (Float) (low_f) – low frequency cutoff to calculate median psd value. Defaults to 5Hz

• (float) (high_f) – high frequency cutoff to calculate median psd value. Defaults to 10Hz

Returns:

wn

Return type:

Float array of white noise levels for each psd passed into argument.

sotodlib.tod_ops.fft_ops.noise_model(f, p)

Noise model for power spectrum with white noise, and 1/f noise.

sotodlib.tod_ops.fft_ops.fit_noise_model(aman, signal=None, f=None, pxx=None, psdargs=None, fwhite=(10, 100), lowf=1, merge_fit=False, f_max=100, merge_name='noise_fit_stats', merge_psd=True)

Fits noise model with white and 1/f noise to the PSD of signal. This uses a MLE method that minimizes a log likelihood. This is better for chi^2 distributed data like the PSD.

Reference: http://keatonb.github.io/archivers/powerspectrumfits

Parameters:

• aman (AxisManager) – Axis manager which has samps axis aligned with signal.

• signal (nparray) – Signal sized ndets x nsamps to fit noise model to. Default is None which corresponds to aman.signal.

• f (nparray) – Frequency of PSD of signal. Default is None which calculates f, pxx from signal.

• pxx (nparray) – PSD sized ndets x len(f) which is fit to with model. Default is None which calculates f, pxx from signal.

• psdargs (dict) – Dictionary of optional argument for scipy.signal.welch

• fwhite (tuple) – Low and high frequency used to estimate white noise for initial guess passed to scipy.signal.curve_fit.

• lowf (tuple) – Frequency below which estimate of 1/f noise index and knee are estimated for initial guess passed to scipy.signal.curve_fit.

• merge_fit (bool) – Merges fit and fit statistics into input axis manager.

• f_max (float) – Maximum frequency to include in the fitting. This is particularly important for lowpass filtered data such as that post demodulation if the data is not downsampled after lowpass filtering.

• merge_name (bool) – If merge_fit is True then addes into axis manager with merge_name.

• merge_psd (bool) – If merg_psd is True then adds fres and Pxx to the axis manager.

Returns:

noise_fit_stats – If merge_fit is False then axis manager with fit and fit statistics is returned otherwise nothing is returned and axis manager is wrapped into input aman.

Return type:

AxisManager