LegendDSP

Digital signal processing for LEGEND waveforms.

The package provides

• DSP chains that reduce tables of waveforms to tables of waveform parameters: dsp_icpc for HPGe detectors, dsp_puls for pulser channels, dsp_aux for auxiliary channels, dsp_sipm for SiPMs and dsp_pmts for PMTs, each with a _compressed variant where compressed data exists;
• grid scans over filter parameters (dsp_*_optimization) that provide the input for choosing optimal filter settings;
• functions for single waveform parameters (get_t0, get_threshold, get_a, get_qdrift, get_lq, get_triggers) and waveform statistics (tailstats, extremestats, thresholdstats, saturation);
• filters for RadiationDetectorDSP.jl waveforms, such as MovingWindowMultiFilter, WeightedSavitzkyGolayFilter and HaarAveragingFilter;
• support-vector-machine quality cuts (get_qc_ml_func, get_qc_classifier);
• the configuration type DSPConfig.

DSPConfig{T <: Real}

Settings of the HPGe DSP chains and filter-parameter scans.

All time-valued fields are Unitful quantities. Create an instance from DSP metadata with DSPConfig(pd::PropDict).

Fields

• enc_pickoff_trap, enc_pickoff_zac, enc_pickoff_cusp: time at which the trapezoidal, ZAC and CUSP filtered waveforms are evaluated in the rise-time scans dsp_trap_rt_optimization, dsp_zac_rt_optimization and dsp_cusp_rt_optimization
• flt_length_cusp, flt_length_zac: total length of the CUSP and ZAC filters
• t0_threshold: threshold in units of the filtered signal for the t0 determination of get_t0
• inTraceCut_std_threshold: threshold of the in-trace pile-up search in standard deviations of the baseline
• sg_flt_degree: polynomial degree of the Savitzky-Golay filters that determine the current
• bl_window: time window for the baseline fit
• tail_window: time window for the decay-time fit and the tail statistics after pole-zero correction
• current_window: time window in which the maximum of the current is searched
• qdrift_int_length: StepRangeLen whose first and last elements are the two integration times of get_qdrift
• lq_windows: vector of [early, gap, late] window lengths of the LQ parameters, see get_lq
• e_grid_rt_trap, e_grid_ft_trap: rise-time and flat-top-time scan ranges of the trapezoidal filter
• e_grid_rt_zac, e_grid_ft_zac: rise-time and flat-top-time scan ranges of the ZAC filter
• e_grid_rt_cusp, e_grid_ft_cusp: rise-time and flat-top-time scan ranges of the CUSP filter
• a_grid_wl_sg: window-length scan range of the Savitzky-Golay current filter
• a_grid_wl_mwa: window-length scan range of the moving-window-average current filter
• a_mwa_windows: fixed window lengths of the moving-window-average current filter; each gives a column a_mwa_<ns>
• a_mwa_passes: number of alternating passes of the moving-window-average current filter, see MovingWindowMultiFilter
• a_sg_windows: fixed window lengths of the Savitzky-Golay current filter; each gives a column a_sg_<ns>
• default_flt_param: PropDict with the default filter parameters: rt and ft for trap, cusp and zac, and the window length for sg and mwa
• kwargs_pars: PropDict with further parameters: fc_bit_depth, t0_flt_pars, t0_mintot, tx_mintot, intrace_mintot, int_interpolation_order, int_interpolation_length, sig_interpolation_order and sig_interpolation_length
• auxpz1_window, auxpz2_window: time windows of the auxiliary pole-zero statistics auxpz1_* and auxpz2_* of dsp_icpc_compressed

Examples

LegendDataManagement is not a dependency of LegendDSP and is only used here to read the metadata.

using LegendDSP
using LegendDataManagement

l200 = LegendData(:l200)
filekey = start_filekey(l200, (:p03, :r000, :cal))
dsp_config = DSPConfig(dataprod_config(l200).dsp(filekey).default)

DSPConfig(pd::PropDict)

Create a DSPConfig{Float64} from DSP metadata.

pd needs one entry per field of DSPConfig, with these differences:

• bl_window, tail_window, current_window, auxpz1_window and auxpz2_window are given as (min, max) entries;
• e_grid_trap, e_grid_zac and e_grid_cusp each contain rt and ft entries with start, step and stop; a_grid_wl_sg and a_grid_wl_mwa have the same three entries. a_grid_wl_mwa is optional and defaults to a_grid_wl_sg;
• qdrift_int_length holds a start and a stop time; the resulting range has a step of 0.1 in the unit of the start time;
• default_flt_param is read from the entry flt_defaults;
• every element of lq_windows must have exactly three entries, otherwise an ArgumentError is thrown.

struct DerivativeFilter{T<:RealQuantity} <: AbstractRadIIRFilter

Backward difference of a signal: y[i] = gain * (x[i] - x[i-1]), and y[1] = gain * (x[2] - x[1]) for the first sample.

The difference is not divided by the sampling step, so the result is the change per sample times gain. The time axis is not changed.

Working example:

using LegendDSP
using RadiationDetectorSignals
using Unitful

signal = rand(100)
t = range(0u"ms", 20u"ms", 100)

wf = RDWaveform(t, signal)
flt = DerivativeFilter()
wf_new = flt(wf)

Constructors:

• DerivativeFilter(; fields...)

Fields:

• gain::Union{Real, Unitful.AbstractQuantity{<:Real}}: Filter gain

FltOptPrep{W, T, B, E, P}(wvfs_pz, t50, bl_stats, signal_estimator, sampling_period)
FltOptPrep(wvfs_pz, t50, bl_stats, signal_estimator, sampling_period)

Baseline-subtracted, pole-zero-corrected waveforms with the quantities derived from them, prepared once for repeated filter evaluations. Create it with dsp_flt_optimization_prepare; it is the input of dsp_energy_rtft and of the FltOptPrep methods of dsp_trap_rt_optimization, dsp_cusp_rt_optimization and dsp_zac_rt_optimization.

Fields

• wvfs_pz: decoded waveforms with the baseline mean in config.bl_window subtracted and pole-zero-corrected
• t50: time at which each pole-zero-corrected waveform reaches 50% of its maximum
• bl_stats: baseline statistics (mean, sigma, slope, offset) of each waveform in config.bl_window, from signalstats
• signal_estimator: SignalEstimator that reads a filtered waveform at a given time
• sampling_period: step of the time axis of the first waveform, assumed to be the same for all waveforms

struct HaarAveragingFilter <: AbstractRadSigFilter{LinearFiltering}

Haar averaging filter: sums each pair of adjacent samples, scales the sum by 1/√2 and keeps every down_sampling_rate-th result.

With down_sampling_rate = 2 the pairs do not overlap, which corresponds to the approximation step of a discrete Haar wavelet transform. The output has ceil(n / down_sampling_rate) samples for n input samples and the output time axis consists of every down_sampling_rate-th element of the input time axis. Sample indices beyond the signal are clamped to the first and last sample.

Constructor:

• HaarAveragingFilter(down_sampling_rate::Int)

Fields:

• down_sampling_rate: distance in samples between the starts of successive pairs
• length: number of averaged samples, always 2

struct IntersectMaximum{T<:RealQuantity} <: Function

Find the up-crossings of a signal through a threshold and the maximum following each of them.

An up-crossing is a position where the signal rises to the threshold and stays at or above it for at least mintot. A crossing at the first sample is not counted. The crossing time is interpolated linearly.

Call an instance as

(f::IntersectMaximum)(signal::SamplesOrWaveform, threshold::RealQuantity)

It returns a named tuple. All entries except multiplicity are vectors with one element per up-crossing:

• x: crossing time
• x_high: time at which the signal falls below the threshold again (interpolated linearly), or the last time of the signal if it does not
• x_tot: x_high - x
• max: maximum of the signal within maxtot after the crossing, interpolated by a parabola if it is an interior sample
• multiplicity: number of up-crossings

mintot and maxtot are converted to numbers of samples using the step of the time axis. For a time axis with a unit they need to be quantities of the same dimension; the defaults of 4 and 100 are plain numbers.

Constructors:

• IntersectMaximum(; fields...)
• IntersectMaximum(mintot, maxtot)

Fields:

• mintot::Union{Real, Unitful.AbstractQuantity{<:Real}}: minimum time-over-threshold

• maxtot::Union{Real, Unitful.AbstractQuantity{<:Real}}: maximum time-over-threshold for max to appear

struct ModifiedSincFilter{T<:RealQuantity} <: AbstractRadFIRFilter

A Modified-Sinc filter.

The kernel is a windowed sinc function with a half-width of m; the degree d determines the number of extrema in the kernel. At the signal boundaries the signal is extended by weighted linear extrapolation. d has to be even and in 2:10, and m has to be at least d/2 + 2 samples; otherwise an AssertionError is thrown.

Working example:

using RadiationDetectorSignals
using Unitful

n = 600
noise = 1.
t = range(0u"μs", 20u"μs", 2*n)
signal = vcat(zeros(n), 10*ones(n)) + (noise*rand(2*n) .- noise/2)
wf = RDWaveform(t, signal)

# define filter parameters and filter
flt = ModifiedSincFilter(d=4, m=1u"μs")

# apply filter to signal
wf_new = flt(wf)

Constructors:

• ModifiedSincFilter(; fields...)
• ModifiedSincFilter(d, m)

Fields:

• d::Int64: degree of the filter determining the number of extrema in the kernel

• m::Union{Real, Unitful.AbstractQuantity{<:Real}}: half-width of the kernel

struct MovingWindowFilter{T<:RealQuantity} <: AbstractRadFIRFilter

