Thomas Jefferson National Accelerator Facility, Newport News, Virginia 23606,
USA

# Insight into Emergence of Hadron Mass from $\boldmath N^{*}$
Electroexcitation Amplitudes

V.I. Mokeevins:jlab  D.S. Carmanins:jlab ins:jlabins:jlab

###### Abstract

The emergence of hadron mass represents one of the most challenging and still
open problems in contemporary hadron physics. The results on the nucleon
resonance electroexcitation amplitudes available from the CLAS data on $\pi N$
and $\pi^{+}\pi^{-}p$ electroproduction analyzed within the continuum
Schwinger method open up a new avenue for gaining insight into the strong
interaction dynamics that are responsible for the generation of the dominant
part of hadron mass. Future prospects of these studies in experiments of the
12-GeV era with CLAS12 and after a potential increase of the CEBAF energy up
to 22 GeV will offer a unique opportunity to explore the full range of
distances where the dominant part of hadron mass and $N^{*}$ structure emerge
from QCD.

##  1 Introduction

An extensive research program is now in progress in Hall B at Jefferson Lab
with the objective to determine the nucleon resonance ($N^{*}$)
electroexcitation amplitudes, also known as the $\gamma_{v}pN^{*}$
electrocouplings, for most states in the mass range $W$ up to 2.0 GeV, for
photon virtualities $Q^{2}<5$ GeV2 from the data on the exclusive meson
electroproduction channels measured with CLAS and analyzed within different
reaction models [1, 2, 3]. This information will be further extended towards
$Q^{2}$ up to 10 GeV2 from the data expected in the measurements with the
CLAS12 detector that are currently in progress [4, 5]. These studies offer
unique insight into many facets of the strong interaction dynamics in the
regime of large (comparable with unity) QCD running coupling $\alpha_{s}/\pi$
seen in the generation of $N^{*}$ states of different quantum numbers with
different structural features [6, 7, 8, 9, 10, 11, 12, 13]. Analyses of the
results on the $Q^{2}$-evolution of the $\gamma_{v}pN^{*}$ electrocouplings
within continuum Schwinger methods (CSM) that provide a connection to the QCD
Lagrangian, open up a promising new avenue in the exploration of the strong
interaction dynamics responsible for the emergence of the dominant part of
hadron mass [14, 15].

Figure 1: Concepts underlying emergent hadron mass (EHM) within the continuum
Schwinger method (CSM) [15, 16]: QCD-driven bare quark and gluon dressing
(left) and the momentum dependencies of the dressed quark and gluon masses
(right) [17].

##  2 CSM Concept for Emergent Hadron Mass

The emergence of hadron mass is the most challenging open problem in hadron
physics. This problem is revealed through the comparison between the observed
nucleon masses and the sum of the masses of their quark constituents. The sum
of the three bare $u$ and $d$ quark masses [18] that enter into the QCD
Lagrangian and generated through the Higgs mechanism accounts for $<$2% of the
measured nucleon masses. This observation suggests that the dominant part of
hadron mass is generated in processes other than the Higgs mechanism.

The QCD gauge gluons of zero mass should be dressed in the processes imposed
by the interaction part of the QCD Lagrangian shown in the bottom left panel
of Fig. EHM_CSM. The gluon self-interaction creates dynamically generated
gluon mass that emerges in the transverse part of the gluon vacuum
polarization tensor $P_{\mu\nu}$ that preserves gauge invariance [15, 16]. The
scalar function for the gluon vacuum polarization tensor $P(p^{2})$ behaves as
$\frac{m_{g}^{2}(p^{2})}{p^{2}}$, where $p^{2}$ is the squared four-momentum
of the dressed gluon. Consequently, the dressed gluon propagator
$\propto\frac{1}{p^{2}(1+P(p^{2}))}=\frac{1}{p^{2}+m_{g}^{2}(p^{2})}$,
demonstrating the emergence of non-zero running mass $m_{g}(p^{2})$ of the
dressed gluon. The gluon self-interaction also causes a rapid increase of the
QCD running coupling $\alpha_{s}/\pi$ in the transition from the perturbative
to the strongly coupled regimes. Furthermore, at gluon momenta $p<0.5$ GeV,
the coupling $\alpha_{s}/\pi$ becomes almost momentum independent owing to the
finite value of the gluon mass as $p\to 0$ [19]. The bare QCD quark should be
dressed in the processes shown in the top left panel in Fig. 1. Its propagator
behaves
$\propto\frac{1}{A(p^{2})p^{2}+B(p^{2})}=\frac{Z}{p^{2}+(B(p^{2})/A(p^{2}))}$,
demonstrating the dynamical generation of the dressed quark mass
$m_{q}^{2}=B(p^{2})/A(p^{2})$ running with $p^{2}$, the dressed quark squared
four-momentum, and the renormalization parameter $Z$ for the quark wave
function. When the QCD running coupling $\alpha_{s}/\pi$ becomes comparable
with unity ($\alpha_{s}/\pi>0.3$), the energy stored in the gluon field is
transformed into a momentum dependence of the dressed quark masses as shown in
the right panel of Fig. 1. The momentum dependencies of the dressed quark and
gluon masses have been computed within CSM as the solution of the QCD
equations of motion for the quark and gluon fields. The CSM predictions on the
dressed quark/gluon mass functions have been confirmed in lQCD studies [20].
The rapid increase seen in the momentum dependence of the dressed quark mass
when the quark momentum decreases underlies the emergence of the dominant part
of the masses of the ground and excited states of the nucleon as bound systems
of three dressed quarks.

##  3 Insight into EHM from Results on $Q^{2}$-evolution of $\gamma_{v}pN^{*}$
Electrocouplings

The ground and excited states of the proton emerge as bound systems of three
dressed quarks. In the processes of $N^{*}$ electroexcitation, the virtual
photon interacts with the dressed quarks. Since the running mass of the
dressed quark is encoded in dressed quark propagator, the information on the
$Q^{2}$ evolution of the $\gamma_{v}pN^{*}$ electrocouplings allows us to map
out the momentum dependence of the dressed quark mass.

Analyses of the CLAS data on $\pi N$, $\eta N$, and $\pi^{+}\pi^{-}p$
electroproduction have provided the only available results on the $Q^{2}$
evolution of the $\gamma_{v}pN^{*}$ electrocouplings for most $N^{*}$ states
in the mass range up to 1.8 GeV for $Q^{2}<5$ GeV2 [1, 2, 21, 22]. Consistent
results on the $\gamma_{v}pN^{*}$ electrocouplings of the resonances in the
mass range up to 1.6 GeV available from independent studies of $\pi N$ and
$\pi^{+}\pi^{-}p$ electroproduction with different non-resonant contributions
demonstrated the capabilities of the reaction models [23, 24] developed by the
CLAS Collaboration for the credible extraction of these quantities.

Figure 2: Description of the CLAS results on the $N\to\Delta$ magnetic
transition form factor (left) [23] and $N(1440)1/2^{+}$ electrocouplings [23,
24, 25] (middle) achieved within CSM [26, 27] by employing the same dressed
quark mass function deduced from the QCD Lagrangian (blue solid lines). The
CSM evaluations with momentum-independent (frozen) quark mass are shown by the
red dashed lines [26, 28]. (Right) CLAS results on the $\Delta(1600)3/2^{+}$
electrocouplings obtained from analyses of $\pi^{+}\pi^{-}p$ electroproduction
cross sections within $W$-intervals as indicated for $Q^{2}$ from 2.0 to 5.0
GeV2 [25]. The CSM predictions are shown by the red line [29].

The experimental results on the $Q^{2}$-evolution of the dominant $N\to\Delta$
magnetic transition form factor and all electrocouplings of the
$N(1440)1/2^{+}$ have been described successfully within CSM [26, 27] starting
from the QCD Lagrangian for $Q^{2}>0.8$ GeV2 for the $\Delta(1232)3/2^{+}$ and
for $Q^{2}>2.0$ GeV2 for the $N(1440)1/2^{+}$. Within these $Q^{2}$ ranges the
core of three dressed quarks becomes the biggest contributor to the structure
of these states, allowing for direct comparison between the data and the CSM
predictions that account for the contribution from only the quark core. A good
description of the electrocouplings of these resonances of different structure
(spin-isospin flip and the first radial excitation of three dressed quarks for
the $\Delta(1232)3/2^{+}$ and $N(1440)1/2^{+}$, respectively) shown in the
left and the middle panels in Fig 2, was achieved with the same dressed quark
mass function deduced from the QCD Lagrangian (see Fig. 1 right), that was
already used in the successful description of the pion and nucleon elastic
electromagnetic form factors. This success offers strong evidence for dressed
quarks with running mass as the active degrees of freedom in the structure of
the pion, ground states of the nucleon, and the $\Delta(1232)3/2^{+}$ and
$N(1440)1/2^{+}$, and demonstrates the capability for mapping out the momentum
dependence of the dressed quark mass from the results on the $Q^{2}$-evolution
of the $\gamma_{v}pN^{*}$ electrocouplings.

The electrocouplings of the $\Delta(1600)3/2^{+}$ were predicted for the
$Q^{2}$-range from 2.0 to 5.0 GeV2 within the CSM approach in 2019 with the
dressed quark mass function described above [29], when experimental results
were not available. The experimental results on the $\Delta(1600)3/2^{+}$
electrocouplings were obtained in 2023 [25] from analyses of the nine
independent one-fold differential $\pi^{+}\pi^{-}p$ electroproduction cross
sections measured with CLAS [30, 31] in three overlapping $W$-intervals for
$Q^{2}$ from 2.0 to 5 0 GeV2 (see Fig. 2 right). The $\Delta(1600)3/2^{+}$
contributes substantially into all three $W$ intervals with different non-
resonant backgrounds. The $\Delta(1600)3/2^{+}$ electrocouplings extracted
from the data in the three $W$-intervals are consistent, demonstrating
credible extraction of these quantities. These experimental results have
confirmed the CSM predictions, which solidifies evidence for insight into the
dressed quark mass function.

The $Q^{2}$-range for the $\gamma_{v}pN^{*}$ electrocouplings covered by the
meson electroproduction data from CLAS spans up to 5 GeV2, allowing us to map
out the range of quark momenta where $<$30% of hadron mass is expected to be
generated (see Fig. 1 right). CLAS12 is the only available facility capable of
extending the $Q^{2}$-coverage of the results on the $\gamma_{v}pN^{*}$
electrocouplings up to 10 GeV2 from measurements of $\pi N$, $KY$, and
$\pi^{+}\pi^{-}p$ electroproduction [1]. The foreseen CLAS12 results will
allow us to explore the range of distances where $\approx$50% of hadron mass
is generated. A further increase of the CEBAF energy up to 22 GeV [32] will
make it possible to explore the $\gamma_{v}pN^{*}$ electrocouplings over a
still broader $Q^{2}$-range. Simulation studies performed at 22 GeV with
CLAS12 suggest that luminosities up to $2-5\times 10^{35}$ cm-2s-1 would be
sufficient in order to obtain the $\gamma_{v}pN^{*}$ electrocouplings for
$Q^{2}$ up to 30 GeV2. These results will offer the only foreseen opportunity
to explore the full range of distances where the dominant part of hadron mass
and $N^{*}$ structure emerge from QCD (see Fig. 1), addressing the key open
problem of the Standard Model on the emergence of hadron mass.

##  4 Conclusions and Outlook

Nucleons and their excited states are the most fundamental three-body systems
in Nature. If we do not understand how QCD builds each state in the complete
$N^{*}$ spectrum, then our understanding of the strongly coupled QCD regime
remains incomplete. High-quality meson electroproduction data from the 6-GeV
era at JLab have allowed for the determination of the electrocouplings of most
$N^{*}$s in the mass range up to 1.8 GeV for $Q^{2}<7.5$ GeV2. A good
description of the $\Delta(1232)3/2^{+}$, $N(1440)1/2^{+}$, and
$\Delta(1600)3/2^{+}$ electrocouplings for $Q^{2}<5$ GeV2 achieved within CSM
with the same dressed quark mass function inferred from the QCD Lagrangian and
used in the successful description of the elastic pion and nucleon
electromagnetic form factors, offers sound evidence for insight into the
momentum dependence of the dressed quark mass. Extension of the results on the
$\gamma_{v}pN^{*}$ electrocouplings into the $Q^{2}$ range from 5 - 30 GeV2
from the CLAS12 measurements in the 12-GeV era and after the increase of the
CEBAF energy up to 22 GeV, will offer the only foreseen opportunity to explore
how the dominant part of hadron mass and $N^{*}$ structure emerge from QCD.

