draco.analysis.delay

Delay space spectrum estimation and filtering.

Functions

`delay_power_spectrum_gibbs`(data, N, Ni, ...)	Estimate the delay power spectrum by Gibbs sampling.
`delay_spectrum_fft`(data, N[, window])	Estimate the delay transform from an input frequency spectrum by IFFT.
`delay_spectrum_gibbs_cross`(data, N, Ni, ...)	Estimate the delay power spectrum by Gibbs sampling.
`delay_spectrum_wiener_filter`(delay_PS, data, ...)	Estimate the delay spectrum from an input frequency spectrum by Wiener filtering.
`flatten_axes`(dset, axes_to_keep[, match_dset])	Move the specified axes of the dataset to the back, and flatten all others.
`fourier_matrix`(N[, fsel])	Generate a Fourier matrix to represent a real to complex FFT.
`fourier_matrix_c2c`(N[, fsel])	Generate a Fourier matrix to represent a complex to complex FFT.
`fourier_matrix_c2r`(N[, fsel])	Generate a Fourier matrix to represent a complex to real FFT.
`fourier_matrix_r2c`(N[, fsel])	Generate a Fourier matrix to represent a real to complex FFT.
`match_axes`(dset1, dset2)	Make sure that dset2 has the same set of axes as dset1.
`null_delay_filter`(freq, delay_cut, mask[, ...])	Take frequency data and null out any delays below some value.

Classes

`DelayCrossPowerSpectrumEstimator`()	A delay cross power spectrum estimator.
`DelayFilter`()	Remove delays less than a given threshold.
`DelayFilterBase`()	Remove delays less than a given threshold.
`DelayPowerSpectrumBase`()	Base class for delay power spectrum estimation (non-functional).
`DelayPowerSpectrumContainerMixin`()	Mixin for creating a delay power spectrum output container.
`DelayPowerSpectrumGeneralEstimator`()	Deprecated.
`DelayPowerSpectrumGibbs`()	Use a Gibbs sampler to estimate the delay power spectrum.
`DelayPowerSpectrumNRML`()	Use a NRML method to estimate the delay power spectrum.
`DelayPowerSpectrumStokesIEstimator`()	Deprecated.
`DelaySpectrumBase`()	Base class for delay spectrum estimation (non-functional).
`DelaySpectrumContainerMixin`()	Mixin for creating a delay transform output container.
`DelaySpectrumFFT`()	Class to measure the delay spectrum of a general container via ifft.
`DelaySpectrumToPowerSpectrum`()	Compute a delay power spectrum from a delay spectrum.
`DelaySpectrumWienerFilter`()	Class to measure delay spectrum of general container via Wiener filtering.
`DelaySpectrumWienerFilterIteratePS`()	Class to estimate the delay spectrum using Wiener filtering.
`DelayTransformBase`()	Base class for transforming from frequency to delay (non-functional).
`GeneralInputContainerMixin`()	Mixin for freq->delay transforms that collapse over several dataset axes.

class draco.analysis.delay.DelayCrossPowerSpectrumEstimator[source]

Bases: DelayPowerSpectrumGibbs, RandomTask

A delay cross power spectrum estimator.

This takes multiple compatible FreqContainer`s as inputs and will return a `DelayCrossSpectrum container with the full pair-wise cross power spectrum.

class draco.analysis.delay.DelayFilter[source]

Bases: SingleTask

Remove delays less than a given threshold.

This is performed by projecting the data onto the null space that is orthogonal to any mode at low delays.

Note that for this task to work best the zero entries in the weights dataset should factorize in frequency-time for each baseline. A mostly optimal masking can be generated using the draco.analysis.flagging.MaskFreq task.

delay_cut

Delay value to filter at in seconds.

Type:: float

za_cut

Sine of the maximum zenith angle included in baseline-dependent delay filtering. Default is 1 which corresponds to the horizon (ie: filters out all zenith angles). Setting to zero turns off baseline dependent cut.

Type:: float

extra_cut

Increase the delay threshold beyond the baseline dependent term.

