# Waveform resampling with LMN method

L. Gerlach  W. Gu  N. Nayak  X. Qian  B. Viren11footnotetext: Corresponding
author.

###### Abstract

Resampling is a common technique applied in digital signal processing. Based
on the Fast Fourier Transformation (FFT), we apply an optimization called here
the LMN method to achieve fast and robust re-sampling. In addition to
performance comparisons with some other popular methods, we illustrate the
effectiveness of this LMN method in a particle physics experiment: re-sampling
of waveforms from Liquid Argon Time Projection Chambers.

##  1 Introduction

Resampling is the process of transforming a signal in the form of samples
collected at discrete points in some domain to a new set of samples collected
at novel points in that domain. The quality of a resampling is typically
judged by how accurately these novel samples represent the original,
underlying signal. One central category of resampling methods is
interpolation. It asserts some model of the underlying signal and is
normalized to preserve the local signal amplitude. Other normalization
policies may be asserted. For example, an application may require the
resampling to leave the integral of the signal invariant.

The computational expense of a particular software implementation of a
resampling method is another important consideration. Various software
packages that are used in the analyses of data from particle physics
experiments provide resampling method implementations including Boost [1],
Eigen [2], SciPy [3] and ROOT [4].

Methods can be placed in two broad categories. The first category includes
those methods that operate in the sampled domain (e.g. interval or time
domain) and include Cardinal B-spline interpolation [5], Barycentric Rational
interpolation [6], Hermite interpolation [7], and Catmull Rom spline
interpolation [8]. Interval domain methods can be fast but tend to not
preserve spectral features.

The second category includes methods that consider the signal spectrum or
otherwise operate in the frequency domain. These include Whittaker-Shannon
interpolation [9] and trigonometric interpolation [10, 11]. Typically, these
methods apply operations in the frequency domain that are efficient while the
discrete Fourier transform (DFT) as implemented with the Fast Fourier
Transform (FFT) algorithm [12]222Gauss developed an algorithm akin to FFT in
1805 for a similar purpose of interpolating orbital measurements [13]. provide
an efficient means to transform a signal between interval and frequency domain
representations. While the FFT is indeed fast, its complexity still grows as
$\mathcal{O}(N\log N)$ in the size $N$ of the input signal. Naive resampling
methods can require inflating the original signal size by a large factor. This
inflation subjects the frequency domain method to ruinously large $N$.

Scaling problems can be exacerbated when the signal is large and multi-
dimensional. Interval-domain resampling can be extended to this case. For
example the bilinear interpolation [14] is a commonly employed method to
interpolate images in the two-dimensional interval or pixel domain. As with
the one-dimensional case, these extensions to two-dimensional intervals tend
to distort spectral features which can be problematic in some applications.
One dimensional methods that resample in the frequency domain may also be
extended to multiple dimensions in a simple manner when those dimensions are
independent. They require application of a series of DFT along each dimension
to reach their frequency domain. An FFT applied to a two-dimensional signal of
size $N\times M$ then costs $\mathcal{O}(MN\log N+NM\log M)$. The inflation
required by naive resampling methods then become yet more costly.

In this paper, we describe a variant of trigonometric interpolation, the LMN
method, that implements an optimized resampling of discrete waveforms in the
frequency domain by avoiding the inflation that comes with naive algorithms
consisting of an integer upsample factor followed by an integer downsample
factor. The quality and computational expense of this method are compared to
others that are used in the field.

This work is motivated by a technical task required in the processing of data
from detectors used by the Deep Underground Neutrino Experiment (DUNE) [15], a
next-generation, long-baseline neutrino oscillation experiment. DUNE’s primary
physics objectives include determining the neutrino mass hierarchy, search for
leptonic CP violation, and search for physics beyond the standard model of
particle physics. DUNE’s primary detector technology is the liquid argon time
projection chamber [16, 17, 18] (LArTPC). To understand the signals from these
detectors, their operation is summarized. When a charged particle traverses
the liquid argon in a LArTPC detector, ionization electrons are produced along
the particle’s trajectory. Under the influence of an external electric field,
the ionization electrons drift at a constant velocity toward a stack of three
parallel “anode planes”. Each anode plane is comprised of an array of parallel
electrodes spanning the width of their plane. The ionized electrons induce a
negative electric current as they approach an electrode and a positive current
as they drift past. This induced current is amplified, filtered, digitized and
sampled and the resulting discrete waveform is recorded. Given the location of
each electrode, these sampled waveforms represent a tomographic measurement of
the distribution of drifting electrons in two dimensions: transversely across
the electrode array and longitudinally along the electron drift direction. No
electron position information is available in the dimension along the length
of the electrodes. The anode planes are oriented so that their arrays of
electrodes are at mutually unique angles. Given a collection of these two-
dimensional tomographic views from all anode planes, a three-dimensional image
of the ionized electrons can be constructed. From that image, it is possible
to reconstruct information about the neutrino-argon interaction. Finally,
neutrino properties can be deduced from the ensemble of many such
reconstructed interactions. Lengthy and complex software developed by DUNE and
the larger LArTPC community is applied to theses recorded waveforms to perform
this extraction of physics information.

DUNE’s far detector underground cavern provides for four independent
cryostats, each containing approximately 17.5 kt of liquid argon. Currently,
two far detector module designs are in the process of being implemented. Each
is based on a variant of common LArTPC design principles. The modules are
known by the direction of their electron drift: the first far-detector module
FD1-HD [19] employs a horizontal drift and the second FD2-VD [20] drifts in
the vertical direction. They also utilize different anode plane and readout
technologies. In particular, the FD2-VD design divides its liquid argon volume
into two regions and supplies anode plane electrodes as traces on large
printed circuit boards installed near the top and bottom of the cryostat.
Supported by the international DUNE collaboration, the design specifies
different types of readout electronics for top and bottom volumes. The top
electronics are based on the design developed for the dual-phase liquid argon
detector technology [21]. The bottom electronics design [22] is in common with
that used for the FD1-HD [19] module. It is chosen to minimize the overall
detector capacitance, which in turn increases the signal-to-noise ratio. This
choice is driven by the necessarily long cable paths from the bottom anode
planes to the cryostat roof. One consequence of employing different readout
electronics technologies in FD2-VD is that a variety of sampling periods are
employed. Specifically, the top electronics employs analog-digital converters
(ADC) with a sampling period of 500ns while the bottom ADCs sample at 512ns.
In addition, the FD1-HD will sample at 512ns. Prior LArTPC experiments
(MicroBooNE [23], ProtoDUNE-SP [24], SBND) for which most of the existing data
processing software was developed sample at 500ns.

The readout electronics apply low-pass filters to their input analog signals
in order to satisfy the Nyquist-Shannon [25] condition in the subsequent
digitization stage. Any sampling period that places the Nyquist frequency
comfortably above the signal bandwidth may be chosen so that the information
in the signal will be faithfully captured. While the exact choice of sampling
period is not critical, employing multiple sampling periods across the
detectors leads to a technical complication for the data-processing software.
It is desirable to resample the ensemble of waveforms across the variety of
detectors to produce waveforms with a common sampling period prior to applying
the large software suite. If the resampling assures the Nyquist frequency
remains above the signal bandwidth and it does not otherwise distort the
waveforms then this resampling will preserve all information relevant to the
underlying physics.

This paper is organized as the following. In section 2, the signal processing
procedure applied to LArTPC data is briefly reviewed as it comprises the first
stage of data processing and it places strong constraints on the quality of
any resampling method. The LMN resampling method is then described in section
3. In section 4, the performance of the LMN method in terms of accuracy and
computational cost is compared with some popular interpolation methods. The
application of the LMN method to realistic LArTPC data is shown in section 5
using simulated data. Finally, a summary and discussions can be found in
section 6.

##  2 LArTPC signal processing

When charged particles traverse the liquid argon medium, electrons are ionized
along the particle trajectory and drift toward the anode planes under the
influence of the external electric field. While drifting, the distribution of
ionization electrons undergoes diffusion and absorption effects [26]. The TPC
signal is then formed with i) the distribution of ionized electrons arriving
at the anode plane, ii) the field response describing the resulting induced
current on the anode electrodes [27] (see Ref. [28] for a recent example of
calculating field response in three dimensions), iii) the electronics
response, which amplifies, shapes and converts current to a voltage signal
that is band limited, and iv) the ADC that samples the continuous voltage
signal to produce discrete-time waveforms.

The goal of the signal processing procedure is to effectively invert this
signal formation process. It reconstructs a measure of the distribution of
ionization electrons arriving at the anode electrodes given the recorded
waveforms. The kernel of the procedure deconvolves a model of the overall
detector response, which includes field and electronics responses, from the
waveforms. A technique to perform this deconvolution as a per-waveform filter
that is applied in the frequency domain via discrete Fourier transform (DFT)
was introduced to the LArTPC field in Ref. [29]. It relied on a one-
dimensional (1D) model that described how the drifting electrons induce
current in an electrode as a function of time. This was later improved by a
new deconvolution that extending the model two dimensions (2D) by adding
variation as a function of distance transverse to an electrode. This new 2D
method was first applied in the MicroBooNE experiment [30, 31] where it
demonstrated improvements compared to the 1D method. The 2D method has since
been utilized in the ProtoDUNE detectors, (DUNE prototype detectors), where it
also achieved good performance [32].

For simplicity and to share the notation of section 3, we describe the 1D
technique that performs a deconvolution over the time dimension. The 2D
deconvolution is conceptually a straight-forward extension. The technique
extracts a finite, sampled sequence $\hat{s}[n]$ that estimates the number of
electrons drifting past a given electrode $s[n]\equiv
s(t_{n})\equiv\int_{t=t_{n}}^{t_{n}+\Delta t}q(t)dt,\ n\in[0,...,N-1]$ over a
sample time period $[t_{n},t_{n}+\Delta t]$ to produce a sequence of $N$
samples. This signal $s[n]$ is not directly observable as it represents the
true distribution of ionization electrons. Instead, a detector provides a
measure $m[n]$ of the signal after the signal has been altered by a detector
response modeled in general as $r[n^{\prime},n]$. In addition, noise is
inescapable and is modeled as $\eta[n]$. In terms of these sampled sequences,
the measure is modeled as

$m[n]=\sum_{n^{\prime}}r[n^{\prime},n]\cdot s[n^{\prime}]+\eta[n].$ (2.1)

The first term is a convolution of the detector response and the signal. For
time-invariant detectors, the response may be replaced with
$r[n^{\prime},n]\to r[n^{\prime}-n]$. This form allows the discrete Fourier
transform (DFT, defined below in eq. 3.1) to be applied to rewrite eq. 2.1 as
$M[k]=R[k]S[k]+\mathcal{N}[k]$ were $k$ enumerates frequency domain samples.
Neglecting noise for the moment allows for an estimate of the signal to be
formed that may be expressed in the frequency domain as,

