# Towards a proof of the Riemann Hypothesis

Guilherme Rocha de Rezende  Federal Institute of Brasilia-Brazil111E-mail:
<EMAIL_ADDRESS>
(August 27, 2024)

###### Abstract

In this article we will prove the Riemann Hypothesis for a infinite number of
choices of the imaginary part of the argument - $\Im(s)=T$.

##  1 Introduction

The Riemann Hypothesis is a famous conjecture made by Bernhard Riemann in his
article on prime numbers. Riemann, as indicated by the title of his article
[1], wanted to know the number of prime numbers in a given interval of the
real line, so he extended a Euler observation and defined a function called
Zeta. Riemann obtained an explicit formula, which depends on the non-trivial
zeros of the Zeta function, for the quantity he was looking for. Along the
way, Riemann mentions that probably all non-trivial zeros of the Zeta function
are, in the now called critical line, that is, when the complex argument
$s=\sigma+IT$ of the Zeta function has a real part equal to one-half. -
$\sigma=\frac{1}{2}$.

Despite being one of the most famous and extensively studied problems in
mathematics, the Riemann Hypothesis remains unsolved. However, progress has
been made in understanding its properties and implications, and many new
techniques and ideas have emerged from the study of this important problem. In
fact, the Riemann Hypothesis has been the driving force behind many
significant developments in mathematics over the past century and a half, and
it continues to be a central focus of research in the field of number theory
today. We will prove that the hypothesis is true for a large number of heights
T.

##  2 The main result

The result is as long as there does not exist a zero in $s=\frac{1}{2}+iT$
then cannot exist a zero in the line $\Im(s)=T$. The reasoning is based on the
argument of Backlund in the book Riemann’s Zeta Function [2] (pages 128 and
129). This occurs because the difference below is at most one so cannot exist
two zeros on the line $\Im(s)=T$.

$\displaystyle N(T+\delta)-N(T-\delta)=$
$\displaystyle\frac{1}{\pi}\vartheta(T+\delta)-\frac{1}{\pi}\vartheta(T-\delta)+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)$
(1)

N(M) is the number of zeros up to M (M excluded), $\vartheta(T)$ is a known
function and can be expanded in terms of T -

$\displaystyle\vartheta(T)=\frac{T}{2}\ln(\frac{T}{2\pi})-\frac{T}{2}-\frac{\pi}{8}+\frac{1}{48T}+\frac{7}{5760T^{3}}+...$
(2)

and the last term of the right hand side of (1) is

$\displaystyle\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)=\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)-\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right).$
(3)

The maximum of the last term of right hand side is one when the integrals
below take the values

$\displaystyle\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)=\frac{1}{2}$
(4)
$\displaystyle\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)=-\frac{1}{2}$
(5)

this occurs because the primitive is the logarithm and the imaginary parts of
the above contour integrals vary in the range $[-\frac{\pi}{2},\frac{\pi}{2}]$
where we are assuming that the real part of the zeta function -$\Re(\zeta(s))$
is never zero on contours \- $C_{1.1}$, $C_{1.2}$ and $C_{2}$ \- see the
graphs of the contours below.

Figure 1: Contour C 1.1

Figure 2: Contour C 1.2

Figure 3: Contour C 2

##  3 Proof

###  3.1 Argument Principle

The argument principle is a theorem in complex analysis that relates the
number of zeros and poles of a complex function to the behavior of its
argument (angle) as one travels around a closed contour. The theorem states
that the difference between the number of zeros and poles of a complex
function inside a closed contour is equal to the change in the argument of the
function as one travels around the contour. More formally, let $f(z)$ be a
complex function that is analytic inside and on a simple, closed contour C,
and let $Z(f)$ and $P(f)$ be the number of zeros and poles of $f(z)$ inside
$C$, respectively. Then:

$\displaystyle\frac{1}{2\pi
i}\ointctrclockwise_{C}\frac{f^{\prime}(z)}{f(z)}dz=Z(f)-P(f)$ (6)

where C denotes the contour integral around C in the counter clockwise
direction. In words, the argument principle states that the difference between
the number of zeros and poles of $f(z)$ inside C is equal to the winding
number of the image of C under $f(z)$ around the origin, divided by $2\pi i$.

Here is a figure to illustrate the use of the argument principle:

Figure 4: Contour C, Zeros and poles

The simple contour C (black), the zeros of f (blue) and the poles of f (red).
Here we have $\frac{1}{2\pi i}\ointctrclockwise_{\partial
R}\frac{f^{\prime}(z)}{f(z)}dz=4-5=-1$.

###  3.2 Backlund formula

Backlund proved using the argument principle in the Riemann’s $\xi(s)$
function $\xi(s)=\frac{1}{2}s(s-1)\pi^{-s/2}\Gamma(s/2)\zeta(s)$ that

$\displaystyle N(T)=\frac{1}{2\pi i}\ointctrclockwise_{\partial
R_{1}}\frac{\xi^{\prime}(s)}{\xi(s)}ds$ (7)

the below formula

$\displaystyle\text{N(T)}=\frac{1}{\pi}\vartheta(T)+1+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)$
(8)

for the number of zeros with multiplicity on critical strip up height - T,
where the region $\partial R_{1}$ is the rectangle
$\left(-\epsilon\leq\Re(s)\leq 1+\epsilon,0\leq\Im(s)\leq T\right)$ and
$C_{1}$ is the portion of $\partial R_{1}$ from $1+\epsilon$ to
$\frac{1}{2}+IT$, where we assume that there are no zeros in the region
$\partial R_{1}$. Backlund observed also that if there are no zeros of the
real part of zeta function on contour $C_{1}$ then the term
$\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)$
has a maximum value $\frac{1}{2}$ and a minimum value -$\frac{1}{2}$. The
contours are:

Figure 5: Contour $\partial R_{1}$

Figure 6: Contour $C_{1}$

####  3.2.1 Avoiding zeros on the contour $R_{1}$

In the case of the above contour $\partial R_{1}$ is sufficient to have no
zeros in the line $\Im(s)=T$ because zeta function is not zero for real part
of argument below $\sigma=0$ and above $\sigma=1$ and there are no zeros of
$\xi(s)$ when $\Im(s)=0$ as we can see in the plot of $\xi(s)$ when $s$ is
real

Figure 7: $\xi(s)$ function for real argument on critical strip

###  3.3 The main idea

The central idea is to improve the approximation made by Backlund in order to
show that the zeros are always on the critical line and they are simple at
least in some cases. We take the difference $N(T+\delta)-N(T-\delta)$ ,
between the number of zeros with multiplicities, where T is the imaginary
height where are the supposed zeros and prove that the difference has an upper
bound equal to one. Therefore at most only one zero with multiplicity one is
possible for each height - T. Take the difference $N(T+\delta)-N(T-\delta)$ we
get

$\displaystyle N(T+\delta)-N(T-\delta)=\frac{1}{2\pi
i}\ointctrclockwise_{\partial
R_{1.1}}\frac{\xi^{\prime}(s)}{\xi(s)}ds-\ointctrclockwise_{\partial
R_{1.2}}\frac{\xi^{\prime}(s)}{\xi(s)}ds.$ (9)

and following Backlund argument [2] we arrive at

$\displaystyle N(T+\delta)-N(T-\delta)=$
$\displaystyle\frac{1}{\pi}\vartheta(T+\delta)-\frac{1}{\pi}\vartheta(T-\delta)+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)-\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)$
(10)

or

$\displaystyle
N(T+\delta)-N(T-\delta)=\frac{1}{\pi}\vartheta(T+\delta)-\frac{1}{\pi}\vartheta(T-\delta)+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right).$
(11)

The contours are

Figure 8: Contour $\partial R_{1.1}$

Figure 9: Contour $\partial R_{1.2}$

Figure 10: Contour $\partial R_{2}$

Figure 11: Contour $C_{1.1}$

Figure 12: Contour $C_{1.2}$

####  3.3.1 The crucial point of the main idea

The crucial point is the observation made my Backlund [2] which we adapt to
this case:if the real part of zeta function is not null on the contours
$C_{1.1}$,$C_{1.2}$ and $C_{2}$ then the error in the last right hand side of
the below equation

$\displaystyle
N(T+\delta)-N(T-\delta)=\frac{1}{\pi}\vartheta(T+\delta)-\frac{1}{\pi}\vartheta(T-\delta)+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right).$
(12)

or

$\displaystyle N(T+\delta)-N(T-\delta)=$
$\displaystyle\frac{1}{\pi}\vartheta(T+\delta)-\frac{1}{\pi}\vartheta(T-\delta)+\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)-\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)$
(13)

is at most one when

$\displaystyle\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.1}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)=\frac{1}{2}$
(14)

and

$\displaystyle\frac{1}{\pi}\Im\left(\ointctrclockwise_{C_{1.2}}\frac{\zeta^{\prime}(s)}{\zeta(s)}ds\right)=-\frac{1}{2}.$
(15)

We can always choose $\delta$ and $\epsilon$ appropriately to avoid zeros on
the contour of line integral. See figure below

Figure 13: Contour $C_{2}$

Suppose we have the nearest zero of the real part of zeta function in the
vicinity of the line segment $\Im(s)=T$. We always can choose $\delta$ and
$\epsilon$ in an convenient manner to avoid the nearest zero of the real part
of zeta function and as a consequence all the null real parts of $\zeta(s)$ if
there is no zero on $s=\frac{1}{2}+iT$. This occurs because we cannot have
accumulation points where zeros are ”glued together”. Therefore the Riemann
hypothesis is true for line segments $\Im(s)=T$ when there is no zero on
$s=\frac{1}{2}+iT$.

##  4 Conclusion

After the efforts of many mathematicians [3] the problem was partially solved,
we have proved that Riemann’s hypothesis is true for each $T$ when there is no
zero in $s=\frac{1}{2}+iT$. We have hope to complete solve the problem with an
appropriate deformation of the contour when there is a zero on
$s=\frac{1}{2}+iT$ for an infinite number of Ts.

##  5 Acknowledgement

I, posthumously, dedicate to my mother Edna Vieira Rocha de Rezende.
I dedicate to my father Rodolpho Antonio de Rezende.
I also dedicate to my sister Gisella Rocha de Rezende.
I also dedicate to my brother Gustavo Rocha de Rezende.

## References

  * [1] B. Riemann. Ueber die Anzahl der Primzahlen unter einer gegebenen GrÃ¶sse. Monatsberichte der Berliner Akademie, 1859. Translation: R. Wilkins, David. On the Number of Prime Numbers less than a Given Quantitiy, 1998.
  * [2] H. M. Edwards. Riemann’s Zeta Function. Dover edition (2001).
  * [3] J. Derbyshire. Prime Obsession: Bernhard Riemann and the Greatest Unsolved Problem in Mathematics. Dover books (2003).


remarkRemark hypothesisHypothesis claimClaim pyFFS: Fast Fourier Series
Computation with PythonE. Bezzam, S. Kashani, P. Hurley, M. Vetterli, and M.
Simeoni ex_supplement