Type:: float

weight_tol

Maximum weight kept in the masked data, as a fraction of the largest weight in the original dataset.

Type:: float

telescope_orientation

Determines if the baseline-dependent delay cut is based on the north-south component, the east-west component or the full baseline length. For cylindrical telescopes oriented in the NS direction (like CHIME) use ‘NS’. The default is ‘NS’.

Type:: one of (‘NS’, ‘EW’, ‘none’)

window

Apply the window function to the data when applying the filter.

Type:: bool

Notes

The delay cut applied is max(za_cut * baseline / c + extra_cut, delay_cut).

Initialize pipeline task.

May be overridden with no arguments. Will be called after any config.Property attributes are set and after ‘input’ and ‘requires’ keys are set up.

process(ss)[source]

Filter out delays from a SiderealStream or TimeStream.

Parameters:: ss (containers.SiderealStream) – Data to filter.
Returns:: ss_filt – Filtered dataset.
Return type:: containers.SiderealStream

setup(telescope)[source]

Set the telescope needed to obtain baselines.

Parameters:: telescope (TransitTelescope) – The telescope object to use

class draco.analysis.delay.DelayFilterBase[source]

Bases: SingleTask

Remove delays less than a given threshold.

This is performed by projecting the data onto the null space that is orthogonal to any mode at low delays.

Note that for this task to work best the zero entries in the weights dataset should factorize in frequency-time for each baseline. A mostly optimal masking can be generated using the draco.analysis.flagging.MaskFreq task.

delay_cut

Delay value to filter at in seconds.

Type:: float

window

Apply the window function to the data when applying the filter.

Type:: bool

axis

The main axis to iterate over. The delay cut can be varied for each element of this axis. If not set, a suitable default is picked for the container type.

Type:: str

dataset

Apply the delay filter to this dataset. If not set, a suitable default is picked for the container type.

Type:: str

Notes

The delay cut applied is max(za_cut * baseline / c + extra_cut, delay_cut).

Initialize pipeline task.

May be overridden with no arguments. Will be called after any config.Property attributes are set and after ‘input’ and ‘requires’ keys are set up.

process(ss: FreqContainerType) → FreqContainerType[source]

Filter out delays from a SiderealStream or TimeStream.

Parameters:: ss – Data to filter.
Returns:: Filtered dataset.
Return type:: ss_filt

setup(telescope: ProductManager | BeamTransfer | TransitTelescope)[source]

Set the telescope needed to obtain baselines.

Parameters:: telescope – The telescope object to use

class draco.analysis.delay.DelayPowerSpectrumBase[source]

Bases: DelayPowerSpectrumContainerMixin, DelayTransformBase

Base class for delay power spectrum estimation (non-functional).

class draco.analysis.delay.DelayPowerSpectrumContainerMixin[source]

Bases: GeneralInputContainerMixin

Mixin for creating a delay power spectrum output container.

nsamp

Number of samples to compute for each power spectrum. Default is 1.

Type:: int

save_samples

When using a sampling-based power spectrum estimator, save out every sample in the chain. Otherwise, only save the final power spectrum. Default is False.

Type:: bool

save_spectrum_mask

Save a mask which flags spectra which have significant error, as determined by the estimator. Default is False.

Type:: bool

class draco.analysis.delay.DelayPowerSpectrumGeneralEstimator[source]

Bases: DelayPowerSpectrumGibbs

Deprecated.

setup(requires=None)[source]: Raise a deprecation warnings.

class draco.analysis.delay.DelayPowerSpectrumGibbs[source]

Bases: DelayPowerSpectrumBase, RandomTask

Use a Gibbs sampler to estimate the delay power spectrum.

The spectrum returned is the median of the final half of the samples calulated.

initial_amplitude

The Gibbs sampler will be initialized with a flat power spectrum with this amplitude. Unused if maxpost=True (flat spectrum is a bad initial guess for the max-likelihood estimator). Default: 10.

Type:: float, optional

class draco.analysis.delay.DelayPowerSpectrumNRML[source]