$S[k]\approx\hat{S}[k]=\frac{M[k]}{R[k]}.$ (2.2)

Naively, the inverse DFT (defined below in eq. 3.4) may be applied to recover
a time domain estimate $\hat{s}[i]$. However, several aspects of real
detectors coincide to make that final step impractical. First, the low-pass
filter that is applied by the detector electronics to the analog input is
designed to vanish at high frequency in order to satisfy the Nyquist-Shannon
condition for the subsequent sampling. Next, the electrodes comprising the
initial anode planes do not collect drifting electrons. This means their net
induced current is zero and their response is bipolar in time resulting in
$R[k\to 0]\to 0$. Finally, real-world noise can not be neglected. The noise to
which LArTPC detectors are subjected fills a broad-band spectrum
$\mathcal{N}[k]$ that can extend up to and beyond the Nyquist frequency while
some sources will introduce noise in the circuit downstream of the band-
limiting low-pass filters. Thus, the presence of noise in the measure $M[k]$
and a response $R[k]$ that is necessarily small or vanishing leads to eq. 2.2
producing a naive signal estimate that is artificially amplified at high
frequency and, for the initial anode planes, will diverge at low frequency.

To combat this practical problem, digital filtering is applied to the measure
in both the frequency and time domains. The time domain filter333The 2D
technique also filters in the transverse or “channel” domain. is introduced by
rewriting eq. 2.2 as

$S[k]\approx\hat{S}[k]=\frac{M[k]}{R[k]}F[k].$ (2.3)

From this result, an estimate $\hat{s}[n]$ in the time domain may be
constructed by applying the inverse DFT. In fact two such estimates are found,
each with their own filter. The first is constructed with a “Wiener-inspired”
filter that maximizes signal-to-noise. The second uses a (truncated) Gaussian-
shaped filter to give an unbiased measure of the signal.

The time-domain digital filtering consists of identifying “signal regions-of-
interest” (ROIs). Each ROI selects a fragment of the signal estimate after the
Wiener-inspired filtered has been applied and that is expected to have a high
signal-to-noise ratio. The algorithms that identify the signal-ROIs are
complex heuristics created by a development process that considered data from
real LArTPC detectors such as MicroBooNE [30, 33]. These ROIs are then applied
to the Gaussian-filtered estimate and a new signal baseline is calculated for
each region based on values at either ends of the ROI. This process
effectively combats the problematic combination of noise with the small or
vanishing portions of the response spectrum.

##  3 LMN resampling method

We review established resampling methods on which the LMN method is based. We
then describe the LMN method, which is essentially identical to the basic
method but adds a constructed, special-case optimization. Finally, a
discussion of normalization, information loss and mitigation of artifacts is
given.

###  3.1 Resampling basics

A continuous signal $g(t)$ may be sampled at times $t_{n}=n\cdot
T_{s},n\in[0,...,N_{s}-1]$ with period $T_{s}$ to produce a time-domain
sequence $\mathbf{g}_{s}\equiv\\{g_{s}[n]\\},\ g_{s}[n]\equiv g(t_{n})$. From
this sequence, another sequence
$\mathbf{g}_{r}\equiv\\{g_{r}[n]\\},n\in[0,...,N_{r}-1]$ may be constructed
that corresponds to a _resampling_ of $g(t)$ with a sampling period $T_{r}\neq
T_{s}$. This correspondence may be exact or aproximate depending on the
resampling method, the bandwidth of the signal and whether an overall
_upsampling_ ($T_{r}<T_{s}$) or _downsampling_ ($T_{r}>T_{s}$) is performed.

One general resampling method known as _fractional resampling_ [34] decomposes
the resampling procedure into two stages. First, the original sequence
$\mathbf{g}_{s}$ is upsampled by an integer factor $L$ to produce an
intermediate sequence $\mathbf{g}_{u}$ with a larger number of samples
$N_{u}=LN_{s}$. Second, this intermediate sequence is downsampled by an
integer factor $M$ to produce the final sequence $\mathbf{g}_{r}$ with a
number of samples $N_{r}=\frac{N_{u}}{M}$. The _simple fraction_
$\mathcal{R}\equiv\frac{L}{M}=\frac{N_{r}}{N_{s}}$ characterizes the
fractional resampling.

If the duration of the original signal is unchanged in the resampling then
$N_{s}T_{s}=N_{r}T_{r}$ and thus $\mathcal{R}=\frac{T_{s}}{T_{r}}$. However,
if $\frac{T_{s}}{T_{r}}$ is irrational then fractional resampling can not be
applied exactly. In such cases, a simple fraction
$\mathcal{R}\approx\frac{T_{s}}{T_{r}}$ may be chosen to approximate the
desired resampling ratio. The choice can be arbitrarily precise at the cost of
potentially large values of $L$ and $M$.

###  3.2 Discrete Fourier transform

The resampling used here is based on a subclass of these integer fractional
methods that operate on the Fourier representation of the sampled signal [35].
The discrete Fourier transform (DFT) is applied to the time-domain sequence
$\mathbf{g}$ of size $N$ to produce the frequency-domain sequence
$\mathbf{G}\equiv\\{G[k]\\}=DFT(\mathbf{g})$ with

$G[k]=\sum_{n=0}^{N-1}g[n]\cdot e^{-i2\pi\frac{nk}{N}},\ k\in\mathbb{I}.$
(3.1)

This sequence is infinite and periodic, $G[k]=G[k+N]$ over size $N$, and has
two points of Hermitian symmetry which are important for the resampling. The
first symmetry is around the “zero frequency” sample $G[0]$ so that
$G[k]=G^{*}[-k]$. The second is around the Nyquist frequency and is expressed
differently if the number of samples is even or odd as

$G[H-k]=\begin{cases}G^{*}[H+k+2],&N\in\mathcal{N}_{even}\\\
G^{*}[H+k+1],&N\in\mathcal{N}_{odd}\end{cases}$ (3.2)

Here, $H$ gives the size the “half spectrum” that consists of the sub-sequence
of frequency-domain samples that are not Hermitian-symmetric with themselves.
It is given by

$H=\begin{cases}\frac{N}{2}-1,&N\in\mathcal{N}_{even}\\\
\frac{N-1}{2},&N\in\mathcal{N}_{odd}\end{cases}$ (3.3)

For $N\in\mathbb{N}_{even}$, the sample $G[H+1]$ is real and called the
“Nyquist bin” as the Nyquist frequency lands at the center of its frequency
band. For $N\in\mathbb{N}_{odd}$, the Nyquist frequency is at the lower edge
of the sample $G[H+1]$.

There are two common ways to represent a finite subsequence of $\mathbf{G}$ of
size $N$. First is the “Nyquist-centered” range where the subset of samples
$\\{G[1],...,G[H]\\}$ are Hermitian-reflected about the Nyquist frequency
starting with $G[H+1]=G^{*}[H]$ for $N\in\mathbb{N}_{odd}$ or with
$G[H+2]=G^{*}[H]$ for $N\in\mathbb{N}_{even}$. The second is the “zero-
centered” range where the same subsequence of size $H$ is instead Hermitian-
reflected about $G[0]$. Influenced by this second choice, the reflected
subsequence is often said to consist of the “negative frequency” samples.

The inverse DFT may be applied to $\mathbf{G}$ to recover the original time-
domain sequence $\mathbf{g}$ as

$g[n]=\frac{1}{N}\sum_{k=0}^{N-1}G[k]\cdot e^{i2\pi\frac{nk}{N}},\
n\in\mathbb{I}.$ (3.4)

The result of the inverse DFT is also an infinite, periodic sequence. Note,
the form of the normalization term $\frac{1}{N}$ shown here is the convention
followed by most software providing the DFT. Additional discussion of
normalization is in Sec. 3.6.

###  3.3 Resampling in the Fourier domain

An upsampling from $N_{s}$ to $N_{u}=LN_{s}$ can be accomplished in the
frequency-domain by inserting $(L-1)N$ samples at both ends of the zero-
centered range (or equivalently in the middle of the Nyquist-centered range)
of $\mathbf{G}_{s}$ to produce an intermediate spectrum sequence
$\mathbf{G}_{u}$ with samples giving,

$G_{u}[k]=\begin{cases}G_{s}[k],&k\in[-H_{s},...,0,...,H_{s}]\\\
0,&\text{otherwise}\end{cases}$ (3.5)

Here, $H_{s}$ represents the “half spectrum” size for the original sampling as
defined in Eq. (3.3) and it is assumed that $N_{s}\in\mathbb{N}_{odd}$. When
$N_{s}\in\mathbb{N}_{even}$, Eq. (3.5) applies with the addition that the
Nyquist bin at $G_{s}[H+1]$ requires special treatment which is described in
Sec. 3.5.

If $\frac{N_{u}}{M}\not\in\mathbb{N}$ then the final $M$-wise downsampling
must be performed in the time domain after applying the inverse DFT of (large)
size $N_{u}=LN_{s}$ to $\mathbf{G}_{u}$. This allows any remainder sequence of
interval domain samples of size $N_{u}\ \mathrm{mod}\ M$ to be ignored and
corresponds to an effective change in the duration of the resampled signal
compared to the original.

In the special case that $N_{u}/M\in\mathbb{N}$, the downsampling may be
performed in the Fourier domain by removing $\frac{N_{u}}{M}(M-1)$ samples
from the middle of the Nyquist-centered range to produce $\mathbf{G}_{r}$ of
size $N_{r}$ with terms (expressed zero-centered),

$G_{r}[k]=G_{u}[k],k\in[-H_{r},...,0,...,H_{r}].$ (3.6)

Here again, $H_{r}$ gives the “half spectrum” size as defined in Eq. (3.3) for
the resampled sequence. Eq. (3.6) is strictly for $N_{r}\in\mathbb{N}_{odd}$
and treatment of the Nyquist bin for the case $N_{r}\in\mathbb{N}_{even}$ is
put aside until Sec. 3.5.

Obviously, in this special case, the upsampling and downsampling stages may be
combined so that the net number of frequency-domain samples is changed by
$\Delta N=N_{r}-N_{s}$. When $\Delta N>0$ a net upsampling is performed by
inserting this number of zero-valued samples similar to Eq. (3.5). Conversely,
$\Delta N<0$ indicates a net downsampling is performed with (the absolute
value of) this number of samples removed as in Eq. (3.6). To finish, the
intermediate spectrum $\mathbf{G}_{r}$ of size $N_{r}$ is formed and the
inverse DFT is applied to produce a resampled time-domain sequence
$\mathbf{g}_{r}$. In this special case, the size of the forward DFT is $N_{s}$
and the size of the inverse DFT is $N_{r}$, both of which are typically many
factors smaller than the inflated $N_{u}$. It is seeking this special case
that leads to the LMN optimization.

