# Sparse Deconvolution Methods for Online Energy Estimation in Calorimeters
Operating in High Luminosity Conditions

Tiago Teixeira  Luciano Andrade  and José Manoel de Seixas

###### Abstract

Energy reconstruction in calorimeters operating in high luminosity particle
colliders has become a remarkable challenge. In this scenario, pulses from a
calorimeter front-end output overlap each other (pile-up effect), compromising
the energy estimation procedure when no preprocessing for signal
disentanglement is accomplished. Recently, methods based on signal
deconvolution have been proposed for both online and offline reconstructions.
For online processing, constraints concerning fast processing, memory
requirements, and cost implementation limit the overall performance. Offline
reconstruction allows the use of Sparse Representation theory to implement
sophisticated Iterative Deconvolution methods. This paper presents Iterative
Deconvolution methods based on Sparse Representation algorithms whose
computational cost is effective for online implementation. Using simulated
data, current techniques were compared to the proposed Sparse Representation
ones for performance validation in the online environments. Analysis has shown
that, despite the higher computational cost, when compared to standard
methods, the performance improvement may justify the use of the proposed
techniques, in particular for the Separable Surrogate Functional, which is
shown to be feasible for implementation in modern FPGAs.

##  1 Introduction

High-Energy Physics (HEP) experiments usually require sophisticated
instrumentation systems, especially when particle colliders are involved [1].
In order to analyze increasingly rare physical phenomena, there is a tendency
of increasing their collision rates and the number of particles per collision,
thus increasing the luminosity achieved by the experiment [2]. This is the
case of the most potent collider ever built, the LHC (Large Hadron Collider)
[3] at CERN (Organisation Européenne pour la Recherche Nucléaire) [4]. The LHC
collides bunches of protons every 25 nanoseconds and has been establishing
records in terms of both collider energy and luminosity [5]. Moreover,
upgrades that are on the way for the next LHC operation periods (Run 3 and
beyond) will increase even more the luminosity [6], towards a projected peak
around 5 to 7.5 $\times 10^{34}cm^{-2}s^{-1}$ [7].

For new physics identification at the LHC, calorimeters [8] play an important
role, as they measure the energy of the incoming particles and are essential
for triggering purposes [9]. The two general-purpose LHC experiments (ATLAS
[10] and CMS [11]) comprise hundred of thousands of calorimeter readout
channels each and employ state-of-art detector technology. However, the LHC’s
enormous luminosity imposes severe restrictions on calorimetry performance as,
despite being fast, the signal generated in the detector’s front-end
electronics requires a certain number of bunches crossings (BCs) to be fully
developed. Therefore, signal pile-up effects may arise when the luminosity is
increased.

The trigger system is often implemented through an online cascade signal
processing chain consisting of two or more levels, where the rate reduction is
proportional to the level of the algorithm complexity for final event
selection. Due to the high event rate, the lower filtering levels are usually
implemented in hardware and rely very much on online energy estimation, which
is addressed in this paper for applications where signal pile-up becomes an
issue.

###  1.1 The Pile-up Effect

Figure 1: Example of the pile-up effect. See text.

An example of out-of-time signal pile-up in collider experiments can be seen
in Figure 1, where signals were digitized by a 40 MHz clock ($25ns$ period, as
in the LHC). The unipolar signal generated in a hypothetical calorimeter cell
by a given collision is referred to as $s_{1}$, and $50ns$ later a new signal
($s_{2}$) from two bunch collisions after reaches that particular cell.
Supposing a $150ns$ readout acquisition length, the $s_{3}$ signal is formed,
which is actually the stacking of both $s_{1}$ and $s_{2}$ signals. This pile-
up effect leads to a sub-optimal operation of the classical energy estimation
algorithms, which are based on Optimal Filtering (OF) theory [12, 13, 14, 15].

The calorimeter signal shape is a linear combination of a signal with
neighboring pileup signals, where the signal amplitude is proportional to the
target energy value [8]. Therefore, the sensor output can be modeled as the
convolution between the target sequence (energy deposition per bunch crossing)
and the calorimeter normalized (unitary amplitude) reference pulse shape. To
recover the target sequence, a deconvolution algorithm may be applied [17].

Different deconvolution algorithms have been proposed for energy estimation in
high-luminosity calorimeter operation. The work in [16] proposes an offline
procedure in two steps, where a signal detection per bunch crossing is
performed first, followed by a multi-amplitude estimation procedure. Still
focusing on offline reconstruction, the work in [18] implements an approach
based on Sparse Representation (SR) [19] of data, presenting remarkable
results, over-performing standard deconvolution methods. Sparse Representation
was originally developed and applied in signal compression and image denoising
routines [20, 21]. The work in [18] inaugurated the Sparse Representation
energy reconstruction applications. However, its cost implementation is
presently forbidden for online data-stream.

Concerning online processing, a deconvolution method based on Finite Impulse
Response (FIR) [22] filters was proposed in [23] and it has shown to be very
attractive due to its lower cost implementation in hardware. However, recently
Sparse Representation theory has developed families of algorithms targeting
cost implementation reduction [24, 25, 26]. Although FIR filter structures are
known to be the most economic approach in terms of hardware resources, the
cost-effectiveness of such modern Sparse Representation techniques has paved
the way for online energy reconstruction implementation. Those approaches have
already been used in several online applications like machine learning [27]
and big data [28].

Therefore, the main goal of this paper is to use modern Sparse Representation
theory in a design that fulfills online requirements. The adaptation of some
of those techniques is proposed envisaging energy estimation in calorimeter
cells, opening a new research area in online energy reconstruction. Since
Sparse Representation has shown to be the most accurate energy reconstruction
procedure under pile-up, the methods proposed here can be used as an
alternative for FIR filter implementation in applications where the hardware
resource minimization is not the main constraint. A comparison among the
aforementioned methods, in terms of cost implementation and energy
reconstruction performance, is also provided, indicating the best situation
for the usage of each approach.

In order to complete the discussion, a possible implementation for one of the
proposed Sparse Representation methods (in this case, the Separable Surrogate
Functional) is outlined at the end of the paper. The circuitry was developed
in Hardware Description Language (HDL) [29] envisaging a straightforward
implementation in Field-Programmable Gate Array (FPGA) technology [30].

##  2 Techniques

In this section, firstly the FIR approach for high-luminosity online energy
estimation is briefly reviewed, as it is being used here for baseline
comparison. In the sequence, the basics of modern sparse signal processing are
presented, as a way to exploit the nonlinear statistics that comes with the
out-of-time signal pileup.

###  2.1 FIR Filter Based Energy Estimation

Figure 2 shows the basic structure for online energy estimation based on a FIR
deconvolutional filter [23]. The signal $r[n]$ represents the measured
information from a single front-end calorimeter sensor, after applying the
free-running (an uninterrupted) Analog-to-Digital Converter (ADC). A cascade
of registers, synchronous with the ADC clock, implements a shift-register
structure [30], which is responsible for storing the latest $M$ digitized
samples. The filter order ($M-1$) depends on both the calorimeter reference
pulse-shape width and the signal pile-up intensity. For deconvolution, the
filter order is often higher then the reference pulse width (two times on
average), since the filter needs tail information of neighbor superimposed
signals to perform the desired source separation. The equation employed in
this deconvolution procedure is:

$x^{\prime}[n]=\stackrel{{\scriptstyle[}}{{l}}=0]{M-1}{\sum}(w[n-l]r[n])$
(2.1)

Figure 2: The FIR filter structure implementing signal deconvolution in online
energy estimation, based on [23].

The weighting vector $\mathbf{w}$ is obtained offline, based on deconvolution
principles. Details on the filter design can be found in the original paper
[23]. The samples are, then, linearly combined at the filter output, which
delivers new information at each clock. A threshold is applied to the filter
output to avoid small artifacts from the deconvolution process in adjacent
bunch crossings, mainly small negative energy estimations.

The system output corresponds to the reconstructed target information sequence
$\hat{x}[n]$, which is a straightforward estimation of the energy deposit per
bunch crossing ($x[n]$). The system starts delivering the output sequence
$\Delta$ clock cycles after the first sample goes through the processing
chain. This delay is necessary because the system needs information from
previous samples in order to perform the deconvolution procedure. Typical
values of $\Delta$ are around half of the filter order.

Despite of its very suitable hardware implementation, the purely linear
operation of FIR filters does not allow access to the Higher Order Statistics
(HOS) [31] present in the embedded pile-up noise, which may be exploited for
further performance improvement. The proposed Sparse Representation based
deconvolution methods access the Higher Order Statistics information
indirectly through the usage of non-linear functions in the deconvolution
process.

###  2.2 The Sparse Representation

The Sparse Representation theory has been applied for energy reconstruction
with great success for offline environments in [18]. On that paper, the
implementation was based on Linear Programming (LP) [33], which is a standard
procedure with high computational cost. Recently, modern Sparse Theory has
raised with more cost-effective algorithms [19]. Modifications on standard
Matching Pursuit (MP) methods [34] have presented an intrinsic relationship
with energy estimation tasks. In those approaches, a signal detection
iterative procedure is executed first, using a bank of linear Matching-Filter
(MF) [35] based filters. After some iterations, all the superimposed signals
within an acquisition window are detected. Their amplitudes are, then,
estimated simultaneously using the Least Square (LS) algorithm [36].

Another Sparse Representation family of methods that is being employed due to
its excellent cost-benefit properties is the Iterative-Shrinkage (IS)
algorithms [37]. These methods combine standard linear Coordinate-Descent (CD)
algorithms [38] with a linear-by-part Shrinkage function, in order to promote
sparsity, and becoming one of the most cost-effective algorithms for Sparse
Representation.

##  3 Proposed Sparse Representation Based on Deconvolution Methods

The Sparse Representation based deconvolution methods are designed by using a
windowed acquisition point of view. The method is guided by the matrix
formulation from the convolution operation of finite size sequences, as
proposed in [18]. Given a normalized calorimeter reference pulse-shape,
represented by the vector $\mathbf{h}$ of size $V$ and a target sequence,
represented by the vector $\mathbf{x}$ (energy per bunch crossing) of size
$G$, the convolution between these two signals produces the measured signal
$\mathbf{r}$; that is, the calorimeter readout channel response. Because there
are $G$ bunch crossings, and the fact that the calorimeter responses extends
$V$ samples beyond the last bunch crossing, the total samples considered for
the vector $\mathbf{r}$ are $G+V-1$ [22]. This convolution process can be
written as:

$\mathbf{r}=\mathbf{H}\mathbf{x}$ (3.1)

where $\mathbf{H}$ is the convolution matrix of size $U\times G$ (columns of
$\mathbf{H}$ contain shifted versions of the vector $\mathbf{h}$, with zero
padding). The windowed deconvolution procedure consists of recovering the
sequence $\mathbf{x}$ when $\mathbf{r}$ is measured and $\mathbf{H}$ is known.
When the size of the vector $\mathbf{r}$ is greater than the one of
$\mathbf{x}$, there are infinite solutions for this system of equations. In
[18], it was shown that, for deconvolution purposes, the sparsest
representation of $\mathbf{x}$ is the correct solution when the target
sequence is characterized by impulsive signals.

According to the Sparse Representation theory [19], the sparsest solution for
Equation 3.1 can be found solving the optimization Problem $P_{\ell}$ for
$0\leq\ell\leq 1$

$(P_{\ell}):\underset{\mathbf{x}}{min}\left\|\mathbf{x}\right\|_{\ell}^{\ell}\>subjected\,to\>\mathbf{r}=\mathbf{Hx}$
(3.2)

where

$\left\|\mathbf{x}\right\|_{\ell}^{\ell}=\underset{b}{\sum}\left|x_{b}\right|^{\ell}$
(3.3)

is the $\ell$-norm of the vector $\mathbf{x}$.

For $\ell=1$, this results in a typical Linear Programming problem [32]. The
use of Linear Programming for Sparse Representation aiming at energy
estimation has been outlined in [18], where this method was shown to
outperform other windowed deconvolution methods. However, the focus of that
paper was on offline energy reconstruction, since Linear Programming presents
a very high computational cost, becoming usually prohibitive for online
execution. In the sequence, two families of Sparse Representation algorithms
are presented in order to solve the $P_{\ell}$ problem with different
perspectives, both focusing on computational cost-effectiveness procedures.

###  3.1 Matching Pursuit

For $\ell=0$, Equation 3.3 results in the number of non-zero elements in
$\mathbf{x}$ and the Sparse Representation solution given by $P_{0}$ becomes
evident. However, $P_{0}$ is not a mathematically tractable problem [19]. This
is a combinatorial problem infeasible to be computationally implemented in
practice.

To avoid the test of all combinatorial possibilities, Matching Pursuit seeks
for similarities between columns of the matrix $\mathbf{H}$ and the vector
$\mathbf{r}$, selecting the smallest number of columns that, when linearly
combined, produce the vector $\mathbf{r}$ added to a small residue vector
$\mathbf{d}$.

Given the support matrix $\mathbf{H}_{s}$, which comprises the selected
columns of $\mathbf{H}$, the residue vector $\mathbf{d}$ is

$\mathbf{d}=\mathbf{H}_{s}\mathbf{x}_{s}-\mathbf{r}$ (3.4)

where the vector $\mathbf{x}_{s}$ has the amplitudes for each signal in
$\mathbf{H}_{s}$ (the columns of $\mathbf{H}_{s}$). The components in
$\mathbf{x}_{s}$ are chosen in order to minimize the square module of the
residue:

$\epsilon^{2}=\mathbf{d}^{T}\mathbf{d}=|\mathbf{H}_{s}\mathbf{x}_{s}-\mathbf{r}|^{2}$
(3.5)

where $T$ indicates a vector transposition.

The vector $\mathbf{x}_{s}$ that minimizes the Equation 3.5 is obtained by the
Least Square equation

$\mathbf{x}_{s}=(\mathbf{H}_{s}^{T}\mathbf{H}_{s})^{-1}\mathbf{H}_{s}^{T}\mathbf{r}$
(3.6)

The Matching Pursuit algorithms select columns of the matrix $\mathbf{H}$ with
the highest potential of minimizing $\epsilon^{2}$ [39]. Depending on the way
this task is performed, a specific variation of the basic Matching Pursuit
algorithm is applied. In this work, the Orthogonal Matching Pursuit (OMP) [40]
and the Least Square Orthogonal Matching Pursuit (LS-OMP) [41] will be
outlined.