Causal moving average of a signal over a window of length.

With the window length l in samples (length divided by the step of the time axis, rounded) and L samples, the output is

y₁ = x₁
yₙ = yₙ₋₁ + (xₙ - x₁)/l,       for n ∈ {2, …, l}
yₙ = yₙ₋₁ + (xₙ - xₙ₋ₗ)/l,     for n ∈ {l+1, …, L}

which equals a moving average over the l latest samples with the signal extended to the left by its first sample. The time axis is not changed.

Working example:

using RadiationDetectorSignals
using Unitful

signal = rand(100)
t = range(0u"ms", 20u"ms", 100)

wf = RDWaveform(t, signal)
flt = MovingWindowFilter(1u"ms")
wf_new = flt(wf)

Constructors:

• MovingWindowFilter(; fields...)

Fields:

• length::Union{Real, Unitful.AbstractQuantity{<:Real}}: length of the moving average window

struct MovingWindowMultiFilter{T<:RealQuantity} <: AbstractRadFIRFilter

Apply passes moving average windows of size length to the signal, alternating forward and backward in time.

With the default three passes the directions are forward, backward, forward. An even number of passes is symmetric in time. Each pass is a MovingWindowFilter. An ArgumentError is thrown for passes < 1.

Working example:

using RadiationDetectorSignals
using Unitful

signal = rand(100)
t = range(0u"ms", 20u"ms", 100)

wf = RDWaveform(t, signal)
flt = MovingWindowMultiFilter(1u"ms")              # three passes: forward, backward, forward
flt = MovingWindowMultiFilter(1u"ms"; passes = 2)  # forward, backward
wf_new = flt(wf)

Constructors:

• MovingWindowMultiFilter(; fields...)
• MovingWindowMultiFilter(length; passes = 3)

Fields:

• length::Union{Real, Unitful.AbstractQuantity{<:Real}}: size of the moving average window

• passes::Int64: number of moving average passes, alternating forward and backward in time

struct MultiIntersect{T<:RealQuantity} <: Function

Find the first times at which a signal reaches given fractions of its maximum.

For each ratio in threshold_ratios, the threshold is that ratio times the maximum of the signal. The first position where the signal stays at or above the threshold for at least mintot is refined: a polynomial of degree d is fitted to the 2n samples around it, evaluated on a grid with sampling_rate times as many points, and searched for the threshold crossing. The thresholds are searched in the order given, each starting at the position found for the previous one, so the ratios should be ascending. An AssertionError is thrown if the fit window would reach beyond the signal.

Call an instance as

(f::MultiIntersect)(signal::SamplesOrWaveform)

It returns a vector with one crossing time per ratio, and an array of zeros for an empty signal.

mintot is converted to a number of samples using the step of the time axis. For a time axis with a unit it needs to be a quantity of the same dimension; the default of 4 is a plain number.

Constructors:

• MultiIntersect(; fields...)

Fields:

• threshold_ratios::Vector{Float64}: ratios which determine the thresholds

• mintot::Union{Real, Unitful.AbstractQuantity{<:Real}}: minimum time-over-threshold

• n::Int64: half window length of polynomial fit

• d::Int64: degree of polynomial

• sampling_rate::Int64: upsampling rate

struct TimeAxisFilter{T<:RealQuantity} <: AbstractRadIIRFilter

Replace the time axis of a waveform by a range with step period that starts at first(time) + offset. The signal is not changed.

The time axis of the waveform needs to be a range, otherwise an ArgumentError is thrown.

Working example:

using LegendDSP
using RadiationDetectorSignals
using Unitful

signal = rand(100)
t = range(0u"ms", 20u"ms", 100)

wf = RDWaveform(t, signal)
flt = TimeAxisFilter(4u"ns")
wf_new = flt(wf)

Constructors:

• TimeAxisFilter(; fields...)
• TimeAxisFilter(period, offset=false*period)

Fields:

• period::Union{Real, Unitful.AbstractQuantity{<:Real}}: Filter period

• offset::Union{Real, Unitful.AbstractQuantity{<:Real}}: Filter offset

struct WeightedSavitzkyGolayFilter{T<:RealQuantity} <: AbstractRadFIRFilter

A Weighted-Savitzky-Golay filter.

A polynomial of degree degree is fitted by weighted least squares in a window of total length length around each sample, and the fitted value is the output. The weights fall off with the distance from the window center according to weightType:

• 0: constant weights, which gives the ordinary Savitzky-Golay filter
• 1: Gaussian-like weights
• 2: Hann weights, cos²
• 3: squared Hann weights, cos⁴
• 4: cubed Hann weights, cos⁶

The window is adapted near the boundaries of the signal. Other values of weightType raise an AssertionError, as do a signal shorter than the window and a degree larger than twice the half-window in samples.

Working example:

using RadiationDetectorSignals
using Unitful

n = 600
noise = 1.
t = range(0u"μs", 20u"μs", 2*n)
signal = vcat(zeros(n), 10*ones(n)) + (noise*rand(2*n) .- noise/2)
wf = RDWaveform(t, signal)

# define filter parameters and filter
flt = WeightedSavitzkyGolayFilter(length=1u"μs", degree=3, weightType=2)

# apply filter to signal
wf_new = flt(wf)

Constructors:

• WeightedSavitzkyGolayFilter(; fields...)
• WeightedSavitzkyGolayFilter(length, degree, weightType)

Fields:

• length::Union{Real, Unitful.AbstractQuantity{<:Real}}: total filter length

• degree::Int64: Polynomial degree

• weightType::Int64: weight function to use

struct WhittakerHendersonFilter <: AbstractRadSigFilter{LinearFiltering}

A Whittaker-Henderson filter.

The output y for an input x minimizes sum((x - y).^2) + λ * sum(Δᵖy.^2), where Δᵖy are the finite differences of order p of y. It solves (I + λ DᵀD) y = x with the difference matrix D by a banded Cholesky decomposition. A larger λ gives a smoother output.

Working example:

using RadiationDetectorSignals
using Unitful

n = 600
noise = 1.
t = range(0u"μs", 20u"μs", 2*n)
signal = vcat(zeros(n), 10*ones(n)) + (noise*rand(2*n) .- noise/2)
wf = RDWaveform(t, signal)

# define filter parameters and filter
flt = WhittakerHendersonFilter(p=3, λ=4)

# apply filter to signal
wf_new = flt(wf)

Constructors:

• WhittakerHendersonFilter(; fields...)

Fields:

• p::Int64: highest differentiation degree

• λ::Real: penalty strengh imposed on non-smooth functions

align_waveform(wvf::RDWaveform, t_align, t_target)

Shift wvf in time so that the feature at t_align lands at t_target. The output has the time axis and length of the input; its value at time t is the input at t + (t_align - t_target), interpolated with a cubic Catmull-Rom spline through the four nearest samples. Positions outside the record repeat the first or last sample.

Arguments

• wvf: waveform to shift
• t_align: time of the feature in wvf, for example t50
• t_target: time at which the feature lies in the returned waveform

Returns

• RDWaveform with the time axis of wvf and floating-point samples

align_waveforms(wvfs::ArrayOfRDWaveforms, t_align::AbstractVector, t_target)

align_waveform for every waveform of wvfs with its own alignment time from t_align, one per waveform. Throws a DimensionMismatch if the axes of wvfs and t_align differ.

aligned_windows(wvfs::ArrayOfRDWaveforms, t_target, n_pre::Integer, n_post::Integer)

Samples from n_pre before to n_post after the sample nearest to t_target, one column per waveform. Throws if the window leaves the record, or if the waveforms do not share one time axis.

Returns

• matrix of size (n_pre + n_post + 1, length(wvfs))

dplms_binned_reference(windows::AbstractMatrix, bins::AbstractVector{<:Integer}, plateau::Real)

Reference constraints per drift-time bin: the superpulse r_b of every bin b in bins, each divided by the common plateau, combined into mat = mean_b r_b r_bᵀ and vec = mean_b r_b. Bins are numbered 1:maximum(bins) and every bin must hold an event. Throws if a bin holds no event, if a bin index is not positive, or if windows and bins differ in the number of events.

Returns

• named tuple (mat, vec, refs) with refs the per-bin references stacked into one column per bin

dplms_default_grid(n_rise::Integer; nm = 1.0, za = 1e3, ft = 1.0, ct = 0.0, zo = 0.0, ff_factors = 1:4, ff = nothing)

Grid of DPLMS weights and flat-top widenings for dsp_dplms_optimization.

The widenings are ff = k·n_rise samples for k in ff_factors, or the sample counts ff if given. The grid holds every combination of zo (a number or a collection) and widening, all at the weights nm, za, ft, ct, ordered by zo, then by ff. Only weight ratios matter, and in units of the reference amplitude any za above about 1 already enforces zero area; za = 0 switches the zero-area condition off. zo = 0 leaves the area of the filter output free, zo = 1e3 makes it zero. See dplms_synthesis for the meaning of the weights.

Arguments

• n_rise: number of samples of the signal rise, used to scale the widenings

Keywords

• nm, za, ft, ct, zo: weights
• ff_factors: multiples of n_rise for the widenings
• ff: explicit widenings in samples; replaces ff_factors

Returns

• vector of NamedTuples (nm, za, ft, ct, zo, ff), all weights Float64 and ff an Int

dplms_flattop_matrix(windows_wide::AbstractMatrix, ff::Integer)