###  3.4 Special-case optimization

To give it a name, creating this special case or simply exploiting it when it
is accidentally met is called here the LMN method. The special case
$\frac{NL}{M}\in\mathbb{N}$ can be met for any given $T_{s}$ and $T_{r}$ that
meet a rationality condition described below and when there is freedom to
choose a value of $N_{s}$ even if it may differ from the size of a given input
signal. It is possible to generate a series of values for $N_{s}$ that satisfy
$N_{s}L/M\in\mathbb{N}$. Constructing the special case begins by considering
that the resampling will shift the sampling frequency by an amount

$\Delta F=\frac{1}{T_{r}}-\frac{1}{T_{s}}=\frac{T_{s}-T_{r}}{T_{s}T_{r}}.$
(3.7)

Next, this shift must span an integer number $\delta N\in\mathbb{N}$ of
original frequency-domain samples and thus must cover a frequency range

$\delta F=\frac{\delta N}{N_{s}T_{s}}.$ (3.8)

Equating $\delta F=\Delta F$ results in the relation

$N_{s}=\frac{\delta N\cdot T_{r}}{T_{s}-T_{r}}\in\mathbb{N}.$ (3.9)

The smallest value of $|\delta N|$ that satisfies Eq. (3.9) results in a
minimal sequence size $N_{s}=N_{min}$ for which the special case applies. Any
multiple of $N_{min}$ will also satisfy the special case. If found, we may
select the multiple of $N_{min}$ that best approximates the size of a given
input signal.

The satisfactory value for $\delta N$ may be found as

$\delta N=\frac{T_{s}-T_{r}}{\gcd(T_{r},\ T_{s}-T_{r})}.$ (3.10)

Here, $\gcd(\cdot,\cdot)$ returns the greatest common divisor of its
arguments. As the arguments are in general real valued, this function may
apply the Euclidean algorithm to calculate its result. Computer floating-point
round-off error will require the algorithm to test for convergence against a
small, non-zero error parameter. Note, a lack of convergence within this error
indicates that the ratio of its arguments is not rational. When this LMN
rationality condition is not met, this form of fractional resampling can not
be applied exactly. This is similar to the rationality condition of the
integer fraction method on which LMN is based.

The choice of the special case size $N_{s}$ of course determines the size of
the forward DFT and (along with the chosen resampling ratio $\mathcal{R}$) the
size of the inverse DFT. The freedom to select a multiple of $\delta N$ may be
exploited for a secondary optimization. The speed of the DFT strongly depends
on the size of prime factors that make up the length of the sequence being
transformed [36]. A DFT size with large prime factors requires substantially
more processing than nearby larger sizes. Thus, all else being equal, an
application of the LMN method may find that an otherwise unnecessary extension
of the original interval domain sequence by an additional $\delta N$ may lead
to an overall faster resampling. The effects of this dependence on the sizes
of the prime factors of a sequence length can be seen in the performance
results presented in section 4.

###  3.5 Treatment of Nyquist bins

In general, when $N\in\mathbb{N}_{even}$ the sequence $\mathbf{G}$ possesses a
Nyquist bin sample. When resampling in the Fourier domain, this sample
requires special treatment in two cases. The first case [35] is when a net
upsampling from size $N_{s}\in\mathbb{N}_{even}$ is performed. Half of the
Nyquist bin sample value $G_{s,nyquist}=G_{s}[H_{s}+1]$ is transferred to each
of the two equivalent samples in the resampled spectrum,
$G_{r}[H_{s}+1]=\frac{1}{2}G_{s,nyquist}$ and
$G_{r}[-H_{s}-1]=\frac{1}{2}G_{s,nyquist}$. Note, neither of these two samples
in the upsampled spectrum $\mathbf{G}_{r}$ are themselves a Nyquist bin. The
second and opposite case is when a net downsampling to size
$N_{r}\in\mathbb{N}_{even}$ is performed. The Nyquist bin in the resampled
spectrum must be set to the absolute value of either of the two corresponding
samples from the original spectrum, for example as
$G_{r}[H_{r}+1]=|G_{s}[H_{r}+1]|$.

###  3.6 Normalization

A “round trip” of a sequence through the forward and inverse DFT requires an
overall normalization of $\frac{1}{N}$ shown in Eq. (3.4). We follow here the
dominant convention to apply this normalization fully as part of the inverse
DFT. Other conventions are possible as long as a total of $\frac{1}{N}$
normalization is applied in the round trip. For example the symmetric
normalization places $\frac{1}{\sqrt{N}}$ on both forward and inverse DFT. All
such conventions produce the same result as long as $N$ is unchanged
throughout the round trip.

However, resampling in the Fourier domain necessarily performs the forward DFT
with size $N_{s}$ and the inverse DFT with size $N_{r}\neq N_{s}$. Care is
needed to select a normalization convention that preserves a desired
invariant. The choice of invariant may depend on how an application interprets
the quantity measured by signal samples.

For each of three interpretations considered below, the nominal $\frac{1}{N}$
applied as part of the inverse DFT is taken for granted and an additional
normalization of the form $g_{r}[n]\to g^{\prime}_{r}[n]=Ag_{r}[n]$ is
provided. It is also assumed that only zero-value frequency samples are
inserted or removed. Later we consider the case where this assumption is not
held.

  * •

Interpolation: An interval-domain sample gives a measure of the instantaneous
value of the underlying continuous signal. This interpretation requires
resampling to preserve the local amplitude of the interval-domain samples and
thus act as an interpolation. The nominal DFT convention applies a
normalization of $\frac{1}{N_{r}}$ and so the signal amplitude will be reduced
in an upsampling and increased in a downsampling compared the a direct
resampling of the continuous signal. Thus, an additional normalization of
$A_{0}=\frac{N_{r}}{N_{s}}$ must be applied to the resampled sequence for
correct interpolation. Alternatively, applying this normalization is
equivalent to applying the basic DFT normalization fully as part of the
forward DFT instead of following the accepted convention to apply it fully as
part of the inverse DFT.

  * •

Integral:  An interval-domain sample represents the integral of the underlying
continuous signal over the sample period. This interpretation must preserve
the sum of the time-domain sequence samples through the resampling process as
expressed by

$\sum_{n=0}^{N_{s}-1}g_{s}[n]=A_{1}A_{0}\sum_{n=0}^{N_{r}-1}g_{r}[n].$ (3.11)

If we were to apply only the interpolation normalization $A_{0}$ it would be
equivalent to integrating each of $N_{r}$ points over the original sample
period $T_{s}$ while integration over the resampled period $T_{r}$ is proper
to preserve the equality between sampled and resampled sums. Thus an
additional normalization of $A_{1}=\frac{T_{r}}{T_{s}}=\frac{N_{s}}{N_{r}}$
must be applied in addition to the interpolation normalization. These two
normalization terms cancel $A_{1}A_{0}=1$ and thus the nominal DFT
normalization convention is directly appropriate for the integral
interpretation.

  * •

Energy:  An interval-domain sample gives a measure of a probability amplitude,
field or other quantity where the sum-of-squares must be conserved through the
resampling. This sum is conserved through the forward DFT as expressed by
Parseval’s theorem which equates an energy measure in the time domain $E_{t}$
to an energy measure in the frequency domain $E_{f}$. In terms of the original
sampling, this may be written as

$E_{t,s}\equiv\sum_{n=0}^{N_{s}-1}\left|g_{s}[n]\right|^{2}=\frac{1}{N_{s}}\sum_{k=0}^{N_{s}-1}\left|G_{s}[k]\right|^{2}\equiv
E_{f,s}.$ (3.12)

The energy after resampling may be written as

$E_{f,r}\equiv\frac{1}{N_{r}}\sum_{k=0}^{N_{r}-1}\left|G_{r}[k]\right|^{2}=\frac{1}{N_{r}}\sum_{k=0}^{N_{s}-1}\left|G_{s}[k]\right|^{2}=\frac{N_{s}}{N_{r}}E_{f,s}.$
(3.13)

The effective replacement $r\to s$ in the sum of right hand term is allowed as
we have assumed only zero-amplitude frequency samples have been added or
removed in the resampling. Thus, in order for Parseval energy to be conserved
through the resampling, an additional normalization
$A_{2}=\sqrt{\frac{N_{r}}{N_{s}}}$ must be multiplied to the resampled time-
domain sequence, $g_{r}[n]\to g^{\prime}_{r}[n]=A_{2}g_{r}[n]$.

In the above discussion, it is assumed that only zero-amplitude frequency-
domain samples are inserted or removed in the resampling. In the case of
upsampling, this may assured by construction. In the case of downsampling it
is only assured if the original signal is band-limited to be below the new,
smaller Nyquist frequency. When this band limit is not assured, a downsampling
effectively applies a perfect low-pass filter with a cut-off at the resampled
Nyquist frequency $\frac{1}{2T_{r}}$. The removal of the non-zero samples
represents actual loss of signal information.

###  3.7 Potential artifacts

The method just described, and indeed any based on a round-trip via DFT,
assumes the signals are periodic and band-limited. Real-world signals are
typically neither. Resampling such signals can induce artifacts.

The original sampling produces a sequence over a finite interval window. It is
typically inescapable that the end of that sequence will not “wrap around” to
form a smooth connection to the beginning of the sequence. Even in a special
case where an underlying signal is constructed to be periodic it typically
includes aperiodic noise.

In the frequency domain, another source of artifacts can arise. While it is
typical to band-limit signals with a low-pass filter prior to the original
sampling, real world low-pass filters merely attenuate and do not fully remove
power above their “cut-off” frequency. Another source of high-frequency
content is noise that is picked up by electronics circuit components that
exist after the low-pass filter. Insertion or removal of frequency samples is
effectively equivalent to the application of a sharp low-pass filter. This
results in some amount of ringing in the interval domain after the inverse DFT
is applied.

The LMN method provides ways to mitigate the endpoint discontinuities of non-
periodic signals. To achieve the special case $N_{s}$ size, the original
signal must either be truncated or extended in the interval domain. When
extended, there is freedom to choose the values applied when padding the
signal to size $N_{s}$. Simple schemes can be employed to allow the endpoints
of the signal to be smoothly, periodically connected. For example, padding can
be filled with an interval-domain linear interpolation between last and first
sampled values of the original signal. Depending on the degree of
discontinuity, the sharp “kinks” this can create may still produce ringing. If
significant, more sophisticated padding models can be applied to form a
smoother connection. If the application allows, these artifacts can also be
mitigated by exploiting the intermediate frequency-domain representation.
There, a low-pass filter spectrum can be multiplied. Both of these forms of
mitigation require minimal additional computational expense.

##  4 Performance comparison with selected methods