##  5 Acknowledgments

This work was supported in part by the U.S. Department of Energy (DOE) under
Contract No. DE-AC05-06OR23177 and Jefferson Science Associates (JSA).

## References

  * [1] Mokeev V. I. Carman D. S., Few Body Syst.63202259.
  * [2] Carman D. S., Joo K. Mokeev V. I., Few Body Syst.61202029.
  * [3] Aznauryan I. G. Burkert V. D., Prog. Part. Nucl. Phys.6720121.
  * [4] Brodsky S. J. et al., Int. J. Mod. Phys. E2920202030006.
  * [5] Burkert V. D., Ann. Rev. Nucl. Part. Sci.682018405.
  * [6] Burkert V. D. Roberts C. D., Rev. Mod. Phys.912019011003.
  * [7] Burkert V. D., EPJ Web Conf.241202001004.
  * [8] Giannini M. M. Santopinto E., Chin. J. Phys.532015020301.
  * [9] Aznauryan I. G. Burkert V., Phys. Rev. C952017065207.
  * [10] Aznauryan I. G. et al., Int. J. Mod. Phys. E2220131330015.
  * [11] Obukhovsky I. T., Faessler A., Fedorov D. K., Gutsche T. Lyubovitskij V. E., Phys. Rev. D842011014004.
  * [12] Obukhovsky I. T., Faessler A., Fedorov D. K., Gutsche T. Lyubovitskij V. E., Phys. Rev. D1002019094013.
  * [13] Lyubovitskij V. E. Schmidt I., Phys. Rev. D1022020094008.
  * [14] Carman D. S., Gothe R. W., Mokeev V. I. Roberts C. D., Particles62023416.
  * [15] Ding M., Roberts C. D. Schmidt S. M., Particles6202357.
  * [16] Roberts C. D., Richards D. G., Horn T. Chang L., Prog. Part. Nucl. Phys.1202021103883.
  * [17] Roberts C. D., AAPPS Bull.3120216.
  * [18] Workman R. L. et al., PTEP20222022083C01.
  * [19] Cui Z.-F., Zhang J.-L., Binosi D., de Soto F., Mezrag C., Papavassiliou J., Roberts C. D., Rodríguez-Quintero J., Segovia J. Zafeiropoulos S., Chin. Phys. C442020083102.
  * [20] Oliveira O., Silva P. J., Skullerud J.-I. Sternbeck A., Phys. Rev. D992019094506.
  * [21] Burkert V. D., Nucleon Resonances and Transition Form Factors, to be published (2024).
  * [22] Hiller Blin A. N. et al., Phys. Rev. C1002019035201.
  * [23] Aznauryan I. G. et al., Phys. Rev. C802009055203.
  * [24] Mokeev V. I. et al., Phys. Rev. C862012035203.
  * [25] Mokeev V. I., Achenbach P., Burkert V. D., Carman D. S., Gothe R. W., Hiller Blin A. N., Isupov E. L., Joo K., Neupane K. Trivedi A., Phys. Rev. C1082023025204.
  * [26] Segovia J., Cloët I. C., Roberts C. D. Schmidt S. M., Few Body Syst.5520141185.
  * [27] Segovia J., El-Bennich B., Rojas E., Cloët I. C., Roberts C. D., Xu S.-S. Zong H.-S., Phys. Rev. Lett.1152015171801.
  * [28] Wilson D. J., Cloet I. C., Chang L. Roberts C. D., Phys. Rev. C852012025205.
  * [29] Lu Y., Chen C., Cui Z.-F., Roberts C. D., Schmidt S. M., Segovia J. Zong H. S., Phys. Rev. D1002019034001.
  * [30] Isupov E. L. et al., Phys. Rev. C962017025209.
  * [31] Trivedi A., Few Body Syst.6020185.
  * [32] Accardi A. et al., Strong Interaction Physics at the Luminosity Frontier with 22 GeV Electrons at Jefferson Lab, to be published (2024).


# Ramanujan subspace pursuit for signal periodic decomposition

Deng Shi-Wen*, Han Ji-Qing Manuscript received May X, 2015; revised XX XX,
2015. *Deng Shi-Wen is with School of Mathematical Sciences, Harbin Normal
University, Harbin, China (e-mail: dengswen@gmail.com).*Han Ji-Qing is with
School of Computer Science and Technology, Harbin Institute of Technology,
Harbin, China (e-mail: jqhan@hit.edu.cn).

###### Abstract

The period estimation and periodic decomposition of a signal are the long-
standing problems in the field of signal processing and biomolecular sequence
analysis. To address such problems, we introduce the Ramanujan subspace
pursuit (RSP) based on the Ramanujan subspace. As a greedy iterative
algorithm, the RSP can uniquely decompose any signal into a sum of exactly
periodic components, by selecting and removing the most dominant periodic
component from the residual signal in each iteration. In the RSP, a novel
periodicity metric is derived based on the energy of the exactly periodic
component obtained by orthogonally projecting the residual signal into the
Ramanujan subspace, and is then used to select the most dominant periodic
component in each iteration. To reduce the computational cost of the RSP, we
also propose the fast RSP (FRSP) based on the relationship between the
periodic subspace and the Ramanujan subspace, and based on the maximum
likelihood estimation of the energy of the periodic component in the periodic
subspace. The fast RSP has a lower computational cost and can decompose a
signal of length $N$ into the sum of $K$ exactly periodic components in
$\mathcal{O}(KN\log N)$. In addition, our results show that the RSP
outperforms the current algorithms for period estimation.

###### Index Terms:

Period estimation, periodic decomposition, period detection, periodic signals,
Ramanujan subspace.

##  I Introduction

There exist many periodic patterns, such as periodicities, repetitions and
regularities, in the Biomolecular sequence [1, 2, 3, 4], speech [5], music
[6], and machine vibration signal [7]. To identify these periodic patterns is
a long-standing problem in many application domains. However, the traditional
methods, such as the discrete Fourier transform (DFT), periodogram, and
autocorrelation, cannot effectively detect these patterns, especially for the
detection of the hidden periodicities [8]. Recently, various approaches have
been proposed to perform the period estimation and the periodic decomposition
of a signal based on subspace decomposition [8, 7, 9, 10], sparse
decomposition [12, 11], and integer period discrete Fourier transform (IPDFT)
[13, 15]. Unfortunately, most of these approaches usually suffer from some
limitations and drawbacks, such as failure to accurately identify periodic
components, limitation to be detected all periods, and higher computational
cost.

By generating a set of so called ‘periodic subspaces’, the algorithms based on
subspace decomposition [8, 7, 9, 10] can model and capture the periodic
components of a signal by projecting the signal into these subspaces. Each
subspace consists of all the possible periodic signals (or components) with a
specific period, and hence it is important for these methods to first define
‘what is the periodic signal?’. With the conventional concept of the
periodicity, i.e., if $x(n)$ is the signal with the period $q$, then it
satisfies that $x(n)=x(n+q)$, Sethares and Staley [8] constructed the periodic
subspaces and proposed the periodicity transforms (PTs) [8]. The PTs can
decompose a signal into a sum of the periodic components by projecting the
signal into these periodic subspaces and updating residuals by removing the
periodicities, iteratively. However, there exist two drawbacks in the PTs. The
first one is that the periodic subspaces in the PTs are not orthogonal to each
other, which results in that a periodic signal (or component) can lie in
different subspaces. For example, the periodic signal $\mathbf{1}$ (the vector
of all ones) with the period $1$ can lie in any periodic subspace. The second
one is that the result of the PTs depends on the order in which the periodic
components are extracted. Due to these shortcomings, the PTs fails to
accurately identify the periodic components in a signal.

To eliminate the ambiguity in the definition of the periodic subspace in the
PTs, Muresan and Parks [7] first introduced the key concept of the exactly
periodic signal. Specifically, if $x(n)$ is an exactly periodic signal, then
it satisfies that $x(n)=x(n+q)$ where $q$ is the smallest integer period.
Based on the definition of the exactly periodic signal, the authors
constructed the exactly periodic subspaces (which are orthogonal to each
other) by calculating the intersection of the periodic subspaces defined in
the PTs [8], and proposed the orthogonal exactly periodic subspace
decomposition (EPSD). Moreover, when the aim is only to detect the periods in
the signal and is not interested in its exactly periodic components, it is
unnecessary to calculate the intersection of the periodic subspaces and the
orthogonal projections in them. The EPSD can perform the period estimation by
directly and efficiently computing the energy of each exactly periodic
component based on the maximum likelihood estimation (MLE) [17], which is a
significant advantage of the EPSD. However, the EPSD has the limitation that
the exactly periodic subspaces are orthogonal to each other only when their
periods are the divisors of the signal length. For example, when the signal
length is equal $12$, only the exactly periodic subspaces associated with the
periods $\\{1,2,3,4,6,12\\}$ are orthogonal to each other, but the exactly
periodic subspaces corresponding to the periods $\\{5,7,8,9,10,11\\}$ are not
orthogonal. Due to this limitation, the EPSD cannot accurately identify all
periods in a signal.

More recently, Vaidyanathan [9, 10] presented the Ramanujan subspace based on
the Ramanujan sums. As an important contribution, the Ramanujan subspace is
the first exactly periodic subspace that can be directly and definitely
constructed. Many interesting properties of the Ramanujan subspace, such as
orthogonality and dimension, are also discussed. As another important
contribution, the relationship between the exactly periodic components and
their frequencies are completely presented. All these works are important to
further explore the periodic patterns of signals. Moreover, based on the
Ramanujan subspace, Vaidyanathan also proposed the Ramanujan periodic
transform (RPT) [10] to decompose a signal into a sum of exactly periodic
components, by orthogonally projecting the signal into a set of Ramanujan
subspaces which are orthogonal to each other. Unfortunately, the RPT also
suffers from the same limitation as the ESPD [7], i.e., to construct
orthogonal Ramanujan subspaces, the signal length must be a multiple of the
least common multiplier of the expected periods.

Based on the framework of the sparse representation of a signal, Vaidyanathan
[12] and Nakashizuka [11] proposed their methods for the periodic
decomposition, respectively, with the predefined dictionary and the algorithms
of sparse decomposition such as Basis Pursuit [18] and LASSO [19]. In [11],
the dictionary is constructed based on the basis vectors defined in [8]. These
basis vectors corresponding to different periods can span a set of periodic
subspaces rather than the exactly periodic subspaces, and hence the method in
[11] cannot accurately represent the exactly periodic components in the
signal. In [12], the frequency dictionary (named Farey dictionary) is
constructed by partitioning the full bandwidth according to the Farey
sequence. In fact, the Farey dictionary is the frequency representation of the
union of all possible Ramanujan subspaces. For example, given the maximum
period $Q=6$, the frequency bins contained in the Farey dictionary are
$\left\\{\frac{0}{1},\frac{1}{6},\frac{1}{5},\frac{1}{4},\frac{1}{3},\frac{2}{5},\frac{1}{2},\frac{3}{5},\frac{2}{3},\frac{3}{4},\frac{4}{5},\frac{5}{6},\frac{1}{1}\right\\}$.
The frequency bins contained in the Ramanujan subspaces,
$\mathcal{S}_{1},\mathcal{S}_{2},\mathcal{S}_{3},\mathcal{S}_{4},\mathcal{S}_{5}$,
and $\mathcal{S}_{6}$, are $\left\\{\frac{1}{1}\right\\}$,
$\left\\{\frac{1}{2}\right\\}$, $\left\\{\frac{1}{3},\frac{2}{3}\right\\}$,
$\left\\{\frac{1}{4},\frac{3}{4}\right\\}$,
$\left\\{\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5}\right\\}$ and
$\left\\{\frac{1}{6},\frac{5}{6}\right\\}$, respectively. Obviously, the
frequencies in the Farey dictionary with the maximum period $6$ are the union
of that contained in all the Ramanujan subspaces with periods from $1$ to $6$
(note that the frequencies $0$ and $1$ are same). Therefore, the sparse
representation of a signal based on the Farey dictionary can reveal its
exactly periodic components. However, all the methods based on the sparse
decomposition suffer from the significantly higher computational costs. For
the method in [12], the size of the Farey dictionary is dramatically increased
with the maximum period $Q$. For example, when $Q=512$, the number of the
atoms in the Farey dictionary is $79852$, which requires a much higher
computational cost to implement the sparse decomposition.