Flat-top matrix for a filter of n = size(windows_wide, 1) - 2ff samples, from aligned pulse shapes (one column each, e.g. events or superpulses) widened by ff samples on each side. It is the mean over columns and over the 2ff shifts of the outer product of the column's first difference, so xᵀ F x is the mean squared slope of the filter output around the pick-off. Throws if ff < 1 or if the windows are too short for the widening.

Returns

• symmetric n × n matrix

dplms_noise_matrix(wvfs::ArrayOfRDWaveforms, bl_window::ClosedInterval, n::Integer)

Noise matrix for a filter of n samples: the mean of b bᵀ over every window b of n consecutive samples inside bl_window of every waveform. The waveforms must be baseline subtracted; a residual offset appears as a rank-one term, not as noise. When the baseline mean subtracted from the waveforms was estimated over the same window, the lowest-frequency noise power is underestimated by about n divided by the number of baseline samples.

Throws if the waveforms do not share one time axis or if bl_window holds fewer than n samples.

Returns

• symmetric n × n matrix in the units of the squared waveform samples

dplms_reference(windows::AbstractMatrix, plateau::AbstractUnitRange)

Superpulse of the aligned windows (one column per event): the row mean divided by its mean over the plateau rows, so the reference has unit amplitude after the rise. Throws if that mean is not positive.

Returns

• named tuple (ref, plateau) with the reference and the divisor, the plateau amplitude of the superpulse

dplms_synthesis(noise, ref, flattop, binned; nm, za, ft, ct = 0, zo = 0)

Solve the DPLMS system

(nm·noise + ref refᵀ + za·𝟙 + ft·flattop + ct·binned.mat + zo·m mᵀ) x = ref + ct·binned.vec

for the filter coefficients x along the window and rescale them so that x ⋅ ref == 1. binned is the result of dplms_binned_reference and is used only if ct is not zero. 𝟙 is the all-ones matrix, so za penalizes the squared coefficient sum. m is the window sample index relative to the window center, divided by the window length n, so zo penalizes (mᵀx)². When the coefficients sum to zero, mᵀx is the area of the filter output for a unit step divided by n: with za > 0, zo > 0 gives an output of zero area, as for a ZAC filter, and a linear baseline slope gives no output. Throws an ArgumentError if the system is singular or the solution is not finite.

Arguments

• noise: noise matrix from dplms_noise_matrix, in units of the reference amplitude squared
• ref: reference from dplms_reference
• flattop: matrix from dplms_flattop_matrix for the same filter length
• binned: reference constraints from dplms_binned_reference

Keywords

• nm, za, ft, ct, zo: weights of the noise, zero-area, flat-top, binned-reference and output-area terms; ct and zo default to 0

Returns

• named tuple with

• coeffs: the impulse response, reverse(x), for ConvolutionFilter(DirectConvolution(), coeffs); the filter output at the last sample of a window equals x ⋅ window,

• noise_out: sqrt(xᵀ noise x), the predicted output noise in units of the reference amplitude,

• zero_area: sum(x),

• output_area: Σᵢ i·xᵢ, the area in samples of the filter output for a unit step when zero_area is zero.

nm, za, ft, ct and zo act only through their ratios: scaling all five by the same factor gives the same filter after normalization.

dsp_aux(data::Table, config::DSPConfig)

DSP chain for auxiliary channels whose waveforms are spikes. Returns a table with one row per waveform.

The config argument is accepted for consistency with the other DSP chains but is not used.

Input columns

waveform, baseline, timestamp, eventnumber, daqenergy

Output columns

• e_max: maximum sample of the waveform minus the FADC baseline
• t_max: time of the maximum sample
• blfc: baseline from the FADC
• timestamp: timestamp from the FADC
• eventID_fadc: event number from the FADC
• e_fc: energy from the FADC

dsp_aux_compressed(data::Table, config::DSPConfig)

Variant of dsp_aux for compressed data. The waveform is read from the column waveform_presummed, and e_max is the maximum sample minus the FADC baseline multiplied by presum_rate. The input also needs the column presum_rate; the output columns are the same.

dsp_cusp_ft_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, rt::Quantity{T}) where T<:Real

Energy of the CUSP filter for the rise time rt and each flat-top time in config.e_grid_ft_cusp.

The waveforms are prepared with dsp_flt_optimization_prepare (decoded, baseline mean in config.bl_window subtracted, pole-zero corrected with the decay time τ). For each flat-top time ft the energy is determined with dsp_energy_rtft from a CUSPChargeFilter with the rise time rt, the flat-top time ft, the length config.flt_length_cusp and the CR filter switched off, read at t50 + config.flt_length_cusp/2. The energies are stored as Float32.

Returns

• Matrix{Union{Missing, Float32}} of size (length(config.e_grid_ft_cusp), length(wvfs)): row f holds the energies for the f-th flat-top time. No entry is missing; an entry is NaN if the pickoff time lies outside the filtered time axis.

dsp_cusp_rt_optimization(prep::FltOptPrep, config::DSPConfig; ft::Quantity{T}=2.0u"µs") where T<:Real
dsp_cusp_rt_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}; ft::Quantity{T}=2.0u"µs") where T<:Real

Values of the CUSP filter at a fixed time for each rise time in config.e_grid_rt_cusp, for estimating the equivalent noise charge (ENC).

The waveforms are decoded, the baseline mean in config.bl_window is subtracted, and the waveforms are pole-zero corrected with the decay time τ; with a FltOptPrep as input this is already done. For each rise time rt the waveforms are filtered with a CUSPChargeFilter with the rise time rt, the flat-top time ft, the length config.flt_length_cusp and the CR filter switched off and the filtered waveform is evaluated at config.enc_pickoff_cusp by the signal estimator of prep. The method with wvfs calls dsp_flt_optimization_prepare and then the method with prep.

Arguments

• prep: prepared waveforms from dsp_flt_optimization_prepare
• wvfs, τ: waveforms and decay time for the preparation; τ and ft must have the same numeric type T

Keywords

• ft: flat-top time

Returns

• Matrix{Float64} of size (length(config.e_grid_rt_cusp), length(wvfs)): row r holds the values for the r-th rise time

dsp_decay_times(wvfs::ArrayOfRDWaveforms, bl_window::ClosedInterval{<:Unitful.Time}, tail_window::ClosedInterval{<:Unitful.Time})
dsp_decay_times(wvfs::ArrayOfRDWaveforms, config::DSPConfig)

Decay time of the exponential tail of each waveform in wvfs.

Each waveform is decoded, the mean in bl_window is subtracted, and tailstats fits a straight line to the logarithm of the samples in tail_window. The second method uses config.bl_window and config.tail_window.

Returns

• vector of decay times τ with unit µs, one entry per waveform. The entry is zero if a sample in tail_window is not positive.

dsp_dplms_energy(wvfs::ArrayOfRDWaveforms, coeffs::AbstractVector, t50::AbstractVector, t_pickoff, current_window::ClosedInterval, signal_estimator::SignalEstimator)

Apply the DPLMS filter coeffs to baseline-subtracted, pole-zero-corrected waveforms and read the energy at t50 + t_pickoff with signal_estimator.

The waveforms are decoded and filtered by ConvolutionFilter(DirectConvolution(), coeffs). Throws if coeffs is longer than a waveform or if the axes of wvfs and t50 differ.

Arguments

• wvfs: waveforms, one time axis per waveform
• coeffs: impulse response of the filter, for example coeffs of dplms_synthesis or coefficients of get_dplms_pars
• t50: alignment time of every waveform
• t_pickoff: read-out time after t50
• current_window: time window of the signal rise
• signal_estimator: SignalEstimator of RadiationDetectorDSP that reads the value at the read-out time

Returns

• named tuple (e, emax): energy at the read-out time, and the maximum of the filter output inside current_window shifted by t_pickoff, the interval that holds the read-out times of all events with t50 inside current_window

dsp_dplms_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity, drift::AbstractVector;
    n_pre = nothing, n_post = nothing, n_noise_windows = 8, ff_level = 0.999, n_bins = 5, plateau_length = 20, flat_tol = 1e-3, grid = dplms_default_grid)

Synthesize DPLMS energy filters on presummed waveforms for every point of grid and return the energies they give. A grid point is a NamedTuple with the weights nm, za, ft, ct, zo and the flat-top widening ff in samples.

Steps: baseline subtraction over config.bl_window; pole-zero correction with InvCRFilter(τ); t50, the time where the waveform reaches 50 % of the maximum of the 10/10 µs trapezoid, found with get_threshold (config.kwargs_pars.tx_mintot) as t50_pre in dsp_icpc_compressed; alignment of every waveform to the sample nearest the median t50; windows of n_pre samples before to n_post after that sample; noise matrix from the baseline window; reference and one reference per drift-time bin, the bins being n_bins quantiles of drift; flat-top matrix from the drift-bin references widened by ff samples. All matrices are in units of the reference amplitude, the mean over the last plateau_length window samples.

Unless given, n_pre and n_post equal (n_bl - n_noise_windows) ÷ 2 with n_bl the number of samples inside config.bl_window: the longest symmetric window that leaves n_noise_windows sliding windows per event for the noise matrix. A window as long as the baseline leaves the constant mode of the noise matrix undetermined, because the baseline mean is subtracted over the same region, so only zero-area filters (za > 0) are meaningful for long windows.

Every grid point yields a filter applied to the unaligned pole-zero-corrected waveforms with the energy read at t50 + t_pickoff, t_pickoff = n_post·Δt. The flat-top term penalizes the output slope of every drift-bin reference shifted by up to ff samples, so the output stays flat while any of these shapes moves by ff samples. n_rise is the number of samples after the alignment point until the slowest drift-bin reference reaches ff_level of the plateau. grid is either a vector of grid points or a function of n_rise returning one, e.g. grid = n -> dplms_default_grid(n; za = 0) for filters without the zero-area condition.