The performance of a resampling algorithm can be evaluated in different ways
and it will vary subject to the shape of the considered signal as well as the
resampling ratio. Therefore, it is important to evaluate various performance
metrics and to do so for a range of signal frequencies, noise levels, and
sampling ratios. In this section, the performance of the LMN method in a
waveform resampling task will be compared with selected interpolation methods
including: i) linear interpolation, ii) polyphase filtering, and iii) cardinal
cubic B-spline.

###  4.1 Considered signals

Figure 1: Examples of generated input signals. The following underlying true
functions were used from top to bottom: single freq, multi freq, gauss pulse.
All three were generated with a frequency parameter of $c_{f}=4$ and a noise
amplitude of $c_{n}=0.1$.

The measure introduced in eq. 2.1 is idealized by a signal generated from an
analytic function $-1\leq f\leq 1$ representing the convolution $s*r$ and a
white-noise term $c_{n}\eta$ formed by scaling samples drawn from the normal
distribution. This idealized measure is then sampled at period
$T_{s}\equiv\frac{1}{N_{s}}$ to produce an initial sampled sequence
$g[n]=f[n]+c_{n}\cdot\eta[n],\ n\in[0,...,N_{s}-1]$. The following analytic
functions were used to generate sampled sequences that will be input to
ressampling.

$\displaystyle f(t)$ $\displaystyle=\sin(c_{f}\cdot t)$ single freq (4.1)
$\displaystyle f(t)$ $\displaystyle=\frac{1}{2}(\cos(c_{f}\cdot
t)+\cos(1.4\cdot c_{f}\cdot t+\frac{\pi}{3}))$ multi freq (4.2) $\displaystyle
f(t)$ $\displaystyle=\exp(-t^{2})\cdot\cos(c_{f}\cdot t^{2})$ gauss pulse
(4.3)

Here, $c_{f}$ is a frequency parameter that allows to generate a similar
signal with a different frequency. Examples of such signals can be seen in
figure 1.

###  4.2 Performance metrics

Figure 2: Comparison of true underlying signal with resampled signal. For each
case, the signal was resampled from $N_{s}=500$ to $N_{r}=512$ points. The
assumed signal here is gauss pulse with a frequency parameter of $c_{f}=5$.

The initial sequences are resampled $g[n]\to g_{r}[n^{\prime}],\
n\in[0,...,N_{s}-1],\ n^{\prime}\in[0,...,N_{r}-1]$. The top of figure 2 shows
an example of a signal of the gauss pulse type (eq. 4.3) with one choice of
frequency and noise parameters. In the middle, the results of resampling this
signal from $N_{s}=500$ to $N_{r}=512$ points while conserving total signal
duration are shown for the considered resampling methods. This is
representative of upsampling a signal with sample period 512ns to one with
500ns. The bottom panel shows the difference between the resampled discrete
noisy signal (“Model”) and the “pure” noise-free signal function $f(t)$
evaluated at the resampled times (“Truth”).

###  4.3 Scaling with signal frequency

In order to evaluate the quality of a resampling across a wide spectrum of
signal frequencies, sequences were generated with the single freq function
(eq. 4.1) for different values of the frequency parameter, $c_{f}$. Each
sequence was then resampled from $N_{s}=500$ to $N_{r}=512$ points using the
considered resampling methods. To quantify how well the resampled signal
matches the underlying true signal, the mean-squared error (MSE) statistic of
eq. 4.4 is calculated.

$\text{MSE}=\frac{1}{N_{r}}\sum_{i=0}^{N_{r}-1}{(f[n]-g_{r}[n])^{2}}.$ (4.4)

figure 3 shows the MSE as a function of the signal frequency parameter $c_{f}$
for the case of no noise (top panel) and a chosen non-zero value of $c_{n}$
(bottom panel). For the noiseless case, it can be seen that the LMN method
results in mean squared errors very close to or identical to zero across all
considered values for $c_{f}$. The other considered resampling methods exhibit
higher MSE with increasing signal frequency parameter. In the case of a noisy
signal, the MSE for signals with a small frequency parameter is lowest for the
linear interpolation. This is not necessarily desired, as the noise level
should in principle remain untouched by resampling. This reduction of noise is
the consequence of an implicit spectral distortion.

Figure 3: Mean squared error as a function of the frequency of the simulated
signal. The single freq signal was used with a noise amplitude of $c_{n}=0$
and $c_{n}=0.3$ in the top and bottom, respectively. For each case, the signal
was upsampled from $N_{s}=500$ to $N_{r}=512$ points.

In order to understand the inaccuracy of some of the resampling methods at
high signal frequencies, frequency-domain spectra of original and resampled
sequences are show in in figure 4 for the case $c_{f}=200$. In the figure, the
original signal spectrum is nearly identical to the resampled spectra from the
LMN and polyphase filtering methods. Resampling with the other methods,
however, introduces artifacts below and above the signal frequency and show
reduced power at the sole frequency of the true signal.

Figure 4: Frequency-domain spectra of original and resampled signals. The
single freq underlying signal was used with a frequency parameter of
$c_{f}=200$ and a noise amplitude of $c_{n}=0$ and $c_{n}=1$ in the top and
bottom, respectively. For each case, the signal was upsampled from $N_{s}=500$
to $N_{r}=512$ points. The amplitude at the peak frequency is given in
parentheses. The original Nyquist frequency is marked with the dashed line.

###  4.4 Processing time scaling with sampling rate

Figure 5: Average processing (run) time vs. resampled size for different
methods. See text for explanation of the Naive Resampler.

Besides being lossless, a good resampling method should also be
computationally efficient. As a way to evaluate this efficiency, the various
resampling algorithms are deployed to resample the same noisy signal 100 times
from which average processing times are determined. The original signal with
$N_{s}=500$ data points is resampled to sizes between 501 and 565 points.
figure 5 shows the average processing time as a function of the resampled
size. In the figure, the Naive Resampler case gives an expectation for run
times when LMN optimization is not used but rather the intermediate signals
are inflated to large sizes in order to enable independent integer upsampling
followed by integer downsampling. All methods utilizing DFT exhibit large
variance in average processing times across the different resampled sizes,
$N_{r}$. Comments on this origin of this speed variance may be found in
subsection 3.4. Depending on the precise number of resampled data points,
either the LMN method or linear interpolation requires the least computation,
whereas the naive approach is typically two orders of magnitude slower.

##  5 LArTPC waveform resampling

The motivation for developing the LMN resampling method was to solve the
practical and technical challenges of processing large amounts of data with
differing sampling periods from the different DUNE LArTPC detectors. Thus it
is imperative to demonstrate the method on waveforms that are representative
of what the LArTPC detectors produce. This section illustrates the resampling
of waveforms with period 512ns to those with period 500ns. The waveforms are
produced by the current state-of-the-art detailed LArTPC signal and noise
simulation algorithms from the Wire-Cell Toolkit [30, 37].

This simulation models the physical processes involved in electron drift[26,
38] through liquid argon, induced current in anode electrodes via field and
electronics response, noise effects and ADC digitization as introduced in
section 2. The kernel of the simulation involves a convolution of the
distribution of ionization electrons after they have drifted through the
volume with the field and electronics responses. The field responses are
calculated with a 2D model of the anode electrodes, the various applied
electrical potentials and a model of electron mobility in liquid argon. An
example of the resulting field responses for one anode plane, calculated with
GARFIELD [39], is shown in figure 6. It is represented as a 2D array giving
the amount of induced electric current in an electrode of interest at
different electron drift path starting locations transverse to (“Pitch”) the
electrode and along (“Time”) the nominal drift path. The longitudinal
dimension is sampled with a period of 100ns. As part of the convolution kernel
operation, the simulation will apply an integer downsampling. With 100ns
responses, a downsampling factor of five achieves the target ADC sample period
of 500ns.

Figure 6: The induced electric current field responses for the first anode
plane of the DUNE horizontal drift detector design. The “pitch” axis samples
transverse to the nominal electron drift direction while the “time” axis
samples at 100ns steps along the detailed drift path. Color gives a
logarithmic scale for the induced current at a given step.

The 100ns sampling can not be downsampled by an integer factor to the target
ADC sample period of 512ns and so a new set of field responses are needed. As
the field calculations are computationally costly, we elect here to apply LMN
resampling to the results at the nominal 100ns to produce results sampled at
64ns. This allows the simulation to downsample by an integer factor of eight
and achieve the desired 512ns ADC sampling period. This resampling of the
field responses also provides an auxiliary test of LMN method.

The induced current in an anode plane electrode that is due to a nearby
electron drift path (e.g. the central row from figure 6) is shown in figure
7). Such a path contributes a large fraction of the total induced current
compared to more distant paths. The $I_{s}$ shows the original response while
$I_{r}$ shows the LMN-resampling of this sequence. Note, this resampling is
normalized following the interpolation interpretation described in subsection
3.6.

Figure 7: One induced electric current field response for an example drift
path near an electrode of interest from the first anode plane in the DUNE
horizontal drift detector design. The left panel shows time domain and the
right shows the magnitude of the frequency domain (half) spectrum. The
originally calculated response is in black and red shows the result of LMN-
resampling from 100ns to 64ns.

High frequency ringing can be seen in the resampled field response. This is
due to the fact that the original field response is not fully band limited at
100ns. This can be seen in the small step down at the original 5MHz Nyquist
frequency. After convolution with the slower electronics response, which is a
low-pass filter, in order to form the overall detector response this ringing
is naturally abated as illustrated in figure 8, The overall detector response
is in units of induced electric charge and not induced electric current. That
is, the integrated normalization introduced in subsection 3.6 is taken by
multiplying the current by the sample period.

Figure 8: One induced charge field response for an example drift path near an
electrode of interest from the first anode plane in the DUNE horizontal drift
detector design. The left panel shows time domain and the right shows the
magnitude of the frequency domain (half) spectrum. The originally calculated
induced current field response convolved with the electronics response is
shown in black. The red shows the induced current field response LMN-resampled
from 100ns to 64ns convolved with electronics response sampled at 64ns.

Both sets of responses, each at a given sampling, are first input to the
simulation in order to output signal-only, voltage-level waveforms. Such
measures are not observable in a real detector but by comparing the native
500ns, 512ns and the $\qty{512}{ns}\to\qty{500}{ns}$ LMN-resampling waveforms
we may observe deviations or artifacts that may otherwise be hidden with
realistic noise and ADC integer truncation. Examples from one representative
channel of each anode plane show very good alignment in the time domain as
seen in figure 9.

Figure 9: Examples of a simulated noise-free, signal-only voltage-level
waveforms from a select channel of each anode plane. Two are from otherwise
identical simulated events with the first sampled at 500ns (blue) and the
second at 512ns (orange). The third waveform (green) is the result of
$\qty{512}{ns}\to\qty{500}{ns}$ LMN resampling.