By considering the limitations of the DFT for period estimation and periodic
decomposition, Epps et al. [13] presented the IPDFT by transforming a signal
according to the period rather than the frequency as done in the DFT. The
IPDFT is defined by

$\displaystyle X(q)=\sum_{n=0}^{N-1}x(n)\exp\left(-j\frac{2\pi
n}{q}\right),q=1,2,\cdots,Q\leq N$

where $Q$ is the maximum period and $N$ is the length of the signal. By
combining the IPDFT and autocorrelation of a signal, the authors further
proposed the autocorrelation-IPDFT method which was later used in [14] for
detecting periodic fragments in DNA sequence. Although the IPDFT improves the
ability of period estimation compared with the DFT, it only represents a
subset of the signals with a special period rather than the whole exactly
periodic subspace. For example, Let $Q=N=6$, the frequency bins in the DFT are
$\\{0,\frac{1}{6},\frac{2}{6},\frac{3}{6},\frac{4}{6},\frac{5}{6}\\}$ and that
in the IPDFT are
$\\{\frac{1}{1},\frac{1}{2},\frac{1}{3},\frac{1}{4},\frac{1}{5},\frac{1}{6}\\}$.
It seems that the IPDFT can represent all the integer periods from $1$ to $6$,
but it is not the case. The frequency bins contained in the exactly periodic
subspace corresponding to the period $q=5$ are
$\\{\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5}\\}$, but the IPDFT only
represents the signals in this subspace with frequency $\\{\frac{1}{5}\\}$.
Therefore, it is not possible for the IPDFT-based methods to accurately
identify the periodic components in a signal.

In this paper, we derive the Ramanujan subspace pursuit (RSP) algorithm, which
can overcome the drawbacks and limitations mentioned above, based on the
Ramanujan subspace. The RSP is a greedy iterative algorithm that can decompose
any signal into a sum of exactly periodic components, by selecting and
removing the most dominant periodic component from the residual signal in each
iteration. Unlike the PTs, the RSP is independent of the order and yields a
unique decomposition of the signal. To accurately select the most dominant
periodic component in each iteration, a novel periodicity metric in the RSP
algorithm is derived based on the energy of the periodic component obtained by
orthogonally projecting the residual signal into the each Ramanujan subspace.
To reduce the computational cost of the RSP, we also propose the fast RSP
based on the analysis of the relationship between the periodic subspace and
the exactly periodic subspace and based on the MLE of the energy of the
periodic component in the periodic subspace. The fast RSP has a lower
computational cost and can decompose a signal of length $N$ into the sum of
$K$ exactly periodic components in $\mathcal{O}(KN\log N)$.

The remainder of the paper is organized as follows. In the next section, the
concepts of exactly periodic signal and Ramanujan subspace are summarized.
Section III introduces the orthogonal projection of a signal with any length
into the Ramanujan subspace. The periodicity metric and the original RSP
algorithm is proposed in Section IV and the fast RSP is introduced in Section
V. Moreover, the evaluations of the proposed RSP algorithm are performed in
Section VI. Finally, Section VII concludes this paper.

The following notations are used throughout:

1) The quantity $\phi(q)$ is the Euler’s totient function. It is equal to the
number of integers $k\in[1,q]$ satisfying $(k,q)=1$.

2) $\mathcal{S}_{q}$ and $\mathcal{P}_{q}$ denote the Ramanujan subspace (or
Exactly $q$-periodic) subspace and $q$-periodic subspace, respectively. The
$\text{dim}(\mathcal{S}_{q})$ or $\text{dim}(\mathcal{P}_{q})$ denotes the
dimension of the subspace $\mathcal{S}_{q}$ or $\mathcal{P}_{q}$.

3) $\mathbf{x}_{q}$ and $\check{\mathbf{x}}_{q}$ denote the Exactly
$q$-periodic component and the $q$-periodic component of a signal
$\mathbf{x}$, respectively.

4) The symbols $\mathbb{Z}$ and $\mathbb{Z}_{+}$ stand for the set of integers
and the set of nonnegative integers, respectively.

5) The uppercase bold letters ($\mathbf{A}$, $\mathbf{P}$, etc.) represent
matrices. $\text{range}(\mathbf{P})$ denotes the range space of $\mathbf{P}$,
which is spanned by the columns of $\mathbf{P}$.

6) The notation $q_{i}|q$ means that $q_{i}$ is a divisor of $q$. And
$q_{i}\nmid q$ represents that $q_{i}$ is not a divisor of $q$.

##  II Exactly periodic signal and Ramanujan subspace

This section gives the definitions of the periodic signal and the exactly
periodic signal, based on which the concepts of the periodic subspace and the
exactly periodic subspace are naturally defined. The exactly periodic subspace
can be achieved by constructing the Ramanujan subspace based on the Ramanujan
sums.

###  II-A Periodic signal and exactly periodic signal

The meaning of the periodicity of a signal seems to be easily understood, but
there exists some ambiguity when we use the concept of the periodicity to
characterize the temporal nature of the signal. It is necessary to formally
give the concept of the periodicity before the further discussion, especially
the concept of the exactly periodic signal which is first introduced in [7].

######  Definition 1 ($q$-periodic)

The signal $x(n)$ of length $N$ ($1\leq n\leq N$) is $q$-periodic if there
exists an integer $q$ such that

$\displaystyle x(n+q)=x(n)$ (1)

for all $n$ and $q$ is called the period of $x(n)$.

However, the definition of the $q$-periodic signal is ambiguous, since any
$q$-periodic signal is also $qM$-periodic for $M\in\mathbb{Z}_{+}$. For
example, the periodic signal $\mathbf{1}$ (the vector of all ones) with the
period $1$ can be considered as any period according to Definition 1. To
eliminate this ambiguity, a more definite definition about the periodicity is
given as follows.

######  Definition 2 (exactly $q$-periodic)

The signal $x(n)$ of length $N$ ($1\leq n\leq N$) is exactly $p$-periodic, if
and only if $q$ is the smallest integer in Eq. (1) for all $n$ and $q$ is
called the exactly period of $x(n)$.

Although the concept of the exactly period is first introduced in [7], it has
been implicitly or definitely used in [9, 10, 20, 11]. According to the
definition of the exactly $q$-periodic signal, the periodic signal
$\mathbf{1}$ with the period $1$ is just the exactly $1$-periodic one without
any ambiguity. Obviously, the exactly period is an important concept for
characterizing the periodicity of the signal.

With the definitions about the periodicity of a signal, we can naturally
define two subspaces. One is referred to as the $q$-periodic subspace denoted
by $\mathcal{P}_{q}$, which is the set of all the $q$-periodic signals. The
other is referred to as the exactly $q$-periodic subspace denoted by
$\mathcal{S}_{q}$, which is the set of all the exactly $q$-periodic signals.
The framework proposed in this paper is mainly based on the exactly periodic
subspace. Although the method for constructing the exactly periodic subspace
was proposed based on the subspace intersection in [7], a more effective and
direct way to generate the subspace is presented by constructing the Ramanujan
subspace based on the Ramanujan sums.

###  II-B Ramanujan sums

For any fixed integer $q$, the Ramanujan sums is a sequence with period $q$
and is defined as follows

$\displaystyle c_{q}(n)=\underset{(k,q)=1}{\sum_{k=1}^{q}}e^{j2\pi kn/q}$ (2)

where $(k,q)$ denotes the greatest common divisor of $k$ and $q$. According to
[9], the Ramanujan sums has the following properties.

Property 1. If $q$ is a prime number, then

