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', rfft=None, **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

• rfft – A RFFTObj for performing the FFTs. Must be valid for the input data and contain both t_forward and t_backward.

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_sine2(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.gain(gain=1.0)

Filter that simply applies a simple gain (which could be complex) to the entire spectrum.

This can be used (with gain=1) as an identity filter for testing the filter preprocessing.

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), or a scalar value, or a string indicating what field 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”, or an AxisManager including the sub-iir_params of each stream_id. In the later case, if the filter parameters of each stream_id is different, raises an error.

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, mask=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.

• mask – If specifed, a boolean array to select which dets are to be considered in the PCA decomp. This is achieved by modifying the cov to make the non-considered dets independent and low significance.

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.

sotodlib.tod_ops.pca.pca_cuts_and_cal(tod, pca_aman, xfac=2, yfac=1.5, calc_good_medianw=False)

Finds the bounds of the pca box using IQR statistics

Parameters:

• tod (AxisManager) – observation axismanagers

• pca_aman (AxisManager) – output pca axismanager from get_pca_model

• xfac (int) – multiplicative factor for the width of the pca box. Default is 2.

• yfac (int) – multiplicative factor for the height of the box. Default is 1.5.

• calc_good_medianw (bool) – If true, the resulting median weight is calculated excluding bad dets. Default is false.

Returns:

pca_relcal – AxisManager pca and relcal information.

Return type:

AxisManager

sotodlib.tod_ops.pca.get_common_mode(tod, signal='signal', method='median', wrap=None, weights=None)

Returns common mode timestream between detectors. This uses method ‘median’ or ‘average’ across detectors as opposed to a principle component analysis to get the common mode.

Parameters:

• tod (axis manager)

• signal (str, optional) – The name of the signal to estimate common mode or ndarray with shape of (n_dets x n_samps). Defaults to ‘signal’.

• method (str) – method of common mode estimation. ‘median’ or ‘average’.

• wrap (str or None.) – If not None, wrap the common mode into tod with this name.

• weights (array with dets axis) – If not None, estimate common mode by taking average with this weights.

Return type:

common mode timestream

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, win_size: int = 20, 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 in blocks of win_size.

• win_size – Window size to downsample by.

• 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, clean: float = 80) → 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.

• clean – Make each jump range a constant height set by this percentile. If this is 0 then it is ignored.

Returns:

Array of jump heights.

Return type:

heights

sotodlib.tod_ops.jumps.twopi_jumps(aman, signal, win_size, nsigma, atol, max_tol, fix: Literal[True], inplace, merge, overwrite, name, ds, **filter_pars) → Tuple[RangesMatrix, csr_array, ndarray[Any, dtype[floating]]]

sotodlib.tod_ops.jumps.twopi_jumps(aman, signal, win_size, nsigma, atol, max_tol, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps_2pi', ds=10, **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.

• max_tol – Upper bound of the nsigma based thresh. atol ignores this.

• 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.

• ds – Downsample factor used when computing noise level, the actual factor used is ds*win_size.

• **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, win_size, thresh, abs_thresh, fix: Literal[True], inplace, merge, overwrite, name, clean, **filter_pars) → Tuple[RangesMatrix, csr_array, ndarray[Any, dtype[floating]]]

sotodlib.tod_ops.jumps.slow_jumps(aman, signal, win_size, thresh, abs_thresh, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps_slow', clean=80, **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.

• clean – Cleaning value to pass to estimate_heights. See that function for details.

• **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, min_sigma, min_size, win_size, exact, fix: Literal[False] = False, inplace=False, merge=True, overwrite=False, name='jumps', ds=10, clean=80, **filter_pars) → Tuple[RangesMatrix, csr_array]

sotodlib.tod_ops.jumps.find_jumps(aman, signal, min_sigma, min_size, win_size, exact, fix: Literal[True], inplace, merge, overwrite, name, ds, clean, **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 when peak finding. Also used for height estimation, should be of order jump width.

• exact – If True search for the exact jump location. If False flag allow some undertainty within the window (cheaper).

• 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.

• ds – Downsample factor used when computing noise level, the actual factor used is ds*win_size.

• clean – Cleaning value to pass to estimate_heights. See that function for details.

• **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='signal', az=None, azrange=None, bins=100, flags=None, scan_flags=None, apodize_edges=True, apodize_edges_samps=1600, apodize_flags=True, apodize_flags_samps=200, apply_prefilt=True, prefilt_cfg=None, prefilt_detrend='linear', method='interpolate', max_mode=None, modes_axis_name='azss_modes', subtract_in_place=False, merge_stats=True, azss_stats_name='azss_stats', merge_model=True, azss_model_name='azss_model', coverage_threshold=0.95, exclude_turnarounds=True, return_det_mask=False)

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.