Waveforms measured by real detectors includes signal and noise added at the
voltage level and then scaled and truncated to integer by the ADC. These
waveforms for native 500ns and 512ns simulation as well as the
$\qty{512}{ns}\to\qty{500}{ns}$ LMN-resampling at ADC level are shown in
figure 10. The detailed noise components in the 512ns waveform are closely
followed in the $\qty{512}{ns}\to\qty{500}{ns}$ resampling. It can also be
seen that the main signal body is faithfully resampled, within the variation
due to differently generated noise between the native 500ns simulation and the
noise shared by the native 512ns simulation and the
$\qty{512}{ns}\to\qty{500}{ns}$ resampling.

Figure 10: The voltage-level signal waveforms of figure 9 after noise is added
and the ADC digitizer model is applied. The 500ns (blue) and 512ns (orange)
sampling cases are shown. The third waveform (green) is the result of
resampling the 512ns ADC-level waveform to 500ns.

The 2D signal processing described in section 2 is sensitive to spectral
distortions. As a final comparison, signal processing using the responses
corresponding to the 500ns sampling is applied to the native 500ns ADC-level
signal and noise simulation and the $\qty{512}{ns}\to\qty{500}{ns}$
resampling. The results from the same channels selected above are shown in
figure 11. The small variations are due to the noise components that differ
between the two waveforms.

Figure 11: The output of the signal processing applied to native 500ns ADC-
level simulation and the 512ns ADC-level simulation after resampling to 500ns.

##  6 Summary

In this paper we describe the LMN method which achieves a fast and robust
waveform resampling based on the Fast Fourier Transformation algorithm. The
principle of the LMN method is introduced. It exploits freedom to resize the
original signal to enable a special case optimization that in turn allows for
upsampling and downsampling by integer factors to be performed simultaneously
in the frequency domain. Intermediate signals are kept to modest sizes and
thus ruinous computational costs are avoided. Normalization for a number of
different signal interpretations is provided. The potential for a resampling
to produce artifacts is raised and mitigation procedures enabled by the method
are outlined. The performance in terms of quality and computational cost of
the LMN method is compared with some popular methods. The LMN method is
generally among the fastest and most accurate. By exercising the simulation
and signal processing procedures used for Liquid Argon Time Projection
Chambers the effectiveness of the LMN method is demonstrated in a particle
physics experiment.

## Acknowledgments

We would like to thank Laura Zambelli for her insightful comments and valuable
suggestions during the DUNE collaboration review of this paper.

## References

  * [1] Boost, “Boost C++ Libraries.” http://www.boost.org/, 2015.
  * [2] G. Guennebaud, B. Jacob et al., “Eigen v3.” https://eigen.tuxfamily.org/, 2010.
  * [3] P. Virtanen et al., _SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python_ , _Nature Methods_ 17 (2020) 261.
  * [4] R. Brun and F. Rademakers, _Root — an object oriented data analysis framework_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 389 (1997) 81.
  * [5] I. Schoenberg, _Cardinal Spline Interpolation_ , vol. 12 of _CBMS-NSF Regional Conference Series in Applied Mathematics_ , SIAM (1973).
  * [6] D. Kong and S. Xiang, _Fast linear barycentric rational interpolation for singular functions via scaled transformations_ , 2021.
  * [7] _Cardinal hermite interpolation_ , in _Cardinal Spline Interpolation_ , pp. 43–57 DOI.
  * [8] E. Catmull and R. Rom, _A class of local interpolating splines_ , _Computer Aided Geometric Design_ 1 (1974) 317.
  * [9] E. Whittaker, _On the interpolatory function in theory of numbers_ , _Proceedings of the Royal Society of Edinburgh_ 35 (1915) 181\.
  * [10] _Trigonometric interpolation_ , in _Theoretical Numerical Analysis_ , Springer New York (2009), DOI.
  * [11] _Trigonometric interpolation_ , in _Numerical Analysis_ , Springer New York (1998), DOI.
  * [12] J.W. Cooley and J.W. Tukey, _An algorithm for the machine calculation of complex fourier series_ , _Mathematics of Computation_ 93 (1965) 297.
  * [13] M.T. Heideman, D.H. Johnson and C.S. Burrus, _Gauss and the history of the fast fourier transform_ , _Archive for History of Exact Sciences_ 34 (1985) 265—277.
  * [14] E.J. Kirkland, _Bilinear interpolation_ , in _Advanced Computing in Electron Microscopy_ , (Boston, MA), pp. 261–263, Springer US (2010), DOI.
  * [15] DUNE collaboration, _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume I Introduction to DUNE_ , _JINST_ 15 (2020) T08008 [2002.02967].
  * [16] C. Rubbia, _The Liquid Argon Time Projection Chamber: A New Concept for Neutrino Detectors_ , 1977.
  * [17] H.H. Chen, P.E. Condon, B.C. Barish and F.J. Sciulli, _A Neutrino Detector Sensitive to Rare Processes: I. A Study of Neutrino Electron Reactions_ , 1976.
  * [18] W.J. Willis and V. Radeka, _Liquid Argon Ionization Chambers as Total Absorption Detectors_ , _Nucl. Instrum. Meth._ 120 (1974) 221.
  * [19] DUNE collaboration, _Deep Underground Neutrino Experiment (DUNE), Far Detector Technical Design Report, Volume IV: Far Detector Single-phase Technology_ , _JINST_ 15 (2020) T08010 [2002.03010].
  * [20] DUNE collaboration, _The DUNE Far Detector Vertical Drift Technology, Technical Design Report_ , 2312.03130.
  * [21] DUNE collaboration, _The DUNE Far Detector Interim Design Report, Volume 3: Dual-Phase Module_ , 1807.10340.
  * [22] V. Radeka et al., _Cold electronics for ’Giant’ Liquid Argon Time Projection Chambers_ , _J. Phys. Conf. Ser._ 308 (2011) 012021.
  * [23] MicroBooNE collaboration, _Design and Construction of the MicroBooNE Detector_ , _JINST_ 12 (2017) P02017 [1612.05824].
  * [24] DUNE collaboration, _Design, construction and operation of the ProtoDUNE-SP Liquid Argon TPC_ , _JINST_ 17 (2022) P01005 [2108.01902].
  * [25] C.E. Shannon, _Communication in the presence of noise_ , _Proceedings of the IRE_ 37 (1949) 10.
  * [26] Y. Li et al., _Measurement of Longitudinal Electron Diffusion in Liquid Argon_ , _Nucl. Instrum. Meth. A_ 816 (2016) 160 [1508.07059].
  * [27] S. Ramo, _Currents induced by electron motion_ , _Proceedings of the IRE_ 27 (1939) 584.
  * [28] S. Martynenko et al., _A hybrid 3D/2D field response calculation for liquid argon detectors with PCB based anode plane_ , _JINST_ 18 (2023) P04033 [2303.10224].
  * [29] B. Baller, _Liquid argon TPC signal formation, signal processing and reconstruction techniques_ , _JINST_ 12 (2017) P07010 [1703.04024].
  * [30] MicroBooNE collaboration, _Ionization electron signal processing in single phase LArTPCs. Part I. Algorithm Description and quantitative evaluation with MicroBooNE simulation_ , _JINST_ 13 (2018) P07006 [1802.08709].
  * [31] MicroBooNE collaboration, _Ionization electron signal processing in single phase LArTPCs. Part II. Data/simulation comparison and performance in MicroBooNE_ , _JINST_ 13 (2018) P07007 [1804.02583].
  * [32] DUNE collaboration, _First results on ProtoDUNE-SP liquid argon time projection chamber performance from a beam test at the CERN Neutrino Platform_ , _JINST_ 15 (2020) P12004 [2007.06722].
  * [33] H. Yu et al., _Augmented signal processing in Liquid Argon Time Projection Chambers with a deep neural network_ , _JINST_ 16 (2021) P01036 [2007.12743].
  * [34] R. Schafer and L. Rabiner, _A digital signal processing approach to interpolation_ , _Proceedings of the IEEE_ 61 (1973) 692.
  * [35] D. Fraser, _Interpolation by the fft revisited – an experimental investigation_ , _IEEE Transactions on Acoustics, Speech, and Signal Processing_ 37 (1989) 665.
  * [36] I.J. Good, _The interaction algorithm and practical fourier analysis_ , _Journal of the Royal Statistical Society. Series B (Methodological)_ 20 (1958) 361.
  * [37] B. Viren, “Wire-cell toolkit web site.” https://wirecell.bnl.gov/, 2024\.
  * [38] Y. Li, C. Bromberg, M. Diwan, S. Kettell, S. Martynenko, X. Qian et al., _Parameterization of electron attachment rate constants for impurities in lartpc detectors_ , _Journal of Instrumentation_ 17 (2022) T11007.
  * [39] R. Veenhof, _Garfield, a drift chamber simulation program_ , _Conf. Proc. C_ 9306149 (1993) 66.


# Scaling of Seismic Memory with Earthquake Size

Zeyu Zheng  Department of Environmental Sciences, Tokyo University of
Information Sciences, Chiba 265-8501, Japan  Kazuko Yamasaki  Department of
Environmental Sciences, Tokyo University of Information Sciences, Chiba
265-8501, Japan  Center for Polymer Studies and Department of Physics, Boston
University, Boston, MA 02215, USA  Joel Tenenbaum  Center for Polymer Studies
and Department of Physics, Boston University, Boston, MA 02215, USA  Boris
Podobnik  Center for Polymer Studies and Department of Physics, Boston
University, Boston, MA 02215, USA  Faculty of Civil Engineering, University of
Rijeka, Rijeka, Croatia  H. Eugene Stanley  Center for Polymer Studies and
Department of Physics, Boston University, Boston, MA 02215, USA

###### Abstract

It has been observed that the earthquake events possess short-term memory,
i.e. that events occurring in a particular location are dependent on the short
history of that location. We conduct an analysis to see whether real-time
earthquake data also possess long-term memory and, if so, whether such
autocorrelations depend on the size of earthquakes within close spatiotemporal
proximity. We analyze the seismic waveform database recorded by 64 stations in
Japan, including the 2011 “Great East Japan Earthquake”, one of the five most
powerful earthquakes ever recorded which resulted in a tsunami and devastating
nuclear accidents. We explore the question of seismic memory through use of
mean conditional intervals and detrended fluctuation analysis (DFA). We find
that the waveform sign series show long-range power-law anticorrelations while
the interval series show long-range power-law correlations. We find size-
dependence in earthquake auto-correlations—as earthquake size increases, both
of these correlation behaviors strengthen. We also find that the DFA scaling
exponent $\alpha$ has no dependence on earthquake hypocenter depth or
epicentral distance.