$\displaystyle c_{q}(n)=\left\\{\begin{array}[]{ll}q-1&\text{if $n$ is
multiple of $q$}\\\ -1&\text{otherwise}\end{array}\right.$ (5)

Property 2. If $q=p^{m}$ for some prime $p$ and integer $m>1$ then

$\displaystyle c_{p^{m}}(n)=\left\\{\begin{array}[]{ll}0&\text{ if
$p^{m-1}\nmid n$}\\\ -p^{m-1}&\text{ if $p^{m-1}|n$ but $p^{m}\nmid n$}\\\
p^{m-1}(p-1)&\text{ if $p^{m}|n$}\end{array}\right.$ (9)

Property 3. if $q_{1}$ and $q_{2}$ are coprime, i.e., $(q_{1},q_{2})=1$, then

$\displaystyle c_{q_{1}q_{2}}(n)=c_{q_{1}}(n)c_{q_{2}}(n)$ (10)

Instead of calculating the Ramanujan sums with recursive computation in [9],
we give a more direct method based on the above three properties and the
fundamental theorem of arithmetic [16]. The procedure of calculating the
Ramanujan sums consists of the following steps.

Firstly, according to the fundamental theorem of arithmetic[16], any positive
integer $q>1$ can be uniquely represented as a product of prime powers

$\displaystyle q=p_{1}^{\alpha_{1}}\cdots
p_{k}^{\alpha_{k}}=\prod_{i=1}^{k}p_{i}^{\alpha_{i}}$ (11)

where $p_{1}<p_{2}<\cdots<p_{k}$ are primes and $\alpha_{i}\in\mathbb{Z}_{+}$.

Secondly, the Ramanujan sums $c_{p_{i}^{\alpha_{i}}}(n)$ of the prime power
$p_{i}^{\alpha_{i}}$ can be calculated by using the Property 1 and Property 2.

Finally, the Ramanujan sums $c_{q}(n)$ of $q$ can be calculated from
$\\{c_{p_{i}^{\alpha_{i}}}(n)\\}_{i=1}^{k}$ by using the Property 3.

With the above three steps, any Ramanujan sums $c_{q}(n)$ of the positive
integer $q$ can be directly calculated. For example, let $q=12$, which can be
represented as

$\displaystyle q=2^{2}\cdot 3$

The Ramanujan sums of $2^{2}$ and $3$ are

$\displaystyle\begin{array}[]{ll}c_{2^{2}}(12)&=2,0,-2,0,2,0,-2,0,2,0,-2,0\\\
c_{3}(12)&=2,-1,-1,2,-1,-1,2,-1,-1,2,-1,-1\end{array}$

Then the Ramanujan sums $c_{12}(12)$ can be obtained, that is

$\displaystyle c_{12}(12)=4,0,2,0,-2,0,-4,0,-2,0,2,0$

###  II-C Ramanujan subspace

Based on the Ramanujan sums $c_{q}(n)$, we can directly construct the
Ramanujan subspace $\mathcal{S}_{q}$ associated with the period $q$. The
Ramanujan subspace $\mathcal{S}_{q}$ is just the exactly $q$-periodic subspace
that will be used to decompose a signal according to its periodic structure in
this paper.

Given the period $q\in\mathbb{Z}+$, let $c_{q}(n)$ be the Ramanujan sums and
$\mathbf{B}_{q}$ be its $q\times q$ integer circulant matrix defined by

$\displaystyle\small\mathbf{B}_{q}=\left[\begin{array}[]{cccc}c_{q}(0)&c_{q}(q-1)&\cdots&c_{q}(1)\\\
c_{q}(1)&c_{q}(0)&\cdots&c_{q}(2)\\\ c_{q}(2)&c_{q}(1)&\cdots&c_{q}(3)\\\
\vdots&\vdots&\ddots&\vdots\\\ c_{q}(q-2)&c_{q}(q-3)&\cdots&c_{q}(q-1)\\\
c_{q}(q-1)&c_{q}(q-2)&\cdots&c_{q}(0)\end{array}\right]$ (18)

Thus, the Ramanujan subspace $\mathcal{S}_{q}\subset\mathcal{R}^{q}$ can be
spanned by the columns of $\mathbf{B}_{q}$, i.e.,
$\mathcal{S}_{q}=\text{Range}(\mathbf{B}_{q})$.

According to Theorem 3 and 4 in [9], we know that both the column rank of
$\mathbf{B}_{q}$ and the dimension of the subspace $\mathcal{S}_{q}$ are equal
to $\phi(q)$, i.e.,
$\text{rank}(\mathbf{B}_{q})=\text{dim}(\mathcal{S}_{q})=\phi(q)$. More
importantly, according to Theorem 10 in [9], the Ramanujan subspace
$\mathcal{S}_{q}$ is the exactly $q$-periodic subspace. Let
$\mathcal{S}_{q}^{\bot}$ denote the orthogonal complement of
$\mathcal{S}_{q}$. The orthogonal projection operator (or projection matrix)
$\mathbf{P}_{q}\in\mathcal{R}^{q\times q}$ along $\mathcal{S}_{q}^{\bot}$ into
$\mathcal{S}_{q}$ is given by

$\displaystyle\mathbf{P}_{q}=\frac{\mathbf{B}_{q}}{q}$ (19)

which satisfies that $\mathbf{P}_{q}^{2}=\mathbf{P}_{q}$ and
$\mathbf{P}_{q}^{T}=\mathbf{P}_{q}$. The Ramanujan subspace $\mathcal{S}_{q}$
is also determined by the range space of $\mathbf{P}_{q}$, i.e.,
$\mathcal{S}_{q}=\text{Range}(\mathbf{P}_{q})$.

The Ramanujan subspace is an important concept used in this paper, based on
which the whole signal space is decomposed into a set of subspaces. Before
doing this, it is required to solve the problem of ‘how to orthogonally
project a signal of any length into the Ramanujan subspace?’ in the following
section.

##  III Orthogonal projection into $\mathcal{S}_{q}$

In this section, we focus on the problem of orthogonally projecting a signal
$\mathbf{x}$ of any length $N$ into the Ramanujan subspace $\mathcal{S}_{q}$.
Although a method for calculating the projection is given in [10] when $q$ is
the divisor of the length $N$, it is not suitable for the signal of any
length, i.e., $q$ is not the divisor of the length $N$. To solve the problem,
we define the orthogonal projection operator into $\mathcal{S}_{q}$ as follows

$\displaystyle\text{Proj}:\mathcal{R}^{N}\mapsto\mathcal{S}_{q}$ (20)

Then the orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ into
$\mathcal{S}_{q}$ can be represented as
$\mathbf{x}_{q}=\text{Proj}(\mathbf{x},\mathcal{S}_{q})$. To implement the
projection operator, we first present a method of projecting a signal of
length $N$ equal to the integer multiple of the period $q$ into
$\mathcal{S}_{q}$, and then extend the method to the general case.

###  III-A Projection of signal with length of multiple periods

Given the signal $\mathbf{x}\in\mathcal{R}^{N}$, let $N$ be equal to the
integer multiple of the period $q$ such that $N=qM$ where
$M\in\mathbb{Z}_{+}$. The Ramanujan subspace
$\mathcal{\widetilde{S}}_{q}\subset\mathcal{R}^{N}$ can be spanned by the
columns of the following matrix

$\displaystyle\mathbf{C}_{qM}=\frac{1}{\sqrt{M}}\left[\begin{array}[]{c}\mathbf{P}_{q}\\\
\mathbf{P}_{q}\\\ \vdots\\\
\mathbf{P}_{q}\end{array}\right]\in\mathcal{R}^{qM\times q}$ (25)

which is the $M$ repetitions of the projection matrix $\mathbf{P}_{q}$ of
$\mathcal{S}_{q}$ defined in Eq. (19), i.e.,
$\mathcal{\widetilde{S}}_{q}=\text{Range}(\mathbf{C}_{qM})$. The subspace
$\mathcal{\widetilde{S}}_{q}$ has the same dimension as $\mathcal{S}_{q}$,
i.e.,
$\text{dim}(\mathcal{\widetilde{S}}_{q})=\text{dim}(\mathcal{S}_{q})=\phi(q)$.
The orthogonal projection matrix $\mathbf{\widetilde{P}}_{q}$ along the
orthogonal complement subspace $\mathcal{\widetilde{S}}_{q}^{\perp}$ into the
subspace $\mathcal{\widetilde{S}}_{q}$ is given by

$\displaystyle\mathbf{\widetilde{P}}_{q}$
$\displaystyle=\mathbf{C}_{qM}(\mathbf{C}_{qM}^{T}\mathbf{C}_{qM})^{-1}\mathbf{C}_{qM}^{T}$
(26)
$\displaystyle=\frac{1}{M}\left[\begin{array}[]{ccc}\mathbf{P}_{q}&\cdots&\mathbf{P}_{q}\\\
\vdots&\ddots&\vdots\\\
\mathbf{P}_{q}&\cdots&\mathbf{P}_{q}\end{array}\right]$ (30)

which is a $qM\times qM$ matrix and satisfies that
$\mathbf{\widetilde{P}}_{q}^{2}=\mathbf{\widetilde{P}}_{q}$ and
$\mathbf{\widetilde{P}}_{q}^{T}=\mathbf{\widetilde{P}}_{q}$. Since the
subspace $\mathcal{\widetilde{S}}_{q}\subset\mathcal{R}^{N}$ is completely
determined by $\mathcal{S}_{q}\subset\mathcal{R}^{q}$, it is also denoted by
$\mathcal{S}_{q}$ to avoid a multitude of notations as done in [9]. Then, the
orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ into
$\mathcal{{S}}_{q}$ can be obtained by

$\displaystyle\mathbf{x}_{q}=\text{Proj}(\mathbf{x},\mathcal{S}_{q})=\mathbf{\widetilde{P}}_{q}\mathbf{x},\text{
}\forall\mathbf{x}\in\mathcal{R}^{N}$ (31)

###  III-B Projecting signal with any length

Instead of modifying the project matrix as done in Eq. (30), we give a more
effective method to calculate the projection by partitioning the $\mathbf{x}$
into a set of blocks and by calculating the projection of the mean of these
blocks into $\mathcal{S}_{q}$. Moreover, the method of calculating mean
projection of blocks can be easily applied to the signal of any length.

###### Theorem 1

Let the signal $\mathbf{x}\in\mathcal{R}^{N}$ where $N=qM$ be partioned into
$M$ blocks

$\displaystyle\mathbf{x}=\left[\begin{array}[]{c}\mathbf{x}^{(1)}\\\
\mathbf{x}^{(2)}\\\ \vdots\\\ \mathbf{x}^{(M)}\end{array}\right]$ (36)

where $\mathbf{x}^{(m)}\in\mathcal{R}^{q}$ is the $m$-th block, $1\leq m\leq
M$. Let $\mathbf{s}\in\mathcal{R}^{q}$ be the mean of the $M$ blocks, defined
by

$\displaystyle\mathbf{s}=\frac{1}{M}\sum_{m=1}^{M}\mathbf{x}^{(m)}$ (37)

Let $\mathbf{\bar{x}}_{q}\in\mathcal{S}_{q}$ be the projection of $\mathbf{s}$
into the Ramanujan subspace $\mathcal{S}_{q}$, which obtained by
$\mathbf{\bar{x}}_{q}=\mathbf{P}_{q}\mathbf{s}$ where $\mathbf{P}_{q}$ is the
orthogonal projection matrix of $\mathcal{S}_{q}$ defined in Eq. (19). Then,
the orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ into the Ramanujan
subspace $\mathcal{S}_{q}$ can be obtained by repeating $M$ times of
$\mathbf{\bar{x}}_{q}$ as follows

$\displaystyle\mathbf{x}_{q}=\text{Proj}(\mathbf{x},\mathcal{S}_{q})=\left[\begin{array}[]{c}\mathbf{\bar{x}}_{q}\\\
\mathbf{\bar{x}}_{q}\\\ \vdots\\\ \mathbf{\bar{x}}_{q}\end{array}\right]$ (42)

###### Proof:

The orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ into the subspace
$\mathcal{S}_{q}$ is given by

$\displaystyle\mathbf{x}_{q}$
$\displaystyle=\mathbf{\widetilde{P}}_{q}\mathbf{x}$
$\displaystyle=\frac{1}{M}\left[\begin{array}[]{ccc}\mathbf{P}_{q}&\cdots&\mathbf{P}_{q}\\\
\mathbf{P}_{q}&\cdots&\mathbf{P}_{q}\\\ \vdots&\ddots&\vdots\\\
\mathbf{P}_{q}&\cdots&\mathbf{P}_{q}\end{array}\right]\left[\begin{array}[]{c}\mathbf{x}^{(1)}\\\
\mathbf{x}^{(2)}\\\ \vdots\\\ \mathbf{x}^{(M)}\end{array}\right]$
$\displaystyle=\frac{1}{M}\left[\begin{array}[]{c}\mathbf{P}_{q}\sum_{k=1}^{M}\mathbf{x}^{(k)}\\\
\mathbf{P}_{q}\sum_{k=1}^{M}\mathbf{x}^{(k)}\\\ \vdots\\\
\mathbf{P}_{q}\sum_{k=1}^{M}\mathbf{x}^{(k)}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{c}\mathbf{\bar{x}}_{q}\\\
\mathbf{\bar{x}}_{q}\\\ \vdots\\\ \mathbf{\bar{x}}_{q}\end{array}\right]$

where $\mathbf{\bar{x}}_{q}=\mathbf{P}_{q}\mathbf{s}$ and
$\mathbf{s}=\frac{1}{M}\sum_{k=1}^{M}\mathbf{x}^{(k)}$.

∎

The theorem 1 reveals that the orthogonal projection $\mathbf{x}_{q}$ of
$\mathbf{x}$ into $\mathcal{S}_{q}$ is only depended on the mean $\mathbf{s}$
of the blocks from $\mathbf{x}$, which provides a way to calculate the
orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ of any length. Consider
the general case where the length $N$ of the signal $\mathbf{x}$ is not the
integer multiple of the period $q$. Let $M=\lceil{N}/{q}\rceil$, where $\lceil
z\rceil$ is the smallest integer greater than or equal to $z$. The signal
$\mathbf{x}$ is padded with $\tilde{N}$ zeros to obtain a new signal
$\tilde{\mathbf{x}}$ of length $qM$, where $\tilde{N}=qM-N$. The signal
$\tilde{\mathbf{x}}$ is partitioned into $M$ blocks
$\\{\mathbf{x}^{(m)}\\}_{m=1}^{M}$, where $\mathbf{x}^{(m)}\in\mathcal{R}^{q}$
is the $m$-th block. The mean of these blocks can be calculated by

$\displaystyle\mathbf{s}=\sum_{m=1}^{M}\mathbf{W}\mathbf{x}^{(m)}$ (43)

where $\mathbf{W}$ is a diagonal matrix defined by

$\displaystyle\mathbf{W}=\text{diag}(\overset{q-\tilde{N}}{\overbrace{\frac{1}{M},\cdots,\frac{1}{M}}},\overset{\tilde{N}}{\overbrace{\frac{1}{M-1},\cdots,\frac{1}{M-1}}})$
(44)

Thus, the orthogonal projection $\tilde{\mathbf{x}}_{q}$ of the signal
$\tilde{\mathbf{x}}$ into $\mathcal{S}_{q}$ can be calculated with theorem 1.
Finally, the orthogonal projection $\mathbf{x}_{q}$ of $\mathbf{x}$ into
$\mathcal{S}_{q}$ can be obtained by simply preserving the first $N$ elements.

##  IV Ramanujan subspace pursuit

A signal can be decomposed into a sum of elementary building blocks so that it
can be represented in a more meaningful way. For example, the Fourier and
wavelet transform can decompose a signal into a sum of frequency components;
the atom decomposition can decompose a signal into a sum of time-frequency
blocks or atom components. Similarly, our purpose is to decompose a signal
into a sum of exactly periodic components to obtain its exactly periodic
representation. The periodic decomposition can be modeled as

$\displaystyle\mathbf{x}=\sum_{k=1}^{K}\mathbf{x}_{q_{k}}+\mathbf{r}$ (45)

where $\mathbf{x}_{q_{k}}=\text{Proj}(\mathbf{x},\mathcal{S}_{q_{k}})$ is the
exactly $q_{k}$-periodic component and $\mathbf{r}$ is the residual component.
Note that, if $\mathbf{r}$ is allowed for $\|\mathbf{r}\|>0$, Eq. (45) is the
problem of the signal approximation, while for $\|\mathbf{r}\|=0$, Eq. (45) is
the problem of the signal representation. If the periodic decomposition can
accurately indicate the periodic components in a signal, then it is useful for
periodic estimation, periodic analysis, and other applications in signal
processing, since it reveals the intrinsic periodic structure of the signal.

To obtain the unique periodic representation, some strategy or constraint is
required, such as greedy strategy as done in Matching Pursuit (MP) [26] or
other methods of subspace pursuit (SP). The greedy strategy for the periodic
decomposition is performed by the following procedures: finding the most
dominant periodic components and removing it from the current signal to obtain
the residual signal, and then repeating the same work for the residual signal
at the next iteration until a stopping criterion is met. Compared with the
methods of MP or SP, the major challenge in the periodic decomposition is how
to select the dominant periodic component in each iteration, due to the
definition of the dominant periodic component in a signal is not clear yet.

Given the maximum expected period $Q$, a set of Ramanujan subspaces
$\\{\mathcal{S}_{q}\\}_{q=1}^{Q}$ are generated according to the period $q$
ranging from 1 up to $Q$ and are used to capture the exactly periodic
components in the signal. Note that, since the dimension of all these
subspaces, $\sum_{q=1}^{Q}\phi(q)$, is far greater than the length $N$ of the
signal, i.e., $\sum_{q=1}^{Q}\phi(q)\gg N$, these Ramanujan subspaces are not
orthogonal to each other (only when $N=lcm(1,\cdots,Q)$, these Ramanujan
subspaces are orthogonal to each other). However, the most dominant periodic
component cannot be simply determined in terms of the value of the energy of
the component captured by the subspace, as done in MP or SP. The reason is
that the dimension $\phi(q)$ of the Ramanujan subspace is directly related to
the period $q\in[1,Q]$ and the higher dimensional subspace tends to capture
more energy from the signal. Therefore, it is necessary to first define a
suitable metric that can be used to select the dominant periodic component at
each iteration.

###  IV-A Periodicity metric

The periodicity metric of the exactly periodic component is first defined
based on the autocorrelation function. Then the relationship of the
periodicity metric with the energy of the exactly periodic component and the
period is revealed.

######  Definition 3 (Periodicity metric)

Given a signal $\mathbf{x}\in\mathcal{R}^{N}$ and the maximum period $Q$,
there are $Q$ exactly periodic components
$\\{\mathbf{x}_{q}=\text{Proj}(\mathbf{x},\mathcal{S}_{q})\\}_{q=1}^{Q}$ from
which the most dominant periodic component will be selected. The periodicity
metric of the exactly periodic component $\mathbf{x}_{q}$ is defined by

$\displaystyle
P(\mathbf{x}_{q},q)\triangleq\sum_{l=0}^{M-1}\varphi_{\mathbf{x}_{q}}(lq)$
(46)

where $M=\lceil\frac{N}{q}\rceil$ and $\varphi_{\mathbf{x}_{q}}(\cdot)$ is the
the autocorrelation function (ACF) of $\mathbf{x}_{q}$ is defined by

$\displaystyle\varphi_{\mathbf{x}_{q}}(k)=\sum_{j=1}^{N-k}\mathbf{x}_{q}(j)\mathbf{x}_{q}(j+k)$
(47)

The periodicity metric $P(\mathbf{x}_{q},q)$ simultaneously measures the
energy $\|\mathbf{x}_{q}\|^{2}=\varphi_{\mathbf{x}_{q}}(0)$ of the exactly
periodic component $\mathbf{x}_{q}$ for $l=0$ and the periodic energies
$\sum_{l=1}^{M-1}\varphi_{\mathbf{x}_{q}}(lq)$ contributed from its
periodicity for $l>0$. Although the definition of the periodicity metric of
$\mathbf{x}_{q}$ depends on its ACF, the following theorem reveals that the
periodicity metric can be calculated based on the energy
$\|\mathbf{x}_{q}\|^{2}$ of $\mathbf{x}_{q}$, which can be obtained by
projecting the signal $\mathbf{x}$ into the Ramanujan subspace
$\mathcal{S}_{q}$.

###### Theorem 2

The periodicity metric of $\mathbf{x}_{q}$ of length $N$ can be directly
calculated from the energy of $\mathbf{x}_{q}$ as follows

$\displaystyle P(\mathbf{x}_{q},q)=\frac{N+q}{2q}\|\mathbf{x}_{q}\|^{2}$ (48)

Moreover, if $N\gg q$, then the equivalent metric can be approximated by

$\displaystyle P(\mathbf{x}_{q},q)\approx\frac{1}{2q}\|\mathbf{{x}}_{q}\|^{2}$
(49)

###### Proof:

Let $\mathbf{x}_{q}=\text{Proj}(\mathbf{x},\mathcal{S}_{q})$ where
$\mathbf{x}\in\mathcal{R}^{N}$, then $\mathbf{x}_{q}$ is partitioned into $M$
blocks

$\displaystyle\mathbf{x}_{q}=\left[\begin{array}[]{c}\mathbf{x}_{q}^{(1)}\\\
\mathbf{x}_{q}^{(2)}\\\ \vdots\\\ \mathbf{x}_{q}^{(M)}\end{array}\right]$

where $M=\lceil\frac{N}{q}\rceil$. The energy $\|\mathbf{x}_{q}\|^{2}$ of
$\mathbf{x}_{q}$ can be represented as

$\displaystyle\|\mathbf{x}_{q}\|^{2}$
$\displaystyle=M\left\|\frac{1}{M}\sum_{m=1}^{M}\mathbf{x}_{q}^{(m)}\right\|^{2}=\frac{1}{M}\left\|\sum_{m=1}^{M}\mathbf{x}_{q}^{(m)}\right\|^{2}$
$\displaystyle=\frac{1}{M}\left\|\left[\begin{array}[]{cccccc}\mathbf{x}_{q}(1)&+&\cdots&+&\mathbf{x}_{q}(1+(M-1)q)\\\
\mathbf{x}_{q}(2)&+&\cdots&+&\mathbf{x}_{q}(2+(M-1)q)\\\
\vdots&\vdots&\cdots&\vdots&\vdots\\\
\mathbf{x}_{q}(q)&+&\cdots&+&\mathbf{x}_{q}(Mq)\end{array}\right]\right\|^{2}$
$\displaystyle=\frac{1}{M}\left(\|\mathbf{x}_{q}\|^{2}+2\sum_{l=1}^{M-1}\varphi_{\mathbf{x}_{q}}(lq)\right)$

where $\varphi_{\mathbf{x}_{q}}(\cdot)$ is the ACF defined in Eq. (47). Thus,
we can obtain that

$\displaystyle
P(\mathbf{x}_{q},q)=\sum_{l=0}^{M-1}\varphi_{\mathbf{x}_{q}}(lq)=\frac{M+1}{2}\|\mathbf{x}_{q}\|^{2}$

And, hence

$\displaystyle P(\mathbf{x}_{q},q)=\frac{N+q}{2q}\|\mathbf{x}_{q}\|^{2}$

Note that, $P(\mathbf{x}_{q},q)$ and $P(\mathbf{x}_{q},q)/N$ are completely
equivalent when they are used to select the dominant periodic component from
the periodic components $\\{\mathbf{x}_{q}\\}_{q=1}^{Q}$, where $Q$ is the
maximum period. Hence, the periodicity metric $P(\mathbf{x}_{q},q)$ is
equivalent to

$\displaystyle P(\mathbf{x}_{q},q)=\frac{1+N/q}{2q}\|\mathbf{{x}}_{q}\|^{2}$

Moreover, if $N\gg q$, then the exactly periodicity metric can be approximated
by

$\displaystyle P(\mathbf{x}_{q},q)\approx\frac{1}{2q}\|\mathbf{{x}}_{q}\|^{2}$

∎

###  IV-B Ramanujan subspace pursuit

The Ramanujan subspace pursuit (RSP) is iteratively performed by approximating
a signal $\mathbf{x}$ with its dominant exactly periodic components and by
calculating its residual. Let $\mathbf{x}$ be a signal of length $N$, $Q$ be
the maximum period, and $\Phi(Q)=\sum_{q=1}^{Q}\phi(q)$. To represent the
signal $\mathbf{x}$ with its exactly periodic components, it requires that
$\Phi(Q)\geq N$ so that the sum of the Ramanujan subspaces
$\\{\mathcal{S}_{q}\\}_{q=1}^{Q}$ forms the complete space $\mathcal{R}^{N}$,
i.e., $\mathcal{S}_{1}+\cdots+\mathcal{S}_{Q}=\mathcal{R}^{N}$. We call
$\Phi(Q)\geq N$ the condition of the periodic representation of the RSP. Then,
the RSP is performed as follows.

Let the initial residual $\mathbf{r}^{0}=\mathbf{x}$. In iteration $k$, given
the maximum period $Q$ and $\Gamma=\\{1,\cdots,Q\\}$, the residual
$\mathbf{r}^{k}$ is calculated by

$\displaystyle\mathbf{r}^{k}=\mathbf{r}^{k-1}-\mathbf{x}_{q_{k}}$ (50)

where $\mathbf{x}_{q_{k}}=\text{Proj}(\mathbf{r}^{k-1},\mathcal{S}_{q_{k}})$.
The period $q_{k}$ is selected so that $\mathbf{x}_{q_{k}}$ has the maximum
value of the periodicity metric among all the periods from $1$ to $Q$, that is

$\displaystyle
q_{k}=\underset{q\in\Gamma}{\text{argmax}}\left\\{P\left(\mathbf{x}_{q}^{(k)},q\right)\right\\}$
(51)

where $\mathbf{x}_{q}^{(k)}=\text{Proj}(\mathbf{r}^{k-1},\mathcal{S}_{q})$.
Since the projection $\mathbf{x}_{q_{k}}$ is orthogonal to $\mathbf{r}^{k-1}$
in Eq. (50), it implies that

$\displaystyle\|\mathbf{r}^{k}\|^{2}=\|\mathbf{r}^{k-1}\|^{2}-\|\mathbf{x}_{q_{k}}\|^{2}$
(52)

Summing Eq. (50) over $k$ between $1$ and $K$ yields

$\displaystyle{\mathbf{x}}=\sum_{k=1}^{K}\mathbf{x}_{q_{k}}+\mathbf{r}^{K}$
(53)

Similarly, summing Eq. (52) over $k$ between $1$ and $K$ gives

$\displaystyle\|\mathbf{x}\|^{2}=\sum_{k=1}^{K}\|\mathbf{x}_{q_{k}}\|^{2}+\|\mathbf{r}^{K}\|^{2}$
(54)

The following theorem proves that the RSP converges exponentially and the norm
of the residual in the RSP converges to zero.

###### Theorem 3

Let $\mathbf{x}$ be a signal of length $N$ and $Q$ be the maximum period
satisfying $\Phi(Q)>N$, where $\Phi(Q)=\underset{q\in\Gamma}{\sum}\phi(q)$ and
$\Gamma=\\{1,\cdots,Q\\}$. The residual $\mathbf{r}^{k}$ computed by the
Ramanujan subspace pursuit in the iteration $k$ satisfies

$\displaystyle\|\mathbf{r}^{k}\|^{2}\leq(1-\mu(\Gamma))^{k}\|\mathbf{x}\|^{2}$
(55)

where $\mu(\Gamma)$ is a constant number and satisfies $\mu(\Gamma)\in(0,1]$.
As a consequence,

$\displaystyle\mathbf{x}=\sum_{k=1}^{+\infty}\mathbf{x}_{q_{k}}$ (56)

and

$\displaystyle\|\mathbf{x}\|^{2}=\sum_{k=1}^{+\infty}\|\mathbf{x}_{q_{k}}\|^{2}$
(57)

###### Proof:

Suppose that the residual $\mathbf{r}^{k}$ is already computed. In the next
iteration, the period $q_{k+1}$ is chosen with the maximum value of the
periodicity metric defined in Eq. (48), which satisfies

$\displaystyle\frac{N+q_{k+1}}{2q_{k+1}}\|\mathbf{x}_{q_{k+1}}\|^{2}\geq\underset{q\in\Gamma}{\text{max}}\frac{N+q}{2q}\|\mathbf{x}_{q}^{(k+1)}\|^{2}$
(58)

or

$\displaystyle\|\mathbf{x}_{q_{k+1}}\|^{2}\geq\underset{q\in\Gamma}{\text{max}}\frac{(N+q)q_{k+1}}{(N+q_{k+1})q}\|\mathbf{x}_{q}^{(k+1)}\|^{2}$
(59)

where $\mathbf{x}_{q_{k+1}}=\text{Proj}(\mathbf{r}^{k},\mathcal{S}_{q_{k+1}})$
and $\mathbf{x}_{q}^{(k+1)}=\text{Proj}(\mathbf{r}^{k},\mathcal{S}_{q})$.
Then, the decomposition at the $k+1$ iteration is

$\displaystyle\mathbf{r}^{k+1}=\mathbf{r}^{k}-\mathbf{x}_{q_{k+1}}$ (60)

and

$\displaystyle\|\mathbf{r}^{k+1}\|^{2}=\|\mathbf{r}^{k}\|^{2}-\|\mathbf{x}_{q_{k+1}}\|^{2}$
(61)

The rate of decay is

$\displaystyle\frac{\|\mathbf{r}^{k+1}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}=1-\frac{\|\mathbf{x}_{q_{k+1}}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}$
(62)

According to Eq. (59) and (62), we have

$\displaystyle\frac{\|\mathbf{r}^{k+1}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}\leq
1-\underset{q\in\Gamma}{\text{max}}\frac{(N+q)q_{k+1}}{(N+q_{k+1})q}\frac{\|\mathbf{x}_{q}^{(k+1)}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}\leq
1$ (63)

Let

$\displaystyle\eta(\mathbf{f},\Gamma)\triangleq\underset{q\in\Gamma}{\text{max}}\frac{(N+q)q_{k+1}}{(N+q_{k+1})q}\frac{\|\text{Proj}(\mathbf{f},\mathcal{S}_{q})\|^{2}}{\|\mathbf{f}\|^{2}}$
(64)

where $\mathbf{f}\in\mathcal{R}^{N}$, then

$\displaystyle\frac{\|\mathbf{r}^{k+1}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}<1-\eta(\mathbf{r}^{k},\Gamma)$
(65)

Note that, for any $q$ and $q_{k+1}$ taken from $\Gamma$, it satisfies that

$\displaystyle\frac{(N+q)q_{k+1}}{(N+q_{k+1})q}\geq\alpha>0$ (66)

where $\alpha=\frac{N/Q+1}{N+1}$, and we have

$\displaystyle\eta(f,\Gamma)\geq\alpha\underset{q\in\Gamma}{\text{max}}\frac{\|\text{Proj}(\mathbf{f},\mathcal{S}_{q})\|^{2}}{\|\mathbf{f}\|^{2}}$
(67)

In terms of the Theorem 12.6 in [27], it shows that

$\displaystyle\underset{\mathbf{f}\in\mathcal{R}^{N}}{\text{inf}}\underset{q\in\Gamma}{\text{max}}\frac{\|\text{Proj}(\mathbf{f},\mathcal{S}_{q})\|^{2}}{\|\mathbf{f}\|^{2}}>0$
(68)

Hence, we have

$\displaystyle\mu(\Gamma)=\underset{\mathbf{f}\in\mathcal{R}^{N}}{\text{inf}}\eta(\mathbf{f},\Gamma)>0$
(69)

Thus, the rate of decay satisfies

$\displaystyle\frac{\|\mathbf{r}^{k+1}\|^{2}}{\|\mathbf{r}^{k}\|^{2}}<1-\mu(\Gamma)$
(70)

Iterating on Eq. (70) proves that

$\displaystyle\|\mathbf{r}_{k}\|^{2}\leq(1-\mu(\Gamma))^{k}\|\mathbf{x}\|^{2}$
(71)

and such that

$\displaystyle\underset{k\rightarrow 0}{\lim}\|\mathbf{r}_{k}\|^{2}=0$ (72)

Hence, both the Eq. (56) and Eq. (57) are established.

∎

Since the RSP requires to calculate all the projections of the signal
$\mathbf{x}$ into each Ramanujan subspace in each iteration, the RSP suffers
from higher computational cost. For a given signal $\mathbf{x}$ of length $N$
and the maximum period $Q$, the total computational cost of the RSP with $K$
iterations is $\mathcal{O}(KQN^{2})$. Fortunately, we will show that the
computational cost of the RSP can be significantly reduced in the next
section.

##  V Fast Ramanujan subspace pursuit and periodic distance

The high cost of the RSP is mainly attributed to the calculation of the
projections of the exactly periodic components of the signal $\mathbf{x}$ into
each Ramanujan subspace. However, it is unnecessary for the RSP to calculate
these projections in each iteration. In this section, we present a fast
algorithm for the RSP based on the following results. Firstly, the
relationship between the periodic subspace $\mathcal{P}_{q}$ and the Ramanujan
subspace $\mathcal{S}_{q}$ reveals that the energy $\|\mathbf{x}_{q}\|^{2}$ of
the orthogonal projection in $\mathcal{S}_{q}$ can be iteratively calculated
from the energy $\|{\check{\mathbf{x}}}_{q}\|^{2}$ of the orthogonal
projection in $\mathcal{P}_{q}$. Secondly, the
$\|{\check{\mathbf{x}}}_{q}\|^{2}$ can be easily estimated from the ACF of the
signal $\mathbf{x}$ based on the MLE. As a result, the computational cost of
the fast RSP can be reduced from $\mathcal{O}(KQN^{2})$ to
$\mathcal{O}(KN^{2})$ and can be further reduced to $\mathcal{O}(KN\log N)$
with the Fast Fourier Transformation (FFT) algorithm.

###  V-A Periodic subspace and exactly periodic subspace

###### Lemma 1

Let $q_{1},\cdots,q_{K}$ be all the divisors of $q$, including $1$ and $q$,
and $\mathcal{S}_{q_{1}},\cdots,\mathcal{S}_{q_{K}}$ be the Ramanujan
subspaces corresponding to these divisors, then

$\displaystyle\mathcal{S}_{q_{1}},\cdots,\mathcal{S}_{q_{K}}\subseteq\mathcal{P}_{q}$
(73)

where the “$=$” is hold only when $q=1$, i.e.,
$\mathcal{S}_{1}=\mathcal{P}_{1}$.

###### Proof:

Let $q_{k}$ be the divisor of $q$ and $M_{k}=q/q_{k}$. With the Definition 1
and 2, any exactly $q_{k}$-periodic signal $\mathbf{x}_{q_{k}}$ must be
$(q_{k}M_{k})$-periodic for $M_{k}$, i.e.,
$\forall\mathbf{x}_{q_{k}}\in\mathcal{S}_{q_{k}}$, then
$\mathbf{x}_{q_{k}}\in\mathcal{P}_{q}$. However, not every $q$-periodic signal
$\mathbf{x}_{q}$ is exactly $q$-periodic, i.e.,
$\exists\mathbf{x}_{q}\in\mathcal{P}_{q}$, but
$\mathbf{x}_{q}\notin\mathcal{S}_{q_{k}}$. Then, we have

$\displaystyle\mathcal{S}_{q_{k}}\subset\mathcal{P}_{q}$

Hence, we have

$\displaystyle\mathcal{S}_{q_{1}},\cdots,\mathcal{S}_{q_{K}}\subseteq\mathcal{P}_{q}$

where the “$=$” is held only when $q=1$, i.e.,
$\mathcal{S}_{1}=\mathcal{P}_{1}$.

∎

To clearly describe the relationship between the Ramanujan subspace
$\mathcal{S}_{q}$ and the periodic subspace $\mathcal{P}_{q}$, we introduce
the following Theorem.

###### Theorem 4

Let $q_{i}$ and $q_{j}$ be any two different divisors of $q$, the
corresponding Ramanujan subspaces $\mathcal{S}_{q_{i}}$ and
$\mathcal{S}_{q_{j}}$ are orthogonal to each other, that is

$\displaystyle\mathcal{S}_{q_{i}}\perp\mathcal{S}_{q_{j}}$ (74)

where $q_{i}\neq q_{j}$.

###### Proof:

Let $\mathbf{\widetilde{P}}_{q_{i}}$ and $\mathbf{\widetilde{P}}_{q_{j}}$ be
the orthogonal projectors of $\mathcal{S}_{q_{i}}$ and $\mathcal{S}_{q_{j}}$
defined in Eq. (26), respectively, where both $\mathbf{\widetilde{P}}_{q_{i}}$
and $\mathbf{\widetilde{P}}_{q_{j}}$ are the $N\times N$ matrices. Suppose
$\mathbf{x}\in\mathcal{S}_{q_{i}}$, then
$\mathbf{\widetilde{P}}_{q_{i}}\mathbf{x}=\mathbf{x}$. The orthogonal
projection of $\mathbf{x}$ into $\mathcal{S}_{q_{j}}$ is

$\displaystyle\mathbf{\widetilde{P}}_{q_{j}}\mathbf{x}=\mathbf{\widetilde{P}}_{q_{j}}\mathbf{\widetilde{P}}_{q_{i}}\mathbf{x}$

According to the orthogonality of the Ramanujan sequence in Eq. (15) in [9],
we know that any row of $\mathbf{\widetilde{P}}_{q_{i}}$ and any column of
$\mathbf{\widetilde{P}}_{q_{j}}$ are orthogonal to each other if $q_{i}$ and
$q_{j}$ are the divisors of $q$ and $q_{i}\neq q_{j}$, that is

$\displaystyle\mathbf{\widetilde{P}}_{q_{j}}\mathbf{\widetilde{P}}_{q_{i}}=\mathbf{0}_{N\times
N}$

Consequently, we have
$\mathbf{\widetilde{P}}_{q_{j}}\mathbf{x}=\mathbf{0}_{N\times 1}$. Similarly,
we can also obtain that
$\mathbf{\widetilde{P}}_{q_{i}}\mathbf{x}=\mathbf{0}_{N\times 1}$,
$\forall\mathbf{x}\in\mathcal{S}_{q_{j}}$. Thus, both Ramanujan subspaces
$\mathcal{S}_{q_{i}}$ and $\mathcal{S}_{q_{j}}$ are orthogonal to each other

$\displaystyle\mathcal{S}_{q_{i}}\bot\mathcal{S}_{q_{j}}$

∎

According to Lemma 1 and Theorem 4, we can obtain the relationship between the
periodic subspace and the Ramanujan subspace, which is introduced in the
following corollary.

###### Corollary 1

Let $q_{1},\cdots,q_{K}$ be all the divisors of $q$, including $1$ and $q$,
then the corresponding Ramanujan subspaces
$\mathcal{S}_{q_{1}},\cdots,\mathcal{S}_{q_{K}}\subset\mathcal{P}_{q}$ are
orthogonal to each other. Moreover, these Ramanujan subspaces form an
orthogonal decomposition of the $q$-periodic space $\mathcal{P}_{q}$, that is,

$\displaystyle\mathcal{P}_{q}=\mathcal{S}_{q_{1}}\oplus\cdots\oplus\mathcal{S}_{q_{K}}$
(75)

and the dimensions of these subspaces satisfy that

$\displaystyle\text{dim}(\mathcal{P}_{q})=\sum_{k=1}^{K}\text{dim}(\mathcal{S}_{q_{k}})=\sum_{k=1}^{K}\phi(q_{k})=q$
(76)

The Corollary 1 provides a way to orthogonally decompose a $q$-periodic signal
into a sum of its exactly periodic components. Let
${\mathbf{x}}\in\mathcal{P}_{q}$ denote the $q$-periodic signal and
$q_{1},\cdots,q_{K}$ be all the possible divisors of $q$ (including $1$ and
$q$). With Corollary 1, ${\mathbf{x}}$ can be orthogonally decomposed into a
sum of the exactly periodic components $\\{\mathbf{x}_{q_{k}}\\}_{k=1}^{K}$
which satisfy that

$\displaystyle{{\mathbf{x}}}=\sum_{k=1}^{K}\mathbf{x}_{q_{k}}$ (77)

and

$\displaystyle\|{{\mathbf{x}}}\|^{2}=\sum_{k=1}^{K}\|\mathbf{x}_{q_{k}}\|^{2}$
(78)

The Corollary 1 indicates the relationship between the energies of the
projections in the Ramanujan subspaces and the energy of the projection in the
$q$-periodic subspace in Eq. (78). The former are directly related to the
calculation of the periodicity metrics of these exactly periodic components in
the corresponding Ramanujan subspaces. In fact, we will show that these
periodicity metrics can be calculated without projecting the signal into each
Ramanujan subspace.

###  V-B Fast Ramanujan subspace pursuit

Figure 1: Demonstration of the fast Ramanujan subspace pursuit algorithm.
Algorithm 1 Estimating the periodic metrics (EPM)

Input: signal $\mathbf{x}\in\mathcal{R}^{N}$
Output: periodicity metrics $P(\mathbf{x}_{1},1),\cdots,P(\mathbf{x}_{Q},Q)$

Initialize: Compute $\|\mathbf{x}_{1}\|^{2}=\|\check{\mathbf{x}}_{1}\|^{2}$
with Eq. (80)

for $q=2$ to $Q$

Estimate $\|\check{\mathbf{x}}_{q}\|^{2}$ with Eq. (80)

Calculate $\|\mathbf{x}_{q}\|^{2}$ with Eq. (79)

Calculate periodicity metric $P(\mathbf{x}_{q},q)$ with Eq. (48) or (49)

end

Given a signal $\mathbf{x}\in\mathcal{R}^{N}$, our purpose is to calculate the
energy $\|\mathbf{x}_{q}\|^{2}$ of the exactly $q$-periodic component in the
Ramanujan subspace $\mathcal{S}_{q}$ without projecting $\mathbf{x}$ into
$\mathcal{S}_{q}$. Suppose that the energy $\|\check{\mathbf{x}}_{q}\|^{2}$ of
the $q$-peroidic component and the energies
$\\{\|\mathbf{x}_{q_{k}}\|^{2}\\}_{k=1}^{K-1}$ of the exactly periodic
components are already computed, where $q_{1},\cdots,q_{K-1}$ are the first
$K-1$ possible divisors (except $q$) of $q$. Then, according to Eq. (78), the
energy $\|\mathbf{x}_{q}\|^{2}$ of the exactly $q$-periodic component can be
calculated by

$\displaystyle\|\mathbf{x}_{q}\|^{2}=\|\check{\mathbf{x}}_{q}\|^{2}-\sum_{k=1}^{K-1}\|\mathbf{x}_{q_{k}}\|^{2}$
(79)

which can be used to calculate the periodicity metric with Eq. (48) or (49).

Although the energy $\|\check{\mathbf{x}}_{q}\|^{2}$ of the $q$-peroidic
component of $\mathbf{x}$ can be obtained by projecting $\mathbf{x}$ into
$\mathcal{P}_{q}$ with the periodic transform [8], it can be estimated by
using the MLE [17, 7], that is

$\displaystyle\|\check{\mathbf{x}}_{q}\|^{2}=\frac{q}{N}\left(\varphi_{\mathbf{x}}(0)+2\sum_{l=1}^{M-1}\varphi_{\mathbf{x}}(lq)\right)$
(80)

where $\varphi_{\mathbf{x}}(\cdot)$ is the ACF of $\mathbf{x}$ defined in Eq.
(47) and $M=\lfloor N/q\rfloor$. Thus, the energy $\|\mathbf{x}_{q}\|^{2}$ of
the exactly $q$-periodic component of $\mathbf{x}$ can be iteratively computed
according to Eq. (79) with the initial condition
$\mathbf{x}_{1}=\check{\mathbf{x}}_{1}$. Base on $\|\mathbf{x}_{q}\|^{2}$, the
periodicity metric of the exactly periodic component $\mathbf{x}_{q}$ can be
computed with Eq. (48) or (49) for $q$ from $1$ up to $Q$. The implementation
of estimating the periodic metrics (EPM) is shown in Algorithm 1. Based on the
Algorithm 1, the fast RSP (FRSP) can be performed as shown in 2. Moreover, the
whole FRSP algorithm is demonstrated in Fig. 1.

Since the periodicity metrics of the exactly periodic components can be
directly estimated from the ACF of the residual in each iteration, the
computational cost of the FRSP can be reduced to $\mathcal{O}(KN^{2})$
compared with the original RSP algorithm. Moreover, when the ACF of the
residual is more effective calculated by using the FFT in each iteration, the
computational cost of the FRSP can be further reduced to $\mathcal{O}(KN\log
N)$.

Algorithm 2 Fast Ramanujan subspace pursuit (FRSP)

Input: signal $\mathbf{x}\in\mathcal{R}^{N}$
Output: Exactly periodic components $\\{\mathbf{x}_{q}\\}_{q=1}^{Q}$

Initialize: Set $\mathbf{r}_{0}=\mathbf{x}$
for $k=1,\cdots,K$

Calculate exactly periodic measures $\\{P(\mathbf{r}^{k},q)\\}_{q=1}^{Q}$ with
Algorithm 1

Select dominant period
$q_{k}=\underset{q\in\Gamma}{\text{argmax}}\left\\{P(\mathbf{r}^{k},q)\right\\}$

Generate Ramanujan subspace $S_{q_{k}}$

Calculate $\mathbf{x}_{q_{k}}=\text{Proj}(\mathbf{r}^{k},\mathcal{S}_{q_{k}})$

Update residual $\mathbf{r}^{k+1}=\mathbf{r}^{k}-\mathbf{x}_{q_{k}}$

end

###  V-C Periodic energy spectrum and periodic distance

With the FRSP algorithm, any signal $\mathbf{x}$ can be decomposed into a sum
of the exactly periodic components. The decomposition can provide a visual
representation named as the periodic energy spectrum (PES) to represent the
periodic structures of the signal. The PES of the signal $\mathbf{x}$ is
defined by

$\displaystyle PES(q)=\|\mathbf{x}_{q}\|^{2},\text{ for }q=1,\cdots,Q$ (81)

where $\mathbf{x}_{q}$ is the exactly $q$-periodic component obtained by using
the FRSP. Unlike the power spectrum which describes the energy of a signal is
distributed over the continuous frequency in the field of signal processing,
the PES describes the energy of the signal is distributed over the discrete
and integer period. For example, Fig. 2 (c) and (d) show the PESs of the
signal $\mathbf{x}$ of length $200$ in Fig. 2 (a) and the signal $\mathbf{y}$
of length $241$ in Fig. 2 (b), respectively. The PESs of the two signals
reveal their similar periodic structures that both of them contain the exactly
periodic components with periods: $7$, $15$, $24$, and $27$.

Figure 2: (a) Signal $\mathbf{x}$ of length $200$; (b) Signal $\mathbf{y}$ of
length $241$; (c) Periodic energy spectrum of the signal $\mathbf{x}$; (d)
Periodic energy spectrum of the signal $\mathbf{y}$.

The PESs of signals provide a way to measure their distance and similarity
based on their periodic structures. Although the concept of the periodic
distance is first introduced in [23], we will redefine it in a more general
form based on the Hellinger distance (or Hellinger divergence) [21, 22], which
can be seen as a symmetric approximation to the Kullback-Leibler (KL)
divergence [24]. The Hellinger distance can be used to measure the distance of
two histograms $\mathbf{h}_{A}$ and $\mathbf{h}_{B}$

$\displaystyle
D(\mathbf{h}_{A},\mathbf{h}_{B})=\frac{1}{2}\sum_{n=1}^{N}\left(\sqrt{\mathbf{h}_{A}(n)}-\sqrt{\mathbf{h}_{B}(n)}\right)^{2}$
(82)

where $\mathbf{h}_{A}=[\mathbf{h}_{A}(1),\cdots,\mathbf{h}_{A}(N)]$ and
$\mathbf{h}_{B}=[\mathbf{h}_{B}(1),\cdots,\mathbf{h}_{B}(N)]$ are two
histograms.

Similar with the PES in Eq. (81), the histogram of the energy distribution of
the periodic components in the signal $\mathbf{x}$, named periodic energy
histogram, can be defined by

$\displaystyle
h_{E}(\mathbf{x})=\frac{1}{\mathbf{\|x\|^{2}}}\left[\|\mathbf{x}_{1}\|^{2},\cdots,\|\mathbf{x}_{Q}\|^{2}\right]$
(83)

where $Q$ is the maximum period and $\\{\mathbf{x}_{q}\\}_{q=1}^{Q}$ are the
exactly $q$-periodic components. Note that only $K$ elements in
$h_{E}(\mathbf{x})$ are non-zero and are obtained by using the FRSP with $K$
iterations. Then, the Hellinger distance between $\mathbf{x}$ and $\mathbf{y}$
(note that they may have different length) is defined by

$\displaystyle
D\left(h_{E}(\mathbf{x}),h_{E}(\mathbf{y})\right)=\frac{1}{2}\sum_{q=1}^{Q}\left(\frac{\|\mathbf{x}_{q}\|}{\|\mathbf{x}\|}-\frac{\|\mathbf{y}_{q}\|}{\|\mathbf{y}\|}\right)^{2}$
(84)

where $\mathbf{x}_{q}$ and $\mathbf{y}_{q}$ are the exactly $q$-periodic
components of $\mathbf{x}$ and $\mathbf{y}$, respectively. For simplicity, the
periodic energy histogram of $\mathbf{x}$ is defined by

$\displaystyle\mathbf{h_{x}}=\frac{1}{\mathbf{\|x\|}}\left[\|\mathbf{x}_{1}\|,\cdots,\|\mathbf{x}_{Q}\|\right]^{T}$
(85)

Then, the Hellinger distance between $\mathbf{x}$ and $\mathbf{y}$ is defined
by

$\displaystyle
D(\mathbf{x},\mathbf{y})=\frac{1}{2}\|\mathbf{h_{x}}-\mathbf{h_{y}}\|^{2}$
(86)

where $\mathbf{h_{x}}$ and $\mathbf{h_{y}}$ are the periodic energy histograms
of $\mathbf{x}$ and $\mathbf{y}$ defined by Eq. (85), respectively. And, the
cosine Hellinger distance between $\mathbf{x}$ and $\mathbf{y}$ is defined by

$\displaystyle
D_{cos}(\mathbf{x},\mathbf{y})=\frac{1}{2}\frac{\|\mathbf{h_{x}}-\mathbf{h_{y}}\|^{2}}{\|\mathbf{h_{x}}\|\|\mathbf{h_{y}}\|}$
(87)

where $D_{cos}(\mathbf{x},\mathbf{y})\in[0,1]$. By the cosine Hellinger
distance, we can define the periodic similarity between $\mathbf{x}$ and
$\mathbf{y}$ as follows

$\displaystyle S(\mathbf{x},\mathbf{y})=1-D_{cos}(\mathbf{x},\mathbf{y})$ (88)

where $S(\mathbf{x},\mathbf{y})\in[0,1]$.

With the periodic similarity defined in Eq. (88), we can compare the
similarity of two signals based on the periodicity. For example, the signals
$\mathbf{x}$ and $\mathbf{y}$, as shown in Fig. 2 (a) and (b), are very
different in their waveforms. However, their periodic similarity
$S(\mathbf{x},\mathbf{y})$ is $0.9703$ and characterizes them as very similar,
because they have the similar periodic structures, i.e., both signals contain
the same periodic components with periods: $7,15,24$ and $27$, as shown in
Fig. 2 (c) and (d).

##  VI Evaluation

To evaluate the effectiveness of the FRSP algorithm in identifying the
periods, the proposed algorithm is applied to a set of synthetic signals and
is compared with several approaches: Ramanujan periodicity transforms (RPT)
[9, 10], Exactly Periodic Subspace Decomposition (EPSD) [7], and four
periodicity transforms (PTs) including Mbest, BestCorrlation, BestFrequence,
Small2Large proposed in [8]. For all these algorithms, we first show their
performances of identifying the periods in a synthetic signal, then evaluate
their capabilities of detecting periods of signals with the different signal
length and the noise level.

###  VI-A Illustration of identifying periods

The synthetic periodic signal is given by

$\displaystyle x[n]=\cos\frac{2\pi n}{q_{1}}+\cos\frac{2\pi
n}{q_{2}}+\cos\frac{2\pi n}{q_{3}}$ (89)

which consists of three cosine signals with the periods $q_{1}=17,q_{2}=36$,
and $q_{3}=45$. Let the length $N$ of the signal be equal to $3060$, which is
the least common multiplier of the periods, as shown in Fig. 3. The range of
the detected periods is from $1$ up to $60$ for these algorithms. Fig. 4
(a)-(g) show the PESs obtained by using these methods. Since the signal in Eq.
(89) consists of integer periods for each periodic components, the FRSP, RPT
EPSD, BestFrequency, and Small2Large can exactly identify three periodic
components. Although there exists a little distortion in the energies of the
exactly periodic components for the BestCorrelation in Fig. 4 (e), it can also
exactly identify the three periods. However, the Mbest in Fig. 4 (d) only
correctly identifies the periodic component at period $45$.

Figure 3: Periodic signal of length $3060$ with three periods: $17$, $36$, and
$45$.

Figure 4: Periodic energy spectrums of the signal in Fig. (3) by using several
methods: (a) FRSP; (b) RPT; (c) EPSD; (d) Mbest; (e) BestCorrelation; (f)
BestFrequency; (g) Small2Large.

Next, consider the performance of identifying the hidden periods of the
synthetic periodic signal of length $511$ as shown in Fig. 5. Since the length
of the signal is not the multiplier of any hidden periods, the signal consists
of non-integer periods for three periodic components of periods: $17,36$, and
$45$. Fig. 6 shows the PESs obtained by using these algorithms. The FRSP in
Fig. 6 (a) can exactly identify three hidden periods when the length of the
signal is reduced from $3060$ to $511$. However, the RPT in Fig. 6 (b) is
completely failed since none of the hidden periods $17,36,45$ is the divisors
of the signal length $N=511$. In fact, two false periods $1$ and $511$ (the
second period peak is not shown in this plot) are detected by the RPT. Fig. 6
(c) show that some false periods are obviously detected by the ESPD due to the
same reason as the RPT. Although several false periods are detected in Fig. 6
(e), the BestCorrlation achieves the similar result when the length of the
signal is changed from $3060$ to $511$. Both the BestFrequency in Fig. 6 (f)
and Small2Large in Fig. 6 (g) can correctly identify two periods $36$ and
$45$, in spite of missing the period $17$. Although several false periods are
detected, the Mbest in Fig. 6 (d) identifies all three periods and obtains a
better result when the length of the signal is changed.

###  VI-B Robustness of identifying hidden periods

Instead of using the signal in Eq. (89) with fixed periods and length, we
generate the test signal in a more general way to investigate the robustness
of these algorithms for identifying hidden periods. The test signal is mixed
by four periodic components, each of which is generated as follows: 1)
Randomly selecting a period $q\in[1,100]$; 2) Generating a random component of
length $q$ and repeating it multiple times.