####  3.1.1 Orthogonal Matching Pursuit

This method selects the column $\mathbf{h}_{j}$ of the matrix $\mathbf{H}$
that best matches the current residue $\mathbf{d}$, minimizing the cost
function

$Q(x_{j})=|x_{j}\mathbf{h}_{j}-\mathbf{d}|^{2}$ (3.7)

The value $x_{j}$ that minimizes this function is the Matching Pursuit
equation

$x_{j}=\frac{\mathbf{h}_{j}^{T}\mathbf{d}}{\mathbf{h}_{j}^{T}\mathbf{h}_{j}}$
(3.8)

This value is computed for all $\mathbf{h}_{j}$ columns not yet inserted in
the support matrix $\mathbf{H}_{s}$. The column with the smallest $Q(x_{j})$
value is selected, the new residue vector is computed and a new iteration is
performed until $\epsilon$ becomes lower than a given threshold
$\epsilon_{0}$. It can be shown [41] that all the columns of $\mathbf{H}_{s}$
are orthogonal to the remaining residue $\mathbf{d}$.

####  3.1.2 Least Squares Orthogonal Matching Pursuit

The LS-OMP method updates temporarily $\mathbf{H}_{s}$ with the current column
$\mathbf{h}_{j}$ in test and computes $\epsilon$ directly (Equation 3.5). The
temporary matrix $\mathbf{H}_{s}$ that gives the smallest $\epsilon$ value is
kept. The process stops when the current residue is lower than $\epsilon_{0}$
as well. Different from OMP, which computes a simple inner product (Matching
Pursuit equation) to each column not yet inserted in $\mathbf{H}_{s}$, LS-OMP
computes the Least Square (Equation 3.6) instead. Therefore, the computational
cost of LS-OMP is higher than OMP.

###  3.2 Iterative-Shrinkage

Iterative-Shrinkage algorithms are based on the Problem $P_{1}$. However, a
relaxation is employed in order to allow a small residual square error
$\epsilon_{0}$ on the constraint.

$(P_{1,\epsilon_{0}}):\underset{\mathbf{x}}{min}\left\|\mathbf{x}\right\|_{1}^{1}\>subjected\,to\>|\mathbf{H}\mathbf{x}-\mathbf{r}|^{2}<\epsilon_{0}$
(3.9)

The problem $P_{1,\epsilon_{0}}$ is a typical Linear Programming problem, as
directly performed in [18] for offline signal reconstruction. However,
focusing on cost-effectiveness, firstly, this problem is transformed to an
optimization routine without constraint by the use of a Lagrange multiplier
$\lambda$ [19]

$(P_{1,\lambda}):\underset{\mathbf{x}}{min}\left\|\mathbf{x}\right\|_{1}^{1}+\lambda|\mathbf{H}\mathbf{x}-\mathbf{r}|^{2}$
(3.10)

where $\epsilon_{0}$ is absorbed by $\lambda$. The Iterative-Shrinkage
algorithms propose solving Problem $P_{1,\lambda}$ by means of iterative
Coordinate-Descent methods: given an initial guess vector $\mathbf{x}^{0}$,
the optimal $\mathbf{x}_{opt}$ value can be recursively inferred in a small
number of steps. However, the use of standard Coordinate-Descent methods, like
Newton-Raphson [42], becomes prohibitive for $P_{1,\lambda}$ because of the
discontinuity of the $\ell$-norm for $\ell=1$ [19].

In [43], the insertion of terms in $P_{1,\lambda}$ is proposed, so that the
minimum coordinate position $\mathbf{x}_{opt}$ is not changed, but the
multivariate problem is split into separated one-dimensional problems, which
can be solved by parts. The resulting equation is named as Surrogate Function,
and the respective Iterative-Shrinkage method is the Separable Surrogate
Functional (SSF). After some algebraic manipulation, the iterative procedure
can be compacted as follows

$\mathbf{x}^{i+1}=\Gamma_{\lambda}\Bigg{(}\mathbf{x}^{i}+\mu[\mathbf{H}^{T}(\mathbf{r}-\mathbf{H}\mathbf{x}^{i})]\Bigg{)}$
(3.11)

where $\mu$ is the step size in the direction of the minimum. The function
$\Gamma_{\lambda}(\theta)$ is the one-dimensional shrinkage function that
should be applied to each component of the vector $\mathbf{x}$. In this paper,
a modification on the shrinkage function is proposed in order to promote
positive energy estimation, as shown in Figure 3. This approach increases the
performance of the method for further iterations, since negative (and wrong)
estimations are discarded along the process. With this modification, the
resulting function can be simply implemented through a subtraction followed by
a threshold ($\lambda$) operation. The argument in Equation 3.11 is identified
as a linear Gradient-Descent (GD) iteration [44] and comprises only addition
and multiplication operations. The parameters $\mu$ and $\lambda$ are both
obtained from simulation.

Figure 3: The original shrinkage function and the proposed modified version to
promote positive energy estimation.

####  3.2.1 Separable Surrogate Functional Initialization

One important prerequisite for online implementation of the Separable
Surrogate Functional algorithm is to avoid a higher number of iterations. This
can be accomplished by choosing a proper initial vector $\mathbf{x}^{0}$. In
this work, a preprocessing routine is proposed, in order to initialize the
solution close to the optimal one. The proposal is to implement the Least
Square (Equation 3.6) for a full-support matrix $\mathbf{H}$:

$\mathbf{x}^{0}=(\mathbf{H}^{T}\mathbf{H})^{-1}\mathbf{H}^{T}\mathbf{r}=\mathbf{H}^{+}\mathbf{r}$
(3.12)

where $\mathbf{H}^{+}=(\mathbf{H}^{T}\mathbf{H})^{-1}\mathbf{H}^{T}$ is the
pseudo-inverse of matrix $\mathbf{H}^{T}$ [45]. This is a constant matrix,
which may be obtained offline.

##  4 Performance Results

In order to validate the proposed algorithms, a toy Monte Carlo simulation was
used for generating a set of consecutive samples, which represent the target
energy values, per bunch crossing, for a single readout channel. Energy
distributions and sensor’s front-end were simulated following the
characteristics found in modern calorimeters [46, 47, 48]. The pile-up effect
was assessed through the calorimeter cell occupancy: null occupancy
corresponds to the absence of hits in a given run and on the other extreme,
100% of occupancy means that an amount of energy is deposited at every
particle collision that happens in a given run.

###  4.1 Physical Model

Concerning hardware implementation for the target online operation, Sparse
Representation and FIR filters comprise very different approaches. The latter
produces its output signal in a sample-by-sample sequence, which is very
suitable and straightforward for free-running environments. On the other hand,
Sparse Representation algorithms operate simultaneously on all the samples
previously buffered within a window. Actually, this non-causal characteristics
is one of the reasons of the increasing in performance, when compared to the
FIR filter structures. Therefore, operating a windowed energy reconstruction
algorithm in a free-running environment is one of the challenges addressed by
the present work, as the online implementation of windowed methods demands
extra digital circuitry for proper data-flow. In this paper, the proposed
Sparse Representation methods are first demonstrated. Then, a hardware
solution for online implementation is outlined.