flat_top and shaping_time describe each filter's output for the full-length drift-bin references. The flat top is the interval around the read-out time where every filtered reference stays within flat_tol of its value at the read-out time; the shaping time runs from the first sample where the filtered superpulse exceeds flat_tol to the start of the flat top.

Throws if t50 lies outside config.current_window for any waveform, if the drift proxies are not finite, or if a window or plateau_length does not fit.

Arguments

• wvfs: presummed waveforms, possibly encoded, with one shared time axis
• config: DSPConfig with bl_window, current_window and kwargs_pars
• τ: decay time for the pole-zero correction
• drift: drift-time proxy of every waveform, for example drift_time

Keywords

• n_pre, n_post: window samples before and after the alignment sample
• n_noise_windows: sliding windows per event reserved for the noise matrix when n_pre and n_post are not given
• ff_level: fraction of the plateau that defines the end of the rise for n_rise
• n_bins: number of drift-time bins
• plateau_length: number of samples at the end of the window that define the reference amplitude
• flat_tol: tolerance relative to the reference amplitude for flat_top and shaping_time
• grid: vector of grid points or function of n_rise returning one

Returns

• e_grid, emax_grid: energies and filter-output maxima, size (length(grid), length(wvfs))
• coeffs: impulse response per grid point; t_pickoff
• noise_out, zero_area: predicted relative output noise and coefficient sum per grid point
• output_area: area of the filter output for a unit step per grid point, in time units
• flat_top, shaping_time: one time per grid point
• grid, ref, ref_bins (one column per bin), t50 (per event), plateau, n_pre, n_post, n_rise
• superpulse: the aligned superpulse in units of the reference amplitude
• superpulse_bins: the aligned superpulse of every drift-time bin, same units
• response: superpulse filtered with coeffs, one RDWaveform per grid point
• noise, flattops: the noise matrix and the flat-top matrix per ff (a Dict), in units of the reference amplitude

All of these are fields of one named tuple.

dsp_energy_rtft(prep::FltOptPrep, config::DSPConfig, rt::Unitful.Time, ft::Unitful.Time, flt::Union{Val, Symbol})

Energy of the prepared waveforms at rise time rt and flat-top time ft. flt is Val(:trap), Val(:cusp), Val(:zac) or the corresponding Symbol.

The trapezoidal energy is read at t50 + rt + ft/2, the CUSP and ZAC energies at t50 + flt_length/2, with flt_length of config.flt_length_cusp or config.flt_length_zac; the trapezoidal filter does not use config. The CUSP and ZAC filters are built with the CR filter switched off. rt and ft are converted to floating point of the type of the sampling period (trapezoidal) or the filter length (CUSP, ZAC). An energy is NaN if the pickoff time lies outside the filtered time axis. A filter that is longer than the waveforms throws an ArgumentError. Only reads prep, so it is safe to call concurrently; each call allocates about two (trapezoidal) or five (CUSP, ZAC) times the size of prep.wvfs_pz.

Arguments

• prep: prepared waveforms from dsp_flt_optimization_prepare
• rt, ft: rise time and flat-top time of the filter

Returns

• vector of energies, one per waveform

dsp_flt_optimization_prepare(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Unitful.Time)

Baseline-subtract and pole-zero-correct wvfs once for repeated filter evaluations.

The waveforms are decoded, the baseline mean in config.bl_window is subtracted, and the waveforms are pole-zero corrected with the decay time τ. t50 is the time at which a corrected waveform reaches 50% of its maximum, determined with config.kwargs_pars.tx_mintot as minimum time over threshold. The signal estimator is a SignalEstimator(PolynomialDNI(order, length)) with config.kwargs_pars.sig_interpolation_order and sig_interpolation_length.

Returns

• a FltOptPrep with one entry per waveform

dsp_icpc(data::Table, config::DSPConfig, τ::Quantity, pars_filter::PropDict; f_evaluate_qc::Union{Function, Missing}=missing)

DSP chain for HPGe (ICPC) detector waveforms. Returns a table with one row per waveform.

The baseline is subtracted, the waveforms are pole-zero corrected with the decay time τ, and timing, energy, current, charge-drift, quality and saturation parameters are determined. Use dsp_icpc_compressed for compressed data.

Arguments

• data: table with the input columns below; encoded columns are decoded with decode_data
• config: DSPConfig. config.kwargs_pars needs the entries fc_bit_depth, t0_flt_pars, t0_mintot, tx_mintot, intrace_mintot, int_interpolation_order, int_interpolation_length, sig_interpolation_order and sig_interpolation_length.
• τ: decay time of the preamplifier signal, used for the pole-zero correction
• pars_filter: optimized filter parameters as a PropDict: trap, cusp and zac with the entries rt and ft, and sg and mwa with the entry wl. Missing entries are taken from config.default_flt_param.

Keywords

• f_evaluate_qc: function for get_qc_classifier. With missing no classification is done and qc_label is -1.

Input columns

waveform, baseline, timestamp, eventnumber, daqenergy

Output columns

Amplitudes have the unit of the waveform samples. Times with a stated unit are converted; the others are in the unit of the waveform time axis.

Baseline and tail:

• blmean, blsigma, blslope, bloffset: statistics of the waveform in config.bl_window
• tailmean, tailsigma, tailslope, tailoffset: the same statistics of the pole-zero-corrected waveform in config.tail_window
• tail_τ, tail_mean, tail_sigma: decay time, and mean and standard deviation of the logarithm of the baseline-subtracted waveform in config.tail_window, from tailstats

Timing, from the pole-zero-corrected waveform, in µs:

• t0: start of the charge rise, from get_t0
• t10, t50, t80, t90, t99: first times at which the waveform reaches 10%, 50%, 80%, 90% and 99% of e_101010, from get_threshold
• t50_current: time at which the Savitzky-Golay derivative reaches half of its maximum
• drift_time: t90 - t0, in ns

Energies, from the pole-zero-corrected waveform:

• e_max, e_min: maximum and minimum of the baseline-subtracted waveform before the pole-zero correction
• e_10410, e_535, e_313, e_101010: maximum after a TrapezoidalChargeFilter with rise time and flat-top time of 10 and 4 µs, 5 and 3 µs, 3 and 1 µs, and 10 and 10 µs
• e_trap, e_cusp, e_zac: value of the trapezoidal, CUSP and ZAC filtered waveform at t50 + rt + ft/2 (trapezoidal) and t50 + flt_length/2 (CUSP and ZAC), determined by a signal estimator. rt and ft are the optimized rise and flat-top time of the respective filter in pars_filter, and flt_length is config.flt_length_cusp or config.flt_length_zac.
• e_trap_max, e_cusp_max, e_zac_max: maximum of the filtered waveform
• t_trap_max, t_cusp_max, t_zac_max: time of that maximum
• e_10410_inv, e_313_inv: maximum after the e_10410 and e_313 filters of the waveform multiplied by -1
• t0_inv: t0 of the waveform multiplied by -1, in µs, determined with the default flt_pars of get_t0

Current:

• a_raw: maximum in config.current_window of the simple derivative (DerivativeFilter(1))
• a_sg, t_a_sg: maximum in config.current_window of the Savitzky-Golay derivative with window length sg.wl and polynomial degree config.sg_flt_degree, and its interpolated time
• a_mwa, t_a_mwa: maximum in config.current_window of the simple derivative after MovingWindowMultiFilter(mwa.wl; passes = config.a_mwa_passes), and its interpolated time
• a_sg_<ns>, a_mwa_<ns>: the same maxima for each fixed window length in config.a_sg_windows and config.a_mwa_windows, with the window length in ns in the column name
• t_amax: interpolated time of the maximum of the simple derivative over the whole waveform

Charge drift:

• qdrift: Q-drift parameter from get_qdrift, using t0 and config.qdrift_int_length
• lq80_<early><gap><late>, lq90_<early><gap><late>, lqamax_<early><gap><late>: LQ parameters from get_lq starting at t80, t90 and t_amax, one column for each [early, gap, late] of config.lq_windows. The window lengths in the name are in units of 0.1 µs with two digits each, so lq80_250025 has an early window of 2.5 µs, no gap and a late window of 2.5 µs.

Quality and saturation:

• qc_label: label from f_evaluate_qc, evaluated on the baseline-subtracted waveform before the pole-zero correction
• inTrace_intersect, inTrace_n: time and number of in-trace pile-up triggers, from get_triggers on the Savitzky-Golay derivative with the threshold config.inTraceCut_std_threshold in standard deviations of its baseline
• n_sat_low, n_sat_high: number of samples of the unsubtracted waveform equal to 0 and to 2^fc_bit_depth - fc_bit_depth, from saturation
• n_sat_low_cons, n_sat_high_cons: longest run of consecutive samples at those values

FADC values:

• blfc, timestamp, eventID_fadc, e_fc: baseline, timestamp, event number and energy from the FADC

dsp_icpc_compressed(data::Table, config::DSPConfig, τ::Quantity, pars_filter::PropDict; f_evaluate_qc::Union{Function, Missing}=missing)

DSP chain for HPGe (ICPC) detector waveforms of compressed data. Returns a table with one row per waveform.

Each event has a presummed waveform, in which presum_rate consecutive samples are summed, and a windowed waveform with the original sampling around the signal rise. Energies, the tail and the baseline are determined from the presummed waveform, and timing, current and charge-drift parameters from the windowed waveform. The baseline mean of the presummed waveform is subtracted from it, and blmean / presum_rate from the windowed waveform. Both waveforms are then pole-zero corrected with τ.