To evaluate the performances of these algorithms for identifying hidden
periods of the test signal of different length, the signal length $N$ is
varied from $100$ to $1000$. For each algorithm, the range of the detected
period is from $1$ up to $120$, and the maximum of the iteration is $10$. For
each signal length $N$, the detection process is repeated $50$ times with
randomly selecting four periods and random components. Than, the periodic
similarities of these algorithms are calculated by using Eq. (88) between the
original test signal and detection results of these algorithms. Since the RPT
can correctly detect any periods which are not the divisors of the length of
the tested signal, the RPT is not used for this test. Fig. 7 shows the plot of
the average periodic similarities of these algorithms as functions of the
signal length $N$. According to the results, the FRSP achieves the best
performance among all these algorithms. Especially, it shows that even if the
length of the tested signal is $200$, which may only contain two periods of
the periodic component with the hidden period $100$, the FRSP also can obtain
nearly $0.8$ of the periodic similarity. Note that, the BestCorrelation and
Mbest achieve better performance than EPSD, BestFrequency and Small2Large.

Figure 5: Periodic signal of length $511$ with three periods: $17$, $36$, and
$45$.

Figure 6: Periodic energy spectrums of the signal in Fig. (5) by using several
methods: (a) FRSP; (b) RPT; (c) EPSD; (d) Mbest; (e) BestCorrelation; (f)
BestFrequency; (g) Small2Large.