Bases: DelayPowerSpectrumBase

Use a NRML method to estimate the delay power spectrum.

maxpost_tol

The convergence tolerance used by scipy.optimize.minimize in the maximum likelihood estimator.

Type:: float, optional

class draco.analysis.delay.DelayPowerSpectrumStokesIEstimator[source]

Bases: DelayPowerSpectrumGibbs

Deprecated.

setup(requires=None)[source]: Raise a deprecation warnings.

class draco.analysis.delay.DelaySpectrumBase[source]

Bases: DelaySpectrumContainerMixin, DelayTransformBase

Base class for delay spectrum estimation (non-functional).

class draco.analysis.delay.DelaySpectrumContainerMixin[source]

Bases: GeneralInputContainerMixin

Mixin for creating a delay transform output container.

save_spectrum_mask

Save a mask which flags spectra which have significant error, as determined by the estimator. Default is False.

Type:: bool

class draco.analysis.delay.DelaySpectrumFFT[source]

Bases: DelaySpectrumBase

Class to measure the delay spectrum of a general container via ifft.

class draco.analysis.delay.DelaySpectrumToPowerSpectrum[source]

Bases: SingleTask

Compute a delay power spectrum from a delay spectrum.

Initialize pipeline task.

May be overridden with no arguments. Will be called after any config.Property attributes are set and after ‘input’ and ‘requires’ keys are set up.

process(dspec: DelayTransform) → DelaySpectrum[source]

Get the delay power spectrum from a delay spectrum.

Parameters:: dspec – Delay spectrum container.
Returns:: Delay power spectrum container.
Return type:: pspec

class draco.analysis.delay.DelaySpectrumWienerFilter[source]

Bases: DelaySpectrumBase

Class to measure delay spectrum of general container via Wiener filtering.

The spectrum is calculated by applying a Wiener filter to the input frequency spectrum, assuming an input model for the delay power spectrum of the signal and that the noise power is described by the weights of the input container. See https://arxiv.org/abs/2202.01242, Eq. A6 for details.

setup(dps=None)[source]

Set the delay power spectrum to use as the signal covariance.

Parameters:: dps (containers.DelaySpectrum) – Delay power spectrum for signal part of Wiener filter.

class draco.analysis.delay.DelaySpectrumWienerFilterIteratePS[source]

Bases: DelaySpectrumWienerFilter

Class to estimate the delay spectrum using Wiener filtering.

This class extends DelaySpectrumWienerFilter by allowing the delay power spectrum (dps) to be updated with each call to process instead of being fixed at setup. The updated dps is used to apply the Wiener filter to the input frequency spectrum.

process(ss, dps)[source]

Estimate the delay spectrum.

Parameters:

ss (containers.FreqContainer) – Data to transform. Must have a frequency axis.
dps (containers.DelaySpectrum) – Delay power spectrum for signal part of Wiener filter.

Returns:

out_cont – Output delay spectrum or delay power spectrum.

Return type:

containers.DelayTransform or containers.DelaySpectrum

class draco.analysis.delay.DelayTransformBase[source]

Bases: SingleTask

Base class for transforming from frequency to delay (non-functional).

freq_zero

The physical frequency (in MHz) of the zero channel. That is the DC channel coming out of the F-engine. If not specified, use the first frequency channel of the stream.

Type:: float, optional

freq_spacing

The spacing between the underlying channels (in MHz). This is conjugate to the length of a frame of time samples that is transformed. If not set, then use the smallest gap found between channels in the dataset.

Type:: float, optional

nfreq

The number of frequency channels in the full set produced by the F-engine. If not set, assume the last included frequency is the last of the full set (or is the penultimate if skip_nyquist is set).

Type:: int, optional

skip_nyquist

Whether the Nyquist frequency is included in the data. This is True by default to align with the output of CASPER PFBs.

Type:: bool, optional

apply_window

Whether to apply apodisation to frequency axis. Default: True.

Type:: bool, optional

window

Apodisation to perform on frequency axis. Default: ‘nuttall’.