###### pacs:

PACS numbers:89.65.Gh, 89.20.-a, 02.50.Ey

##  I Introduction

Many complex physical systems exhibit complex dynamics in which subunits of
the system interact at widely varying scales of time and space MFS ; Bunde .
These complex interactions often generate very noisy output signals which
still exhibit scale-invariant structure. Such complex systems span areas
studied in physiology Ashk20 , finance Engle82 , and seismology Bak ; Corral04
; Corral05 ; Lippiello07 ; Lippiello08 ; Bottiglieri ; Kagan ; Lennartz10 ;
Livina05 ; Lennartz08 .

In seismology the study of seismic waves is both scientifically interesting
and of practical concern, particularly in such applied areas as engineering. A
better understanding of seismic waves is immediately applicable in the design
of structures for earthquake-prone areas Mayeda ; Hisada ; Hisada2 . It also
allows scientists to better understand the underlying mechanisms that drive
earthquakes Bensen ; Xu2 ; Okada ; Ide ; Shapiro ; Campillo . In seismology,
temporal and spatial clustering are considered important properties of seismic
occurrences and, together with the Omori law (dictating aftershock timing) and
the Gutenberg-Richter law (specifying the distribution of earthquake size),
comprise the main starting requirements to be fulfilled in any reasonable
seismic probabilistic model. Analyzing the timing of individual earthquakes,
Ref. Bak  introduces the scaling concept to statistical seismology. The
recurrence times are defined as the time intervals between consecutive events,
$\tau_{i}=t_{i}-t_{i-1}$. In the case of stationary seismicity, the
probability density $P(\tau)$ of the occurrence times was found to follow a
universal scaling law

$P(\tau)=Rf(R\tau)$ (1)

where $f$ is a scaling function and $R$ is the rate of seismic occurrence,
defined as the mean number of events with $M\geq M_{c}$ Corral04 . Reference
Corral05 ; Lippiello07  has demonstrated how the structure of seismic
occurrence in time and magnitude can be treated within the framework of
critical phenomena.

Recently, a few papers have analyzed the existence of correlations between
magnitudes of subsequent earthquakes Corral05 ; Lippiello07 . Analyzing
earthquakes with $\tau$ greater than 30 minutes, Ref. Corral05  reported
possible magnitude correlations in the Southern California catalog. Magnitude
correlations have often been interpreted as a spurious effect due to so called
short-term aftershock incompleteness (STAI) Kagan . This hypothesis assumes
that some aftershocks, especially small events, are not reported in the
experimental catalogs, which is in agreement with the standard approach that
assumes interdependence of earthquake magnitudes implying no memory in
earthquakes.

However, recent work has also challenged this interpretation. Reference
Lippiello08  reports the existence of magnitude clustering in which
earthquakes of a given magnitude are more likely to occur close in time and
space to other events of similar magnitude. They find that a subsequent
earthquake tends to have a magnitude similar to but smaller than the previous
earthquake. Reference Lippiello07  also reports the existence of magnitude
correlations and additionally demonstrates the structure of these correlations
and their relationship to $\Delta t$ and $\Delta r$, where the latter
represents the distance between subsequent epicenters. Reference Lennartz10
creates a model to explain these magnitude correlations. They note that the
Omori law and “background tectonic cycles” are responsible for clustering in
interoccurrence times. Additionally, Refs. Livina05  and Lennartz08  find that
the distribution of recurrence times strongly depends on the previous
recurrence time such that small and large recurrence times tend to cluster in
time. This dependence on the past is reflected in both the conditional mean
recurrence time and the conditional mean residual time until the next
earthquake.

Since it is our hypothesis that long-range autocorrelations exist in seismic
waves, we first note that long-range power-law autocorrelations are quite
common in a large number of natural phenomena ranging from weather Yamasaki3 ;
Gozolchiani ; BP05 , and physiological systems Lennartz ; Ashk20 ;
Kantelhardt2 ; Stanley3 ; Karasik , to financial markets Mantegna ; Yamasaki ;
Wang2 ; Wang3 ; Stanley ; BPPnas09 ; BPPnas10 .

In addition to analyzing the raw waveform, it is also common to analyze
related time series, such the time series generated by taking the sign or
magnitude of the waveform Ashk20 . Reference Ashk20  reports an empirical
approximate relation at small time scales for the scaling exponents calculated
for sign, magnitude, and the original time series, $\alpha_{\rm
sign}=1/2(\alpha_{\rm magnitude}+\alpha_{\rm original})$, in physiology. The
study of magnitude and sign time series is important in physiology because the
magnitude time series exhibits weaker autocorrelations and a scaling exponent
closer to the exponent of an uncorrelated series found when a subject is
unhealthy Ashk20 . Diagnostic power in physiology has been confirmed for sign
time series as well—the sign time series of heart failure subjects exhibit
scaling behavior similar to that observed in the original time series, but
significantly different that of healthy subjects Ashk20 . Understanding the
correlation properties of these three time series allows us to also understand
the underlying processes generating them.

Our investigation and discussion is organized as follows. First, we study the
autocorrelations of interval series by using the mean conditional technique.
Second, we employ detrended fluctuation analysis (DFA) CKP ; Hu ; Chen  and
find long-range power-law autocorrelations in the sign and interval time
series. For the interval time series we find a positive regression between the
DFA scaling exponent $\alpha$ and earthquake size (measured by the Richter
magnitude scale $M$ or seismic moment $M_{0}$), while for the sign time series
we find an inverted regression between $\alpha$ and earthquake magnitude. Thus
we report that the observed autocorrelation depends on earthquake size, both
in the sign and interval time series. We also find that the scaling exponent
$\alpha$ has no dependence on hypocenter depth or epicentral distance.

##  II Data

Seismic waves are unique in that they have non-stationarities of a much larger
order than those of any other known natural signal. Large earthquakes are
characterized by a maximum amplitude that is often $>100$ times larger than
the mean amplitude [see Fig. 1(a)]. This is a limitation that makes seismic
waves difficult to analyze using traditional analysis. Although we might want
to use detrended fluctuation analysis (DFA) CKP ; Hu ; Chen ; Chen2 ,
originally proposed to study the correlations in a time series in the presence
of non-stationarities commonly observed in natural phenomena, the level of
non-stationarity in earthquakes is so large that DFA is inappropriate
regardless of the order of the polynomial fit applied Chen . Thus, due to lack
of methods for highly non-stationary signals, we do not analyze correlations
in the series of magnitudes, but instead analyze the correlations in the sign
series [Fig. 1(c)] and interval series [Fig. 1(d)]. For our data, we use the
seismic waveform database from the National Research Institute for Earth
Science and Disaster Prevention (NIED) F-net (Full Range Seismograph Network
of Japan), which records continuous seismic waveform data $w_{t}$ by using
broadband sensors in 64 stations in Japan [see Fig. 1(a)]. In our study we
select 46 stations (ADM, AOG, ASI, HID, HJO, HRO, IGK, IMG, INN, IYG, IZH,
KGM, KMU, KNM, KNP, KNY, KSK, KSN, KSR, KYK, MMA, NKG, NOK, NOP, NRW, NSK,
OSW, SAG, SHR, SIB, TAS, TGA, TGW, TKO, TMC, TSA, TYM, TYS, UMJ, WTR, YAS,
YNG, YSI , YTY, YZK, ZMM), based on locations and integrity of data series.
Seismic signals are recorded in three directions: (1) U (up-down with up
positive), N (north-south with north positive), and E (east-west with east
positive) Okada . In this paper, we report results from the vertical dimension
only (U data), since the results for the horizontal data (N and E) data are
very similar. Sampling intervals have five recording frequencies: 80Hz, 20Hz,
1Hz, 0.1Hz, and 0.01Hz. We study earthquake coda wave data with 1Hz sampling
interval for the year 2003, together with selected earthquake coda wave data
from 11 March 2011. We note that, because of the interaction between
earthquakes, not all earthquakes can be employed in our analysis (see Appendix
A). The data from 11 March 2011 is selected because it contains the notable
2011 Tohoku earthquake (“Great East Japan Earthquake”) which resulted in the
tsunami that caused a number of nuclear accidents. We also add two large
earthquakes ($M=7.3$ and $M=7.6$) to our study, which also occurred the same
day as aftershocks.

We employ the following procedure to create our time series:

  * (i)

For each selected earthquake (see Appendix A) we create a new time series, the
normalized waveform denoted by $w_{t}$ out of the raw seismic acceleration
waveform data

$w_{norm}\equiv(w_{t}-\overline{w})/\sqrt{\overline{w_{t}^{2}}-\overline{w}^{2}}.$
(2)

  * (ii)

From the time series $w_{norm}$ we define a new sub-series $w_{t}^{\prime}$,
starting at time coordinate where maximum $w_{t}$ occurs and terminating at
the end of the normalized waveform $w_{t}^{\prime}$ (see inset in Fig1(a)).

  * (iii)

Let the time series $t_{i}$ denote the points in time when $w_{t}^{\prime}$
changes sign, with $t_{i}<t_{i+1}$. We define (see Fig 1(c)) the interval
series by

$\tau_{i}\equiv t_{i}-t_{i-1}.$ (3)

  * (iv)

The sign series (see Fig. 1(d)) is defined by

$s_{t}\equiv sgn(w_{t}^{\prime})$ (4)

Note that our definition of interval is different than that recently defined
in several papers, where the return intervals $\tau$ have studied between
consecutive fluctuations above a volatility threshold $q$ in different complex
systems. The probability density function (pdf) of return intervals
$P_{q}(\tau)$ scales with the mean return interval as

$P_{q}(\tau)=\overline{\tau}^{-1}f(\tau/\overline{\tau})$ (5)

where $f()$ is a stretched exponential Yamasaki ; Wang2 ; Wang3 . Since, on
average, there is one volatility above the threshold $q$ for every
$\overline{\tau}_{q}$ volatilities, then it holds that BPPnas09

$1/\overline{\tau}_{q}\approx\int_{q}^{\infty}P(|R|)d|R|=P(|R|>q)\sim
q^{-\alpha}.$ (6)

For the time intervals $\tau_{q}$ between events given by fluctuations $R$
where $R>q$ Ref. BPPnas09  derived that $\overline{\tau_{q}}$, the average of
$\tau_{q}$, obeys a scaling law,

$\overline{\tau_{q}}=q^{\alpha}$ (7)

where by $\alpha$ denotes our estimate of the tail exponent probability
density function, $P(|R|^{1+\alpha})$. Similarly, if $P(|R|)$ follows an
exponential function $P(|R|)\propto\exp(-\beta|R|)$, then employing Eq. (6) we
easily derive

$\overline{\tau}_{q}\propto\exp{(\beta q)}.$ (8)

Eq. (8) can be used as a new method for estimation of the exponential
parameter $\beta$.