Arguments

• data: table with the input columns below; encoded columns are decoded with decode_data. presum_rate has to be the same for all events.
• config: DSPConfig. config.kwargs_pars needs the entries fc_bit_depth, t0_flt_pars, t0_mintot, tx_mintot, int_interpolation_order, int_interpolation_length, sig_interpolation_order and sig_interpolation_length.
• τ: decay time of the preamplifier signal, used for the pole-zero correction
• pars_filter: optimized filter parameters as a PropDict, as for dsp_icpc

Keywords

• f_evaluate_qc: function for get_qc_classifier_compressed. With missing no classification is done and qc_label is -1.

Input columns

waveform_presummed, waveform_windowed, presum_rate, baseline, timestamp, eventnumber, daqenergy, t_sat_lo, t_sat_hi, deadtime

Output columns

Amplitudes have the unit of the waveform samples; those from the presummed waveform are sums over presum_rate samples. Times with a stated unit are converted; the others are in the unit of the waveform time axis.

FADC values and saturation:

• blfc, timestamp, eventID_fadc, e_fc, deadtime, t_sat_lo, t_sat_hi: copied from the input
• n_sat_low, n_sat_high, n_sat_low_cons, n_sat_high_cons: number of samples and longest run of consecutive samples of the unsubtracted presummed waveform equal to 0 and to (2^fc_bit_depth - fc_bit_depth) * presum_rate, from saturation

Baseline and tail, from the presummed waveform:

• blmean, blsigma, blslope, bloffset, bl_slope_sigma: statistics in config.bl_window, with bl_slope_sigma the standard deviation of the residuals of the slope fit
• tailmean, tailsigma, tailslope, tailoffset: the same statistics of the pole-zero-corrected waveform in config.tail_window
• tail_τ, tail_mean, tail_sigma: decay time, and mean and standard deviation of the logarithm of the baseline-subtracted waveform in config.tail_window, from tailstats
• auxpz1_mean, auxpz1_sigma, auxpz1_slope_sigma, auxpz2_mean, auxpz2_sigma, auxpz2_slope_sigma: mean, standard deviation and slope residual of the pole-zero-corrected waveform in config.auxpz1_window and config.auxpz2_window

Extrema:

• e_max, t_max: maximum and its time from the baseline-subtracted presummed waveform
• e_min, t_min: minimum and its time from the baseline-subtracted windowed waveform
• e_max_pre, e_min_pre, t_max_pre, t_min_pre: maximum, minimum and their times from the baseline-subtracted presummed waveform

Timing, in µs, from the pole-zero-corrected waveforms:

• t0: start of the charge rise of the windowed waveform, from get_t0
• t10, t50, t80, t90, t99: first times at which the windowed waveform reaches 10%, 50%, 80%, 90% and 99% of e_101010 / presum_rate
• t50_pre: first time at which the presummed waveform reaches 50% of e_101010
• drift_time: t90 - t0, in ns
• t50_current: time at which the Savitzky-Golay derivative of the presummed waveform reaches half of its maximum. The filter window length is sg.wl * presum_rate / 2.
• t_amax: interpolated time of the maximum of the simple derivative of the windowed waveform

Energies, from the presummed waveform:

• e_10410, e_535, e_313, e_101010: maximum after a TrapezoidalChargeFilter with rise time and flat-top time of 10 and 4 µs, 5 and 3 µs, 3 and 1 µs, and 10 and 10 µs
• t_10410, t_535, t_313, t_101010: times of those maxima
• e_trap, e_cusp, e_zac: value of the trapezoidal, CUSP and ZAC filtered waveform at t50_pre + rt + ft/2 (trapezoidal) and t50_pre + flt_length/2 (CUSP and ZAC), determined by a signal estimator, with the parameters as for dsp_icpc
• e_trap_max, e_cusp_max, e_zac_max, t_trap_max, t_cusp_max, t_zac_max: maximum of the filtered waveform and its time
• e_dplms: value of the presummed waveform filtered with the DPLMS filter of pars_filter.dplms at t50_pre + t_pickoff, determined by the signal estimator. Without pars_filter.dplms, the trapezoidal filter of config.default_flt_param.dplms is used, see get_dplms_pars.
• e_dplms_max: maximum of the DPLMS filtered waveform in config.current_window shifted by t_pickoff
• e_10410_inv, e_313_inv: maximum after the e_10410 and e_313 filters of the presummed waveform multiplied by -1
• t0_inv: t0 of the windowed waveform multiplied by -1, in µs, determined with the default flt_pars of get_t0

Current, from the windowed waveform:

• a_raw, a_sg, t_a_sg, a_mwa, t_a_mwa, a_sg_<ns>, a_mwa_<ns>: as for dsp_icpc

Charge drift, from the windowed waveform:

• qdrift: Q-drift parameter from get_qdrift, using t0 and config.qdrift_int_length
• lq80_<early><gap><late>, lq90_<early><gap><late>, lqamax_<early><gap><late>: LQ parameters as for dsp_icpc

Quality:

• qc_label: label from f_evaluate_qc, evaluated on the baseline-subtracted presummed waveform before the pole-zero correction
• trigger_position, trigger_multiplicity: time and number of in-trace triggers, from get_triggers on the presummed waveform filtered with TrapezoidalChargeFilter(2u"µs", 1u"µs"), with the threshold config.inTraceCut_std_threshold in standard deviations of its baseline and a minimum time over threshold of 2 µs

dsp_mwa_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, pars_filter::PropDict; f_evaluate_qc::Union{Function, Missing}=missing) where T<:Real

Energy and A/E for each window length in config.a_grid_wl_mwa of the moving-window-average current filter.

Same as dsp_sg_optimization, except for the current: the pole-zero-corrected waveforms are filtered with DerivativeFilter(1) and then with MovingWindowMultiFilter(wl; passes = config.a_mwa_passes) for each window length wl, and the maximum in config.current_window is the current. This is the chain of the column a_mwa of dsp_icpc.

dsp_mwa_optimization_compressed(wvfs_wdw::ArrayOfRDWaveforms, wvfs_pre::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, pars_filter::PropDict; presum_rate::Real=T(8), f_evaluate_qc::Union{Function, Missing}=missing) where T<:Real

Variant of dsp_mwa_optimization for compressed data.

The arguments and the split of the work between wvfs_pre and wvfs_wdw are the same as for dsp_sg_optimization_compressed.

dsp_pmts(data::Table, config::PropDict)

DSP chain for PMT waveforms. Returns a table with one row per waveform.

The time axis is replaced using TimeAxisFilter, the baseline is subtracted, and the extrema, triggers, saturation and the extrema of a smoothed waveform are determined.

Arguments

config needs the entries

• time_axis_step_length: sampling period of the new time axis
• baseline_window_start, baseline_window_end: baseline window
• min_tot_intersect, max_tot_intersect, intersect_threshold: minimum time over threshold, maximum time over threshold for the maximum, and threshold of the trigger search with IntersectMaximum
• wsg_window_length, wsg_flt_degree, wsg_weight: window length, polynomial degree and weight type of the smoothing filter. For wsg_weight == 0 this is a SavitzkyGolayFilter without derivative, otherwise a WeightedSavitzkyGolayFilter
• saturation_limit_low, saturation_limit_high: sample values counted as saturated

Input columns

waveform, timestamp, channel, eventnumber, daqenergy

Output columns

• timestamp, eventID_fadc, e_fc, channel: from the input
• raw_pulse_height, raw_pulse_low: maximum and minimum of the baseline-subtracted waveform
• raw_t0_hi, raw_t0_low: times of that maximum and minimum
• trig_max, trig_t: maxima and times of the triggers above intersect_threshold, one vector per waveform
• trig_mult: number of triggers
• sat_low, sat_high: number of samples of the unsubtracted waveform equal to saturation_limit_low and saturation_limit_high
• pulse_height, pulse_low, t0_hi, t0_low: the same as the raw_ values, determined on the smoothed waveform
• bl_mean, bl_sigma, bl_slope: baseline statistics in the baseline window

dsp_puls(data::Table, config::DSPConfig)

DSP chain for pulser waveforms. Returns a table with one row per waveform.

The baseline is the fit in config.bl_window. It is subtracted before the maximum, the time at 50% of the maximum and the energy of a trapezoidal filter are determined.

Input columns

waveform, baseline, timestamp, eventnumber, daqenergy

Output columns

• blmean, blsigma, blslope, bloffset: baseline statistics in config.bl_window
• t50: time at which the waveform reaches 50% of its maximum, in µs
• e_max: maximum of the baseline-subtracted waveform
• e_10410: maximum after a TrapezoidalChargeFilter(10u"µs", 4u"µs") (rise time 10 µs, flat-top time 4 µs)
• blfc: baseline from the FADC
• timestamp: timestamp from the FADC
• eventID_fadc: event number from the FADC
• e_fc: energy from the FADC

dsp_puls_compressed(data::Table, config::DSPConfig)

Variant of dsp_puls for compressed data. The waveform is read from the column waveform_presummed instead of waveform; the other input and all output columns are the same.

dsp_qc_flt_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Unitful.Time, f_evaluate_qc::Function)
dsp_qc_flt_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Unitful.Time, f_evaluate_qc::Missing)

Energy, baseline and quality-cut label of each waveform with the default trapezoidal filter, used to select waveforms for the filter-parameter scans.