• azrange (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 (int or sequence of scalars) – If bins is an int, it defines the number of equal-width bins in the given azrange (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 azrange.

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

• scan_flags (str or Ranges, optional) – Subtract in the scan/time region specified by flags. Typically flags.left_scan or flags.right_scan will be used. If we estimate and subtract azss in left going scan only, then flags should be a “union of glitch_flags and flags.right_scan (or ~flags.left_scan)”, and scan_flags should be flags.left_scan.

• apodize_edges (bool, optional) – If True, applies an apodization window to the edges of the signal. Defaults to True.

• apodize_edges_samps (int, optional) – The number of samples over which to apply the edge apodization window. Defaults to 1600.

• apodize_flags (bool, optional) – If True, applies an apodization window based on the flags. Defaults to True.

• apodize_flags_samps (int, optional) – The number of samples over which to apply the flags apodization window. Defaults to 200.

• apply_prefilt (bool, optional) – If True, applies a pre-filter to the signal before azss extraction. Defaults to True.

• prefilt_cfg (dict, optional) – Configuration for the pre-filter. Defaults to {‘type’: ‘sine2’, ‘cutoff’: 0.005, ‘trans_width’: 0.005}.

• prefilt_detrend (str, optional) – Method for detrending before filtering. Defaults to ‘linear’.

• 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, optional) – The number of Legendre modes to use for azss when method is ‘fit’. Required when method is ‘fit’.

• modes_axis_name (string, optional) – The name assigned to the LabelAxis of azss legendre modes when method is ‘fit’. Set a unique name when fitting azss with different max_mode values to avoid axis name conflicts. Defaults to ‘azss_modes’.

• 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’.

• coverage_threshold (float, optional) – Minimum azimuth coverage fraction required. Default 0.8

• exclude_turnarounds (bool, optional) – Exclude turnarounds when checking coverage. Default True

• return_det_mask (bool, optional) – If True, return detector mask along with model. Default False

Returns:

• azss_stats: core.AxisManager

  • azss statistics including: azumith bin centers, bin counts, 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.

  • If return_det_mask=True then also includes: bad_dets, coverages.

• 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, azss_stats, signal='signal', method='interpolate', max_mode=None, modes_axis_name='azss_modes', azrange=None, scan_flags=None, subtract_name='azss_remove', in_place=False, remove_template=False)

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

Parameters:

• aman (AxisManager) – The axis manager containing the signal and the azss template.

• azss_stats (AxisManager) – Contains AxisManager from get_azss.

• signal (str, optional) – The name of the field in the axis manager containing the signal to be processed. Defaults to ‘signal’.

• method (str) – The method to use for azss modeling. Options are ‘interpolate’ and ‘fit’.

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

• azrange – Azimuth range used for ‘fit’ method.

• scan_flags (str or Ranges, optional) – Subtract in the scan/time region specified by flags. Typically flags.left_scan or flags.right_scan will be used.

• subtract_name (str, optional) – The name of the field in the axis manager that will store the azss-subtracted signal. Only used if in_place is False. Defaults to ‘azss_remove’.

• in_place (bool, optional) – If True, the subtraction is done in place, modifying the original signal in the axis manager. If False, the result is stored in a new field specified by subtract_name. Defaults to False.

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, weight_for_signal=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 ((str or RangesMatrix or Ranges), optional) – Flag indicating whether to exclude flagged samples when binning the signal. If provided by a string, aman.flags.get(flags) is used for the flags. Default is no mask applied.

• weight_for_signal (array-like, optional) – Array of weights for the signal values. If None, all weights are assumed to be 1. You can get a apodizing window by ‘sotodlib.tod_ops.apodize.get_apodize_window_for_ends’ or ‘get_apodize_window_from_flags’.

Returns:

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

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

• bin_counts (dict key): counts of binned samples.

• 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_name='signal', exclude_turnarounds=False, mask=None, exclusive=True, method='legendre', 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_name (string, optional) – The name of 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. 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.

• exclusive – Optional. If True, the mask is used to exclude data from fitting. If False, the mask is used to include data for fitting. Default is True.