# pyFFS: A Python Library for Fast Fourier Series Computation and
Interpolation with GPU Acceleration††thanks: Published on August 18, 2022.
This work was in part funded by the Swiss National Science Foundation (SNSF)
under grants 200021_181978/1 SESAM –- Sensing and Sampling: Theory and
Algorithms (first, second and fifth authors) and CRSII5 180232 FemtoLippmann –
Digital twin for multispectral imaging (second and fifth authors). For part of
this work, the second, third, and fifth authors were also with the Foundations
of Cognitive Solutions group in the IBM Research Laboratory of Zurich.

Eric Bezzam Audiovisual Communications Laboratory, École Polytechnique
Fédérale de Lausanne, Switzerland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Sepand Kashani22footnotemark: 2
Paul Hurley Western Sydney University, Australia ().
<EMAIL_ADDRESS>Martin Vetterli22footnotemark: 2 Matthieu
Simeoni Centre for Imaging, École Polytechnique Fédérale de Lausanne,
Switzerland<EMAIL_ADDRESS>
###### Abstract

Fourier transforms are an often necessary component in many computational
tasks, and can be computed efficiently through the fast Fourier transform
(FFT) algorithm. However, many applications involve an underlying continuous
signal, and a more natural choice would be to work with e.g. the Fourier
series (FS) coefficients in order to avoid the additional overhead of
translating between the analog and discrete domains. Unfortunately, there
exists very little literature and tools for the manipulation of FS
coefficients from discrete samples. This paper introduces a Python library
called pyFFS for efficient FS coefficient computation, convolution, and
interpolation. While the libraries SciPy and NumPy provide efficient routines
for discrete Fourier transform coefficients via the FFT algorithm, pyFFS
addresses the computation of FS coefficients through what we call the fast
Fourier series (FFS). Moreover, pyFFS includes an FS interpolation method
based on the chirp Z-transform that can make it more than an order of
magnitude faster than the SciPy equivalent when one wishes to perform
distortionless bandlimited interpolation. GPU support through CuPy is readily
available, and allows for further acceleration: an order of magnitude faster
for computing the 2-D FS coefficients of $1000\times 1000$ samples and nearly
two orders of magnitude faster for 2-D interpolation. As an application, we
discuss the use of pyFFS in Fourier optics. pyFFS is available as an open
source package at https://github.com/imagingofthings/pyFFS, with documentation
at https://pyffs.readthedocs.io.

###### keywords:

fast Fourier series, bandlimited interpolation, chirp Z-transform, numerical
library, Python, GPU

65T40, 97N80, 97N50, 42B05

##  1 Introduction

Discretization is an inevitable part of digital signal processing. Although
the universe around us moves and shakes with infinite precision, our computers
can only handle so much. However, with a useful model and the appropriate
analog and digital processing, we can faithfully simulate and process
continuous-domain processes. This idea is the essence of the Nyquist-Shannon
sampling theorem for bandlimited signals [19] and more generally for signals
that have a finite rate of innovation [3]. Simply put: a function $f(t)$ with
finite degrees-of-freedom can be completely determined by a finite number of
its samples. In his paper [19], Shannon even admits that the theorem so often
associated with him is a fact “which is common knowledge in the communication
art.” However, without the proper formulation, common knowledge may never seep
into common practice.

The inception of the fast Fourier transform (FFT) algorithm [5] has a similar
story. Arguably one of the most influential algorithms in computational and
natural sciences, the FFT was nearly not published, as one of the original
authors, John Tukey, felt that it was a simple observation that was probably
already known and unworthy of attention.111As history would later tell, Tukey
was indeed right as Gauss had written a paper about an interpolation technique
with essentially the same idea [8]. With the digital revolution that ensued
afterwards, the FFT quickly became an indispensable algorithm in digital
signal processing. Once again, we see that without the proper formulation and
implementation, common knowledge may not establish itself in common practice.

Continuing along these lines, we formulate old Fourier analysis tricks with a
new perspective and propose new tools to efficiently compute and interpolate
Fourier series (FS) coefficients. As is well known, for discrete samples of a
periodic, bandlimited signal, we can perfectly recover the FS coefficients of
the underlying function when sampling is done according to the Nyquist-Shannon
sampling theorem [21]. Moreover, continuous-domain operations, such as
convolution and interpolation, can be implemented in terms of distortionless
discrete operations on the FS coefficients. These mathematical facts can be
used in practice to design an algorithm for fast FS computation and
interpolation, which we call the _fast Fourier series_ (FFS) algorithm. While
deceptively simple, the FFS algorithm is surprisingly neither described in
signal processing textbooks nor implemented in numerical computing libraries.
For example, NumPy [11] and SciPy [22] in the Python ecosystem focus mainly on
the discrete Fourier transform (DFT) and related operations. We aim to change
this by providing an efficient and easy-to-use interface to the FFS algorithm
via a Python package called pyFFS. Furthermore, distortionless discrete
operations with pyFFS are not limited to periodic signals. In fact, one is
often interested in working with signals of compact-support, which can be cast
as periodic by repeating the signal with a period larger than or equal to its
support.

So why use pyFFS rather than the FFT routines from NumPy/SciPy? One reason is
convenience when working with continuous-domain signals. The philosophy of
pyFFS is to retain the continuous-domain perspective, often neglected when
using numerical libraries such as NumPy and SciPy, which allows for much
clearer code as seen in LABEL:lst:prop. This can prevent common pitfalls due
to an invalid conversion between discrete and continuous domains. Moreover, FS
coefficients are an important component of the Non Uniform Fast Fourier
Transform (NUFFT), extensively used in MRI imaging and other computational
imaging modalities [2]. Another reason is efficiency. We benchmark pyFFS with
equivalent functions in SciPy, observing scenarios in which the proposed
library is more than an order of magnitude faster, e.g. for interpolation.
Moreover, GPU support has been seamlessly incorporated for an even faster
implementation. Just as the FFT implementation via NumPy and SciPy can be
readily used for an efficient $\mathcal{O}(N\log N)$ analysis and synthesis of
discrete sequences, pyFFS offers the same ease-of-use and performance
capabilities for discrete representations of continuous-domain signals.

The paper is organized as follows. In Section 2, we go over the theory behind
pyFFS and motivate why working with FS coefficients may be preferable. We also
present theorems showing how FS coefficients of periodic, bandlimited signals
can be computed and interpolated exactly and efficiently. In Section 3, the
Python user interface of pyFFS is discussed.222A more extensive documentation
can be found at https://pyffs.readthedocs.io. Section 4 presents benchmarking
tests showing the efficiency of pyFFS against SciPy and the gains of GPU
acceleration. In Section 5, we discuss an application of pyFFS in Fourier
optics. Another application of pyFFS can be found in [7], where the authors
make use of the FFS algorithm to efficiently compute multidimensional periodic
splines. Section 6 concludes the paper.

##  2 Theory

In this section, we give an overview of the theory behind pyFFS. This can be
useful to those who wish to understand what is happening under the hood and
sheds light as to when using FS coefficients over DFT ones may be of interest.

###  2.1 Numerically compute FS coefficients

At first glance, it seems unlikely that FS coefficients could be of practical
use in computational scenarios, as their very definition necessitates the
integral of a continuous function [21]

(1)
$X_{k}^{\text{FS}}=\frac{1}{T}\int_{T_{c}-\frac{T}{2}}^{T_{c}+\frac{T}{2}}x(t)\exp\left(-j\frac{2\pi}{T}kt\right)dt,$

where $x:\mathbb{R}\to\mathbb{C}$ is a $T$-periodic function and $T_{c}$ is
any period mid-point.

The synthesis equation, on the other hand, expresses this function in terms of
a discrete, albeit infinite, set of samples, namely the FS coefficients

(2)
$x(t)=\sum_{k\in\mathbb{Z}}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}kt\right).$

For a bandlimited signal, we can write Eq. 2 as

(3)
$x(t)=\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}kt\right),\quad
N\geq 0,$

where $x(t)$ is said to have a bandwidth of $N_{\text{FS}}=2N+1$.333This is
equivalent to a maximal frequency of $N/T$ Hz. Moreover, by taking uniform
samples with a sampling period of $T_{s}=T/N_{\text{FS}}$, we obtain

$\displaystyle x(nT_{s})=\left[{\bf{x}}\right]_{n}$
$\displaystyle=\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}k(nT_{s})\right)$
$\displaystyle\overset{(a)}{=}\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{N_{\text{FS}}}kn\right)$
(4)
$\displaystyle\overset{(b)}{=}\sum_{k=-N}^{N}X_{k}^{\text{FS}}W_{N_{\text{FS}}}^{-nk},$

where $(a)$ uses $(T_{s}/T)=1/N_{\text{FS}}$ and $(b)$ uses
$W_{N}=\exp\left(-j\frac{2\pi}{N}\right)$. Uniform sampling is equivalent to
convolution with a Dirac stream in the FS coefficient domain. Therefore the
sequence $\left\\{X_{k}^{\text{FS}}\in\mathbb{C},\;k\in\mathbb{Z}\right\\}$ is
$N_{\text{FS}}$-periodic, as shown in Figure 1. Moreover, choosing $T_{s}\leq
T/N_{\text{FS}}$ ensures that the FS coefficients remain “intact”, i.e. there
is no distortion due to aliasing from the overlapping spectrum replicas.

In Section 2.1 we have an expression that is equivalent to what is commonly
referred to as the discrete Fourier series (DFS) [16]

(5)
$\left[{\bf{x}}\right]_{n}=\sum_{k=0}^{N_{\text{FS}}-1}\left[{\bf{X}}\right]_{k}\exp\left(j\frac{2\pi}{N_{\text{FS}}}kn\right),\quad
n\in\mathbb{Z},$

which is similar to the inverse discrete Fourier transform (IDFT), except that
the time index $n$ spans all integers. As the coefficients of the sampled
function are $N_{\text{FS}}$-periodic, shifting the summation in Section 2.1
to $[0,N_{\text{FS}}-1]$, as in Eq. 5, is equivalent to the summation over
$[-N,N]$.

Figure 1: Visualizing the Fourier series (FS) coefficients of a sampled,
bandlimited signal: (left) original FS coefficients of a signal whose
bandwidth is $N_{\text{FS}}=2N+1$; (right) sampled FS coefficients where
$N_{\text{s}}\geq N_{\text{FS}}$ such that the sampling period is
$T_{\text{s}}=T/N_{\text{s}}$. When performing synthesis, i.e. interpolating
coefficients for time-domain values, our formulation Section 2.1 sums over
coefficients centered around zero (dashed box), while the discrete Fourier
series Eq. 5 sums over positive coefficients (dotted box). This is equivalent
as the sampled Fourier series coefficients are a periodic sequence.

With Section 2.1, we can numerically compute the FS coefficients of a
bandlimited, periodic signal by setting up a system of equations with at least
as many samples as FS coefficients:

$\begin{bmatrix}\tiny\left[{\bf{x}}\right]_{-N}\\\
\left[{\bf{x}}\right]_{-N+1}\\\ \vdots\\\
\left[{\bf{x}}\right]_{N}\end{bmatrix}=\scalebox{0.9}{
$\begin{bmatrix}W_{N_{FS}}^{-(-N)(-N)}&W_{N_{FS}}^{-(-N)(-N+1)}&\cdots&W_{N_{FS}}^{-(-N)(N)}\\\\[5.0pt]
W_{N_{FS}}^{-(-N+1)(-N)}&W_{N_{FS}}^{-(-N+1)(-N+1)}&\cdots&W_{N_{FS}}^{-(-N+1)(N)}\\\
\vdots&\vdots&\ddots&\vdots\\\\[5.0pt]
W_{N_{FS}}^{-(-N)(N)}&W_{N_{FS}}^{-(-N+1)(N)}&\cdots&W_{N_{FS}}^{-(N)(N)}\\\
\end{bmatrix}$}\begin{bmatrix}X_{-N}^{\text{FS}}\\\ X^{\text{FS}}_{-N+1}\\\
\vdots\\\ X^{\text{FS}}_{N}\end{bmatrix}.$

A naive approach could solve this system of equations for the FS coefficients
by multiplying the left-hand side with the inverse of the matrix, requiring
$\mathcal{O}(N_{FS}^{3})$ operations. However, the above matrix is unitary,
allowing us to multiply both sides by the complex conjugate of the matrix
rather than having to invert it. This gives a complexity of
$\mathcal{O}(N_{FS}^{2})$. With some manipulation, the matrix-vector
relationship can be turned into one where the matrix is the DFT matrix. We can
then exploit the FFT algorithm to further reduce the complexity to
$\mathcal{O}(N_{FS}\log N_{FS})$ [5]. This result is given in the following
theorem, with the proof provided in Appendix A.

######  Theorem 2.1 (Fast Fourier series).

Let $x:\mathbb{R}\to\mathbb{C}$ be a $T$-periodic function of bandwidth
$N_{\text{FS}}=2N+1$, with $T_{c}\in\mathbb{R}$ the mid-point of any period.
Let the zero-padding amount $Q\in 2\mathbb{N}$ be an arbitrary even integer
such that the number of samples $N_{\text{s}}=N_{\text{FS}}+Q$. Then

(6)
$\displaystyle{\bf{x}}=N_{\text{s}}\operatorname{IDFT}_{N_{s}}\left({\bf{X}}^{\text{FS}}\odot
B_{1}^{{\bf{E}}_{1}}\right)\odot B_{2}^{N{\bf{E}}_{2}},$ (7)
$\displaystyle{\bf{X}}^{\text{FS}}=\frac{1}{N_{\text{s}}}\operatorname{DFT}_{N_{s}}\left({\bf{x}}\odot
B_{2}^{-N{\bf{E}}_{2}}\right)\odot B_{1}^{-{\bf{E}}_{1}},$

where

(8)
$\displaystyle{\bf{x}}=\left[{\bf{x}}(t_{0}),\ldots,{\bf{x}}(t_{M}),{\bf{x}}(t_{-M}),\ldots,{\bf{x}}(t_{-1})\right]\in\mathbb{C}^{N_{\text{s}}},$
(9) $\displaystyle t_{n}=T_{c}+\frac{T}{N_{\text{s}}}n,\quad n\in\mathbb{Z},$
$\displaystyle M=(N_{\text{s}}-1)/2,$
$\displaystyle{\bf{X}}^{\text{FS}}=\left[X_{-N}^{\text{FS}},\ldots,X_{N}^{\text{FS}},{\bf{0}}_{Q}\right]\in\mathbb{C}^{N_{\text{s}}},$

and  $\displaystyle B_{1}=\exp\left(j\frac{2\pi}{T}T_{c}\right)\in\mathbb{C},$
$\displaystyle{\bf{E}}_{1}=\left[-N,\ldots,N,{\bf{0}}_{Q}\right]\in\mathbb{Z}^{N_{\text{s}}},$
$\displaystyle
B_{2}=\exp\left(-j\frac{2\pi}{N_{\text{s}}}\right)\in\mathbb{C},$
$\displaystyle{\bf{E}}_{2}=\left[0,\ldots,M,-M,\ldots,-1\right]\in\mathbb{Z}^{N_{\text{s}}}.$

$\operatorname{DFT}_{N}$ and $\operatorname{IDFT}_{N}$ denote the length-$N$
DFT and IDFT respectively. The operation $\odot$ is an element-wise
multiplication, with which in Eqs. 6 and 7 we modulate the input and the
output of the standard DFT and IDFT operations. This modulation shifts the
summation bounds in Section 2.1 to those of the DFT, allowing direct use of
the FFT for an efficient computation. The zero-padding $Q$ can be used to set
the FFT length to a highly composite value for faster computation.

The above theorem assumes an odd-length $N_{\text{s}}$. For an even-length
sequence, a slight modification of the sample locations $t_{n}$ and modulation
terms $B_{1},B_{2}$ is necessary. As the general idea is the same, we refer
the reader to Appendix A for the presentation and proof for the odd- and the
even-length cases.

###  2.2 Efficient interpolation of FS coefficients

There is nothing particularly unique about the FS coefficients given by
Theorem 2.1 in Eq. 7. They are essentially DFT coefficients modulated,
reordered, and scaled so that the resulting sequence is directly the FS
coefficients within $[-N,N]$. However, their intrinsic link to a continuous
signal makes them all the more convenient when we consider interpolation at
arbitrary sample locations, rather than griding or sub-griding at the discrete
level. From a practitioner’s perspective, it may also be more intuitive to
think and work directly with the continuous-domain sample locations rather
than the corresponding sample (or subsample) values.

Equation 3 could be used to compute the value of the underlying bandlimited
function at arbitrary sample locations, requiring
$\mathcal{O}(MN_{\text{FS}})$ operations for $M$ sample points. Below we show
one way these FS coefficients can be interpolated at $M$ regularly-spaced
samples in a more efficient manner. This is done by making use of the chirp
Z-transform (CZT).

######  Definition 2.2 (Chirp Z-transform).

Let ${\bf{x}}\in\mathbb{C}^{N}$. The length-$M$ _chirp Z-transform_
$\operatorname{CZT}_{N}^{M}({\bf{x}})\in\mathbb{C}^{M}$ of parameters
$A,W\in\mathbb{C}^{*}$ is defined as [17]

(10)
$\left[\operatorname{CZT}_{N}^{M}({\bf{x}})\right]_{k}=\sum_{n=0}^{N-1}\left[{\bf{x}}\right]_{n}A^{-n}W^{nk},\quad
k\in\left\\{0,\ldots,M-1\right\\},$

where $A$ is the complex starting point and $W$ is the complex ratio between
points along a logarithmic spiral contour. The CZT is a generalization of the
DFT which samples the $Z$ plane at uniformly-spaced points along the unit
circle.

$\operatorname{CZT}_{N}^{M}$ can be efficiently computed using the
$\operatorname{DFT}$ and $\operatorname{IDFT}$ in $\mathcal{O}(L\log L)$
operations, where $L\geq N+M-1$. This is done via Bluestein’s algorithm [4]. A
similar formulation of the CZT and its efficient computation is known as the
fast fractional FT algorithm [1].

The theorem below makes use of the CZT and Bluestein’s algorithm to perform
efficient interpolation at regularly-spaced values.

######  Theorem 2.3 (FS interpolation, $x:\mathbb{R}\to\mathbb{C}$).

Let $x:\mathbb{R}\to\mathbb{C}$ be a $T$-periodic function of bandwidth
$N_{\text{FS}}=2N+1$. Let $a<b\in\mathbb{R}$ be the end-points of an interval
on which we want to evaluate $M$ equi-spaced samples of $x$. Then

(11)
${\bf{x}}=A^{N}\operatorname{CZT}_{N_{\text{FS}}}^{M}({\bf{X}}^{\text{FS}})\odot
W^{-N{\bf{E}}},$

where

$\displaystyle{\bf{x}}=\left[x(t_{0}),\ldots,x(t_{M-1})\right]\in\mathbb{C}^{M},$
$\displaystyle{\bf{X}}^{\text{FS}}=\left[X_{-N}^{\text{FS}},\ldots,X_{N}^{\text{FS}}\right]\in\mathbb{C}^{N_{\text{FS}}},$
$\displaystyle t_{n}=a+\frac{b-a}{M-1}n,\quad n\in\mathbb{Z},$

$\displaystyle A=\exp\left(-j\frac{2\pi}{T}a\right),$ $\displaystyle
W=\exp\left(j\frac{2\pi}{T}\frac{b-a}{M-1}\right),$
$\displaystyle{\bf{E}}=\left[0,\ldots,M-1\right]\in\mathbb{N}^{M},$

and $\operatorname{CZT}$ (of parameters $A,W$) is as defined in Definition
2.2.

The proof can be found in Appendix B. A similar interpolation technique with
the fractional FT is presented in [1].

With Theorem 2.3, one can interpolate sub-sections of a period efficiently.
Moreover, it is possible to perform DFTs of a smaller length than what would
be normally required with a more standard IDFT interpolation approach, namely
zero-padding the DFT coefficients and taking a longer IDFT for an increase in
temporal resolution across the _entire_ period.

#### Complexity comparison

To get an idea of when it is preferable to use the proposed technique over
zero-padding DFT coefficients, we present a rough complexity analysis. For a
$T$-periodic function that we would like to evaluate at steps of $\Delta t$,
we would need $N_{\text{target}}=\lceil T/\Delta t\rceil$ samples within a
single period. If we had $N<N_{\text{target}}$ samples,444With
$N>N_{\text{FS}}$ so that we can have ideal reconstruction according to the
Nyquist-Shannon sampling theorem. we would need to pad the DFT coefficients
with $(N_{\text{target}}-N)$ zeros in order to get this temporal resolution,
resulting in an IDFT with complexity $\mathcal{O}(N_{\text{target}}\log
N_{\text{target}})$. With Theorem 2.3, interpolating at steps of $\Delta t$
over the entire period $T$ would lend to a computational complexity of
$\mathcal{O}\big{(}(N_{\text{target}}+N_{\text{FS}})\log(N_{\text{target}}+N_{\text{FS}})\big{)}$,
which is certainly not advantageous. The benefits arise when we wish to
interpolate over a smaller region within the period, i.e. when zooming in on a
section of length $mT$ with $m\in(0,1)$. In such a scenario, the number of
interpolation points is $M=\lceil mT/\Delta t\rceil$. As the complexity of
Theorem 2.3 is
$\mathcal{O}\big{(}(M+N_{\text{FS}})\log(M+N_{\text{FS}})\big{)}$, it can be
more efficient than zero-padding the DFT coefficients when $(M+N_{\text{FS}})$
is less than the padded IDFT length $N_{\text{target}}=\lceil T/\Delta
t\rceil$.

Both the proposed approach in Theorem 2.3 and interpolation by zero-padding
the DFT coefficients are bandlimited interpolation techniques. Given all the
FS or DFT coefficients of a bandlimited signal, both yield a distortionless
interpolation. A comparison between these two approaches is done in Section
4.2.

Another approach for bandlimited interpolation is sinc interpolation with
Dirichlet apodization [20], and it allows one to focus on a specific region
much like Theorem 2.3. However, its complexity is given by
$\mathcal{O}(MN_{\text{FS}})$, which quickly becomes prohibitive as the number
of interpolation points or the bandwidth increases. So we do not consider this
approach in our comparison.

###  2.3 When to use FS coefficients

In our discussion above, we already came across two requirements on the input
in order to numerically compute its FS coefficients from a discrete set of
samples without distortion:

  * •

_Periodic_ as the FS coefficients are defined for such signals.

  * •