##  III  Memory of interval time series

Returning to waveform data, we begin analyzing the series by studying the
conditional mean

$\langle\tau|\tau_{0}\rangle/\overline{\tau}$ (9)

which gives the mean value of $\tau$ (see Eq. (3)) immediately following a
given term $\tau_{0}$, normalized in units of $\overline{\tau}$. The
conditional mean gives evidence of whether seismic memory exists in the
intervals in the form of correlations or anticorrelations. For example, should
correlations exist, one would expect the mean interval to be shorter in the
window immediately following a small interval.

Indeed, Fig. 2 shows that the large intervals $\tau$ tend to follow large
initial $\tau_{0}$ and small $\tau$ follow small $\tau_{0}$ indicating the
existence of (positive) correlations in the interval time series. We also note
that the autocorrelations tend to be stronger for the subset associated with
larger earthquakes than for those associated with smaller earthquakes.

To expand on this we also extend our investigation to longer range effects. We
investigate the mean interval after a cluster of $n$ consecutive intervals
that are either entirely above the series mean or entirely below it. We denote
clusters that are entirely above the series mean with a “$+$” and clusters
below the series mean with a “$-$”. Fig. 3 shows the mean interval $\tau$ that
follows a $\tau_{0}(n)$ defined as a cluster size of $n$. We find that for
“$+$” clusters—shown by open symbols—the mean interval increases with the size
of the cluster $n$. This is the opposite of what we find for “$-$”
clusters—shown as closed symbols. The results indicate the existence of at
least short-term memory in the interval time series. Furthermore, we find that
the mean interval increases with the seismic magnitude. However, this
relationship breaks at the high end of the Richter magnitude scale $M>6.5$.

##  IV Detrended fluctuation analysis

Many physical, physiological, biological, and social systems are characterized
by complex interactions between a large number of individual components, which
manifest in scale-invariant correlations MFS ; Bunde ; Tak ; Kob82 . Since the
resulting observable at each moment is the product of a magnitude and a sign,
many recent investigations have focused on the study of correlations in
magnitude and sign time series Ashk20 ; Kantelhardt2 ; Kant2002 ; Hu ;
Plamen04itt ; Livina ; Engle82 . For example, the time series of changes
$\delta\tau_{i}$ of heartbeat intervals Ashk20 ; Kant2002 ; Kantelhardt2 ,
physical activity levels Hu , intratrading times in the stock market
Plamen04itt , and river flux values Livina  all exhibit power-law
anticorrelations, while their magnitudes $|\delta\tau_{i}|$ are positively
correlated. A common means of finding autocorrelations hidden within a noisy
non-stationary time series is detrended fluctuation analysis (DFA)CKP ; Hu ;
Chen . In the DFA method, the time series is partitioned into pieces of equal
size $n$. For each piece, the local trend is subtracted and the resulting
standard deviation over the entire series is obtained. In general, the
standard deviation $F(n)$ of the detrended fluctuations depends on $n$, with
smaller $n$ resulting in trends that more closely match the data. The
dependence of $F$ on $n$ can generally be represented as a power law such that

$F(n)\propto n^{\alpha},$ (10)

where $\alpha$ is the scaling exponent—sometimes referred to as the Hurst
exponent—to be obtained empirically. DFA therefore can conceptually be
understood as characterizing the motion of a random walker whose steps are
given by the time series. $F(n)$ gives the walker’s deviation from the local
trend as a function of the trend window. Because the root mean square
displacement of a walker with no correlations between his steps scales like
$\sqrt{(}n)$, we can expect a time series with no autocorrelations to yield an
$\alpha$ of 0.5. Similarly, long-range power-law correlations in the signal
(i.e. large terms follow large terms and small terms follow small terms)
manifest as $\alpha>0.5$. Power-law anticorrelations within a signal will
result in $\alpha<0.5$. Additionally, DFA can be related to the
autocorrelation as follows: if the autocorrelation function $C(L)$ can be
approximated by a power law with exponent $\gamma$ such that

$C(L)\propto L^{-\gamma},$ (11)

then $\gamma$ is related to $\alpha$ by CKP

$\alpha\approx 1-\gamma/2.$ (12)

Another reason we employ the DFA method is that it is appropriate for sign
time series Kantelhardt2 . Other techniques for the detection of correlations
in non-stationary time series are not appropriate for sign time series. Also,
because the sign and interval time series have affine relations, the analysis
of sign will be helpful in understanding the intervals. However, the DFA gives
biased estimates for the power-law exponent in analysis of anticorrelated
series Hu , and so in order to improve the accuracy of analysis, we integrate
the time series before we employ the standard DFA procedure.

For the 2011 Tohoku earthquake, also known as the “Great East Japan
Earthquake”, we present the fluctuation function $F(n)$ of the coda wave,
measured at KSN station, as typical examples of sign and intervals time series
(Fig. 4). By using DFA, we find, for most coda waves after earthquakes, that
the time series of the intervals are consistent with a power-law correlated
behavior $\alpha=0.69$, while the sign time series of Eq. (4) are consistent
with a power-law anti-correlated behavior ($\alpha=0.32$). The results
therefore indicate that for the interval series large increments are more
likely to be followed by large increments and small increments by small
increments. These results are in agreement with the results of the correlation
analysis reported in Section 3. In contrast, anticorrelations in the sign time
series indicate that positive increments are more likely to be followed by
negative increments and vice versa.

For the entire set of sign time series comprising our sample we calculate the
average DFA scaling exponent $\overline{\alpha}=0.34\pm 0.09$ indicating
anticorrelations, and for the interval time series we calculate the average
DFA scaling exponent $\overline{\alpha}=0.58\pm 0.08$ indicating correlations.
For the different stations measuring the 2011 Tohoku earthquake we find that
for the sign time series, $\overline{\alpha}=0.29\pm 0.05$ and for the
interval time series, $\overline{\alpha}=0.66\pm 0.07$.

##  V Relation between earthquake moments and scaling exponents of sign and
interval series

Because large earthquake events release such extraordinary amounts of energy,
it is reasonable to ask whether their occurrence influences local wave
dynamics. To this end, we study interval time series of coda waves from
earthquakes occurring in 2003, also including the particularly large events of
11 March 2011, when three events $M>7$ occurred in the same day. Fig. 5(a)
shows the DFA scaling exponent of the sign series versus seismic moment, where
seismic moment is a quantity used to measure the size of an earthquake. We
find a decreasing functional dependence between the DFA exponent of the sign
series and the seismic moment of the proximal earthquake with slope
$\gamma=-0.028\pm 0.002$, indicating that the DFA exponent decreases
approximately with seismic moment. Note that because most of the exponents are
$<0.5$, this indicates the presence of ever stronger anticorrelations in the
time series as earthquake magnitude increases. Note, however, that the data
break with this trend for very large earthquakes (Richter magnitude scale
$>6.6$ or seismic moment $>10^{19}$).

We also find similar results in the interval series, the difference being that
the anticorrelations become correlations. Fig. 5(b) shows that the DFA
interval exponent and seismic moment exhibit a positive functional dependence
with slope $\gamma=0.025\pm 0.002$ so that the DFA exponent increases with
increasing seismic moment. Because most of the exponents for the interval
series are $>0.5$, this indicates that the series show stronger correlations
for increasing seismic moment. Again, as with the sign series, we find a
deviation from this trend for very large earthquakes.

Having observed the influence of seismic moment on autocorrelations, we now
investigate whether other readily observable factors such hypocenter depth and
epicentral distance (the distance from the event to the recording station)
also contribute. Specifically, we would like to explore whether there is
evidence that such long-term memory is affected by the spreading process as
seismic waves disseminate outward from their epicenter to a recording station
or whether the memory observed is strictly due to the seismic activity. Fig. 6
shows that the DFA exponent for both interval and sign series are independent
of both hypocenter depth and epicentral distance. From these results we
speculate that the DFA exponent is mainly a result of the characteristics of
the hypocenter rather than the process by which the seismic waves are spread.

For moderately large earthquakes ($M_{0}=10^{14}\sim 10^{19}$), we approximate
the relation between the DFA scaling exponent and seismic moment through the
empirical formula

$\alpha\approx a~{}log_{10}(M_{0})+c$ (13)

where $a=-0.028$, $c=0.797$ for the sign time series and where $a=0.025$,
$c=0.174$ for the interval time series. Since

$M=(log(M_{0})-9.1)/1.5,$ (14)

we can also write

$\alpha\approx a(1.5M+9.1)+c=a^{\prime}~{}M+c^{\prime},$ (15)

where $a^{\prime}=-0.042$, $c^{\prime}=0.542$ for the sign series, and
$a^{\prime}=0.037$, $c^{\prime}=0.398$ for the interval series.

We note that similar size dependence in Hurst exponent was found in Ref.
Eisler  where Hurst exponents of financial time series increase
logarithmically with company size.

##  VI Summary

We analyze seismic coda waves during earthquakes, finding long-range power-law
autocorrelations in both the interval and sign time series. The sign series
generally display power-law anticorrelated behavior, with anticorrelations
becoming stronger with larger earthquake events, while the interval series
generally display power-law correlated behavior, with correlations also
becoming stronger with larger earthquake events. We also show that while the
DFA autocorrelation exponent is influenced by the size of the earthquake
seismic moment, it is unaffected by earthquake depth or epicentral distance.
Our findings are in contrast with a standard approach which assumes
independence in earthquake signals and thus have strong implications on the
ongoing debate about earthquake predictability SornettePNAS .

##  VII Acknowledgements

We thank S. Havlin for his constructive suggestions, and thank JSPS for grant
of ”Research project for a sustainable development of economic and social
structure dependent on the environment of the eastern coast of Asia” that made
it possible to complete this study. We also thank the National Science
Foundation and the Ministry of Science of Croatia for financial support.

##  VIII Appendix: The Selection of Earthquakes

In some regions it is common for multiple earthquakes to occur in short
succession. In many cases, because the interoccurrence times are so short, the
coda waves can be derived from more than one earthquake. This is especially
true for large earthquakes with many aftershocks Utu . In order to make sure
that the coda waves we study are the effects of only one earthquake, we need a
way of determining which earthquakes are independent. We use the following two
functions to determine the sphere of influence and duration of each earthquake
by using the Richter magnitude scale M Utu . We select only those earthquakes
that have no larger earthquake in their spatiotemporal sphere of influence,

$t\approx 10^{(M-4.71)/1.67}$ (16)

and

$R\approx 2\times 10^{(M+1)/2.7},$ (17)

where $t$ is the duration and $R$ is the sphere radius of influence. The two
functions are empirical formulas based on an analysis of earthquakes in Japan
Utu . The $10^{M+1.0}/2.7$ is an empirical formula that indicates the maximum
radius that a human can feel an earthquake, especially for the earthquakes in
Japan.