Actually, the adaptation of such a non-causal method to online operation is
possible thanks to some characteristics of modern colliders, as is the case of
the LHC, which operates in a train-of-bunches [49, 50]. Due to particle
injection stability issues, the accelerator needs small periods of empty-
bunches intercalated with periods of bunch-trains [51]. At Run 2, for
instance, typical physics runs were operated in a configuration of 48 filled
bunches followed by 7 (or more) empty bunches [52]. It turns out that this
issue encourages the use of windowed estimation methods as the ones proposed
here, since such an operational approach eliminates two important design
difficulties:

  1. 1.

Window size determination \- in adapting windowed methods to free-running
environments, the window size is to be determined. In a train-of-bunches
configuration, it is natural to choose the window size in order to fit an
entire train within the processing window.

  2. 2.

Border effect \- A problem found in the windowing signal processing is that
signals on the window border are split into neighbor windows. As a
consequence, the estimation of bunch crossings closer to the window edge tends
to reduce the processing efficiency. However, due to the train-of-bunches
configuration, providing that the calorimeter reference signal width is
smaller than the empty-bunches period, the border effect does not occur.

In this paper, bunch-trains of 48 bunches followed by 7 empty bunches were
used. Thus, window lengths of 55 bunches (48 central plus the guard comprising
three bunches on each side of the window and an extra bunch between windows)
were employed to avoid the border effect.

Without loss in generality, the cell calorimeter shaped pulse was considered
as unipolar [53, 54, 55], Gaussian-like, with an acquisition length of
$150~{}ns$. The collision rate was of $40~{}MHz$ (as in the LHC), resulting in
a seven-sample pulse representation where the pulse peak corresponded to the
central sample ($4^{th}$ sample). It is worth mentioning that the proposed
method might be applied to any pulse shape, provided it represents the impulse
response ($h[n]$) to be deconvolved.

For optimal parameter determinations and for performance evaluations, the RMS
of the error between the target ($x[n]$) and the estimated ($\hat{x}[n]$)
sequences was employed (Equation 4.1)

$RMS~{}Error=\sqrt{\frac{\sum_{n=1}^{N}(\hat{x}[n]-x[n])^{2}}{N}}$ (4.1)

where $N$ was the total number of samples in the test set (10,010).

###  4.2 Parameters Choice

The simulation framework employed here was also used in recent related works
[18, 23], where more details can be found. The signal characteristics were:

  * •

The cell occupancy range: from 1% to 90%, in steps of approximately 10%.

  * •

Energy deposition: from an exponential distribution of 360 MeV.

  * •

Phase shift: an uniformly random pulse-phase in the range $\pm 1$ ns is
employed to each generated signal.

  * •

Electronic noise: zero-mean white Gaussian noise with $\sigma=15$ MeV.

For general filter design and performance comparisons, signal data from 1,820
windows (with 55 bunches each) were considered, forming a data stream with
100,100 samples. The simulation data were split into ten subsets of ten
thousand and ten events each, in order to evaluate the statistical
fluctuations. The error bars were computed through a cross-validation
procedure [56]. In this procedure, one subset was kept for testing and did not
participate in the filter design, which used the remaining nine subsets. The
testing subset was shifted for each fold, so that each subset could play once
the role of the testing subset along the cross-validation procedure. The
average performance and the Root Mean Square (RMS) value from the ten models
developed along such a procedure were computed as the estimate of the expected
performance and the corresponding error bar, respectively.

####  4.2.1 FIR Filter Design

The baseline implementation was proposed in [23] (see Figure 2). For filter
design, two parameters must be tuned: the filter order and the threshold value
for reducing the estimation error.

The energy estimation performance as a function of the filter order is shown
in Figure 4. In this figure, the output threshold was set to null aiming at
avoiding negative energy estimations. Based on this figure, a filter order of
22 was kept fixed for further analysis (hence, the filter delay was
$\Delta=11$ samples), regardless of the occupancy value, since the RMS error
did not decrease for higher filter order values. For optimal threshold value
determination, the estimation error as a function of the threshold value is
plotted in Figure 5 for each occupancy value. One can conclude that a
threshold around 30 MeV might be used regardless of the occupancy level. As
expected, the estimation error increased with the pile-up level.

Figure 4: The estimation error as a function of the FIR filter order
(threshold set as zero). The black dots point out the FIR Filter order values
for which the minimum RMS errors were achieved.

Figure 5: The estimation error as a function of the threshold value determined
in the output of the FIR filter. The black dots point out the threshold values
for which the minimum RMS errors were achieved.

####  4.2.2 Matching Pursuit Parameters

The unique parameter to be determined for the Matching Pursuit methods is the
squared value for the residue limit ($\epsilon_{0}^{2}$). Figure 6 shows such
error limits as a function of this parameter for both the OMP and the LS-OMP
methods. Since this parameter is measured in MeV2, its squared root value was
employed for the sake of physical interpretation in terms of energy values in
MeV. As expected, the optimal residue value decreased as the occupancy level
was increased, requiring more signals in the support matrix $\mathbf{H}_{s}$.
Although both methods achieved similar performance, it is worth mentioning
that the amount of operations for LS-OMP implementation is much higher, as it
executes a complete estimation with Least Square for each iteration whereas
the OMP method applies the linear regression only at the final stage.

(a)

(b)

Figure 6: The RMS error as a function of the residue energy for both OMP (a)
and LS-OMP (b) methods. The black dots point out the threshold values for
which the minimum RMS errors were achieved.

####  4.2.3 Separable Surrogate Functional Parameters

For Separable Surrogate Functional, two parameters require tuning: the step-
size $\mu$ for the Coordinate-Descent iterations and the $\lambda$ value for
the shrinkage function. For faster convergence, the $\mu$ value may be
determined dynamically, depending on the current residue $\mathbf{d}$. Since
the argument in Equation 3.11 is a linear Gradient-Descent iteration step, the
optimal step-size $\mu$ for the $i$-th iteration is [57]

$\mu=\frac{\mathbf{d}^{iT}\mathbf{d}^{i}}{\mathbf{d}^{iT}(\mathbf{H}^{T}\mathbf{H})\mathbf{d}^{i}}$
(4.2)

where $\mathbf{d}^{i}=(\mathbf{H}\mathbf{x}^{i}-\mathbf{r}$) is the residue
vector at the $i$-th iteration.

Figure 7 shows the evolution of the RMS error as a function of the iteration
index, for four different $\mu$ values, when the occupancy level was kept
fixed at 30%. When evaluating the dependency with respect to the parameter
$\mu$, the other parameter ($\lambda$) was set to zero. For these plots, the
entire set (100,100 samples) was used to increase the statistics and no error
bar is shown.

The four evaluated values for $\mu$ were: the dynamic one (according to
Equation 4.2) and three fixed values, in order to avoid the computation of
Equation 4.2 in hardware. The three fixed values were chosen within the range
of the values obtained through the dynamic approach, so that no multiplication
circuitry would be necessary in a fixed-point numerical representation for any
of these values in hardware (bit-shifting operation only) [58]. One may notice
that $\mu=0.5$ was too high and the algorithm turned out to be unstable. On
the other extreme, for $\mu=0.125$ the convergence was too slow. The best
solution was found for $\mu=0.25$, which produced an evaluation performance
similar to the dynamical one (see Figure 7). Figure 8 summarises the
performance for the Separable Surrogate Functional method as a function of the
number of iterations for different occupancy levels, when $\mu=0.25$ was kept
fixed.