_Bandlimited_ so that we can perfectly recover the FS coefficients from a
finite set of discrete samples.

In fact, the DFT makes similar assumptions on periodicity. As we saw with the
DFS, when we consider indices outside of the sequence’s finite support, we
observe a periodic structure. Consequently, when one wishes to interpolate
outside the support of a finite sequence, it is common to apply a tapering
window in order to avoid sharp discontinuities at the boundaries. Similarly,
in order to apply Theorem 2.1, we can consider the samples of an arbitrary
continuous function as a truncation that we periodize and optionally taper to
remove any discontinuities at the boundaries. For example, in Figure 2 we
compute the FS coefficients of a truncated sinusoid with a large discontinuity
at the border. Applying a tapering window attenuates the boundaries to zero,
thereby removing the discontinuity as seen in the periodization (bottom of
Figure 2(a)) and suppressing higher frequency coefficients that arise due to
this discontinuity (Figure 2(b)).

(a)

(b)

Figure 2: Example to motivate tapering window before Fourier analysis. Notice
in Fig. 2(b), the higher frequency coefficients that appear in the Fourier
analysis due to the discontinuity (top), and how they are suppressed when
applying a tapering window (below).

Although not necessarily required for the DFT, bandlimitness of the input
sequence is desired prior to sampling so that distortion due to aliasing is
minimized. However, after sampling there is not much that can be done to
remove such aliasing. The same is true for computing the FS coefficients.
However, if one has control on the acquisition process, the sampling rate can
be set to avoid spectral overlap or the signal can be bandlimited
appropriately. As we saw in Eq. 3, bandlimiting corresponds to truncating the
FS summation.

FS coefficients coupled with knowledge about the underlying continuous signal
therefore allows us to better understand, and perhaps control, where errors
could arise during our analysis and processing:

  * •

From aliasing if the underlying function could not be bandlimited or sampled
appropriately.

  * •

From truncation of the FS coefficients to ensure bandlimitedness. The mean-
squared error of this truncation can be bounded [9].

  * •

From a tapering window to remove border discontinuities if a function was
simply periodized by truncation and repetition.

Applying the DFT has similar constraints and consequences that tend to get
ignored in practice, as they may not be critical for the application or cannot
be corrected with the discrete samples alone.

In practice, discrete samples are bandlimited due to the sampling operation,
although there might be aliasing if the proper filter could not be used
beforehand, e.g. in medical imaging [20]. In such scenarios, we can only
expect to obtain the FS coefficients of the underlying aliased signal.

###  2.4 Multidimensional

The above methods are not limited to a single dimension. As the DFT plays a
central role, the FS coefficient and interpolation tools can be extended to
multidimensional in a similar fashion as the DFT. By leveraging the property
that the multidimensional DFT can be computed as a composition of one-
dimensional DFTs along each dimension, the FS computation and interpolation
methods above, namely Theorems 2.1 and 2.3, can be similarly computed along
each dimension of the samples and the FS coefficients respectively.

This application of the 1-D FFT along each dimension results in a complexity
of $\mathcal{O}(N\log N)$, where $N=N_{1}\cdot N_{2}\cdot\cdots N_{D}$ for a
general N-D signal. More efficient algorithms for performing the
multidimensional DFT do exist [6], which we make use of in our implementations
of Theorems 2.1 and 2.3 in multiple dimensions.

##  3 pyFFS overview and usage

pyFFS is a Python library for performing efficient and distortionless FS
coefficient computation, convolution, and interpolation for periodic,
bandlimited signals. The goal is to provide an intuitive tool to numerically
work with such signals of any dimension $D$. Just as the FFT functions from
NumPy [11] and SciPy [22] can be used without too much thought about the
internal details, we have created a user interface for FS computations that
can also be used out-of-the-box for the appropriate scenario.

###  3.1 Fourier series analysis and synthesis

The user interface for 1-D functions is shown below. Note that the samples
provided to ffs must be in the same order as specified in Eq. 8, which is not
in chronological order. The method ffs_sample returns the timestamps and
indices necessary for ensuring the samples provided to ffs are in the expected
order.

⬇

# determine appropriate timestamps and indices for rearranging input

sample_points, idx = pyffs.ffs_sample(

T, # function period

N_FS, # function bandwidth, i.e. number of FS coefficients (odd)

T_c, # function center

N_s # number of samples

)

# sample a known function at the correctly ordered timestamps

x = pyffs.func.dirichlet(sample_points, T, T_c, N_FS)

# OR rearrange ordered samples using ‘idx‘

# x = x[idx]

# compute FS coefficients

x_FS = pyffs.ffs(x, T, T_c, N_FS)

# back to samples with inverse transform

x_r = pyffs.iffs(x_FS, T, T_c, N_FS) # equivalent to x

Listing 1: 1-D fast Fourier series analysis and synthesis.

The user interface for the general N-D case is shown below, with the specific
example of 2-D. As in the 1-D case, samples provided to ffsn are not in
increasing order of the input variables. The method ffsn_sample returns the
locations and indices necessary for making sure the samples provided to ffsn
are in the expected order. Alternatively, the method ffs_shift can be used to
reorder the samples.

⬇

T = [T_x, T_y] # list of periods for each dimension

T_c = [T_cx, T_cy] # list of function centers for each dimension

N_FS = [N_FSx, N_FSy] # list of function bandwidths for each dimension

N_s = [N_sx, N_sy] # number of samples per dimension

# determine appropriate timestamps and indices for rearranging input

sample_points, idx = pyffs.ffsn_sample(T=T, N_FS=N_FS, T_c=T_c, N_s=N_s)

# sample a known function at the correctly ordered timestamps

x = pyffs.func.dirichlet_2D(sample_points, T, T_c, N_FS)

# OR rearrange ordered samples

# x = pyffs.ffs_shift(x)

# compute FS coefficients

x_FS = pyffs.ffsn(x, T=T, T_c=T_c, N_FS=N_FS)

# go back to samples

x_r = pyffs.iffsn(x_FS, T=T, T_c=T_c, N_FS=N_FS) # equivalent to x

Listing 2: 2-D fast Fourier series analysis and synthesis.

###  3.2 Circular convolution

The user interface for N-D functions is shown below. f and h are function
values at the sampling points specified by ffsn_sample, namely they must be
functions of the same period T, with the same period center T_c, and of the
same bandwidth N_FS.

Samples can be provided in their natural order or in the order expected by
ffsn. By default, the argument reorder is set to True, such that samples are
expected in their natural order and are reordered internally. The output
samples are returned in the same order as the inputs.

⬇

out = pyffs.convolve(

f=f, # samples of one function in the convolution

h=h, # samples of the other function in the convolution

T=T, # period(s) of both functions along all dimensions

T_c=T_c, # period center(s) of both functions along all dimensions

N_FS=N_FS, # number of FS coefficients for both functions along all dimensions

reorder=True # whether input samples should be reordered into expected order
for ffsn

)

Listing 3: Circular convolution of two N-D functions through Fourier series
coefficients.

###  3.3 Interpolation

The user interface for 1-D functions is shown below.

⬇

x_interp = pyffs.fs_interp(

x_FS, # FS coefficients in increasing order of index

T, # period

a, # start time

b, # stop stop

M # number of points

)

Listing 4: 1-D fast Fourier series interpolation.

The user interface for the general N-D case is shown below, with the specific
example of 2-D.

⬇

x_interp = pyffs.fs_interpn(

x_FS, # multidimensional FS coefficients

T=[T_x, T_y], # list of periods for each dimension

a=[a_x, a_y], # list of start points for each dimension

b=[b_x, b_y], # list of stop points for each dimension

M=[M_x, M_y] # number of samples per dimension

)

Listing 5: 2-D fast Fourier series interpolation.

In both cases, the provided FS coefficients must be ordered such that the
indices are in increasing order, as returned by ffs and ffsn.

###  3.4 GPU support

GPU usage can lead to a significant reduction in computation time if a task
consists of many operations that can be done in parallel. The FFT algorithm
can be parallelized and could therefore benefit from such a reduction in
computation time. As the fast FS algorithm presented in Theorem 2.1 makes use
of the DFT and IDFT, it can directly benefit from this speed-up, and so can
the interpolation of FS coefficients described in Theorem 2.3, as Bluestein’s
algorithm for the CZT employs the DFT and IDFT.

GPU support is available through the CuPy library [15]. If the appropriate
version of CuPy is installed,555See installation guide:
https://docs.cupy.dev/en/stable/install.html nearly all array operations will
take place on the GPU if the provided input is a CuPy array, as shown below.
NumPy arrays can be passed if one wishes to still perform operations on the
CPU.

⬇

import cupy as cp

x_cp = cp.array(x) # convert existing ‘numpy‘ array to ‘cupy‘ array

# apply functions like before, array operations take place on GPU

x_FS = pyffs.ffs(x_cp, T, T_c, N_FS) # compute FS coefficients

x_r = pyffs.iffs(x_FS, T, T_c, N_FS) # back to samples

y = pyffs.convolve(x_cp, x_cp, T, T_c, N_FS) # convolve

x_interp = pyffs.fs_interp(x_FS, T, a, b, M) # interpolate

Listing 6: GPU support through CuPy.

Note that converting between CuPy and NumPy requires data transfer between the
CPU and GPU, which could be costly for large arrays. Therefore, if passing
CuPy arrays to pyFFS, it is recommended to perform as much pre-processing and
post-processing as possible on the GPU in order to limit such data transfer.

###  3.5 Summary

Table 1 compares equivalent functions between pyFFS and SciPy. For a more
extensive documentation and the latest information, we refer to
pyffs.readthedocs.io, and example scripts can be found in the examples folder
of the repository: github.com/imagingofthings/pyFFS/tree/master/examples.

| pyFFS | SciPy
---|---|---
1-D Fourier analysis | pyffs.ffs | scipy.fft.fft
1-D Fourier synthesis | pyffs.iffs | scipy.fft.ifft
N-D Fourier analysis | pyffs.ffsn | scipy.fft.fftn
N-D Fourier synthesis | pyffs.iffsn | scipy.fft.ifftn
N-D convolution | pyffs.convolve | scipy.signal.fftconvolve
1-D bandlimited interpolation | pyffs.fs_interp | scipy.signal.resample
N-D bandlimited interpolation | pyffs.fs_interpn | -

Table 1: Functionality comparison between pyFFS and SciPy [22]. SciPy’s
convolution zero-pads the inputs in order to approximate a linear convolution,
while pyFFS performs a circular convolution. Within SciPy, circular
convolution is only supported for 2-D by calling scipy.signal.convolve2d with
the parameter boundary=‘wrap’. For N-D bandlimited interpolation with SciPy,
it is possible to use scipy.signal.resample along each dimension. However,
there is no one-shot function.

##  4 Benchmarking

Computational efficiency is a primary objective for pyFFS. In this section, we
present several benchmarking results to compare the computational speed
between pyFFS and SciPy for convolution and interpolation, and to demonstrate
the benefits of GPU acceleration.666The scripts to reproduce these results can
be found in the in the profile folder of the repository:
github.com/imagingofthings/pyFFS/tree/master/profile. All benchmarking is
performed on a Lenovo ThinkPad P15 Gen 1 laptop, with an Intel i7-10850H six-
core processor and an NVIDIA Quadro RTX 3000 GPU (when applicable).