## References

  * (1) M. F. Shlesinger, Ann. NY Acad. Sci. 504, 214 (1987).
  * (2) A. Bunde and S. Havlin (Editors), Fractals in Science (Springer, Berlin, 1994).
  * (3) Y. Ashkenazy, P. C. Ivanov, S. Havlin, C. K. Peng, A. L. Goldberger, and H. E. Stanley, Phys. Rev. Lett. 86, 1900 (2001).
  * (4) R. F. Engle, Econometrica 50, 987 (1982).
  * (5) P. Bak, K. Christensen, L. Danon, and T. Scanlon, Phys. Rev. Lett. 88, 178501 (2002).
  * (6) A. Corral, Phys. Rev. Lett. 92, 108501 (2004).
  * (7) A. Corral, Phys. Rev. Lett. 95, 028501 (2005).
  * (8) E. Lippiello, C. Godano, and L. de Arcangelis, Phys. Rev. Lett. 98, 098501 (2007).
  * (9) Y. Y. Kagan, Bull. Seismol. Soc. Am. 94, 1207 98, 098501 (2004).
  * (10) E. Lippiello, L. de Arcangelis, and C. Godano, Phys. Rev. Lett. 100, 038501 (2008).
  * (11) M. Bottiglieri, L. de Arcangelis, C. Godano, and E. Lippiello, Phys. Rev. Lett. 104, 158501 (2010).
  * (12) S. Lennartz, A. Bunde, and D. L. Turcotte, Geophys. J. Int. 184 1214 (2010).
  * (13) V. N. Livina, S. Havlin, and A. Bunde, Phys. Rev. Lett. 95, 208501 (2005).
  * (14) S. Lennartz, V. N. Livia, A. Bunde, and S. Havlin, Europhys, Lett. 81, 69001 (2008).
  * (15) K. Mayeda and W. Walter, Journal of Geophysical Research 101, 11195 (1996).
  * (16) Y. Hisada, Bulletin of the Seismological Society of America 90 387 (2000).
  * (17) Y. Hisada, A. Shibaya, and M.R. Ghayamghamian, Bulletin Earthquake Research Institute University of Tokyo 79 81 (2004).
  * (18) G. D. Bensen, M. H. Ritzwoller, M. P. Barmin, A. L. Levshin, F. Lin, M. P. Moschetti, N. M. Shapiro, and Y. Yang, Geophys. J. Int. 169, 1239 (2007).
  * (19) Z. J. Xu and X. Song, Proc. Natl. Acad. Sci. USA 106, 14207 (2009).
  * (20) Y. Okada, K. Kasahara, S. Hori, K. Obara, S. Sekiguchi, H. Fujiwara, and A. Yamamoto, Earth Planets Space 56, xv-xxviii (2004).
  * (21) S. Ide, G. C. Beroza, D. R. Shelly, and T. Uchide, Nature 447, 76 (2007).
  * (22) N. M. Shapiro, M. Campillo, L. Stehly, and M. H. Ritzwoller, Science 307, 1615 (2005).
  * (23) M. Campillo and A. Paul, Science 299 547 (2003).
  * (24) K. Yamasaki, A. Gozolchiani, S. Havlin, Phys. Rev. Lett. 100, 228501 (2008).
  * (25) A. Gozolchiani, K. Yamasaki, O. Gazit, S. Havlin, Europhys. Lett. 83, 28005 (2008).
  * (26) B. Podobnik, P. Ch. Ivanov, V. Jazbinsek, Z. Trontelj, H. E. Stanley, and I. Grosse, Phys. Rev. E Rapid Communication 71, 025104(R) (2005).
  * (27) S. Lennartz, V. N. Livina, A. Bunde, and S. Havlin, Europhys. Lett. 81, 69001 (2008).
  * (28) J. W. Kantelhardt, Y. Ashkenazy, P. C. Ivanov, A. Bunde, S. Havlin, T. Penzel, J. H. Peter, and H. E. Stanley, Phys. Rev. E 65, 051908 (2002).
  * (29) H. E. Stanley, L. A. N. Amaral, A. L. Goldberger, S. Havlin, P. C. Ivanov, and C.K. Peng, Physica A 270, 309 (1999).
  * (30) R. Karasik, N. Sapir, Y. Ashkenazy, P. C. Ivanov, I. Dvir, P. Lavie, and S. Havlin, Phys. Rev. 66, 62902 (2002).
  * (31) R. Mantegna and H. E. Stanley, Nature 376, 46 (1995).
  * (32) K. Yamasaki, L. Muchnik, S. Havlin, A. Bunde, and H. E. Stanley, Proc. Natl. Acad. Sci. USA 102, 9424 (2005).
  * (33) F. Wang, K. Yamasaki, S. Havlin, and H. E. Stanley, Phys. Rev. E 73, 026117 (2006).
  * (34) F. Wang, K. Yamasaki, S. Havlin, and H. E. Stanley, Phys. Rev. E 77, 016109 (2008).
  * (35) H. E. Stanley, V. Plerou, and X. Gabaix, Physica A 387, 3967 (2008).
  * (36) B. Podobnik, D. Horvatic, A. M. Petersen, and H. E. Stanley, Proc. Natl. Acad. Sci. USA 106, 22079 (2009).
  * (37) B. Podobnik, D. Horvatic, A. M. Petersen, B. Urosevic, and H. E. Stanley, Proc. Natl. Acad. Sci. USA 107, 18325 (2010).
  * (38) C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E. Stanley, A. L. Goldberger, Phys. Rev. E 49, 1685 (1994).
  * (39) K. Hu, P. C. Ivanov, Z. Chen, P. Carpena, H. E. Stanley, Phys. Rev. E 64, 011114 (2001).
  * (40) Z. Chen, P. C. Ivanov, K. Hu, and H.E. Stanley, Phys. Rev. E 65, 041107 (2002).
  * (41) Z. Chen, K. Hu, P. Carpena, P. Bernaola-Galvan, H. E. Stanley, and P. Ch. Ivanov, Phys. Rev. E 71, 011104 (2005).
  * (42) H. Takayasu, Fractals in the Physical Sciences (Manchaster U. Press, Manchaster, 1997).
  * (43) M. Kobayashi and T. Musha, IEEE Trans. Biomed. Eng. 29, 456 (1982).
  * (44) P. Ch. Ivanov, A. Yuen, B. Podobnik, Y. Lee, Phys. Rev. E 69, 056107 (2004).
  * (45) V. N. Livina, Y. Ashkenazy, P. Braun, R. Monetti, A. Bunde, S. Havlin, Phys. Rev. E 67, 042101 (2003).
  * (46) J. W. Kantelhardt, S. A. Zschiegner, E. Koscielny-Bunde, S. Havlin, A. Bunde, and H. E. Stanley, Physica A 316, 87 (2002).
  * (47) Z. Eisler and J. Kertesz, Eur. Phys. J. B 51, 145 (2006).
  * (48) D: Sornette, Proc. Natl. Acad. Sci. USA 99, 2522 (2002).
  * (49) T. Utu, Seismicity studies: a comprehensive review (University of Tokyo Press, Tokyo, 1999).

Figure 1: (a) Location map for the 46 broadband stations of Full Range
Seismograph Network of Japan (F-net) (red snow marks). Inset: An example of a
record of a seismic wave (Up-Down component). (b) A part of the coda wave
series indicated in inset of (a), as an example. (c) An example sign time
series where the positive sign (+1) represents a positive waveform, and the
negative sign (-1) represents a negative waveform in coda wave series of
seismic wave. (d) Interval time series ($\tau$) of the coda wave series for a
subset of the record shown in (b).
Figure 2:  Scaled mean conditional interval
$\langle\tau|\tau_{0}\rangle/\overline{\tau}$ vs $\tau_{0}/\overline{\tau}$ .
Five groups, one with no proximal earthquake and earthquakes with Richter
magnitude scale $M<4.5$, $M=4.5\sim 5.5$, $M=6.5\sim 6.5$, $M>6.5$. An
increasing trend implies a short-range correlation in the interval series.
Figure 3:  Long-range memory in interval clusters. $\tau_{0}$ signifies a
cluster of intervals, consisting of $n$ consecutive values that all are above
(denote as ”$+$”) or below (denote as ”$-$”) the median of the entire interval
records. Plots display the scaled mean interval conditioned on a cluster,
$\langle\tau|\tau_{0}\rangle/\overline{\tau}$ vs the size $n$ of the cluster
for five group intervals. The upper part (overplotted) of curves is for ”$+$”
clusters while the lower part is for ”$-$” clusters. The plots show that ”$+$”
clusters are likely to be followed by large intervals and ”$-$” clusters by
small intervals, consistent with long-term correlations in interval records.
Similar to Fig.2, the long-term correlation increases with earthquake size,
with exceptions for very large earthquakes.
Figure 4:  DFA fluctuation function $F(n)$ of 2011 Tohoku earthquake as a
function of time scale $n$ ($F(n)\propto n^{\alpha}$) for (a) sign time series
($\alpha+1=1.32$, ($\alpha<0.5$), indicates anticorrelations) and (b) interval
time series ($\alpha+1=1.69$, ($\alpha>0.5$), indicates correlations).
Figure 5:  Scaling exponent $\alpha$ vs seismic moment (Richter magnitude
scale) for (a) sign time series (correlation coefficient $Cor=-0.3604$), and
(b) interval time series (correlation coefficient $Cor=0.3602$). The values of
$\gamma$ show negative slope in the regression $\alpha$ vs seismic moment of
the sign series, and positive slope in the regression of the interval series.
Triangular symbols show the mean of exponent within each bin ( bins:
$<1e+15,1e+15\sim 1e+16,1e+16\sim 1e+17,1e+17\sim 1e+18,1e+18\sim
1e+19,1e+19\sim 1e+20,>1e+21$), the error bar shows the $\pm$ standard
deviation. The plots show a linear relationship between logarithmic earthquake
moment and scaling exponent $\alpha$ in the sign and interval series, with
exceptions for very large earthquakes.
Figure 6:  Scaling exponent $\alpha$ vs hypocenter depth for events where
Richter magnitude scale $M<5$ for (a) sign time series (b) interval time
series. Inset: scaling exponent $\alpha$ vs hypocenter depth requiring that
Richter magnitude scale $M>5$. (c) and (d) show Scaling exponent $\alpha$ vs
epicentral distance for events where Richter magnitude scale $M<5$. Inset:
scaling exponent $\alpha$ vs epicentral distance requiring that Richter
magnitude scale $M>5$. (c) sign time series, (d) interval series. All absolute
values of correlation coefficient are smaller than $0.1$, showing that
$\alpha$ is uncorrelated with both hypocenter depth and epicentral distance.