Figure 7: The RMS error as a function of the number of iterations. A cell with
30% of occupancy level was employed.

Figure 8: The RMS error as a function of the number of Separable Surrogate
Functional iterations for several occupancy levels and $\mu=0.25$.

Concerning the tuning for the parameter $\lambda$, Figure 9 shows how the
estimation performance behaved around the optimal value for several occupancy
levels. It can be seen that the value of the parameter $\lambda$ tended to
zero while the pile-up level increased. Therefore, for sake of hardware
simplicity, the parameter $\lambda$ was set to zero, regardless of the
occupancy level.

Figure 9: The RMS error as a function of the parameter $\lambda$ in the
Separable Surrogate Functional method, considering several occupancy levels
and $\mu=0.25$. A total of 330 iterations were applied in order to guarantee
full convergence of the algorithm. The black dots refer to the optimal values.

The results in Figure 8 were obtained by initializing the vector
$\mathbf{x}^{0}$ with the 48 central elements of the measured vector
$\mathbf{r}$. This was shown to be faster than initializing $\mathbf{x}^{0}$
as a null vector. However, as can be seen from this figure, more than 200
iterations were still necessary for algorithm convergence for most of the
occupancy levels. Thus, in this work, the initialization of $\mathbf{x}$ is
proposed to follow Equation 3.12 and Figure 10 compares the evaluation of the
estimation error as the number of iterations increases, with and without this
preprocessing scheme, for a cell with $30\%$ of occupancy level. One can
notice that the convergence proved to be faster when the proposed
initialization was applied.

Figure 10: The RMS error as a function of the number of Separable Surrogate
Functional iterations with and without the preprocessing step for
$\mathbf{x}^{0}$ determination. A cell with $30\%$ of occupancy level was
employed.

###  4.3 Performance

Once the parameters for the aforementioned methods were all set to their
design tuned values, a performance comparison was carried out. Table 1
summarizes the final performance. In this table, Linear Programming (as
described in [18]) corresponds to a direct implementation of the problem
$P_{1,\epsilon_{0}}$ (Equation 3.9), and FIR corresponds to the method
proposed in [23]. The last three columns refer to the methods proposed in this
paper.

Table 1: Performance estimation in terms of the RMS error values (in MeV).
Best performance for each occupancy level is shown in boldface.

| Occupancy
---
level
| Linear
---
Programming
| FIR
---
Filter
OMP | LS-OMP | SSF∗
1 % | 4.44$\pm 1.07$ | 15.22$\pm 2.33$ | 3.27$\pm 0.75$ | 2.86$\pm 0.49$ | 3.72$\pm 0.33$
5 % | 10.13$\pm 0.67$ | 18.48$\pm 0.91$ | 8.77$\pm 0.56$ | 7.29$\pm 0.35$ | 9.21$\pm 0.62$
10 % | 14.9$\pm 0.93$ | 22.3$\pm 0.74$ | 13.95$\pm 0.9$ | 11.8$\pm 0.91$ | 13.28$\pm 0.84$
20 % | 22.12$\pm 0.83$ | 29.76$\pm 0.68$ | 23.27$\pm 0.64$ | 19.36$\pm 0.79$ | 20.51$\pm 0.92$
30 % | 28.47$\pm 0.9$ | 35.33$\pm 0.69$ | 30.79$\pm 1.03$ | 27.12$\pm 0.73$ | 26.26$\pm 0.82$
40 % | 34.32$\pm 0.92$ | 40.27$\pm 0.65$ | 37.26$\pm 0.86$ | 33.45$\pm 0.93$ | 31.67$\pm 0.9$
50 % | 39.03$\pm 0.91$ | 44.12$\pm 0.7$ | 43.25$\pm 1.09$ | 39.91$\pm 0.83$ | 36.57$\pm 0.68$
60 % | 45.6$\pm 0.97$ | 47.72$\pm 0.8$ | 48.17$\pm 0.77$ | 45.67$\pm 0.97$ | 41.01$\pm 0.98$
90 % | 59.21$\pm 0.87$ | 54.15$\pm 0.76$ | 57.68$\pm 0.88$ | 57.68$\pm 0.88$ | 52.63$\pm 0.65$

∗Separable Surrogate Functional (SSF) parameters set as $\mu=0.25$ and
$\lambda=0$

For lower occupancy levels, the target sequence ($x[n]$) presented high
sparsity and the Sparse Representation methods outperformed the FIR filter
considerably. As the occupancy level was increased, the sparsity assumption
gradually lost its strength. Therefore, the performance of the Sparse
Representation methods got closer to the one from FIR filters for occupancy
values above 30%. Among the Matching Pursuit methods, LS-OMP was slightly
better. However, the Matching Pursuit methods were less efficient with respect
to Separable Surrogate Functional for higher occupancy levels. The first step
of the Matching Pursuit algorithms, the vector support determination,
presented lower efficiency when the sparsity of the target sequence was not so
high. Concerning the offline method used for comparison (Linear Programming),
the achieved results were similar to the ones from LS-OMP.

In general, the occupancy levels shown in Table 1 may be considered practical
for different designs in modern collider’s experiments like ATLAS and CMS, in
LHC, and for technologies which are being considered for future colliders.
Apart from the nominal collider luminosity, the occupancy level is also
strongly correlated to the detector geometry parameters, such as pseudo-
rapidity and radial layer. For forward calorimeters (with very high pseudo-
rapidity), occupancy values at 90% level are already found in the LHC
detectors, while in low pseudo-rapidity regions, the occupancy near 1% is the
one mostly found. However, as the luminosity of the colliders is foreseen to
increase continuously, the necessity of using methods more resilient to pile-
up effects, such the one proposed here, becomes more evident. In the Run 3,
for example, the luminosity is expected to increase considerably.

####  4.3.1 Detection Efficiency

In this section, performance evaluation in terms of detection efficiency is
outlined. The analysis was performed through the Receiver Operating
Characteristic (ROC) curve [59], which shows how the Probability of Detection
(PD) for a given signal varies when a detection threshold is applied in the
energy estimation value and the corresponding impact in the False Alarm (FA)
probability. Here, PD comprised the percentage of bunches crossings with hits
that were correctly detected (estimated energy above the threshold) and FA was
the percentage of bunches crossings with no hits that have been wrongly
selected. Figure 11 presents the ROC curves for two different occupancy levels
(1% and 30%). The two tested Matching Pursuit methods (OMP and LS-OMP)
presented similar results. Thus, only LS-OMP is shown in this analysis. For 1%
of occupancy, both FIR filter and Separable Surrogate Functional methods
presented similar detection efficiency. As the occupancy was increased, the
advantage of the Sparse Representation methods, when compared to the FIR
filter, became more evident. The Separable Surrogate Functional method
remained as one with the highest detection efficiency, regardless of the
occupancy level (1% or 30%).

(a) 1% of occupancy

(b) 30% of occupancy

Figure 11: The ROC curves for low (a) and high (b) occupancy values.

####  4.3.2 Analysis Summary