###  4.1 Convolution

Before presenting the benchmarking results, we show a toy example of 1-D
convolution, to compare the outputs of pyFFS and SciPy.

(a)

(b)

Figure 3: (a) Comparing 1-D convolution between pyFFS and SciPy. pyFFS
performs a circular convolution as is expected by Fourier series coefficients,
while SciPy’s fftconvolve approximates a linear convolution by zero-padding
DFT coefficients. (b) Profiling 2-D bandlimited circular convolution between
pyFFS and SciPy.

Both DFT and FS analysis assume periodic functions, discrete and continuous
ones respectively [21]. Moreover, the convolution of two periodic functions of
the same period results in a circular convolution, such that the output is
also periodic. In certain scenarios, one may be interested in the convolution
of two signals that are not necessarily periodic, e.g. a speech recording and
a room impulse response. By zero-padding the two inputs, a linear convolution
can be computed from the circular convolution of two discrete inputs. This is
precisely what the convolve function of SciPy does. However, SciPy offers no
function for the 1-D circular convolution.777It can be manually done by taking
the inverse DFT of the product of the DFTs of the two inputs. pyFFS’
convolution function performs the circular convolution, as expected by FS
convolution theory. This explains the noticeable difference between the two
approaches in Figure 3(a), in particular at the leftmost boundary. The SciPy
approach tapers off due to zero-padding as it approximates a linear
convolution, while pyFFS performs a faithful circular convolution.

#### Benchmark

As SciPy’s function for 1-D convolution yields a different result, we do not
benchmark pyFFS against it. For 2-D convolution, SciPy has a function
scipy.signal.convolve2d that can perform a circular convolution when the
argument boundary=‘wrap’. In Figure 3(b) we compare pyffs.convolve2d against
scipy.signal.convolve2d, as these two lead to the same output. As the number
of samples per dimension grows, using pyFFS is noticeably faster than the
equivalent function in SciPy. At around $100$ samples per dimension, pyFFS is
already two orders of magnitude faster than SciPy.

###  4.2 Interpolation

Likewise, before presenting the benchmarking results, we show toy examples of
interpolation in 1-D and 2-D, comparing Theorem 2.3 and interpolation by zero-
padding DFT coefficients.

Figures 4(a) and 4(b) show the interpolation of a section of a 1-D Dirichlet
function

(12)
$\displaystyle\phi(t)=\sum_{k=-N}^{N}\exp\Big{(}j\frac{2\pi}{T}k(t-T_{c})\Big{)},$

whose bandwidth is given by $N_{\text{FS}}=2N+1$. Theorem 2.3, which makes use
of pyffs.fs_interp, and interpolating zero-padded DFT coefficients, which
employs scipy.signal.resample, perfectly match the ground truth function, as
shown in Figure 4(b). This is expected as both approaches are bandlimited
interpolation techniques and the target function is bandlimited.

(a) 1-D Dirichlet.

(b) Interpolated section of (a).

(c) 2-D Dirichlet kernel.

(d) Interpolated region of (c).

(e) Cross-section at $y=0.3$.

Figure 4: (Top) Bandlimited interpolation of a 1-D Dirichlet kernel, Equation
12 with $N_{\text{FS}}=51$, $T=1$, $T_{c}=0$, and 64 samples. (Bottom)
Bandlimited interpolation of a 2-D Dirichlet kernel, Equation 13 with
$N_{\text{FS},x}=N_{\text{FS},y}=31$, $T_{x}=T_{y}=1$, $T_{c,x}=T_{c,y}=0$,
and 256 samples along each dimension. The dashed box in (c) indicates the 2-D
interpolated section shown in (d). The dashed line in (d) indicates the cross-
section shown in (e).

Figures 4(c), 4(d), and 4(e) show the interpolation of a section of a 2-D
Dirichlet function, which is essentially the outer product of two 1-D
Dirichlet functions along the $x$\- and $y$\- dimensions

(13) $\displaystyle\phi(x,y)$
$\displaystyle=\sum_{k_{x}=-N_{x}}^{N_{x}}\sum_{k_{y}=-N_{y}}^{N_{y}}\exp\left(j\frac{2\pi}{T_{x}}k_{x}(x-T_{c,x})\right)\exp\left(j\frac{2\pi}{T_{y}}k_{y}(y-T_{c,y})\right),$

whose bandwidth is given by $N_{\text{FS},x}=2N_{x}+1$ and
$N_{\text{FS},y}=2N_{y}+1$ in the $x$\- and $y$\- dimensions respectively. We
can draw similar conclusions as for the 1-D case: the proposed technique and
interpolating zero-padded DFT coefficients perfectly match the ground truth
function as it is bandlimited.

#### Benchmark

(a)

(b)

(c)

(d)

Figure 5: (Top) Benchmarking 1-D bandlimited interpolation. $128$ samples are
taken of a signal that has bandwidth $N_{\text{FS}}=127$. Each point in the
curves is averaged over $10$ trials. $10\%$ of the period is equivalent to the
dotted black region in LABEL:{fig:interp1d_input}. (Bottom) Benchmarking 2-D
bandlimited interpolation. $256\times 256$ samples are taken of a signal that
has bandwidth $N_{\text{FS}}=[255\times 255]$. Each point in the curves is
averaged over $10$ trials. $2\%$ of the period is equivalent to the dashed red
box in LABEL:{fig:interp_2d_input}.

As a reminder, the zero-padding technique results in a resampling across the
entire period, whereas the proposed technique yields finer resolution only in
the selected region through use of the CZT. As a result, the computational
load of both approaches can be very different with respect to the width of the
selected region and the number of interpolation points in that region.

In Figures 5(a) and 5(b), we compare the two bandlimited interpolation
techniques in 1-D. As the number of interpolation points increases as shown in
Figure 5(a), it becomes slightly advantageous to use the CZT-based approach
rather than interpolation by zero-padding DFT coefficients. The clear
advantage of the CZT approach is when we focus in on smaller and smaller
regions as shown in Figure 5(b), as this necessitates a very large IDFT for
the zero-padding approach. Moreover, for a fixed number of interpolation
points, the computational load of the CZT approach is independent of the
interpolation region size. As shown in Figure 5(b), this is beneficial when we
interpolate very small regions but detrimental when interpolating over larger
sections. In the latter scenario, it is better to use the traditional zero-
padding approach.

In Figures 5(c) and 5(d), we compare the two bandlimited interpolation
techniques in 2-D. For a modest-sized input ($256\times 256$ samples), it is
noticeably advantageous to use the CZT-based approach as the number of
interpolation points increases (Figure 5(c)), and very beneficial when we
focus in on smaller regions (Figure 5(d)). The latter result makes the CZT-
based approach an attractive choice when desiring to zoom into bandlimited
images in a distortionless fashion.

###  4.3 GPU acceleration

(a) 1-D

(b) 2-D

(c)

(d)

(e)

(f)

Figure 6: Profiling GPU acceleration; each point is averaged over $10$ trials.
(Top) Fourier series computation; 1-D (middle) and 2-D (bottom) Fourier series
interpolation.

We now quantify the speed-up provided by a GPU for FS computation and
interpolation. There are two important considerations when using a GPU.
Firstly, if the application permits, it is recommended to work with float32 /
complex64 arrays for less memory consumption and potentially faster
computation. By default, NumPy and CuPy create float64 / complex128 arrays,
e.g. when initializing an array with np.zeros, so casting the arrays
accordingly is recommended. In the benchmarking tests below, we use float32 /
complex64 arrays. Secondly, the benefits of using a GPU typically emerge when
the processed arrays are larger than the CPU cache. So the crossover between
CPU and GPU performance can be very hardware dependent.

Figures 6(a) and 6(b) compares the processing time between a CPU and a GPU for
computing an increasing number of FS coefficients. In 1-D, for more than
$1^{\prime}000$ coefficients it starts to become beneficial to use a GPU, and
at around $10^{\prime}000$ coefficients it is an order of magnitude faster to
use a GPU. In 2-D, the crossover point is at around $100$ coefficients per
dimension, and at around $1^{\prime}000$ coefficients per dimension it is more
than an order of magnitude faster to use a GPU. From the 1-D and 2-D cases, it
is clear that using a GPU scales well as the input increases in size. When
considering a 2-D or even a 3-D object, where input sizes quickly grow, it is
attractive to make use of a GPU for even modest input sizes.

Figures 6(c) and 6(d) profiles the processing time for 1-D FS interpolation.
The following three approaches are compared:

  * •

Directly evaluating the bandlimited Fourier synthesis expression Eq. 3 at each
timestamp on a CPU.

  * •

Applying Theorem 2.3 on a CPU.

  * •

Applying Theorem 2.3 on a GPU.

As we vary both the number of FS coefficients and the number of interpolation
points, using Theorem 2.3 greatly reduces computational cost, as we observe an
order of magnitude reduction for the interpolation of $300$ FS coefficients.
The difference between “direct” and “CPU” is essentially the gains we get from
the FFT algorithm, namely $\mathcal{O}(N\log N)$ complexity instead of
$\mathcal{O}(N^{3})$. Using a GPU becomes more attractive as the number of
coefficients and number of samples exceeds $300$. As mentioned earlier, this
is probably the point when the arrays and computation can no longer fit on the
CPU cache.

Figures 6(e) and 6(f) compares the CPU and GPU approaches in 2-D function.
Using a GPU consistently provides two orders of magnitude faster computation
for a varying number of FS coefficients and varying number of interpolation
points per dimension. The direct approach is not even considered as it is much
too slow. The benefits of using a GPU are even more prominent in 2-D as input
sizes quickly grow when considering multidimensional scenarios.

At the time of writing, CuPy has not implemented an equivalent of SciPy’s
resample function to perform interpolation comparisons as in Figure 5.

##  5 Example application in Fourier optics

(a)

(b)

(c)

Figure 7: Visualization of optical wave propagation setup. (7(a)) When
simulating optical wave propagation for holography, one often considers the
propagation along the $z$-axis between two parallel planes, one being the
source - $f(x,y)$ \- and the other being the target plane - $g(x,y)$. (7(b))
In practice we have a finite input and output region, which determines the
maximum angle and therefore maximum spatial frequency we can observe. Focusing
on a single axis $x$, this maximum frequency is given by
$\sin\theta_{\text{max}}/\lambda$, where $\lambda$ is the optical wavelength.
(7(c)) Lower resolution source wavefront for holographic tiling.

In Fourier optics, we are often interested in the propagation of light between
two planes, i.e. a source plane and a target plane as shown in Figure 7(a).
Given an aperture function or phase pattern at the source plane, we would like
to determine the pattern at the target plane, as predicted by the Rayleigh-
Sommerfeld diffraction formula. This propagation is often modeled with one of
three approaches that make use of the FFT for an efficient simulation:
Fraunhofer approximation, Fresnel approximation, or the angular spectrum
method [10]. The choice between these three approaches typically depends on
the requirements of the application, e.g. the distance between the two planes
and the size of input and output regions [18]. For all approaches, we again
find ourselves with a continuous-domain phenomenon that can be considered
bandlimited and periodic. Bandlimited as in practice we consider finite input
and output regions, lending to a restricted set of angles and therefore a
bandlimited spatial frequency response between the source and target planes.
This restriction of angles is shown in Figure 7(b). Even though our input may
not be bandlimited, the resulting output is bandlimited after convolution with
such a response [12]. Finally, we can frame the optical simulation as periodic
as the input and output regions have a compact support and can thus be
replicated to form periodic signals.