The waveforms are decoded, the baseline mean in config.bl_window is subtracted, and the waveforms are pole-zero corrected with the decay time τ. The energy is the value of a TrapezoidalChargeFilter(rt, ft) with rt and ft of config.default_flt_param.trap, determined by a signal estimator at t50 + rt + ft/2. t50 is the time at which the waveform reaches 50% of its maximum. The signal estimator and t50 are configured with config.kwargs_pars (sig_interpolation_order, sig_interpolation_length and tx_mintot).

The label is get_qc_classifier(wvfs, f_evaluate_qc) of the baseline-subtracted and pole-zero-corrected waveforms. With missing as f_evaluate_qc no classification is done and qc_label is -1.

Returns

• table with the columns energy, blmean, blsigma, blslope (per unit of the waveform time axis), t50 and t50_pre (in µs) and qc_label, one row per waveform. t50_pre is the time at which the waveform reaches 50% of the maximum of a TrapezoidalChargeFilter(10u"µs", 10u"µs"), as t50_pre in dsp_icpc_compressed and t50 in dsp_dplms_optimization.

dsp_qc_flt_optimization_compressed(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Unitful.Time, f_evaluate_qc::Function)
dsp_qc_flt_optimization_compressed(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Unitful.Time, f_evaluate_qc::Missing)

Variant of dsp_qc_flt_optimization for presummed waveforms. The label is determined with get_qc_classifier_compressed; everything else is the same.

dsp_qdrift_flt_optimization(wvfs::ArrayOfRDWaveforms, blmean::Vector{<:Real}, config::DSPConfig, τ::Unitful.Time)

Q-drift parameter of each waveform, for scanning the integration times config.qdrift_int_length.

The waveforms are decoded, blmean is subtracted, and the waveforms are pole-zero corrected with τ. Then t0 is determined with get_t0 (threshold config.t0_threshold, config.kwargs_pars.t0_flt_pars and t0_mintot) and the parameter with get_qdrift (integration times config.qdrift_int_length, config.kwargs_pars.int_interpolation_order and int_interpolation_length).

Arguments

• blmean: baseline mean of each waveform

Returns

• vector of Q-drift values, one per waveform

dsp_sg_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, pars_filter::PropDict; f_evaluate_qc::Union{Function, Missing}=missing) where T<:Real

Energy and A/E for each window length in config.a_grid_wl_sg of the Savitzky-Golay current filter.

The waveforms are decoded and the baseline mean in config.bl_window is subtracted. The quality-cut label is determined from these waveforms. The waveforms are then pole-zero corrected with τ. The energy is the value of a TrapezoidalChargeFilter(rt, ft) at t50 + rt + ft/2, determined by a signal estimator, where t50 is the time at which the waveform reaches 50% of its maximum. For each window length wl the current is the maximum in config.current_window of the waveform filtered with SavitzkyGolayFilter(wl, config.sg_flt_degree, 1). A/E is the current without its unit divided by the energy.

Arguments

• pars_filter: needs pars_filter.trap.rt and pars_filter.trap.ft, the rise and flat-top time of the energy filter

Keywords

• f_evaluate_qc: function for get_qc_classifier. With missing no classification is done and qc_label is -1.

Returns

• table with one row per waveform and the columns aoe (A/E for all window lengths of the grid, as a vector per waveform), energy, blmean, blslope, t50 (in µs) and qc_label

dsp_sg_optimization_compressed(wvfs_wdw::ArrayOfRDWaveforms, wvfs_pre::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, pars_filter::PropDict; presum_rate::Real=T(8), f_evaluate_qc::Union{Function, Missing}=missing) where T<:Real

Variant of dsp_sg_optimization for compressed data.

The baseline, quality-cut label, t50 and energy are determined from the presummed waveforms wvfs_pre. The current is determined from the windowed waveforms wvfs_wdw, from which blmean / presum_rate is subtracted. The output table has the same columns; t50 is the time of the presummed waveform.

Keywords

• presum_rate: number of samples summed in the presummed waveforms
• f_evaluate_qc: function for get_qc_classifier_compressed

dsp_sg_sipm_optimization_compressed(wvfs::ArrayOfRDWaveforms, dsp_config::PropDict, optimization_config::PropDict)
dsp_sg_sipm_optimization_compressed(n_max_wvfs::Int, wvfs::ArrayOfRDWaveforms, dsp_config::PropDict, optimization_config::PropDict)

Savitzky-Golay filter window lengths scanned for the SiPM trigger threshold and the maxima of the triggers found with it.

For each window length in optimization_config.e_grid_wl the waveforms are filtered with a SavitzkyGolayFilter(wl, dsp_config.sg_flt_degree, 1). The threshold is dsp_config.n_σ_threshold times thresholdstats of the samples of the first optimization_config.threshold.n_wvfs waveforms (all waveforms if there are fewer) between optimization_config.threshold.min_cut and optimization_config.threshold.max_cut. The triggers are found with IntersectMaximum(dsp_config.min_tot_intersect, dsp_config.max_tot_intersect).

Arguments

• n_max_wvfs: number of waveforms per chunk. The method with this argument is meant to process the waveforms in chunks of this size, to merge the maxima of the chunks per window length and to take the smallest threshold of the chunks. It throws a MethodError.
• wvfs: SiPM waveforms
• dsp_config: PropDict with the entries min_tot_intersect, max_tot_intersect, n_σ_threshold and sg_flt_degree
• optimization_config: PropDict with the entries e_grid_wl and threshold (with min_cut, max_cut and n_wvfs)

Returns

• table with one row per window length and the columns trig_max_grid (maxima of all triggers of all waveforms) and thresholds_grid (threshold)

dsp_sg_sipm_thresholds_compressed(wvfs::ArrayOfRDWaveforms, sg_window_length::Unitful.Time, config::PropDict)

Samples of the Savitzky-Golay derivative and of its integral for all waveforms, as input for choosing the SiPM trigger thresholds.

The waveforms are decoded, converted to Float64 and filtered with a SavitzkyGolayFilter(sg_window_length, config.sg_flt_degree, 1). The derivative is then integrated with IntegratorFilter(gain=1).

Arguments

• wvfs: SiPM waveforms
• sg_window_length: window length of the Savitzky-Golay filter
• config: PropDict with the entry sg_flt_degree

Returns

• table with the columns bsl_deriv (samples of the derivative), bsl (samples of the integrated derivative) and bsl_flipped (-bsl). Each column contains the samples of all waveforms in one vector.

dsp_sipm(data::Table, config::PropDict, pars_optimization::PropDict)

DSP chain for SiPM waveforms. Returns a table with one row per waveform.

Input columns

waveform, baseline, timestamp, eventnumber, daqenergy

Arguments

• data: table of waveforms; encoded columns are decoded with decode_data
• config: PropDict with the entries sg_flt_degree, t0_hpge_window (start and stop of the window in which the HPGe signal is expected, used for the _lar parameters) and filters. filters.sg holds the settings of the Savitzky-Golay pipeline and filters.trap those of the trapezoidal pipeline. Both contain min_tot_intersect, max_tot_intersect, min_threshold, max_threshold, n_σ_threshold, min_dc_threshold, max_dc_threshold and n_σ_dc_threshold; filters.trap also contains rt, ft and pz_tau. filters.xt_tag holds the settings of the cross-talk tag pipeline: window_length, n_σ_threshold, min_tot_intersect and max_tot_intersect.
• pars_optimization: PropDict with the optimized Savitzky-Golay window length pars_optimization.sg.wl

Output columns

Times are in the unit of the waveform time axis unless noted.

• blfc, timestamp, eventID_fadc, e_fc: baseline, timestamp, event number and energy from the FADC
• t_max, t_min: times of the waveform maximum and minimum, in µs
• t_max_lar, t_min_lar: the same within t0_hpge_window, in µs
• e_max, e_min: waveform maximum and minimum
• e_max_lar, e_min_lar: the same within t0_hpge_window
• blmean, blsigma, blslope, bloffset: statistics of the integrated Savitzky-Golay derivative from the start of the waveform to the first Savitzky-Golay trigger (to the end of the waveform if there is none, but at least three samples)
• wfmean, wfsigma, wfslope, wfoffset: the same statistics over the whole integrated waveform
• threshold_sg: trigger threshold of the Savitzky-Golay pipeline
• threshold_DC_sg: discharge threshold of the Savitzky-Golay pipeline
• trig_pos_sg, trig_max_sg: times and maxima of the triggers, one vector per waveform
• trig_pos_DC_sg, trig_max_DC_sg: times and maxima of the discharge triggers
• threshold_trap, threshold_DC_trap: trigger and discharge thresholds of the trapezoidal pipeline
• trig_pos_trap, trig_pos_high_trap, trig_pos_tot_trap, trig_max_trap: time, falling-edge time, time over threshold and maximum of the trapezoidal triggers
• trig_pos_DC_trap, trig_pos_high_DC_trap, trig_pos_tot_DC_trap, trig_max_DC_trap: the same for the trapezoidal discharge triggers
• threshold_xt_tag: threshold of the cross-talk tag pipeline
• trig_pos_xt_tag, trig_pos_high_xt_tag, trig_pos_tot_xt_tag, trig_max_xt_tag: start time, end time, time over threshold and maximum of each region of the long-window integrated waveform that is above the threshold

Processing