Type:: window available in draco.util.tools.window_generalised(), optional

complex_timedomain

Whether to assume the original time samples that were channelized into a frequency spectrum were purely real (False) or complex (True). If True, freq_zero, nfreq, and skip_nyquist are ignored. Default: False.

Type:: bool, optional

use_average_weights

Use noise weights averaged over time samples. This means that only a single covariance matrix needs to be created for unique power spectrum (i.e., for each baseline), but it is only valid if the frequency masking is constant in time. Default: True.

Type:: bool, optional

weight_boost

Multiply weights in the input container by this factor. This causes the task to assume the noise power in the data is weight_boost times lower, which is useful if you want the “true” noise to not be downweighted by the Wiener filter, or have it included in the Gibbs sampler. Default: 1.0.

Type:: float, optional

freq_frac: The threshold for the fraction of time samples present in a frequency for it to be retained. Must be strictly greater than this value, so the default value 0, retains any channel with at least one sample. A value of 0.01 would retain any frequency that has > 1% of time samples unmasked.

time_frac: The threshold for the fraction of frequency samples required to retain a time sample. Must be strictly greater than this value. The default value (-1) means that all time samples are kept. A value of 0.01 would keep any time sample with >1% of frequencies unmasked.

remove_mean: Subtract the mean in time of each frequency channel. This is done after time samples are pruned by the time_frac threshold.

scale_freq: Scale each frequency by its standard deviation to flatten the fluctuations across the band. Applied before any apodisation is done.

Initialize pipeline task.

May be overridden with no arguments. Will be called after any config.Property attributes are set and after ‘input’ and ‘requires’ keys are set up.

process(ss)[source]

Estimate the delay spectrum or power spectrum.

Parameters:: ss (containers.FreqContainer) – Data to transform. Must have a frequency axis.
Returns:: out_cont – Output delay spectrum or delay power spectrum.
Return type:: containers.DelayTransform or containers.DelaySpectrum

class draco.analysis.delay.GeneralInputContainerMixin[source]

Bases: object

Mixin for freq->delay transforms that collapse over several dataset axes.

The delay spectrum or power spectrum output is indexed by a baseline axis. This axis is the composite axis of all the axes in the container except the frequency axis or the sample_axis. These constituent axes are included in the index map, and their order is given by the baseline_axes attribute.

dataset

Calculate the delay spectrum of this dataset (e.g., “vis”, “map”, “beam”). If not set, assume the input is a DataWeightContainer and use the main data dataset.

Type:: str, optional

sample_axis

Assume that every sample over this axis is drawn from the same power spectrum.

Type:: str

draco.analysis.delay.delay_power_spectrum_gibbs(data, N, Ni, initial_S, window='nuttall', fsel=None, niter=20, rng=None, complex_timedomain=False)[source]

Estimate the delay power spectrum by Gibbs sampling.

This routine estimates the spectrum at the N delay samples conjugate to an input frequency spectrum with N/2 + 1 channels (if the delay spectrum is assumed real) or N channels (if the delay spectrum is assumed complex). A subset of these channels can be specified using the fsel argument.

Parameters:

data (np.ndarray[:, freq]) – Data to estimate the delay spectrum of.
N (int) – The length of the output delay spectrum. There are assumed to be N/2 + 1 total frequency channels if assuming a real delay spectrum, or N channels for a complex delay spectrum.
Ni (np.ndarray[freq]) – Inverse noise variance.
initial_S (np.ndarray[delay]) – The initial delay power spectrum guess.
window (one of {'nuttall', 'blackman_nuttall', 'blackman_harris', None}, optional) – Apply an apodisation function. Default: ‘nuttall’.
fsel (np.ndarray[freq], optional) – Indices of channels that we have data at. By default assume all channels.
niter (int, optional) – Number of Gibbs samples to generate.
rng (np.random.Generator, optional) – A generator to use to produce the random samples.
complex_timedomain (bool, optional) – If True, assume input data arose from a complex timestream. If False, assume input data arose from a real timestream, such that the first and last frequency channels have purely real values. Default: False.