The application of the CZT, or equivalently the fractional FT, for
interpolation has already found its use in Fourier optics to resample the
output plane outside of the grid defined by the FFT [13, 14, 23], as
demonstrated with pyFFS in Figs. 4(c) and 4(d).

Below we show how the pyFFS interface can be used in optical wave propagation
for efficient simulation and interpolation.

⬇

# pad input and reorder

f_pad = numpy.pad(f, pad_width=pad_width)

f_pad_reorder = pyffs.ffs_shift(f_pad)

# compute FS coefficients of input

F = pyffs.ffsn(f_pad_reorder, T, T_c, N_FS)

# convolution in frequency domain with free space transfer function

G = F * H

# interpolate at the desired location and resolution

# a and b specify the region while N_out specifies the resolution

g = pyffs.fs_interpn(G, T, a, b, N_out)

Listing 7: Optical free space propagation with pyFFS.

The free space propagation transfer function H in the above code listing can
be obtained by evaluating the analytic expression for the Fresnel
approximation or the angular spectrum method transfer functions at the
appropriate frequency values [10], or by measuring this response and computing
its FS coefficients with pyffs.ffsn.

One may wish to simulate an output window with the same size as the input but
at a finer resolution. In order to circumvent the much larger FFT that this
may require, an approach known as rectangular tiling [13], as shown in Figure
7(c), can be used to split the output window into tiles. In its original
proposition, the tiles were simulated sequentially, but with a GPU they could
be computed in parallel for a significantly shorter simulation time: pyFFS’s
GPU support enables this possibility. Moreover, rectangular tiling in its
original proposition requires that each tile has the same number of samples as
the input window. This restriction is removed by the interpolation approach of
pyFFS.

##  6 Conclusion

In this paper we have presented pyFFS, a Python library for efficient Fourier
series (FS) coefficient computation, convolution, and interpolation. The
intended use of this package is when working with discrete samples that arise
from a continuous-domain signal. When the underlying signal is periodic (or
has finite support and can be periodized) and bandlimited, its FS coefficients
can be computed and interpolated in a straightforward and distortionless
fashion with pyFFS. If either periodicity or bandlimitedness is not met, the
same workarounds as when applying the discrete Fourier transform can be used,
namely windowing to taper discontinuous boundaries or bandlimiting by FS
coefficient truncation.

As computation is posed in the continuous-domain, accuracy loss that may arise
from switching between the discrete- and the continuous-domain can be
minimized. Moreover, this package serves as a handy continuous-domain
complement to the functionalities already available in SciPy [22]. We also
provide functionality not available in SciPy, namely N-D circular convolution,
N-D bandlimited interpolation, and a bandlimited interpolation technique based
on the chirp Z-transform. As shown in our benchmarking results, the latter can
be more than an order of magnitude faster when interpolating sub-regions of a
1-D or 2-D periodic function. Similar results can be expected for a general
N-D function. Furthermore, GPU support has been seamlessly integrated through
the CuPy package [15], offering more than an order of magnitude reduction when
computing and interpolating a large number of FS coefficients.

In summary, pyFFS offers researchers and engineers a convenient and efficient
interface for working with FS coefficients. The source code is made available
on GitHub888github.com/imagingofthings/pyFFS and can be easily installed for
Python through PyPi.999pip install pyffs More extensive and up-to-date
documentation can be found at pyffs.readthedocs.io.

##  Appendix A Fast Fourier series computation for even-length sequences and
proofs

Theorem 2.1 addresses the fast Fourier series computation for odd-length
sequences. For even-length sequences, there is a slight modification in the
timestamps and modulation terms.

######  Theorem A.1 (Fast Fourier series, $N_{\text{s}}\in 2\mathbb{N}$).

Let $x:\mathbb{R}\to\mathbb{C}$ be a $T$-periodic function of bandwidth
$N_{\text{FS}}=2N+1$, with $T_{c}\in\mathbb{R}$ the mid-point of any period.
Let $Q\in 2\mathbb{N}+1$ be an arbitrary odd integer such that
$N_{\text{s}}=N_{\text{FS}}+Q$. Then

(14)
$\displaystyle{\bf{x}}=N_{\text{s}}\operatorname{IDFT}_{N_{s}}\left({\bf{X}}^{\text{FS}}\odot
B_{1}^{{\bf{E}}_{1}}\right)\odot B_{2}^{N{\bf{E}}_{2}},$ (15)
$\displaystyle{\bf{X}}^{\text{FS}}=\frac{1}{N_{\text{s}}}\operatorname{DFT}_{N_{s}}\left({\bf{x}}\odot
B_{2}^{-N{\bf{E}}_{2}}\right)\odot B_{1}^{-{\bf{E}}_{1}},$

where

$\displaystyle{\bf{x}}=\left[x(t_{0}),\ldots,x(t_{M-1}),x(t_{-M}),\ldots,x(t_{-1})\right]\in\mathbb{C}^{N_{\text{s}}},$
$\displaystyle{\bf{X}}^{\text{FS}}=\left[X_{-N}^{\text{FS}},\ldots,X_{N}^{\text{FS}},{\bf{0}}_{Q}\right]\in\mathbb{C}^{N_{\text{s}}},$
(16) $\displaystyle
t_{n}=T_{c}+\frac{T}{N_{\text{s}}}\left(\frac{1}{2}+n\right),\quad
n\in\mathbb{Z},$ $\displaystyle M=N_{\text{s}}/2,$

and

$\displaystyle
B_{1}=\exp\left(j\frac{2\pi}{T}\left[T_{c}+\frac{T}{2N_{\text{s}}}\right]\right)\in\mathbb{C},$
$\displaystyle{\bf{E}}_{1}=\left[-N,\ldots,N,{\bf{0}}_{Q}\right]\in\mathbb{Z}^{N_{\text{s}}},$

$\displaystyle
B_{2}=\exp\left(-j\frac{2\pi}{N_{\text{s}}}\right)\in\mathbb{C},$
$\displaystyle{\bf{E}}_{2}=\left[0,\ldots,M-1,-M,\ldots,-1\right]\in\mathbb{Z}^{N_{\text{s}}}.$

######  Proof A.2.

Theorems 2.1 and A.1

Starting with the Fourier series (FS) synthesis expression for a bandlimited
signal Eq. 3, we plug in $t_{n}$ from Eqs. 9 and 16

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}kt_{n}\right)$
$\displaystyle=\begin{dcases}\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}k\left[T_{c}+\frac{T}{N_{\text{s}}}n\right]\right),\quad\text{\lx@cref{creftype~refnum}{thm:ffs_odd}}\\\
\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}k\left[T_{c}+\frac{T}{N_{\text{s}}}\left(\frac{1}{2}+n\right)\right]\right),\quad\text{\lx@cref{creftype~refnum}{thm:ffs_even}}.\end{dcases}$

With $B_{1}=\exp\left(j\frac{2\pi}{T}T_{c}\right)$ for Theorem 2.1 and
$B_{1}=\exp\left(j\frac{2\pi}{T}\left[T_{c}+\frac{T}{2N_{\text{s}}}\right]\right)$
for Theorem A.1, and $B_{2}=\exp\left(-j\frac{2\pi}{N_{\text{s}}}\right)$ for
both cases, we can write

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=-N}^{N}X_{k}^{\text{FS}}B_{1}^{k}B_{2}^{-nk}.$

We can then shift the summation terms so they are similar to that of an
inverse discrete Fourier transform (IDFT), i.e. summation starting at $k=0$

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=0}^{2N}X_{k-N}^{\text{FS}}B_{1}^{k-N}B_{2}^{-n(k-N)}=B_{2}^{nN}\sum_{k=0}^{N_{\text{s}}-1}X_{k-N}^{\text{FS}}B_{1}^{k-N}B_{2}^{-nk},$

where $N_{\text{s}}=2N+1$.

We now introduce notation to index the $k$-th element of a vector as
$x_{k}=[{\bf{x}}]_{k}$ in order to write

$\displaystyle x(t_{n})$
$\displaystyle=B_{2}^{nN}\sum_{k=0}^{N_{\text{s}}-1}\left[{\bf{X}}^{\text{FS}}\odot
B_{1}^{{\bf{E}}_{1}}\right]_{k}B_{2}^{-nk},$

where ${\bf{E}}_{1}$ is as defined in Theorems 2.1 and A.1, and
${\bf{X}}^{\text{FS}}$ contains the ordered FS coefficients of $x(t)$ within
$[-N,N]$.

As $B_{2}=\exp\left(-j\frac{2\pi}{N_{\text{s}}}\right)$, we can write the
above expression as an IDFT of modulated FS coefficients

$\displaystyle
x(t_{n})=N_{\text{s}}\left[\operatorname{IDFT}_{N_{\text{s}}}({\bf{X}}^{\text{FS}}\odot
B_{1}^{{\bf{E}}_{1}})\right]_{n}B_{2}^{nN}.$

Using the periodicity of $x(t)$, namely $x(t_{n})=x(t_{n+qN_{s}})$ for
$q\in\mathbb{Z}$, we can write

Theorem 2.1 $\displaystyle{\bf{x}}$
$\displaystyle=\left[x(t_{0}),\ldots,x(t_{M}),x(t_{M+1}),\ldots,x(t_{2M})\right]$
$\displaystyle=\left[x(t_{0}),\ldots,x(t_{M}),x(t_{-M}),\ldots,x(t_{-1})\right].$
Theorem A.1 $\displaystyle{\bf{x}}$
$\displaystyle=\left[x(t_{0}),\ldots,x(t_{M-1}),x(t_{M}),\ldots,x(t_{2M-1})\right]$
$\displaystyle=\left[x(t_{0}),\ldots,x(t_{M-1}),x(t_{-M}),\ldots,x(t_{-1})\right].$

We can write this set of reordered samples as the IDFT of modulated FS
coefficients

$\displaystyle{\bf{x}}=N_{\text{s}}\operatorname{IDFT}_{N_{\text{s}}}({\bf{X}}^{\text{FS}}\odot
B_{1}^{{\bf{E}}_{1}})B_{2}^{N{\bf{E}}_{2}},$

where ${\bf{E}}_{2},$ as defined in Theorems 2.1 and A.1, encapsulates this
reordering of samples. We have therefore shown the relations Eqs. 6 and 14 for
Theorems 2.1 and A.1 respectively. In order to obtain Eqs. 7 and 15, which
express the ordered FS coefficients as a function of the reordered samples, we
simply have to invert Eqs. 6 and 14 respectively.

##  Appendix B Fast Fourier series interpolation proof

######  Proof B.1.

Theorem 2.3

Starting with the Fourier series (FS) synthesis expression for a bandlimited
signal Eq. 3, we plug in $t_{n}$ from Theorem 2.3

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}kt_{n}\right)=\sum_{k=-N}^{N}X_{k}^{\text{FS}}\exp\left(j\frac{2\pi}{T}k\Big{(}a+\frac{b-a}{M-1}n\Big{)}\right).$