Similarly, the performances of these algorithms for identifying hidden periods
in noise are evaluated by adding the white noise to the tested signal of fixed
length $500$. The white noise, which comes from NOISEX-92 database [25], is
added to the tested signal at the signal-to-noise ratio (SNR) varied from
$-20$ to $40$ dB. For each algorithm, the range of the detected period is from
$1$ up to $120$, and the maximum of the iteration is $10$. For each noise
level (SNR), the detection process is repeated $50$ times with randomly
selecting four periods and random components. Fig. 8 shows the plot of the
average periodic similarities of these algorithms as functions of SNR of the
tested signal. According to the results, the FRSP achieves the best
performance among all these algorithms. Especially, it shows that even when
the SNR is equal to $0$ dB, the FRSP also can obtain nearly $0.8$ of the
periodic similarity. Note that, the BestCorrelation and Mbest achieve better
performance than EPSD, BestFrequency and Small2Large.

##  VII Conclusion

In this paper, we presented a novel algorithm named Ramanujan subspace pursuit
(RSP) for period estimation and periodic decomposition of a signal based on
the Ramanujan subspace. The RSP is a greedy iterative algorithm and is
implemented by selecting and removing the most dominant periodic component
from the residual signal in the current iteration. The RSP can identify any
periodic components corresponding to the periods from $1$ the maximum period
$Q$ and hence successfully overcomes the limitation in ESPD and RPT, which
only the periods being the divisors of the length of the signal can be
accurately detected. In addition to the period estimation, the RSP also
provides a way to decompose the signal into the sum of the exactly periodic
components, just like the method of the atomic decomposition used in signal
processing. Moreover, the fast RSP, based on the MLE of the energy of the
periodic component in the periodic subspace, has a lower computational cost
and can decompose a signal of length $N$ in $\mathcal{O}(KN\log N)$. Our
results show that the RSP outperforms the current algorithms such as RPT, ESPD
and PTs.