• The waveform is converted to Float64.
• Savitzky-Golay pipeline: the waveform is filtered with a SavitzkyGolayFilter(pars_optimization.sg.wl, config.sg_flt_degree, 1), which gives the smoothed derivative. The threshold is n_σ_threshold times thresholdstats_mad of the samples between min_threshold and max_threshold, and the triggers are found with IntersectMaximum.
• The derivative is integrated with IntegratorFilter(gain=1), which gives a smoothed waveform, and the baseline statistics are determined on it.
• Discharge detection: the median of the integrated waveform is taken as its baseline and subtracted, and the result is multiplied by -1. The thresholds are computed in the same way as above from the dc settings of filters.sg, and the triggers are found with the Savitzky-Golay IntersectMaximum.
• Trapezoidal pipeline: the integrated waveform is pole-zero corrected with pz_tau and filtered with a TrapezoidalChargeFilter(rt, ft). Thresholds and triggers are computed in the same way as for the Savitzky-Golay pipeline, with the settings of filters.trap.
• Trapezoidal discharge detection: the thresholds are computed from the flipped, baseline-subtracted integrated waveform with the dc settings of filters.trap, and the triggers are found with the IntersectMaximum of the trapezoidal settings (filters.trap.min_tot_intersect and max_tot_intersect).
• Cross-talk tag pipeline: the waveform is filtered with a SavitzkyGolayFilter(filters.xt_tag.window_length, config.sg_flt_degree, 1) and integrated, and the median of the result is subtracted. The threshold is filters.xt_tag.n_σ_threshold times thresholdstats_mad of all samples, without lower and upper bounds. The regions above the threshold are found with an IntersectMaximum(filters.xt_tag.min_tot_intersect, filters.xt_tag.max_tot_intersect); the waveform is not flipped.

dsp_sipm_compressed(data::Table, config::PropDict, pars_optimization::PropDict)

Variant of dsp_sipm for compressed data. The waveform is read from the column waveform_bit_drop instead of waveform; the arguments, the processing and the output columns, including the cross-talk tag columns, are the same as for dsp_sipm.

Input columns

waveform_bit_drop, baseline, timestamp, eventnumber, daqenergy

dsp_trap_ft_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, rt::Quantity{T}) where T<:Real

Energy of the trapezoidal filter for the rise time rt and each flat-top time in config.e_grid_ft_trap.

The waveforms are prepared with dsp_flt_optimization_prepare (decoded, baseline mean in config.bl_window subtracted, pole-zero corrected with the decay time τ). For each flat-top time ft the energy is determined with dsp_energy_rtft from a TrapezoidalChargeFilter(rt, ft), read at t50 + rt + ft/2. The energies are stored as Float32.

Returns

• Matrix{Union{Missing, Float32}} of size (length(config.e_grid_ft_trap), length(wvfs)): row f holds the energies for the f-th flat-top time. No entry is missing; an entry is NaN if the pickoff time lies outside the filtered time axis.

dsp_trap_rt_optimization(prep::FltOptPrep, config::DSPConfig; ft::Quantity{T}=2.0u"µs") where T<:Real
dsp_trap_rt_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}; ft::Quantity{T}=2.0u"µs") where T<:Real

Values of the trapezoidal filter at a fixed time for each rise time in config.e_grid_rt_trap, for estimating the equivalent noise charge (ENC).

The waveforms are decoded, the baseline mean in config.bl_window is subtracted, and the waveforms are pole-zero corrected with the decay time τ; with a FltOptPrep as input this is already done. For each rise time rt the waveforms are filtered with a TrapezoidalChargeFilter(rt, ft) and the filtered waveform is evaluated at config.enc_pickoff_trap by the signal estimator of prep. The method with wvfs calls dsp_flt_optimization_prepare and then the method with prep.

Arguments

• prep: prepared waveforms from dsp_flt_optimization_prepare
• wvfs, τ: waveforms and decay time for the preparation; τ and ft must have the same numeric type T

Keywords

• ft: flat-top time

Returns

• Matrix{Float64} of size (length(config.e_grid_rt_trap), length(wvfs)): row r holds the values for the r-th rise time

dsp_zac_ft_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}, rt::Quantity{T}) where T<:Real

Energy of the ZAC filter for the rise time rt and each flat-top time in config.e_grid_ft_zac.

The waveforms are prepared with dsp_flt_optimization_prepare (decoded, baseline mean in config.bl_window subtracted, pole-zero corrected with the decay time τ). For each flat-top time ft the energy is determined with dsp_energy_rtft from a ZACChargeFilter with the rise time rt, the flat-top time ft, the length config.flt_length_zac and the CR filter switched off, read at t50 + config.flt_length_zac/2. The energies are stored as Float32.

Returns

• Matrix{Union{Missing, Float32}} of size (length(config.e_grid_ft_zac), length(wvfs)): row f holds the energies for the f-th flat-top time. No entry is missing; an entry is NaN if the pickoff time lies outside the filtered time axis.

dsp_zac_rt_optimization(prep::FltOptPrep, config::DSPConfig; ft::Quantity{T}=2.0u"µs") where T<:Real
dsp_zac_rt_optimization(wvfs::ArrayOfRDWaveforms, config::DSPConfig, τ::Quantity{T}; ft::Quantity{T}=2.0u"µs") where T<:Real

Values of the ZAC filter at a fixed time for each rise time in config.e_grid_rt_zac, for estimating the equivalent noise charge (ENC).

The waveforms are decoded, the baseline mean in config.bl_window is subtracted, and the waveforms are pole-zero corrected with the decay time τ; with a FltOptPrep as input this is already done. For each rise time rt the waveforms are filtered with a ZACChargeFilter with the rise time rt, the flat-top time ft, the length config.flt_length_zac and the CR filter switched off and the filtered waveform is evaluated at config.enc_pickoff_zac by the signal estimator of prep. The method with wvfs calls dsp_flt_optimization_prepare and then the method with prep.

Arguments

• prep: prepared waveforms from dsp_flt_optimization_prepare
• wvfs, τ: waveforms and decay time for the preparation; τ and ft must have the same numeric type T

Keywords

• ft: flat-top time

Returns

• Matrix{Float64} of size (length(config.e_grid_rt_zac), length(wvfs)): row r holds the values for the r-th rise time

extremestats(signal::SamplesOrWaveform, start::RealQuantity, stop::RealQuantity)
extremestats(signal::RDWaveform)
extremestats(signal::AbstractSamples)

Minimum, maximum and their time positions of signal between start and stop.

The interval includes the samples nearest to start and stop. The methods without an interval use the whole signal.

Returns

• named tuple (min, max, tmin, tmax)

get_a(wvfs::ArrayOfRDWaveforms, current_window::ClosedInterval{<:Unitful.Time{<:Real}})

Maximum of the current waveforms wvfs in current_window, quadratically interpolated around the maximum sample (get_wvf_maximum). Returns the maxima max and their timepoints t.

get_dplms_pars(pd::PropDict, config::DSPConfig, Δt::Unitful.Time)

DPLMS filter for dsp_dplms_energy as (; coefficients, t_pickoff), with t_pickoff the read-out time after t50.

Returns the entries coefficients and t_pickoff of pd.dplms if present. Otherwise builds the trapezoidal filter of config.default_flt_param.dplms (rt, ft) sampled at Δt, read out at rt + ft/2; it gives the same energies as TrapezoidalChargeFilter(rt, ft). Throws an ArgumentError if pd.dplms lacks coefficients or t_pickoff, if there is neither pd.dplms nor config.default_flt_param.dplms, or if the default rise time is shorter than Δt.

Arguments

• pd: filter parameters of the detector, as pars_filter of dsp_icpc_compressed
• config: DSPConfig
• Δt: sampling interval of the presummed waveforms

get_lq(wvfs_int::ArrayOfRDWaveforms, t_start::Array{<:Unitful.Time}, early::Unitful.Time, gap::Unitful.Time, late::Unitful.Time; pol_power::Int=3, sign_est_length::Unitful.Time=100u"ns")

LQ parameter: charge in the late window [t_start + early + gap, t_start + early + gap + late] minus the charge in the early window [t_start, t_start + early].

wvfs_int are the integrated waveforms (IntegratorFilter(1).(wvfs)), so that the integration is done once for all LQ variants. The charges are determined with a SignalEstimator(PolynomialDNI(pol_power, sign_est_length)).

get_qdrift(wvfs, t, Δt) equals get_lq(wvfs_int, t, first(Δt), 0u"µs", last(Δt) - first(Δt)).

get_qc_classifier(wvfs::ArrayOfRDWaveforms, f_evaluate_qc::Function)
get_qc_classifier(wvfs::ArrayOfRDWaveforms, f_evaluate_qc::Function, config::DSPConfig)

Quality-cut labels of the waveforms wvfs, evaluated by the support vector machine function f_evaluate_qc.

Each waveform is filtered five times with HaarAveragingFilter(2), which reduces the number of samples by a factor of 32, and divided by the larger of the absolute values of its minimum and maximum (waveforms that are zero everywhere are not scaled). The filtered signals form the columns of a matrix that is passed to f_evaluate_qc. The method with config subtracts the baseline mean in config.bl_window first; the other method uses wvfs as given.

Arguments

• f_evaluate_qc: function that returns the tuple (labels, decision_values) for a matrix with one column per waveform, as returned by get_qc_ml_func

Returns

• vector of labels, one per waveform

get_qc_classifier_compressed(wvfs::ArrayOfRDWaveforms, f_evaluate_qc::Function)
get_qc_classifier_compressed(wvfs::ArrayOfRDWaveforms, f_evaluate_qc::Function, config::DSPConfig)

Variant of get_qc_classifier for presummed waveforms: each waveform is filtered twice instead of five times with HaarAveragingFilter(2), which reduces the number of samples by a factor of 4.

get_qc_ml_func(dwts_norm::Matrix{<:Real}, dc_labels::Vector{<:Real}, hyperparams::PropDict)