Using $A=\exp\left(-j\frac{2\pi}{T}a\right)$ and
$W=\exp\left(j\frac{2\pi}{T}\frac{b-a}{M-1}\right)$ from Theorem 2.3, we can
write

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=-N}^{N}X_{k}^{\text{FS}}A^{-k}W^{nk}.$

We can then shift the summation terms so they are similar to that of the chirp
Z-transform (CZT) from Definition 2.2, i.e. summation starting at $k=0$

$\displaystyle x(t_{n})$
$\displaystyle=\sum_{k=0}^{2N}X_{k-N}^{\text{FS}}A^{-(k-N)}W^{n(k-N)}=A^{N}W^{-Nn}\sum_{k=0}^{N_{\text{FS}}-1}X_{k-N}^{\text{FS}}A^{-k}W^{nk},$

where $N_{\text{FS}}=2N+1$.

We can then use Definition 2.2 of the CZT to writex

$\displaystyle
x(t_{n})=A^{N}W^{-Nn}\left[\operatorname{CZT}_{N_{\text{FS}}}^{M}({\bf{X}}^{FS})\right]_{n}.$

Rearranging $x(t_{n})$ into vector form concludes the proof.

## References

  * [1] D. H. Bailey and P. N. Swarztrauber, Fractional Fourier transform and applications, SIAM Rev., 33 (1991), pp. 389–404, https://doi.org/10.1137/1033097.
  * [2] A. H. Barnett, J. Magland, and L. af Klinteberg, A parallel nonuniform fast fourier transform library based on an “exponential of semicircle” kernel, SIAM Journal on Scientific Computing, 41 (2019), pp. C479–C504.
  * [3] T. Blu, P. L. Dragotti, M. Vetterli, P. Marziliano, and L. Coulot, Sparse sampling of signal innovations, IEEE Signal Process. Mag., 25 (2008), pp. 31–40, https://doi.org/10.1109/MSP.2007.914998.
  * [4] L. I. Bluestein, A linear filtering approach to the computation of discrete Fourier transform, IEEE Trans. Audio Electroacoust., 18 (1970), pp. 451–455, https://doi.org/10.1109/TAU.1970.1162132.
  * [5] J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comput., 19 (1965), p. 297, https://doi.org/10.2307/2003354.
  * [6] P. Duhamel and M. Vetterli, Fast fourier transforms: A tutorial review and a state of the art, Signal Processing, 19 (1990), pp. 259–299, https://doi.org/10.1016/0165-1684(90)90158-U.
  * [7] J. Fageot and M. Simeoni, Tv-based reconstruction of periodic functions, Inverse Problems, 36 (2020), p. 115015.
  * [8] C. F. Gauss, Theoria interpolationis methodo nova tractata, vol. Band 3, 1866.
  * [9] C. R. Giardina and P. M. Chirlian, Bounds on the truncation error of periodic signals, (1972), pp. 206–207.
  * [10] J. Goodman, Introduction to Fourier optics, 2005.
  * [11] C. R. Harris et al., Array programming with NumPy, Nature, 585 (2020), p. 357, https://doi.org/10.1038/s41586-020-2649-2, https://doi.org/10.1038/s41586-020-2649-2.
  * [12] K. Matsushima and T. Shimobaba, Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields, Optics Express, 17 (2009), p. 19662, https://doi.org/10.1364/oe.17.019662.
  * [13] R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, Shifted Fresnel diffraction for computational holography, Optics Express, 15 (2007), p. 5631, https://doi.org/10.1364/oe.15.005631.
  * [14] V. Nascov and P. C. Logofǎtu, Fast computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution, Appl. Opt., 48 (2009), pp. 4310–4319, https://doi.org/10.1364/AO.48.004310.
  * [15] R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis, CuPy: A NumPy-compatible library for NVIDIA GPU calculations, tech. report, 2017, http://learningsys.org/nips17/assets/papers/paper_16.pdf.
  * [16] P. Prandoni and M. Vetterli, Signal processing for communications, EPFL Press, aug 2008, https://doi.org/10.1201/9781439808009, https://www.taylorfrancis.com/books/9781439808009.
  * [17] L. Rabiner, R. W. Schafer, and C. Rader, The chirp z-transform algorithm, IEEE transactions on audio and electroacoustics, 17 (1969), pp. 86–92.
  * [18] J. D. Schmidt, Numerical simulation of optical wave propagation with examples in MATLAB, SPIE, jul 2010, https://doi.org/10.1117/3.866274, https://spiedigitallibrary.org/ebooks/PM/Numerical-Simulation-of-Optical-Wave-Propagation-with-Examples-in-MATLAB/eISBN-9780819483270/10.1117/3.866274.
  * [19] C. E. Shannon, Communication theory in the presence of noise, Proc. IRE, 37 (1949), pp. 10–21.
  * [20] P. Thévenaz, T. Blu, and M. Unser, Image interpolation and resampling, Handb. Med. Image Process. Anal., (2009), pp. 465–493, https://doi.org/10.1016/B978-012373904-9.50037-4.
  * [21] M. Vetterli, J. Kovačević, and V. K. Goyal, Foundations of signal processing, Cambridge University Press, 2014.
  * [22] P. Virtanen et al., SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nat. Methods, 17 (2020), pp. 261–272, https://doi.org/10.1038/s41592-019-0686-2, https://doi.org/10.1038/s41592-019-0686-2, https://arxiv.org/abs/1907.10121.
  * [23] X. Yu, T. Xiahui, Q. Yingxiong, P. Hao, and W. Wei, Band-limited angular spectrum numerical propagation method with selective scaling of observation window size and sample number, J. Opt. Soc. Am. A, 29 (2012), p. 2415, https://doi.org/10.1364/josaa.29.002415.


Primary Cosmic Ray Energy Spectrum and Mean Mass Composition Using Data from
the TAIGA Astrophysical Complex

V. Prosin1$\star$ I. Astapov6 P. Bezyazeekov2 E. Bonvech2 A. Borodin7 A.
Bulan1 A. Chiavassa4 D. Chernov1 A. Dyachok2 A. Gafarov2 A. Garmash11,9 V.
Grebenyuk7,8 O. Gress2 E. Gress2 T. Gress2 A. Grinyuk7 O. Grishin2 A. D.
Ivanova2 A. L. Ivanova9,2 N. Kalmykov1 V. Kindin6 S. Kiryuhin2 R. Kokoulin6 K.
Komponiets6 E. Korosteleva1 V. Kozhin1 E. Kravchenko9,11 A. Kryukov1 L.
Kuzmichev1 A. Lagutin10 M. Lavrova7 Y. Lemeshev2 B. Lubsandorzhiev3 N.
Lubsandorzhiev1 A. Lukanov3 D. Lukyantsev2 S. Malakhov2 R. Mirgazov2 R.
Monkhoev2 E. Okuneva6 E. Osipova1 A. Pakhorukov2 A. Pan7 L. Panasenko11 L.
Pankov2 A. D. Panov1 A. Petrukhin6 I. Poddubny2 D. Podgrudkov1 V. Poleschuk2
V. Ponomareva2 E. Popova1 E. Postnikov1 V. Ptuskin5 A. Pushnin2 R. Raikin10 A.
Razumov1 G. Rubtsov3 E. Ryabov2 Y. Sagan7,8 V. Samoliga2 A. Silaev1 A.
Silaev(junior)1 A. Sidorenkov3 A. Skurikhin1 A. Sokolov9,11 L. Sveshnikova1 V.
Tabolenko2 A. Tanaev2 B. Tarashchansky2 M. Y. Ternovoy2 L. Tkachev7 R. Togoo12
N. Ushakov3 A. Vaidyanathan11 P. Volchugov1 N. Volkov10 D. Voronin3 A.
Zagorodnikov2 A. Zhaglova2 D. Zhurov2,13 I. Yashin6

1 Skobeltsyn Institute of Nuclear Physics, Moscow State University, Moscow,
Russia

2 Applied Physics Institute of Irkutsk State University, Irkutsk, Russia

3 Institute for Nuclear Research of the RAS, Troytsk, Moscow, Russia

4 Dipartimento di Fisica Generale Universiteta di Torino and INFN, Turin,
Italy

5 Puskov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave
Propagation of the RAS, Troitsk, Moscow. Russia

6 National Research Nuclear University MEPhI, Moscow, Russia

7 Joint Institute for Nuclear Research, Dubna, Moscow Region, Russia

8 DUBNA University, Dubna, Moscow Region, Russia

9 Budker Institute of Nuclear Physics SB RAS, Novosibirsk, Russia

10 Altai State Univeristy, Barnaul, Russia

11 Novosibirsk State University, Novosibirsk, Russia

12 Institute of Physics and Technology Mongolian Academy of Sciences,
Ulaanbaatar, Mongolia

13 Irkutsk National Research Technical University, Irkutsk, Russia

*<EMAIL_ADDRESS>
21st International Symposium on Very High Energy Cosmic Ray Interactions
(ISVHECRI 2022)
Online, 23-27 May 2022
10.21468/SciPostPhysProc.?

## Abstract

The corrected dependence of the mean depth of the EAS maximum $X_{max}$ on the
energy was obtained from the data of the Tunka-133 array for 7 years and the
TAIGA-HiSCORE array for 2 years. The parameter $\langle\ln A\rangle$,
characterizing the mean mass compositon was derived from these results. The
differential energy spectrum of primary cosmic rays in the energy range of
$2\cdot 10^{14}$ – $2\cdot 10^{16}$ eV was reconstructed using the new
parameter $Q_{100}$ the Cherenkov light flux at the core distance 100 m.
Change of the parameter for the energy reconstuction for the TAIGA-HiSCORE
from $Q_{200}$ to $Q_{100}$ provides a decreasing energy threshold for the
spectrum to about 200 TeV.

###### Contents

  1. 1 Introduction
  2. 2 Brief description of the arrays
  3. 3 Reconstruction of the EAS parameters
  4. 4 Mean experimental depth of the EAS maximum
  5. 5 Cherenkov light flux at a core distance 100 m as a new estimator of energy
  6. 6 Experimental energy spectrum by the data of TAIGA-HiSCORE
  7. 7 Conclusions

##  1 Introduction

The energy spectrum and mass composition of primary cosmic rays are the main
characteristics that can be obtained by studying extensive air showers (EAS).
The total flux of Cherenkov light is proportional to the total energy
scattered by the shower in the atmosphere. The lateral distribution function
(LDF) of the EAS Cherenkov light reflects the position of the shower
development maximum, which in turn characterizes the mass of the primary
particle. The main aim of the present work is to achieve a lower energy
threshold by the new method of energy reconstruction, descibed below.

##  2 Brief description of the arrays