From the above analysis, one can conclude that, for lower occupancy
conditions, if only the detection efficiency is taken into account (triggering
purpose only), the FIR filter structure is a competitive choice, as it can be
seen from the ROC curves. On the other hand, if the estimation accuracy is
also an important issue, Matching Pursuit methods may be employed. For high
pile-up levels, the Separable Surrogate Functional method may be chosen when
both estimation accuracy and detection efficiency are of main concern.
Besides, even for lower occupancy values, the Separable Surrogate Functional
method was competitive with respect to general energy estimation performance.

The Separable Surrogate Functional method has also the lowest computational
cost among the Sparse Representation one. Therefore, the next section
evaluates an online implementation of the Separable Surrogate Functional
method in FPGA technology and does a comparison with the FIR filter
implementation in terms of technology resource requirements.

##  5 Hardware Implementation and Cost Analysis

Aiming at optimizing the hardware resource utilization, the argument in
Equation 3.11 is rewritten according to:

$\mathbf{x}^{i}+\mu(\mathbf{H}^{T}\mathbf{r}-\mathbf{H}^{T}\mathbf{H}\mathbf{x}^{i})$
(5.1)

where $\mathbf{H}^{T}\mathbf{H}$ is a constant Toeplitz matrix [60], which is
computed offline. Concerning logic circuitry allocation, the Toeplitz matrix
presents a certain symmetry that allows to reduce the cost through parallel
implementation. In the text and in the figures that follow, this matrix is
identified as $\mathbf{Y}$ for sake of simplicity. For the same reason, the
vector $\mathbf{H}^{T}\mathbf{r}$ is also previously computed in hardware,
just before the first iteration, and it will be indicated as a constant vector
$\mathbf{f}$.

Other important simplification concerns the non-linearity imposed by the
shrinkage function (Figure 3). For $\lambda=0$, this function is implemented
by a simple logic comparison: if the value is lower than zero, the component
of the vector $\mathbf{x}$ is replaced by zero before the next iteration.

Figure 12 shows the block diagram for Separable Surrogate Functional algorithm
implementation. The blocks in gray color execute matrix computations
altogether, in parallel (distributed combinational circuitry) and in fixed-
point arithmetic [61]. In this figure, it is possible to identify a system
with three pipeline stages [62], split into input and output register banks,
where:

Figure 12: The block diagram of the proposed circuitry. It operates
synchronously with the $r[n]$ input signal.

  1. 1.

The first pipeline stage comprises the SIPO (Serial-Input Parallel-Output)
circuitry, the block which computes the $\mathbf{f}=\mathbf{H}^{T}\mathbf{r}$
operation and the preprocessing for $\mathbf{x}^{0}$ computation. The SIPO
assemblies a window of 55 samples (the vector $\mathbf{r}$) buffering the most
recent consecutive $r[n]$ samples. When the SIPO internal registers are fully
occupied, both vectors $\mathbf{f}$ and $\mathbf{x}^{0}$ are computed. In the
next clock rising edge, the 48 resulting components, for both vectors, are
stored in input register banks. This process is repeated to each window (55
clock cycles) uninterruptedly.

  2. 2.

The second pipeline stage comprises the blocks between the input and the
output register banks. The ones in gray color are responsible to execute one
Separable Surrogate Functional iteration. This stage has a feedback path,
synchronous to the data clock (40 MHz, if LHC would be the target
application), in order to compute several iterations before the next window is
delivered by the first pipeline stage. The resulting reconstructed vector
($\mathbf{x}^{i}$) is updated in the output register bank to each clock.

  3. 3.

The third and last stage is executed by the PISO (Parallel-Input Serial-
Output) block. This circuitry serialises the reconstructed vector present in
the output register bank, delivering the free-running $\hat{x}[n-\Delta]$
information, synchronous with the ADC clock (40 MHz). In this way, the
windowed processing of the proposed system becomes transparent.

###  5.1 Firmware Implementation

As described before, the system operates in a three stage pipeline
architecture. Thus, three consecutive windows are evaluated simultaneously,
one at each stage. The first stage buffers the most recent window and sends
the data, in parallel, to the next stage.

The Separable Surrogate Functional iterations are implemented by the second
pipeline stage. In the first clock, the multiplexer is switched to the input
register. In this way, the initial value of the algorithm is $\mathbf{x}^{0}$,
as expected. This signal propagates through the blocks that compute one
Separable Surrogate Functional iteration. For the consecutive 55 clock cycles,
the multiplexer points to the output register, which updates the reconstructed
vector, in a feedback data-flow, at each clock cycle. In this stage, if only
one Separable Surrogate Functional iteration block is implemented, it is
possible to compute 55 iterations before the first stage sends a new fresh
window. Therefore, a cascade of $z$ blocks allows the execution of $55\times
z$ iterations.

The last stage simply serializes the most recent reconstructed window
delivered by the second stage. As a consequence of this proposed architecture,
the response delay $\Delta$ is kept fixed in 2 times the window size ($2\times
55=110$), regardless the number of iterations (which is a multiple of 55).

###  5.2 Cost Analysis

The proposed circuitry was implemented in FPGA. The Software Vivado [63] was
used in the design, and Verilog [64] was the Hardware Description Language
(HDL) [29]. For performance tests, the chip XCZU7EV-FFVC1156-2-E from the Zynq
UltraScale+ family [65] was chosen, since its evaluation is free and can be
crosschecked without costs. For representation, 17 bits (7 bits for range and
10 bits for precision) were selected for fixed-point arithmetic, for all the
parts of the circuitry. This configuration was chosen empirically after
exhaustive simulation with several different numerical representations. Table
2 summarizes the implementation performance. As for comparison, the results
for a FIR filter structure are also outlined, and the details about the
hardware implementation of the FIR filter may be found in [23]. In this table,
it is shown the performance results for up to three cascaded iteration blocks.
For more than this, the operation frequency is below 40 MHz. Among the
Separable Surrogate Functional implementations, the one with only one
iteration block together with the preprocessing circuit is an attractive
solution in terms of cost-efficiency.

Table 2: Implementation resources in terms of different hardware elements and
frequency operation.

Metrics | FIR | SSF 1 Block | SSF 2 Blocks | SSF 3 Blocks
---|---|---|---|---
Pp | nPp | Pp | nPp | Pp | nPp
LUT** | 1,447 | 31,378 | 33,132 | 59,927 | 63,966 | 87,160 | 94,870
FF** | 731 | 2,097 | 2,625 | 2,098 | 2,625 | 2,097 | 2,625
DSP** | 125 | 435 | 308 | 503 | 376 | 571 | 444
Freq (MHz) | 56.76 | 71.20 | 79.17 | 50.34 | 57.39 | 41.47 | 42.22
RMS (MeV) | 40 | 31 | 43 | 29 | 32 | 28 | 29

Pp = with preprocessing nPp = without preprocessing
**Number of: LUT = Look-up-tables, FF = Flip-Flops, DSP = Digital Signal
Processor

As it can be seen from previous analysis, the Separable Surrogate Functional
implementation requires $15\times$ more resources, on average, than FIR filter
structure. However, the proposed method presents a higher estimation
performance, which may justify its implementation, depending on the available
resources. It is worth to mention that the amount of resources for Separable
Surrogate Functional implementation is proportional to the square of the
window length. The value of 55 used here is just one example of direct
implementation. Ongoing analyses are investigating whether splitting the
processing window in several smaller windows may improve the implementation
requirements. For practical considerations, Table 3 presents the resources
used to fit the Separable Surrogate Functional method in three different FPGA
technologies, where the feasibility of implementation becomes evident.