Train a support vector machine with a radial basis function kernel and return a function that classifies waveforms.

Arguments

• dwts_norm: training features, one column per waveform
• dc_labels: training label of each waveform
• hyperparams: PropDict with the LIBSVM settings cost, gamma, weights (class weights, keyed by the class label), probability, cache_size, coef0, shrinking and tolerance

Returns

• function f(x) for a feature matrix x with one column per waveform. It returns the tuple (labels, decision_values) of LIBSVM.svmpredict. Pass it as f_evaluate_qc to the DSP chains and to get_qc_classifier.

get_qdrift(wvfs::ArrayOfRDWaveforms, t_start::Array{<:Unitful.Time}, Δt::StepRangeLen{<:Unitful.Time}; pol_power::Int=3, sign_est_length::Unitful.Time=100u"ns")

Q-drift parameter of each waveform in wvfs: the difference of two charges collected after t_start.

The waveforms are integrated with IntegratorFilter(1). Values of the integrated waveform at t_start, t_start + first(Δt) and t_start + last(Δt) are determined with a SignalEstimator(PolynomialDNI(pol_power, sign_est_length)). The result is the charge collected in [t_start + first(Δt), t_start + last(Δt)] minus the charge collected in [t_start, t_start + first(Δt)].

Arguments

• t_start: start time for each waveform, typically t0
• Δt: only the first and the last element are used, as the two integration times relative to t_start

Keywords

• pol_power: polynomial degree of the signal estimator
• sign_est_length: length of the signal estimator window

get_t0(wvfs_pz::ArrayOfRDWaveforms, t0_threshold::Real; flt_pars::Vector{<:Unitful.Time}=[40u"ns", 100u"ns", 2000u"ns"], mintot::Unitful.Time=1500u"ns")

Start time t0 of the charge rise of each pole-zero-corrected waveform in wvfs_pz.

Each waveform is filtered with an asymmetric TrapezoidalChargeFilter(flt_pars...). t0 is the interpolated time at which the filtered waveform first rises above t0_threshold and stays above it for at least mintot.

Keywords

• flt_pars: pre-rise averaging time, gap time and post-rise averaging time of the trapezoidal filter
• mintot: minimum time over threshold

Returns

• vector of t0 with unit µs, one entry per waveform. The entry is zero if the threshold is not crossed.

get_threshold(wvfs::ArrayOfRDWaveforms, threshold::Union{Array{T}, T}; mintot::Unitful.Time=1000.0u"ns") where T<:Real

Time at which each waveform in wvfs first rises above threshold.

threshold is either one value for all waveforms or an array with one value per waveform. The time is interpolated linearly, and the waveform has to stay above the threshold for at least mintot.

Returns

• vector of times with unit µs, one entry per waveform. The entry is zero if the threshold is not crossed.

get_triggers(wvfs::ArrayOfRDWaveforms, sigma_threshold::Real, bl_window::ClosedInterval{<:Unitful.Time}; mintot::Unitful.Time=100.0u"ns")

Trigger positions and multiplicities found on the waveforms.

The threshold of each waveform is sigma_threshold times its standard deviation in the baseline window, computed on the time axis of the first waveform from first(wvfs[1].time) + leftendpoint(bl_window) to rightendpoint(bl_window). A threshold of zero is replaced by one. The waveforms are searched with Intersect(mintot).

Returns

• named tuple (intersect, n): intersect holds the time of the first crossing (zero if there is none), n the number of crossings of at least mintot.

get_triggers_reversed(wvfs::ArrayOfRDWaveforms, sigma_threshold::Real, bl_window::ClosedInterval{<:Unitful.Time}; mintot::Unitful.Time=100.0u"ns")

Trigger positions and multiplicities found on the time-reversed waveforms.

The threshold of each waveform is sigma_threshold times its standard deviation in the baseline window, computed on the time axis of the first waveform from first(wvfs[1].time) + leftendpoint(bl_window) to rightendpoint(bl_window). A threshold of zero is replaced by one. The waveforms are reversed in time and searched with Intersect(mintot).

Returns

• named tuple (intersect, n): intersect holds the time of the first crossing in the original time axis, n the number of crossings of at least mintot.

get_wvf_maximum(signal::SamplesOrWaveform, start::RealQuantity, stop::RealQuantity)

Maximum of signal between start and stop and its time, from a parabola through the maximum sample and its two neighbors.

The interval includes the samples nearest to start and stop. If the maximum sample is the first or last sample of the interval, its value and time are returned without interpolation. Throws a DimensionMismatch if the time axis and the signal have different axes and an ArgumentError if the interval leaves the record.

Returns

• named tuple (max, t)

realized_flattop(ft::Unitful.Time, flt_length::Unitful.Time, sampling_period::Unitful.Time, flt::Union{Val, Symbol})

Flat-top time that the filter flt (Val(:trap), Val(:cusp), Val(:zac) or the corresponding Symbol) applies when built with flat-top time ft on waveforms sampled with sampling_period. flt_length is the CUSP/ZAC filter length and is unused for :trap.

RadiationDetectorDSP rounds filter times to whole samples. The trapezoidal gap is round(ft / sampling_period) samples. For CUSP and ZAC the flat top is lengthened by one sample if the filter length minus the flat top, in samples, is odd. The result is a multiple of sampling_period.

Returns

• the realized flat-top time, a quantity in the unit of sampling_period

saturation(signal::SamplesOrWaveform, low::RealQuantity, high::RealQuantity)
saturation(signal::SamplesOrWaveform, start::RealQuantity, stop::RealQuantity, low::RealQuantity, high::RealQuantity)

Count the samples of signal that are equal to the lower saturation value low or the upper saturation value high.

The first method uses the whole signal, the second the interval between start and stop, including the samples nearest to both.

Returns

• named tuple (low, high, max_cons_low, max_cons_high): number of samples equal to low and high, and the length of the longest run of consecutive samples equal to low and to high.

tailstats(signal::SamplesOrWaveform, start::RealQuantity, stop::RealQuantity)

Fit a straight line to the logarithm of signal between start and stop, as a model of an exponential decay.

The interval includes the samples nearest to start and stop.

Returns

• named tuple (mean, sigma, τ): mean and standard deviation of the logarithm of the samples in the interval, and the decay time τ, the negative inverse slope of the fit. If any sample in the interval is not positive, all three values are zero.

thresholdstats(signal::SamplesOrWaveform, min::RealQuantity=-Inf, max::RealQuantity=Inf)

Standard deviation of all samples of signal that lie between the lower bound min and the upper bound max.

The default bounds include all samples. The result is NaN if no sample lies between the bounds.

thresholdstats_mad(signal::SamplesOrWaveform, min::RealQuantity=-Inf, max::RealQuantity=Inf)

Median absolute deviation (MAD) of all samples of signal that lie between the lower bound min and the upper bound max, scaled by 1.4826 to estimate the standard deviation of normally distributed data.

The default bounds include all samples. The result is zero if no sample lies between the bounds.

_create_dsp_config(dsp_metadata::PropDicts.PropDict)

Create a DSPConfig from a PropDict of DSP metadata.

Arguments

• dsp_metadata::PropDicts.PropDict: DSP metadata

Returns

• dsp_config::DSPConfig: DSP configuration

_lq_column_name(t_start::Symbol, windows)

Column name of an LQ parameter: lq<start>_<early><gap><late> with the start timepoint without its leading t/t_ and the three window lengths in units of 0.1 µs with two digits, e.g. _lq_column_name(:t80, [2.5, 0.0, 2.5]u"µs") == :lq80_250025.

banded_cholesky_solve!(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)

Solve A Aᵀ x = y and store the result in x, where A is the lower Cholesky factor of a banded symmetric matrix in the storage of choleskyL!: A[1, i] is the diagonal, A[2, i] the first subdiagonal, and so on.

choleskyL!(b::AbstractMatrix{T})

In-place Cholesky decomposition of a banded symmetric matrix in the storage of makeIPlusLambdaDprimeD: b[1, i] holds the diagonal, b[2, i] the first subdiagonal, and so on.

Throws an AssertionError if the matrix is not positive definite.

extrema3points(y1::T, y2::T, y3::T)::T where {T<:RealQuantity}

Value at the extremum of the parabola through three equidistant samples y1, y2 and y3.

get_fltpars(pd::PropDict, flt::Symbol, dspconfig::DSPConfig)

Filter parameters for the filter flt from pd, with missing entries taken from dspconfig.default_flt_param.

Returns the window length pd.sg.wl or pd.mwa.wl for flt equal to :sg or :mwa, and the tuple (rt, ft) of rise time and flat-top time for :trap, :cusp and :zac.

makeFitWeights(d::Integer, m::Int)

Weights of the linear fits that extend the signal beyond its boundaries in ModifiedSincFilter.

The weights are the square root of a cosine that falls to zero at a distance determined by the half-width m of the kernel and the degree d of the filter.

makeIPlusLambdaDprimeD(λ::T, p::Int, N::Int) where {T}

Banded storage of the centrosymmetric N × N matrix I + λ DᵀD, where D is the finite-difference matrix of order p.

Row d + 1 of the returned (p + 1) × N matrix holds the d-th subdiagonal. Throws an AssertionError unless N > p.

weighted_linear_reg(w::AbstractVector, y::AbstractVector)

Weighted linear regression of the data y with weights w. The abscissa is assumed to be 0, 1, …, length(y) - 1.

Returns the tuple (offset, slope). The slope is zero if it evaluates to NaN, for example for a single point, and NaN if the sum of the weights is not positive.