Returns:

spec – List of spectrum samples.

Return type:

list

draco.analysis.delay.delay_spectrum_fft(data, N, window='nuttall')[source]

Estimate the delay transform from an input frequency spectrum by IFFT.

This routine makes no attempt to account for data masking, and only supports complex to complex fft.

Parameters:

data (np.ndarray[nsample, freq]) – Data to estimate the delay spectrum of.
N (int) – The length of the output delay spectrum. There are assumed to be N/2 + 1 total frequency channels if assuming a real delay spectrum, or N channels for a complex delay spectrum.
window (one of {'nuttall', 'blackman_nuttall', 'blackman_harris', None}, optional) – Apply an apodisation function. Default: ‘nuttall’.

Returns:

y_spec – Delay spectrum for each element of the sample axis.

Return type:

np.ndarray[nsample, ndelay]

draco.analysis.delay.delay_spectrum_gibbs_cross(data: ndarray, N: int, Ni: ndarray, initial_S: ndarray, window: str = 'nuttall', fsel: ndarray | None = None, niter: int = 20, rng: Generator | None = None) → list[ndarray][source]

Estimate the delay power spectrum by Gibbs sampling.

This routine estimates the spectrum at the N delay samples conjugate to an input frequency spectrum with N/2 + 1 channels (if the delay spectrum is assumed real) or N channels (if the delay spectrum is assumed complex). A subset of these channels can be specified using the fsel argument.

Parameters:

data – A 3D array of [dataset, sample, freq]. The delay cross-power spectrum of these will be calculated.
N – The length of the output delay spectrum. There are assumed to be N/2 + 1 total frequency channels if assuming a real delay spectrum, or N channels for a complex delay spectrum.
Ni – Inverse noise variance as a 3D [dataset, sample, freq] array.
initial_S – The initial delay cross-power spectrum guess. A 3D array of [data1, data2, delay].
window (one of {'nuttall', 'blackman_nuttall', 'blackman_harris', None}, optional) – Apply an apodisation function. Default: ‘nuttall’.
fsel – Indices of channels that we have data at. By default assume all channels.
niter – Number of Gibbs samples to generate.
rng – A generator to use to produce the random samples.

Returns:

spec – List of cross-power spectrum samples.

Return type:

list

draco.analysis.delay.delay_spectrum_wiener_filter(delay_PS, data, N, Ni, window='nuttall', fsel=None, complex_timedomain=False)[source]

Estimate the delay spectrum from an input frequency spectrum by Wiener filtering.

This routine estimates the spectrum at the N delay samples conjugate to an input frequency spectrum with N/2 + 1 channels (if the delay spectrum is assumed real) or N channels (if the delay spectrum is assumed complex). A subset of these channels can be specified using the fsel argument.

Parameters:

delay_PS (np.ndarray[ndelay]) – Delay power spectrum to use for the signal covariance in the Wiener filter.
data (np.ndarray[nsample, freq]) – Data to estimate the delay spectrum of.
N (int) – The length of the output delay spectrum. There are assumed to be N/2 + 1 total frequency channels if assuming a real delay spectrum, or N channels for a complex delay spectrum.
Ni (np.ndarray[freq]) – Inverse noise variance.
fsel (np.ndarray[freq], optional) – Indices of channels that we have data at. By default assume all channels.
window (one of {'nuttall', 'blackman_nuttall', 'blackman_harris', None}, optional) – Apply an apodisation function. Default: ‘nuttall’.
complex_timedomain (bool, optional) – If True, assume input data arose from a complex timestream. If False, assume input data arose from a real timestream, such that the first and last frequency channels have purely real values. Default: False.

Returns:

y_spec – Delay spectrum for each element of the sample axis.

Return type:

np.ndarray[nsample, ndelay]

draco.analysis.delay.flatten_axes(dset: MemDatasetDistributed, axes_to_keep: list[str], match_dset: MemDatasetDistributed | None = None) → tuple[MPIArray, list[str]][source]