Table 3: Resources spent in commercial FPGA devices.

Metrics | SSF 3 Blocks1 | XCZU7EV2 | XC7VX1140T3 | XCVU4404
---|---|---|---|---
nPp | low cost | middle cost | high cost
LUT | 94,870 | 41.17% (230,400) | 13,32% (712,000) | 3,74% (2,532,960)
FF | 2,625 | 0.56% (460,800) | 0,18% (1,424,000) | 0,05% (5,065,920)
DSP | 444 | 25.69% (1,728) | 13,21% (3,360) | 15,41% (2,880)

1SSF 3 Blocks without preprocessing
2XCZU7EV-FFVC1156-2-E 3XC7VX1140T-1FLG1930C-ND 4XCVU440-2FLGA2892E

##  6 Conclusion

This paper presented an iterative method for a matrix based deconvolution
process, which can be implemented for online operation. The proposed method is
based on established sparse theory, over-performing the previously proposed
FIR filter based deconvolution method. The method based on Iterative Shrinkage
has shown to be a good choice for implementation and a FPGA circuitry was
proposed. Compared with a simpler FIR filter design, despite the increase in
cost implementation, it is feasible when modern FPGA technology is considered.
The gain in estimation accuracy is evident, with improvement between 3% and
300% (depending on the occupancy level), when compared with FIR filter
implementation. Therefore, it may justify its usage in experiments with
intense pile-up levels.

## Acknowledgment

The authors would like to thank CNPq, CAPES, RENAFAE-MCTI, FAPERJ, FAPEMIG
(Brazil) and CERN (Switzerland) for their support to this work.