Several arrays that detect the EAS Cherenkov light were successively
constructed in the Tunka Valley. The most productive arrays were the Tunka-25
(2000 - 2005) [1], which consisted of 25 detectors with a sensitive area of
0.1 $m^{2}$ each, covering a total area of approximately 0.1 $km^{2}$, and the
Tunka-133 (2009 - 2017) [2], consisted finally of 175 detectors with a
sensitive area of 0.03 $m^{2}$, covering an area of approximately 3 $km^{2}$.
The experimetal data was accumulated over 350 clear moonless nights. The total
time of the data acquisition was 2175h. Their modern successor is the TAIGA-
HiSCORE [3] array, a part of the TAIGA experimental complex [4]. TAIGA-HiSCORE
single station has a sensitive area of 0.5 $m^{2}$. Every station has it’s own
trigger. The stations are merged to an EAS event in case of $\geq 3$ station
hits coincident inside a time gate of 3 $\mu$s. This work presents the TAIGA-
HiSCORE data that was obtained using 67 stations (two first clusters) during
135 clear moonless nights in the seasons of 2019 - 2020 and 2020 - 2021. The
total data acquisition time was 327 h.

Figure 1: Mean experimental $X_{max}$ |  Figure 2: Mean $\langle\ln A\rangle$
---|---

##  3 Reconstruction of the EAS parameters

The reconstruction of the EAS parameters for the Tunka-133 array is described
in [2]. The same algorithms and the fitting functions are used for the TAIGA-
HiSCORE data processing [3]. We use the ratio $P=Q(80)/Q(200)$ as a
quantitative parameter of LDF steepness. Here, $Q(R)$ is the Cherenkov light
flux at a distance, R, in meters. One has to control that there are
measurements of light flux at core distances $R_{c}\geq 200$ m and $R_{c}\geq
80$ m. The first of these conditions is applied to the events for the primary
energy $E_{0}\geq 10^{16}$ eV for Tunka-133 and $E_{0}\geq 10^{15}$ eV for the
TAIGA-HiSCORE. CORSIKA simulation [5] confirmed that the Cherenkov light LDF
steepness is determined solely by the thickness of the atmosphere between the
array and the depth of the EAS maximum ($\Delta X_{max}=X_{0}/sec\theta-
X_{max}$). Here, $X_{0}$ is the total depth of the atmosphere. The calculated
connection between $P$ and $\Delta X_{max}$, inside the limited range of
parameter $P$ from 2.5 to 9, can be fitted inside the following expression
[5]:

$\Delta X_{max}=\left\\{\begin{array}[]{lc}\\!\ 929-103\cdot
P,&\mbox{if}\hskip 14.22636ptP\leq 3.9\\\ \\!\ 882-91\cdot P,&\mbox{if}\hskip
14.22636ptP>3.9\\\ \end{array}\right.$ (1)

The standard deviation of simulated points from the fitting line for this
range is approximately 15 $g/cm^{2}$ [5].

##  4 Mean experimental depth of the EAS maximum

The above described parameter of the LDF steepness $P$ was applied to analyze
the data of both the Tunka-133 and TAIGA-HiSCORE arrays. The depth of the
maximum is found by the formula:

$X_{max}=965/sec\theta-\Delta X_{max}$ (2)

where 965 $g/cm^{2}$ is the total depth of the atmosphere at the location of
the arrays in the Tunka Valley. To obtain undistorted estimations of the depth
of the maximum, showers are selected for the zenith angle $\theta\leq
30^{\circ}$ and the energy above $10^{16}$ eV for the Tunka-133 array and
above $10^{15}$ eV for TAIGA-HiSCORE. We have 69000 events for 7 years of
operation of the Tunka-133 and 380000 events for 2 years of operation of
TAIGA-HiSCORE. The expermental results are shown in Fig. 1. The data of both
arrays, despite the difference in their geometry, agree well with each other,
providing a wide energy range from $10^{15}$ to $3\cdot 10^{17}$ eV. Our
experimental data are compared with the direct measurements of the depth of
the maximum obtained by observing the fluorescent EAS light at the Pierre
Auger Observatoy (PAO) [10] and Telescope Array (TA) [11]. A close agreement
of our data with the PAO data is observed at an energy of $\sim 3\cdot
10^{17}$ eV. All the experimental results are compared with theoretical curves
calculated using the QGSJET-II-04 model [12] for primary protons and iron
nuclei.

Fig. 2 shows the results of recalculation from the mean depth of the maximum
to the parameter $\langle\ln A\rangle$ the average logarithm of the atomic
number using the QGSJET-II-04 model. Qualitatively, the behavior of the mean
mass composition repeats what was published in our previous studies [13]
becoming heavier in the energy range of $3\cdot 10^{15}$ \- $3\cdot 10^{16}$
eV and lighter with a further increase in energy. However, the mean
composition in the entire energy range under consideration is estimated as
mostly light.

##  5 Cherenkov light flux at a core distance 100 m as a new estimator of
energy

The TAIGA-HiSCORE array structure is a square net of stations with a step of
about 100 m. So the minimal core distance for which light flux can be
reconstructed for almost all the events is about 100 m. Our previous parameter
for the energy reconstruction was light flux at a core distance 200 m
($Q_{200}$) [3]. When the EAS zenith angle $\theta$ changes from $0^{\circ}$
to $45^{\circ}$, $Q_{200}$ changes by less than 10%. Therefore, it was assumed
in [3], that $Q_{200}$ does not depend on $\theta$ for a fixed particle
energy. It was found by the new CORSIKA simulation that light flux $Q_{100}$
depends on the zenith angle $\theta$ sigificantly more, changing by about 2.5
times with the same change in $\theta$. So first one needs to recalculate from
the measured light flux to the $\theta=0^{\circ}$ using the new CORSIKA
results:

$\log_{10}(Q_{100}(0))=\log_{10}(Q_{100}(\theta))+(\sec\theta-1)\cdot(1.25-0.083\cdot\log_{10}(Q_{100}(\theta))$
(3)

Then $Q_{100}(0)$ can be recalculated to the primary energy $E_{0}$ using the
result of the new CORSIKA simulation:

$\log_{10}(E_{0}/GeV)=0.88\cdot\log_{10}(Q_{100}(0))+5.14$ (4)

Figure 3: Differential primary CR spectrum by TAIGA-HiSCORE data

##  6 Experimental energy spectrum by the data of TAIGA-HiSCORE

The experimental energy estimation differs from that described at the previous
section because the real atmosphere light absorbtion is different from night
to night in contradiction with standard absorption assumed in CORSIKA
simulations. So first we obtain the integral energy spectrum for the single
night using the expression (4). Then we normalize this spectrum to the
reference energy spectrum measured by the QUEST experiment [6]. The mean
difference of normalization constant from that in expression (4) is 0.03. The
differential energy spectrum obtained from the data of TAIGA-HiSCORE array is
shown in Fig. 3. Efficiency of the events at the first left point (starting
from the energy $2\cdot 10^{14}$ eV) is more than 95%. Points for the lower
energy obtained from the events with lower efficiency are removed. The low
energy points of our spectrum are in good agreement with direct balloon [7],
satellite [8] and mountain [9] measurements.

##  7 Conclusions

The new estimations of $X_{max}$, derived from the steepness parameter
$P=Q(80)/Q(200)$ provides good agreement between both the results of our
arrays Tunka-133, TAIGA-HiSCORE and Tunka-133 and the direct measurements of
$X_{max}$ at the Pierre Auger Observatory (PAO).

The primary composition, derived from $X_{max}$ is lighter than it seemed in
our previous publications. It is mostly light (p+He) over the whole energy
range.

The observed increase of $\langle\ln A\rangle$ in the energy range $10^{16}$
\- $10^{17}$ eV demands a new theoretical explanation.

Change of the parameter for the energy reconstruction for the TAIGA-HiSCORE
from $Q_{200}$ to $Q_{100}$ provides a decreasing energy threshold for the
spectrum to about 200 TeV.

The all particle energy spectrum over the energy range 200 TeV - 3 PeV follows
a pure power law with index $2.71\pm 0.01$.

## Acknowledgements

The work was performed at the UNU “Astrophysical Complex of MSU-ISU”
(agreement EB 075-15-2021-675). The work is supported by RFBR (grants
19-52-44002, 19-32-60003), the RSF(grants 19-72-20067(Section 2)), the Russian
Federation Ministry of Science and High Education (projects FZZE-2020-0017,
FZZE-2020-0024, and FSUS-2020-0039).

## References

  * [1] N. Budnev et al, Tunka-25 Air Shower Cherenkov array The main results, Astropart. Phys. 50, 18 (2013), 10.1016/j.astropartphys.2013.09.006.
  * [2] N. M. Budnev et al, The primary cosmic-ray energy spectrum measured with the Tunka-133 array, Astropart. Phys. 117, 102406 (2020), 10.1016/j.astropartphys.2019.102406.
  * [3] V. Prosin et al, Cosmic Ray Energy Spectrum derived from the Data of EAS Cherenkov Light Arrays in the Tunka Valley, EPJ Web of Conf. 210, 01003 (2019), 10.1051/epjconf/201921001003.
  * [4] L. Kuzmichev et al, Experimental Complex TAIGA, Phys. At. Nucl. 83, 1375 (2020), 10.1134/S1063778820090161.
  * [5] V. Prosin et al, Energy Spectrum and Mass Composition of Cosmic Rays from the Data of the Astrophysical Complex TAIGA, Physics of Atomic Nulei, 84, 9, 1563 (2021), 10.1134/S1063778821090283.
  * [6] E. E. Korosteleva et al, Measurement of cosmic ray primary energy with the atmospheric Cherenkov light technique in extensive air showers, Nucl. Phys. B Proc. Suppl. 165, 74 (2007), 10.1016/j.nuclphysbps.2006.11.012.
  * [7] A. D. Panov et al, Energy Spectra of Abundant Nuclei of Primary Cosmic Rays from the Data of ATIC-2 Experiment Final Results, Bull. Russ. Acad. Sci.: Phys. 73, 564 (2009), 10.3103/S1062873809050098.
  * [8] A. N. Turundaevskiy et al, Main Results from the NUCLEON Experiment, Bull. Russ. Acad. Sci.: Phys. 85, 353 (2021), 10.3103/S1062873821040377.
  * [9] R. Alfaro et al (HAWC Collab.), All-particle cosmic ray energy spectrum measured by the HAWC experiment from 10 to 500 TeV, Phys. Rev. D: Part. Fields 96, 122001 (2017), 10.1103/PhysRevD.96.122001.
  * [10] A. Yushkov (for P. A. Collab.), Mass composition of cosmic rays with energies above 1017.2 eV from the hybrid data of the Pierre Auger Observatory, in Proceedings of the 36th ICRC, Madison, USA, PoS(ICRC2019), 482 (2019).
  * [11] R. U. Abbasi et al. (T. A. Collab.), Depth of Ultra High Energy Cosmic Ray Induced Air Shower Maxima Measured by the Telescope Array Black Rock and Long Ridge FADC Fluorescence Detectors and Surface Array in Hybrid Mode, Astrophys. J. 858, 76 (2018), 10.48550/arXiv.1801.09784.
  * [12] S. Ostapchenko and M. Bleicher, Constraining pion interactions at very high energies by cosmic ray data, Phys. Rev. D: Part. Fields 93, 051501 (2016); arXiv: 1601.06567.; 10.48550/arXiv.1601.06567. EPJ Web of Conf. 208, 02002 (2019), 10.1051/epjconf/201920802002.
  * [13] V. V. Prosin et al, Results from Tunka-133 (5 years observation) and from the Tunka-HiSCORE prototype, EPJ Web of Conf. 121, 03004 (2016), 10.1051/epjconf/201612103004.