Move the specified axes of the dataset to the back, and flatten all others.

Optionally this will add length-1 axes to match the axes of another dataset.

Parameters:

dset – The dataset to reshape.
axes_to_keep – The names of the axes to keep.
match_dset – An optional dataset to match the shape of.

Returns:

flat_array – The MPIArray representing the re-arranged dataset. Distributed along the flattened axis.
flat_axes – The names of the flattened axes from slowest to fastest varying.

draco.analysis.delay.fourier_matrix(N: int, fsel: ndarray | None = None) → ndarray[source]

Generate a Fourier matrix to represent a real to complex FFT.

Parameters:

N (integer) – Length of timestream that we are transforming to. Must be even.
fsel (array_like, optional) – Indexes of the frequency channels to include in the transformation matrix. By default, assume all channels.

Returns:

Fr – An array performing the Fourier transform from a real time series to frequencies packed as alternating real and imaginary elements,

Return type:

np.ndarray

draco.analysis.delay.fourier_matrix_c2c(N, fsel=None)[source]

Generate a Fourier matrix to represent a complex to complex FFT.

These Fourier conventions match numpy.fft.fft().

Parameters:

N (integer) – Length of timestream that we are transforming to.
fsel (array_like, optional) – Indices of the frequency channels to include in the transformation matrix. By default, assume all channels.

Returns:

F – An array performing the Fourier transform from a complex time series to frequencies, with both input and output packed as alternating real and imaginary elements.

Return type:

np.ndarray

draco.analysis.delay.fourier_matrix_c2r(N, fsel=None)[source]

Generate a Fourier matrix to represent a complex to real FFT.

Parameters:

N (integer) – Length of timestream that we are transforming to. Must be even.
fsel (array_like, optional) – Indexes of the frequency channels to include in the transformation matrix. By default, assume all channels.

Returns:

Fr – An array performing the Fourier transform from frequencies packed as alternating real and imaginary elements, to the real time series.

Return type:

np.ndarray

draco.analysis.delay.fourier_matrix_r2c(N, fsel=None)[source]

Generate a Fourier matrix to represent a real to complex FFT.

Parameters:

N (integer) – Length of timestream that we are transforming to. Must be even.
fsel (array_like, optional) – Indexes of the frequency channels to include in the transformation matrix. By default, assume all channels.

Returns:

Fr – An array performing the Fourier transform from a real time series to frequencies packed as alternating real and imaginary elements,

Return type:

np.ndarray

draco.analysis.delay.match_axes(dset1, dset2)[source]

Make sure that dset2 has the same set of axes as dset1.

Sometimes the weights are missing axes (usually where the entries would all be the same), we need to map these into one another and expand the weights to the same size as the visibilities. This assumes that the vis/weight axes are in the same order when present

Parameters:

dset1 – The dataset with more axes.
dset2 – The dataset with a subset of axes. For the moment these are assumed to be in the same order.

Returns:

A view of dset2 with length-1 axes inserted to match the axes missing from dset1.

Return type:

dset2_view

draco.analysis.delay.null_delay_filter(freq, delay_cut, mask, num_delay=200, tol=1e-08, window=True, type_='high', lapack_driver='gesvd')[source]

Take frequency data and null out any delays below some value.

Parameters:

freq (np.ndarray[freq]) – Frequencies we have data at.
delay_cut (float) – Delay cut to apply.
mask (np.ndarray[freq]) – Frequencies to mask out.
num_delay (int, optional) – Number of delay values to use.
tol (float, optional) – Cut off value for singular values.
window (bool, optional) – Apply a window function to the data while filtering.
type (str, optional) – Whether to apply a high-pass or low-pass filter. Options are high or low. Default is high.
lapack_driver (str, optional) – Which lapack driver to use in the SVD. Options are ‘gesvd’ or ‘gesdd’. ‘gesdd’ is generally faster, but seems to experience convergence issues. Default is ‘gesvd’.

Returns:

filter – The filter as a 2D matrix.

Return type:

np.ndarray[freq, freq]