## References

  * [1] W. Barletta, M. Battaglia, M. Klute, M. Mangano, S. Prestemon, L. Rossi et al., _Future Hadron Colliders: from physics perspectives to technology R &D_, _Nucl. Instrum. Methods Phys. Res., A_ 764 (Jan, 2014) 352–368. 17 p.
  * [2] O. Brüning and L. Rossi, _The High Luminosity Large Hadron Collider_.  WORLD SCIENTIFIC, Geneve, 2015, 10.1142/9581.
  * [3] L. Evans and P. Bryant, _Lhc machine_ , _Journal of Instrumentation_ 3 (2008) S08001.
  * [4] K. Anthony, _Celebrating the first of a kind_ , _CERN Document Server_ (sep, 2014) .
  * [5] C. Pralavorio, _The new annual report is available_ , _CERN Document Server_ (jun, 2017) .
  * [6] M. Cepeda, S. Gori, P. Ilten and Collaboration, _Report from Working Group 2: Higgs Physics at the HL-LHC and HE-LHC_ , _CERN Yellow Rep. Monogr._ 7 (Dec, 2018) 221–584. 364 p.
  * [7] ATLAS Collaboration collaboration, W. Walkowiak, _ATLAS Plans for the High-Luminosity LHC_ , Tech. Rep. ATL-PHYS-PROC-2018-048, CERN, Geneva, Jun, 2018.  10.22323/1.326.0055.
  * [8] R. Wigmans, _Calorimetry_.  Oxford University Press, 2nd edition ed., Jan, 2018, 10.1093/oso/9780198786351.001.0001.
  * [9] M. Jeitler, _Trigger systems of LHC experiments_ , _Journal of Instrumentation_ 12 (may, 2017) C05012–C05012.
  * [10] ATLAS collaboration, G. Aad et al., _The atlas experiment at the cern large hadron collider_ , _Journal of Instrumentation_ 3 (2008) S08003.
  * [11] CMS collaboration, S. Chatrchyan et al., _The cms experiment at the cern lhc_ , _Journal of Instrumentation_ 3 (2008) S08004.
  * [12] E. Fullana, J. Castelo, V. Castillo, C. Cuenca, A. Ferrer, E. Higon-Rodriguez et al., _Digital signal reconstruction in the atlas hadronic tile calorimeter_ , _Nuclear Science, IEEE Transactions on_ 53 (09, 2006) 2139 – 2143.
  * [13] V. Radeka and S. Rescia, _Speed and noise limits in ionization chamber calorimeters_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 265 (Mar, 1988) 228–242.
  * [14] G. Bertuccio, E. Gatti, M. Sampietro, P. Rehak and S. Rescia, _Sampling and optimum data processing of detector signals_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 322 (Nov, 1992) 271–279.
  * [15] W. Cleland and E. Stern, _Signal processing considerations for liquid ionization calorimeters in a high rate environment_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 338 (Jan, 1994) 467–497.
  * [16] L. M. de Andrade Filho, B. S. Peralva, J. M. de Seixas and A. S. Cerqueira, _Calorimeter response deconvolution for energy estimation in high-luminosity conditions_ , _IEEE Transactions on Nuclear Science_ 62 (Dec, 2015) 3265–3273.
  * [17] S. Haykin, _Adaptative Filter Theory_.  Prentice Hall, NJ, USA, 5th edition ed., 2013.
  * [18] D. P. Barbosa, L. M. d. A. Filho, B. S. Peralva, A. S. Cerqueira and J. de Seixas, _Sparse representation for signal reconstruction in calorimeters operating in high luminosity_ , _IEEE Transactions on Nuclear Science_ 64 (Jul, 2017) 1942–1949.
  * [19] M. Elad, _Sparse and Redundant Representations_.  Springer New York, New York, 2010, 10.1007/978-1-4419-7011-4.
  * [20] L. R. M. Silva, L. M. de Andrade Filho and C. A. Duque, _Sparse representation algorithm applied to power systems signal compression_ , _International Transactions on Electrical Energy Systems_ 29 (Aug, 2018) e2693.
  * [21] J.-L. Starck, F. Murtagh and J. M. Fadili, _Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity_.  Cambridge University Press, New York, 2010.
  * [22] S. K. Mitra, _Digital Signal Processing : A Computer-Based Approach_.  McGraw-Hill Higher Education, New York, 4th edition ed., 2010.
  * [23] J. Duarte, L. de Andrade Filho, E. de Simas Filho, P. Farias and J. de Seixas, _Online energy reconstruction for calorimeters under high pile-up conditions using deconvolutional techniques_ , _Journal of Instrumentation_ 14 (dec, 2019) P12017–P12017.
  * [24] S. H. Lim, J. Won Choi and B. Shim, _Greedy sparse channel estimation for millimeter wave communications_ , in _TENCON 2018 - 2018 IEEE Region 10 Conference_ , pp. 1628–1632, Oct, 2018.  DOI.
  * [25] A. Zaki, A. Venkitaraman, S. Chatterjee and L. K. Rasmussen, _Greedy sparse learning over network_ , _IEEE Transactions on Signal and Information Processing over Networks_ 4 (Sep, 2018) 424–435.
  * [26] M. Niu, W. Li, Z. Liu, Y. Zhang and J. Yang, _An accelerated iterative shrinkage-thresholding algorithm for real-beam scanning radar super-resolution_ , in _2019 IEEE Radar Conference (RadarConf)_ , pp. 1–4, Apr, 2019.  DOI.
  * [27] Y. Ma and T. Zheng, _Stabilized sparse online learning for sparse data_ , 2017\.
  * [28] P. Zhao, D. Wang, P. Wu and S. C. H. Hoi, _A unified framework for sparse online learning_ , _ACM Trans. Knowl. Discov. Data_ 14 (Aug., 2020) .
  * [29] N. Botros, _Hdl with Digital Design_.  Mercury Learning & Information, 2015.
  * [30] U. Meyer-Baese, _Digital Signal Processing with Field Programmable Gate Arrays_.  Springer Berlin Heidelberg, New York, 3 ed., 2007, 10.1007/978-3-540-72613-5.
  * [31] G. Arce, _Nonlinear signal processing : a statistical approach_.  Wiley-Interscience, Hoboken, N.J, 2005.
  * [32] E. Candes and T. Tao, _The dantzig selector: Statistical estimation when p is much larger than n_ , _Ann. Statist._ 35 (12, 2007) 2313–2351.
  * [33] D. G. Luenberger and Y. Ye, _Linear and Nonlinear Programming_.  Springer US, New York, 2008.
  * [34] S. G. Mallat and Zhifeng Zhang, _Matching pursuits with time-frequency dictionaries_ , _IEEE Transactions on Signal Processing_ 41 (Dec, 1993) 3397–3415.
  * [35] D. R. Hadley, _Digital filtering performance in the ATLAS level-1 calorimeter trigger_ , in _Real Time Conference (RT), 2010 17th IEEE-NPSS_ , pp. 1–6, IEEE, 2010.
  * [36] S. M. Kay, _Fundamentals of Statistical Signal Processing, Volume I_.  Prentice Hall, NJ, USA, 1993.
  * [37] S. Chen and M. Yang, _An adaptive fast iterative shrinkage threshold algorithm_ , in _2017 29th Chinese Control And Decision Conference (CCDC)_ , pp. 2190–2194, May, 2017.  DOI.
  * [38] M. Takac, _Randomized Coordinate Descent Methods for Big Data Optimization_.  University of Edinburgh, Edinburgh, 2014.
  * [39] H. Wu and S. Wang, _Adaptive sparsity matching pursuit algorithm for sparse reconstruction_ , _IEEE Signal Processing Letters_ 19 (2012) 471–474.
  * [40] T. T. Cai and L. Wang, _Orthogonal matching pursuit for sparse signal recovery with noise_ , _IEEE Transactions on Information Theory_ 57 (2011) 4680–4688.
  * [41] B. Dumitrescu, _Sparse total least squares: Analysis and greedy algorithms_ , _Linear Algebra and its Applications_ 438 (Mar, 2013) 2661–2674.
  * [42] C. T. Kelley, _Solving Nonlinear Equations with Newtonś Method_.  Society for Industrial and Applied Mathematics, Jan, 2003, 10.1137/1.9780898718898.
  * [43] I. Daubechies, M. Defrise and C. De Mol, _An iterative thresholding algorithm for linear inverse problems with a sparsity constraint_ , _Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences_ 57 (2004) 1413–1457.
  * [44] S. Yun, _A Coordinate Gradient Descent Method for Structured Nonsmooth Optimization: Theory and Applications_.  VDM Verlag Dr. Müller, Dec, 2010.
  * [45] S. Roman, _Advanced linear algebra_.  Springer, New York, 2008.
  * [46] I. Hrivnacova, D. Adamova, V. Berejnoi, R. Brun, F. Carminati, A. Fasso et al., _The Virtual Monte Carlo_ , _CERN Document Server_ (May, 2003) 8 p.
  * [47] J. Chapman, _ATLAS simulation computing performance and pile-up simulation in ATLAS_ , _LPCC_ (2011) .
  * [48] CMS Collaboration collaboration, S. Banerjee, _CMS Simulation Software_ , Tech. Rep. CMS-CR-2012-146, CERN, Geneva, Jun, 2012.  10.1088/1742-6596/396/2/022003.
  * [49] Tevatron Beam-Beam Team collaboration, V. Shiltsev, _Overview of beam-beam effects in the Tevatron_ , in _Proceedings, ICFA Mini-Workshop on Beam-Beam Effects in Hadron Colliders (BB2013): CERN, Geneva, Switzerland, March 18-22 2013_ , pp. 11–18, 2014.  1410.3326.  DOI.
  * [50] W. Herr, _Dynamic behaviour of nominal and pacman bunches for different lhc crossing schemes_ , Tech. Rep. LHC-Project-Report-856. CERN-LHC-Project-Report-856, CERN, 01, 2005.
  * [51] J. Fox, H. Bartosik, E. Bjørsvik, S. De Santis, J. Dusatko, W. Höfle et al., _Control of Intra-Bunch Vertical Motion in the SPS with GHz Bandwidth Feedback_ , _J. Phys. : Conf. Ser._ 1067 (2018) WEPAL079. 6 p.
  * [52] M. Hostettler, K. Fuchsberger, G. Papotti, Y. Papaphilippou and T. Pieloni, _Luminosity Scans for Beam Diagnostics_ , _Phys. Rev. Accel. Beams_ 21 (Apr, 2018) 102801. 10 p.
  * [53] B. S. M. Peralva, L. M. A. Filho, A. S. Cerqueira and J. M. Seixas, _The TileCal Energy Reconstruction for Collision Data Using the Matched Filter_ , in _2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (2013 NSS/MIC)_ , pp. 1–6, 2013.  DOI.
  * [54] Y. Sugiyama, J. Xu, M. Tecchio, N. Whallon, D. McFarland, J. Ma et al., _The data acquisition system for the koto experiment_ , _IEEE Transactions on Nuclear Science_ 62 (Jun, 2015) 1115–1121.
  * [55] J. Keil, _Operational experience with hera_ , in _2007 IEEE Particle Accelerator Conference (PAC)_ , pp. 1932–1934, Jun, 2007.  DOI.
  * [56] R. Duda, P. Hart and D. G.Stork, _Pattern Classification_ , vol. xx, p. 654.  John Wiley and Sons, inc, 01, 2001.
  * [57] D. P. Bertsekas, _Nonlinear Programming_.  Athena Scientific, Sep, 1999.
  * [58] W. E. Milne, _Numerical Calculus (Princeton Legacy Library)_.  Princeton University Press, Apr, 2016.
  * [59] T. Fawcett, _An introduction to ROC analysis_ , _Pattern Recognition Letters_ 27 (Jun, 2006) 861–874.
  * [60] A. Böttcher and S. M. Grudsky, _Toeplitz Matrices, Asymptotic Linear Algebra and Functional Analysis_.  Hindustan Book Agency, 2000, 10.1007/978-93-86279-04-0.
  * [61] W. T. Padgett and D. V. Anderson, _Fixed-point signal processing_ , _Synthesis Lectures on Signal Processing_ 4 (Jan, 2009) 1–133.
  * [62] Y. Oda and T. Oyama, _Fast calculation using parallel processing and pipeline processing in power system analysis_ , _Electrical Engineering in Japan - ELEC ENG JPN_ 116 (Jan, 1996) 85–96.
  * [63] S. Churiwala, _Designing with Xilinx FPGAs: Using Vivado_.  Springer Publishing Company, Incorporated, 1st ed., 2016.
  * [64] S. Monk, _Programming FPGAs : getting started with Verilog_.  McGraw Hill Education, New York, NY, 2017.
  * [65] Xilinx, _Zynq ultrascale $+{}$ device packaging and pinouts_, 2020.