Figure 7: Average periodic similarity v.s. signal length $N$.

Figure 8: Average periodic similarity v.s. $SNR$.

## Appendix

## Acknowledgment

This work was supported in part by the Major Research plan of the National
Natural Science Foundation of China (No. 91120303), National Natural Science
Foundation of China (No. 91220301), Natural Science Foundation of Heilongjiang
Province of China (No. F2015012), and Academic Core Funding of Young Projects
of Harbin Normal University of China (No. KGB201225).

## References

  * [1] D. Anastassiou, “Genomic signal processing,” IEEE Signal Processing Mag., vol. 18, no. 4, pp. 8–20, 2001.
  * [2] D. Sharma, B. Isaac, G. P. S. Raghav and R. Ramaswamy, “Spectral Repeat Finder (SRF): identification of repetitive sequences using Fourier transformation,” Bioinformatics, vol. 20, no. 2, pp. 1405–1412, 2004.
  * [3] R. Arora, W. A. Sethares, J. A. Bucklew, “Latent Periodicities in Genome Sequences,” IEEE Journal of Selected Topics in Signal Processing, vol. 2, no. 3, pp. 332–342, 2008.
  * [4] S. Adalbjornsson, J. Sward, J. Wallin, and A. Jakobsson, “Estimating periodicities in symbolic sequences using sparse modeling,” IEEE Trans. on Signal Process., vol. 63, no. 8, pp. 2142–2150, 2015.
  * [5] J. L. Roux, “Single and multiple F0 contour estimation through parametric spectrogram modeling of speech in noisy environments,” IEEE Trans. Audio, Speech Language Process., vol. 15, no. 4, pp. 1135–1145, 2007.
  * [6] W. A. Sethares, R. D. Morris, J. C. Sethares, “Beat tracking of musical performances using low-level audio features,” IEEE Trans. Speech and Audio Process., vol. 13, no. 2, pp. 438–452, 2005.
  * [7] D. D. Muresan and T. W. Parks, “Orthogonal, exactly periodic subspace decomposition,” IEEE Trans. on Signal Process., vol. 51, no. 9, pp. 2270–2279, 2003.
  * [8] W. A. Sethares and T. W. Staley, “Periodicity transforms,” IEEE Trans. on Signal Process., vol. 47, no. 11, pp. 2953–2964, 1999.
  * [9] P. Vaidyanathan, “Ramanujan sums in the context of signal processing – Part I: fundamentals,” IEEE Trans. on Signal Process., vol. 62, no. 16, pp. 4145–4157, 2014.
  * [10] P. Vaidyanathan, “Ramanujan sums in the context of signal processing – Part II: FIR representations and applications,” IEEE Trans. on Signal Process., vol. 62, no. 12, pp. 4158–4172, 2014.
  * [11] M. Nakashizuka, H. Okumura, and Y. Iiguni, “A sparse periodic decomposition and its application to speech representation,” in Proc. 16th European Signal Processing Conference, 2008, pp. 1–5.
  * [12] P. Vaidyanathan and P. Pal, “The Farey-dictionary for sparse representation of periodic signals,” in Proc. IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP), 2014, pp. 360–364.
  * [13] J. Epps, E. Ambikairajah, and M. Akhtar, “An integer period DFT for biological sequence processing,” in Proc. IEEE International Workshop on Genomic Signal Processing and Statistics, 2008, pp. 1–4.
  * [14] J. Epps, H. Ying, and G. A. Huttley, “Statistical methods for detecting periodic fragments in DNA sequence data,” Biology Direct, vol. 6, no. 21, pp. 1–16, 2011.
  * [15] J. Epps, “A hybrid technique for the periodicity characterization of genomic sequence data,” Eurasip Journal on Bioinformatics & Systems Biology, vol. 2009, no. 1, pp. 1–7, 2009.
  * [16] I. Niven, H. S. Zuckerman, and H. L. Montgomery, “An introduction to the theory of numbers,” John Wiley & Sons Ltd., 2008.
  * [17] J. Wise, J. Caprio, and T. W. Parks, “Maximum likelihood pitch estimation,” IEEE Trans. Acoustics, Speech and Signal Process., vol. 24, no. 5, pp. 418–423, 1976.
  * [18] S. S. Chen, D. L. Donoho and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM Journal on Scientific Computing, vol. 20, no. 1, pp. 33–61, 1998.
  * [19] Tibshirani R. “Regression shrinkage and selection via the LASSO,” Journal of the Royal Statistical Society. Series B (Methodological), vol. 58, no. 1, pp. 267–288, 1996.
  * [20] S. V. Tenneti and P. Vaidyanathan, “Arbitrarily shaped periods in multidimensional discrete time periodicity,” IEEE Signal Process. Lett., vol. 22, no. 10, pp. 1748–1751, 2015.
  * [21] T. Jebara, R. Kondor, A. Howard, K. Bennett, and N. Cesa-Bianchi, “Probability product kernels,” Journal of Machine Learning Research, vol. 5, no. 3, pp. 819–844, 2004.
  * [22] H. D. Tran and H. Li, “Sound event recognition with probabilistic distance SVMs,” IEEE Trans. Audio, Speech Language Process., vol. 19, no. 6, pp. 1556–1568, 2011.
  * [23] M. Vlachos, S. Y. Philip, and V. Castelli, “On Periodicity Detection and Structural Periodic Similarity,” in Proc. SIAM International Conference on Data Mining, 2005, pp. 449–460.
  * [24] T. Jaakkola, M. Meila, and T. Jebara, “Maximum entropy discrimination,” Proc. Advances in Neural Information Processing Systems, vol. 12, 1999.
  * [25] The Rice University, NOISEX-92 database, available at http://www.spib.rice.edu/spib/select-noise.html.
  * [26] S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. on Signal Process., vol. 41, no. 12, pp. 3397–3415, 1993.
  * [27] S. G. Mallat, “A wavelet tour of signal processing: the sparse way,” Academic press, 2008.