• method (str) – Optioal. Method to model the baseline of TOD. In legendre method, baseline model is constructed using orthonormality of Legendre function. In polyfit method, numpy.polyfit is used. legendre is faster. Default is legendre.

• 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=None, rn_range=None, si_range=None, phase_to_pW=None, 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. If None, no flags are not applied from P_SAT.

• rn_range (Tuple) – Tuple (lower_bound, upper_bound) for r_n det selection.

• si_range (Tuple) – Tuple (lower_bound, upper_bound) for s_i det selection.

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

• 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, merge_subscans=True, turnarounds_in_subscan=False)

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.

• merge_subscans (bool) – (Optional). Also merge an AxisManager with subscan information.

• turnarounds_in_subscan (bool) – (Optional). Turnarounds are included as part of a subscan.

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, subscan=False)

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).

• subscan (bool) – If True, compute the glitch threshold on a per-subscan basis. Includes turnarounds.

Returns:

flag – RangesMatrix object containing glitch mask.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_trending_flags(aman, max_trend=1.2, t_piece=500, 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.

• t_piece (float) – Duration in seconds of each 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.

sotodlib.tod_ops.flags.get_dark_dets(aman, merge=True, overwrite=True, dark_flags_name='darks')

Identify and flag dark detectors in the given aman object.

Parameters:

• aman (AxisManager) – The tod.

• merge (bool, optional) – If True, merge the dark detector flags into the aman.flags. Default is True.

• overwrite (bool, optional) – If True, overwrite existing flags with the same name. Default is True.

• dark_flags_name (str, optional) – The name to use for the dark detector flags in aman.flags. Default is ‘darks’.

Returns:

mskdarks – A matrix of ranges indicating the dark detectors.

Return type:

RangesMatrix

Raises:

ValueError – If merge is True and dark_flags_name already exists in aman.flags and overwrite is False.

sotodlib.tod_ops.flags.get_ptp_flags(aman, signal_name='signal', kurtosis_threshold=5, merge=False, overwrite=False, ptp_flag_name='ptp_flag', outlier_range=(0.5, 2.0))

Returns a ranges matrix that indicates if the peak-to-peak (ptp) of the tod is valid based on the kurtosis of the distribution of ptps. The threshold is set by kurtosis_threshold.

Parameters:

• aman (AxisManager) – The tod

• signal_name (str) – Signal to estimate flags off of. Default is signal.

• kurtosis_threshold (float) – Maximum allowable kurtosis of the distribution of peak-to-peaks. Default is 5.

• merge (bool) – Merge RangesMatrix into aman.flags. Default is False.

• overwrite (bool) – Whether to write over any existing data in aman.flags[ptp_flag_name] if merge is True. Default is False.

• ptp_flag_name (str) – Field name used when merge is True. Default is ptp_flag.

• outlier_range (tuple) – (lower, upper) bound of the initial cut before estimating the kurtosis.

Returns:

mskptps – RangesMatrix of detectors with acceptable peak-to-peaks. All ones if the detector should be cut.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_inv_var_flags(aman, signal_name='signal', nsigma=5, merge=False, overwrite=False, inv_var_flag_name='inv_var_flag')

Returns a ranges matrix that indicates if the inverse variance (inv_var) of the tod is greater than nsigma away from the median.

Parameters:

• aman (AxisManager) – The tod

• signal_name (str) – Signal to estimate flags off of. Default is signal.

• nsigma (float) – Maximum allowable deviation from the median inverse variance. Default is 5.

• merge (bool) – Merge RangesMatrix into aman.flags. Default is False.

• overwrite (bool) – Whether to write over any existing data in aman.flags[inv_var_flag_name] if merge is True. Default is False.

• inv_var_flag_name (str) – Field name used when merge is True. Default is inv_var_flag.

Returns:

mskptps – RangesMatrix of detectors with acceptable inverse variance. All ones if the detector should be cut.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_subscans(aman, flags=None, merge=True, include_turnarounds=False, overwrite=True)

Returns an axis manager with information about subscans. This includes direction and a ranges matrix (subscans samps) True inside each subscan.

Parameters:

• aman (AxisManager) – Input AxisManager.

• flags (FlagManager) – AxisManager or FlagManager containing (samps,) Ranges left_scan, right_scan, and turnarounds. Defaults to aman.flags.

• merge (bool) – Merge into aman as ‘subscan_info’

• include_turnarounds (bool) – Include turnarounds in the subscan ranges

• overwrite (bool) – If true, write over subscan_info.

Returns:

subscan_aman – AxisManager containing information about the subscans. “direction” is a (subscans,) array of strings ‘left’ or ‘right’ “subscan_flags” is a (subscans, samps) RangesMatrix; True inside the subscan.

Return type:

AxisManager

sotodlib.tod_ops.flags.get_subscan_signal(aman, arr, isub=None, trim=False)

Split an array into subscans.

Parameters:

• aman (AxisManager) – Input AxisManager.

• arr (Array) – Input array.

• isub (int) – Index of the desired subscan. May also be a list of indices. If None, all are used.

• trim (bool) – Do not include size-zero arrays from empty subscans in the output.

Returns:

out – If isub is a scalar, return an Array of arr cut on the samps axis to the given subscan. If isub is a list or None, return a list of such Arrays.

Return type:

list

sotodlib.tod_ops.flags.apply_rng(arr, rng)

Apply a Ranges object to an array. rng should be True on the samples you want to keep.

Parameters:

• arr (Array) – Array containing the signal. Should have one axis of len (samps).

• rng (Ranges) – Ranges object of len (samps) selecting the desired range.

sotodlib.tod_ops.flags.wrap_stats(aman, info_aman_name, info, info_names, merge=True)

Wrap multiple stats into a new aman, checking for subscan information. Stats can be (dets,) or (dets, subscans).

Parameters:

• aman (AxisManager) – Input AxisManager.

• info_aman_name (str) – Name for info_aman when wrapped into input.

• info (Array) – (stats, dets,) or (stats, dets, subscans) containing the information you want to wrap.

• info_names (list) – List of str names for each entry in the new aman.

• merge (bool) – If True merge info_aman into aman.

Returns:

info_aman – (dets,) or (dets, subscans) aman with a field for each item in info_names.

Return type:

AxisManager

sotodlib.tod_ops.flags.get_stats(aman, signal, stat_names, split_subscans=False, mask=None, name='stats', merge=False)

Calculate basic statistics on a TOD or power spectrum.

The statistics currently implemented are: - 'mean' - 'median' - 'ptp' (peak to peak) - 'std' (standard deviation) - 'var' (variance) - 'kurtosis' - 'skew'

Parameters:

• aman (AxisManager) – Input AxisManager.

• signal (Array) – Input signal. Statistics will be computed over axis 1.

• stat_names (list) – List of strings identifying which statistics to run.

• split_subscans (bool) – If True statistics will be computed on subscans. Assumes aman.subscan_info exists already.

• mask (Array) – Mask to apply before computation. 1d array for advanced indexing (keep True), or a slice object.

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

sotodlib.tod_ops.flags.get_focalplane_flags(aman, merge=True, overwrite=True, invalid_flags_name='fp_flags')

Generate flags for invalid detectors in the focal plane.

The tod.

mergebool

If true, merges the generated flag into aman.

overwritebool

If true, write over flag. If false, don’t.

invalid_flags_namestr

Name of flag to add to aman.flags if merge is True.

Returns:

msk_invalid_fp – RangesMatrix of invalid detectors in the focal plane.

Return type:

RangesMatrix

sotodlib.tod_ops.flags.noise_fit_flags(aman, low_wn, high_wn, high_fk)

Evaluate white noise and fknee cuts based on provided boundaries.

Parameters: aman : object

System Message: ERROR/3 (/home/docs/checkouts/readthedocs.org/user_builds/sotodlib/checkouts/stable/sotodlib/tod_ops/flags.py:docstring of sotodlib.tod_ops.flags.noise_fit_flags, line 5)

Unexpected indentation.

An object containing noise fit statistics and noise model coefficients.

System Message: WARNING/2 (/home/docs/checkouts/readthedocs.org/user_builds/sotodlib/checkouts/stable/sotodlib/tod_ops/flags.py:docstring of sotodlib.tod_ops.flags.noise_fit_flags, line 6)

Block quote ends without a blank line; unexpected unindent.

low_wnfloat or None

The lower boundary for white noise. If None, white noise flagging is skipped.

high_wnfloat or None

The upper boundary for white noise. If None, white noise flagging is skipped.

high_fkfloat or None

The upper boundary for fknee. If None, fknee flagging is skipped.

Returns: tuple or None

System Message: ERROR/3 (/home/docs/checkouts/readthedocs.org/user_builds/sotodlib/checkouts/stable/sotodlib/tod_ops/flags.py:docstring of sotodlib.tod_ops.flags.noise_fit_flags, line 15)

Unexpected indentation.

A tuple containing flags for valid white noise and fknee if both boundaries are provided. If only one boundary is provided, returns the corresponding flag. If no boundaries are provided, returns None.

sotodlib.tod_ops.flags.get_noisy_subscan_flags(aman, subscan_stats, nstd_lim=None, ptp_lim=None, kurt_lim=None, skew_lim=None, noisy_detector_lim=0.9, merge=False, overwrite=False, name='noisy_subscan')

Identify and flag bad subscans based on various statistical thresholds. aman : AxisManager

System Message: ERROR/3 (/home/docs/checkouts/readthedocs.org/user_builds/sotodlib/checkouts/stable/sotodlib/tod_ops/flags.py:docstring of sotodlib.tod_ops.flags.get_noisy_subscan_flags, line 3)

Unexpected indentation.

The tod.

System Message: WARNING/2 (/home/docs/checkouts/readthedocs.org/user_builds/sotodlib/checkouts/stable/sotodlib/tod_ops/flags.py:docstring of sotodlib.tod_ops.flags.get_noisy_subscan_flags, line 4)

Block quote ends without a blank line; unexpected unindent.

subscan_statsdict

Dictionary containing statistical metrics for subscans. Keys should include ‘std’, ‘ptp’, ‘kurtosis’, and ‘skew’.

nstd_limfloat [optional]

Threshold for standard deviation.

ptp_limfloat [optional]

Threshold for peak-to-peak values in pW.

kurt_limfloat [optional]

Threshold for kurtosis.

skew_limfloat [optional]

Threshold for skewness.

noisy_detector_limfloat [optional]

If > noisy_detector_lim fraction of samps are in flagged subscans then the detector is added to noisy_detector_flags.

mergebool

If true, merges the generated flag into aman.

overwritebool

If true, write over flag. If false, don’t.

namestr

Name of flag to add to aman.flags if merge is True.

Returns:

• noisy_subscan_flags (RangesMatrix) – RangesMatrix of bad subscans. (dets, samps), True for flagged (bad) samples.

• noisy_detector_flags (ndarray) – Array indicating detectors that are too noisy for more than noisy_detector_lim fraction of the obs duration. (dets,). We return ~noisy_detector_flags which is True for good detectors.

sotodlib.tod_ops.flags.expand_smurfgaps_flags(aman, buffer=200, name='smurfgaps', merge=True)

smurfgaps flags indicates the samples of each stream_id where the lost frames are filled in the bookbinding process. See sotodlib.io.bookbinder.bind.

This function expands smurfgaps flags of each stream_id to all detectors.

Parameters:

• aman (AxisManager) – Input AxisManager

• buffer (int) – Amount of buffer to apply on smurfgaps

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

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

Returns:

smurfgaps – smurfgaps flag with ‘dets’ and ‘samps’ axis

Return type:

RangesMatrix

sotodlib.tod_ops.flags.get_good_distribution_flags(aman, param_name='wn_signal', outlier_range=(0.5, 2.0), kurtosis_threshold=2.0, blame_max=False, blame_min=False)

Get detector cuts based on the detectors which fit within a gaussian distribution. Outlier and kurtosis thresholds set the cuts and param name determines which statistic’s distribution is evaluated to generate the cuts. This either cuts one detector at a time at the max or min of the distribution (blame max/min) or it determines if the distribution is screwed high or low and cuts max if high and min if low this is the default behavior.

Parameters:

• aman (AxisManager) – Input axis manager.

• param_name (str) – Field in aman where to evaluate distribution of.

• outlier_range (tuple) – Min and max of distribution allowed

• kurtosis_threshold (float) – Maximum allowed kurtosis of final distribution

• blame_max (bool) – If true, iteratively removes the maximum value from the distribution until the kurtosis threshold is met. blame_min cannot be True when blame_max is True. Default is False.

• blame_min (bool) – Same as blame_max but removes the minimum value.

Returns:

det_mask – Boolean mask for cuts True for keep and False for cut.

Return type:

list(bool)

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

class sotodlib.tod_ops.fft_ops.RFFTObj(n_det: int, n: int, a: ~numpy.ndarray[~typing.Any, ~numpy.dtype[~numpy.float32]], b: ~numpy.ndarray[~typing.Any, ~numpy.dtype[~numpy.complex64]], t_forward: ~typing.Callable = functools.partial(<function _t_null>, 'forward'), t_backward: ~typing.Callable = functools.partial(<function _t_null>, 'backward'))

Dataclass to store information needed for rfft.

n_det

Number of detectors this object was built for.

Type:

int

n

Number of samples this object was built for.

Type:

int

a

Buffer for the real part of the FFT.

Type:

numpy.ndarray[Any, numpy.dtype[numpy.float32]]

b

Buffer for the complex part of the FFT.

Type:

numpy.ndarray[Any, numpy.dtype[numpy.complex64]]

t_forward

Function for performing the forward FFT.

Type:

Callable

t_backward

Function for performing the backward FFT.

Type:

Callable

classmethod for_shape(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:

An instance of RFFTObj

Return type:

rfft_obj

classmethod for_tod(tod, direction='FFTW_FORWARD', **kwargs)

Wrapper around for_shape that pulls the shape from a 2d array

Parameters:

• tod – 2D array to build rfft object for.

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

• kwargs – additional arguments to pass to pyfftw.FFTW

Returns:

An instance of RFFTObj

Return type:

rfft_obj

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, merge_suffix=None, overwrite=True, subscan=False, full_output=False, label_axis='dets', **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.

• merge_suffix (str, optional) – Suffix to append to the Pxx field name in aman. Defaults to None (merged as Pxx).

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

• subscan (bool) – if True, compute psd on subscans.

• full_output – if True this also outputs nseg, the number of segments used for welch, for correcting bias of median white noise estimation by calc_wn.

• label_axis (str) – The name of LabelAxis in the input aman. Default is dets.

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

Returns:

array of frequencies corresponding to PSD calculated from welch. Pxx: array of PSD values. nseg: number of segments used for welch. this is returned if full_output is True.

Return type:

freqs

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

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) (nseg (Int or 1d Int) – Psd information to calculate white noise. Defaults to aman.pxx

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

• array) – number of segmnents used for welch. Defaults to aman.nseg. This is necessary for debiasing median white noise estimation. welch PSD with non-overlapping n segments follows chi square distribution with 2 * nseg degrees of freedom. The median of chi square distribution is biased from its average.

• (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

• (str) (method) – median or mean to compute white noise. Default is median.

Returns:

wn

Return type:

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

sotodlib.tod_ops.fft_ops.noise_ratio(aman, pxx, freqs, f_sig=(0.04, 0.14), f_wn=(0.6, 1.0), subscan=False)

Compute the ratio of the mean PSD in two frequency regions to evaluate the noise.

Parameters:

• aman (AxisManager) – Only used for matching the dets and subscans dimensions in the output aman.

• pxx (np.ndarray[float]) – Input PSD. Can be [dets, nufreq] or [dets, nufreq, subscans] (NOT just [nufreq]).

• freqs (np.ndarray[float]) – frequency information related to the psd.

• f_sig (tuple) – 2-tuple of frequencies giving the range of the numerator (“signal band”)

• f_wn (tuple) – 2-tuple of frequencies giving the range of the denominator (“white noise”)

• subscan (bool) – True if the PSD is split by subscans

Returns:

calc_aman – Axis manager with fields “rdets” giving the per-detector ratio and “rmean” giving the ratio for the mean (over detectors) PSD.

Return type:

AxisManager

sotodlib.tod_ops.fft_ops.noise_model(f, params, **fixed_param)

Noise model for power spectrum with white noise, and 1/f noise. If any fixed param is handed, that parameter is fixed in the fit. ‘alpha’ or ‘wn’ can be fixed. params = [wn, fknee, alpha]

sotodlib.tod_ops.fft_ops.get_psd_mask(aman, psd_mask=None, f=None, mask_hwpss=True, hwp_freq=None, max_hwpss_mode=10, hwpss_width=((-0.4, 0.6), (-0.2, 0.2)), mask_peak=False, peak_freq=None, peak_width=(-0.002, 0.002), merge=True, overwrite=True)

Function to get masks for hwpss or single peak in PSD.

Parameters:

• aman (AxisManager) – Axis manager which has ‘nusamps’ axis.

• psd_mask (numpy.ndarray or Ranges) – Existing psd_mask to be updated. If None, a new mask is created.

• f (nparray) – Frequency of PSD of signal. If None, aman.freqs are used.

• mask_hwpss (bool) – If True, hwpss are masked. Defaults to True.

• hwp_freq (float) – HWP frequency. If None, calculated based on aman.hwp_angle

• max_hwpss_mode (int) – Maximum hwpss mode to subtract.

• hwpss_width (array-like) – If given in float, nf-hwpss will be masked like (n * hwp_freq - width/2) < f < (n * hwp_freq + width/2). If given in array like [1.0, 0.4], 1f hwpss is masked like (hwp_freq - 0.5) < f < (hwp_freq + 0.5) and nf hwpss are masked like (n * hwp_freq - 0.2) < f < (n * hwp_freq + 0.2). If given in array like [[-0.4, 0.6], [-0.2, 0.3]], 1f hwpss is masked like (hwp_freq - 0.4) < f < (hwp_freq + 0.6) and nf are masked like (n * hwp_freq - 0.2) < f < (n * hwp_freq + 0.3).

• mask_peak (bool) – If True, single peak is masked.

• peak_freq (float) – Center frequency of the mask.

• peak_width (tuple) – Range to mask signal. The default masked range will be (peak_freq - 0.002) < f < (peak_freq + 0.002).

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

• mode (str) – if “replace”, existing PSD mask is replaced to new mask. If “add”, new mask range is added to the existing one.

• overwrite (bool) – if true will overwrite aman.psd_mask.

Returns:

psd_mask (nusamps)

Return type:

Ranges array. If merge == True, “psd_mask” is added to the aman.

sotodlib.tod_ops.fft_ops.get_binned_psd(aman, f=None, pxx=None, unbinned_mode=3, base=1.05, merge=False, overwrite=True)

Function that masks hwpss in PSD.

Parameters:

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

• f (nparray) – Frequency of PSD of signal.

• pxx (nparray) – PSD of signal.

• mask (bool) – if True calculate binned psd with mask.

• unbinned_mode (int) – First Fourier modes up to this number are left un-binned.

• base (float (> 1)) – Base of the logspace bins.

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

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

Returns:

f_binned, pxx_binned, bin_size

Return type:

binned frequency and PSD.

sotodlib.tod_ops.fft_ops.fit_noise_model(aman, signal=None, f=None, pxx=None, psdargs={}, fknee_est=1, wn_est=4e-05, alpha_est=3.4, merge_fit=False, lowf=None, f_max=100, merge_name='noise_fit_stats', merge_psd=True, mask=False, fixed_param=None, binning=False, unbinned_mode=3, base=1.05, freq_spacing=None, subscan=False, label_axis='dets', bounds=None, fit_method='likelihood', curve_fit_kwargs=None)

Fits noise model with white and 1/f noise to the PSD of binned signal. This uses a least square method since the distrubition of binned points is close to the standard distribution.

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

• fknee_est (float or nparray) – Initial estimation value of knee frequency. If it is given in array, initial values are applied for each detector.

• wn_est (float or nparray) – Initial estimation value of white noise. If it is given in array, initial values are applied for each detector.

• alpha_est (float or nparray) – Initial estimation value of alpha. If it is given in array, initial values are applied for each detector.

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

• lowf (float) – Minimum frequency to include in the fitting. Default is None which selects lowf as the second index of f.

• 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 merge_psd is True then adds fres and Pxx to the axis manager.

• mask (bool) – If True, (nusamps,) mask is taken from aman.psd_mask, or calculated on the fly. Can also be a 1d (nusamps,) bool array True for good samples to keep. Can also be a Ranges in which case it will be inverted before application.

• fixed_param (str) – This accepts ‘wn’ or ‘alpha’ or None. If ‘wn’ (‘alpha’) is given, white noise level (alpha) is fixed to the wn_est (alpha_est).

• binning (bool) – True to bin the psd before fitting. The binning is determined by the ‘unbinned_mode’ and ‘base’ params.

• unbinned_mode (int) – First Fourier modes up to this number are left un-binned.

• base (float (> 1)) – Base of the logspace bins.

• freq_spacing (float) – The approximate desired frequency spacing of the PSD. Passed to calc_psd.

• subscan (bool) – If True, fit noise on subscans.

• label_axis (str) – The name of LabelAxis in the input aman. Default is dets.

• bounds (list of tuple) – List length 2 with number of non-fixed params each list entry is a tuple first is lower bounds second is upper bounds to constrain fit.

• fit_method (str) – Selects optimizer. “likelihood” (default) uses the Nelder-Mead negative log-likelihood fit through scipy.optimize.minimize. “log_curve_fit” performs a log-space non-linear least square fit via scipy.optimize.curve_fit.

• curve_fit_kwargs (dict) – Extra options forwarded to curve_fit when fit_method is “log_curve_fit”. Ignored otherwise.

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

sotodlib.tod_ops.fft_ops.get_mask_for_hwpss(freq, hwp_freq, max_mode=10, width=((-0.4, 0.6), (-0.2, 0.2)))

Function that returns boolean array to mask hwpss in PSD.

Parameters:

• freq (nparray) – Frequency of PSD of signal.

• hwp_freq (float) – HWP frequency.

• max_mode (int) – Maximum hwpss mode to subtract.

• width (float/tuple/list/array) – If given in float, hwpss will be masked like (hwp_freq - width/2) < f < (hwp_freq + width/2). If given in array like [1.0, 0.4], 1f hwpss is masked like (hwp_freq - 0.5) < f < (hwp_freq + 0.5) and Nf hwpss are masked like (hwp_freq - 0.2) < f < (hwp_freq + 0.2). If given in array like [[-0.4, 0.6], [-0.2, 0.3]], 1f hwpss is masked like (hwp_freq - 0.4) < f < (hwp_freq + 0.6) and Nf are masked like (hwp_freq - 0.2) < f < (hwp_freq + 0.3). Usually 1f hwpss distrubites wider than other hwpss.

Returns:

mask – True in this array stands for the index of hwpss to mask.

Return type:

Boolean array to mask frequency and power of the given PSD.

sotodlib.tod_ops.fft_ops.get_mask_for_single_peak(f, peak_freq, peak_width=(-0.002, 0.002))

Function that returns boolean array to masks single peak (e.g. scan synchronous signal) in PSD.

Parameters:

• f (nparray) – Frequency of PSD of signal.

• peak_freq (float) – Center frequency of the mask.

• peak_width (tuple) – Range to mask signal. The default masked range will be (peak_freq - 0.002) < f < (peak_freq + 0.002).

Returns:

mask – True in this array stands for the index of the single peak to mask.

Return type:

Boolean array to mask the given PSD.

sotodlib.tod_ops.fft_ops.log_binning(psd, unbinned_mode=3, base=1.05, mask=None, return_bin_size=False, drop_nan=False)

Function to bin PSD or frequency. First several Fourier modes are left un-binned. Fourier modes higher than that are averaged into logspace bins.

Parameters:

• psd (numpy.ndarray) – PSD (or frequency) to be binned. Can be a 1D or 2D array.

• unbinned_mode (int, optional) – First Fourier modes up to this number are left un-binned. Defaults to 3.

• base (float, optional) – Base of the logspace bins. Must be greater than 1. Defaults to 1.05.

• mask (numpy.ndarray, optional) – Mask for psd. If all values in a bin are masked, the value becomes np.nan. Should be a 1D array.

• return_bin_size (bool, optional) – If True, the number of data points in the bins are returned. Defaults to False.

• drop_nan (bool, optional) – If True, drop the indices where psd is NaN. Defaults to False.

Returns:

• binned_psd (numpy.ndarray) – The binned PSD. If the input is 2D, the output will also be 2D with the same number of rows.

• bin_size (numpy.ndarray, optional) – The number of data points in each bin, only returned if return_bin_size is True.

sotodlib.tod_ops.fft_ops.build_filt_params_dict(filter_name, noise_fit=None, filter_params=None)

Build the filter parameter dictionary from a provided dictionary or from noise fit results.

Parameters:

• filter_name (str) – Name of the filter to build the parameter dict for.

• noise_fit (AxisManager) – AxisManager containing the result of the noise model fit sized nparams x ndets.

• filter_params (dict) – Filter parameters dictionary to complement parameters derived from the noise fit (or to be used if noise fit is None).

Returns:

filter_params – Returns a dictionary of the median values of the noise model fit parameters if noise_fit is not None, otherwise return the provided filter_params.

Return type:

dict

sotodlib.tod_ops.fft_ops.get_common_noise_params(aman, signal=None, merge=False, merge_name='common_mode', **fit_args)

Fits noise model with white and 1/f noise to the PSD of common noise. This can be applied to any 1D data with same size as aman.samps.

Parameters:

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

• signal (nparray) – The pre-computed common mode sized 1 x nsamps, or the signal sized ndets x nsamps. When signal is given in 2D shape, the common mode across detectors is calculated with median method.

• merge (bool) – If True, aman_common is added into axis manager with merge_name.

• merge_name (str) – aman_common is added with this name.

Returns:

noise_fit_stats – The fit parameters of given common mode.

Return type:

AxisManager