|  Deng Shi-Wen Shi-Wen Deng received the B.E degree in the Institute of
Technology from Jia Mu Si University, JiaMuSi, China, in 1997, the M.E in The
school of Computer Science from Harbin Normal University, Harbin, China, in
2005, and the Ph.D degree in the school of Computer Science from Harbin
Institute of Technology in 2012. Currently, he is with the School of
Mathematical Sciences, Harbin Normal University, Harbin, China. His research
interests are in the area of speech and audio signal processing, including
content-based audio analysis, noise suppression, speech/audio
classification/detection.
---|---
Han Ji-Qing Ji-Qing Han received the B.S., M.S. in electrical engineering, and
Ph.D. degrees in computer science from the Harbin Institute of Technology,
Harbin, China, in 1987, 1990, and 1998, respectively. Currently, he is the
associate dean of the school of Computer Science and Technology, Harbin
Institute of Technology. He is a member of IEEE, member of the editorial board
of Journal of Chinese Information Processing, and member of the editorial
board of the Journal of Data Acquisition and Processing. Prof. Han is
undertaking several projects from the National Natural Science Foundation,
863Hi-tech Program, National Basic Research Program. He has won three Second
Prize and two Third Prize awards of Science and Technology of
Ministry/Province. He has published more than 100 papers and 2 books. His
research fields include speech signal processing and audio information
processing.
---