of 48
J. Fluid Mech.
(2024),
v
ol
.
985, A42, doi:10.1017/jfm.2024.70
Spectral proper orthogonal decomposition of
harmonically forced turbulent flows
Liam Heidt
1
,
and Tim Colonius
2
1
Graduate Aerospace Laboratories of the California Institute of Technology, California Institute of
Technology, Pasadena, CA 91101, USA
2
Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, CA
91101, USA
(Received 10 May 2023; revised 11 October 2023; accepted 19 December 2023)
Many turbulent flows exhibit time-periodic statistics. These include turbomachinery flows,
flows with external harmonic
forc
ing and the wakes of bluff bodies. Many existing
techniques for identifying turbulent coherent structures, however, assume the statistics are
statistically stationary. In this paper, we leverage cyclostationary analysis, an extension of
the statistically stationary framework to processes with periodically varying statistics, to
generalize the spectral proper orthogonal decomposition (SPOD) to the cyclostationary
case. The resulting properties of the cyclostationary SPOD (CS-SPOD for short) are
explored, a theoretical connection between CS-SPOD and the harmonic resolvent analysis
is provided, simplifications for the low and high forcing frequency limits are discussed,
and an efficient algorithm to compute CS-SPOD with SPOD-like cost is presented.
We illustrate the utility of CS-SPOD using two example problems: a modified complex
linearized
Ginzburg–Landau model and a high-Reynolds-number turbulent jet.
Key words:
computational methods, low-dimensional models, turbulent flows
1. Introduction
Periodic and quasi-periodic forced turbulent flows are ubiquitous in engineering and
nature. Such flows include those in turbomachinery, weather and climate, flow control with
harmonic actuation, or any other flow that exhibits periodic/quasi-periodic modulation of
the turbulence/statistics. In cases where the forcing is slow compared to the turbulence
time scales, the statistics may be
mod
elled as quasi-stationary (comprising a series of
stationary states). However, in many cases, the forcing is at frequencies commensurate
with the turbulence, and the turbulence structure is not only modulated by, but also altered
Email address for correspondence:
lheidt@caltech.edu
© The Author(s), 2024. Published by Cambridge University Press
985
A42-1
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
L. Heidt and T. Colonius
by, the forcing. In such, a key goal is to identify coherent structures that can be compared
with and contrasted to their occurrence in similar but unforced flows.
The most commonly used technique to identify coherent structures in turbulence is
proper orthogonal decomposition (Lumley
1967
,
1970
; Aubry
et al.
1988
; Sirovich
1989
;
Aubry
1991
), which represents flow data as mutually orthogonal modes whose amplitudes
optimally reconstruct the correlation tensor. When applied in its typical space-only form,
the modes are not coherent in time, leading many researchers to apply
dynamic
mode
decom
po
si
tion and its variants (Rowley
et al.
2009
;Schmid
2010
;Schmid
et al.
2011
).
However, for statistically stationary flows, spectral
proper
orthog
o
nal
decom
po
si
tion
(SPOD) (Lumley
1967
,
1970
; Citriniti & George
2000
; Picard & Delville
2000
; Towne,
Schmidt & Colonius
2018
) leads to an optimal reconstruction of the
space–time statistics
and results in modes that oscillate at a single frequency. A fundamental assumption
required in both space-only proper orthogonal decomposition (POD) and SPOD is
statistical stationarity, meaning that the statistics are
time
invari
ant. This assumption is
appropriate for many unforced flows. However, when forced, this fundamental assumption
is no longer valid as the
flow and its
statis
tics are now correlated to the forcing. To
clarify, by forcing, we mean any system that exhibits a periodic modulation of the statistics
(including by actuators, vortex-shedding in bluff body flows, rotation in turbomachinery,
etc
.). Several works have developed extensions to SPOD to study forced turbulent flows.
Franceschini
et al.
(
2022
) studied flows where a high-frequency turbulent component
develops on a low-frequency periodic motion. Subsequently, a quasi-steady assumption
is
made and conditionally fixed coherent structures at each phase are determined. Glezer,
Kadioglu & Pearlstein (
1989
) developed an extended POD method for flows with periodic
statistics by summing an ensemble of time series. However, since this method is based on
POD, it still contains the shortcomings present in POD. Heidt
et al.
(
2021
) applied SPOD
to the residual component of the triply decomposed fields (Hussain & Reynolds
1970
,
1972
) to isolate the impact of the forcing on the turbulence but still required a stationary
assumption. Clearly, SPOD and the aforementioned extensions are not sufficient to study
forced turbulent flows. This motivates an extension of SPOD to these flows, which is the
primary focus of this paper which we achieve by leveraging cyclostationary analysis.
Cyclostationary analysis is an extension to statistically stationary analysis to processes
with periodic statistics that
have been applied in a range of fields (Gardner
2018
), from
economics to physics and mechanics. Initially developed by Gudzenko (
1959
), Lebedev
(
1959
)
and Gladyshev (
1963
), it was then extensively studied and popularized
by Hurd
(
1969
) and Gardner (
1972
). The theory of second-order cyclostationary processes was
further developed by Boyles & Gardner (
1983
) and Gardner (
1986
b
), while Brown (
1987
)
and Gardner (
1986
c
) furthered the theory of complex-valued processes. Cyclostationary
analysis provides a robust statistical theory to study these processes, and tools analogous
to those used to study stationary processes (e.g. the mean, cross-correlation, cross-spectral
density, etc
.) have been developed which naturally collapse back to their stationary
counterparts when
analysing a stationary process.
Kim, North & Huang (
1996
) developed cyclostationary empirical orthogonal-functions
(CSEOFs) that essentially extends SPOD to cyclostationary processes for one-dimensional
data. Kim & North (
1997
) modified this technique to include multi-dimensional data
by reducing the computational cost through several approximations. However, due to
a lack of clarity in the literature regarding the derivation, properties,
inter
pre
ta
tion
and computation of these techniques, their use has been limited. Furthermore, despite
the aforementioned approximations, both formulations are computationally intractable
for high-dimensional data. In this paper, we extend SPOD to flows with time-periodic
985
A42-2
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
Decomposition of harmonically forced flows
statistics through an extension to the exact form of CSEOFs (Kim
et al.
1996
) to include
large multi-dimensional data. We hereafter refer to this method as cyclostationary SPOD
(CS-SPOD for short).
Methods used to model coherent structures are also considered. Specifically, we consider
resolvent analysis (also known as input/output analysis), where one seeks forcing modes
that give rise to the most amplified response modes with respect to their energetic gain.
When applied to turbulent fluid flows, the nonlinear modal interactions are regarded as
forcing terms to the linearized time-averaged turbulent mean (McKeon & Sharma
2010
).
Resolvent analysis has been used to study a wide range of transitional and turbulent flows
(Cossu, Pujals & Depardon
2009
; McKeon & Sharma
2010
; Meliga, Pujals & Serre
2012
;
Sharma & McKeon
2013
; Oberleithner, Paschereit & Wygnanski
2014
; Jeun, Nichols
& Jovanovi
́
c
2016
;Schmidt
et al.
2018
), amongst others. Towne
et al.
(
2018
) provided
a theoretical connection between SPOD and resolvent, showing that resolvent output
modes equal SPOD modes when the resolvent forcing modes are mutually uncorrelated.
This provides a theoretical basis to use resolvent analysis to develop models of the
space–time statistics of a turbulent flow (Moarref
et al.
2013
; Towne, Lozano-Durán &
Yang
2020
; Amaral
et al.
2021
) and the development of various methods (Morra
et al.
2019
;Pickering
et al.
2021
a
) to help whiten the forcing coefficients, thereby improving
these models. Resolvent analysis was extended to flows with a time-periodic mean flow
by
Padovan, Otto & Rowley (
2020
) and Padovan & Rowley (
2022
)
, and is termed harmonic
resolvent analysis. This leads to a system of
fre
quency
-
cou
pled equations that provide
the ability to study the first-order triadic interactions present in these time-periodic
flows. Analogous to the relationship between SPOD and resolvent analysis, in the present
paper, we establish a theoretical connection between CS-SPOD and harmonic resolvent
analysis.
The remainder of the paper is organized as follows. Section
2
introduces the theory
of cyclostationary processes and reviews an algorithm to compute their statistics. In §
3
,
CS-SPOD is derived, its properties
explored and an efficient computational algorithm is
proposed. After validating the method in §
4
, we demonstrate
the utility of CS-SPOD in
§
5
.In§
6
, we explore the relationship between CS-SPOD and the harmonic resolvent
analysis. Finally, in §
7
, we investigate low- and high-frequency forcing limits. Section
8
concludes the manuscript and summarizes the main points.
2. Cyclostationary theory
This section provides an overview of the theory of cyclostationary analysis and the tools
used to study them, with a focus on fluid dynamics. Comprehensive reviews can be found
from Gardner, Napolitano & Paura (
2006
), Antoni (
2009
)
and Napolitano (
2019
).
A complex-valued scalar process
q
(
t
)
at time
t
is cyclostationary in the wide sense if its
mean and autocorrelation functions are periodic with period
T
0
(Gardner
1986
b
), giving
E
{
q
(
t
)
}=
E
{
q
(
t
+
T
0
)
}
,
(2.1
a
)
R
(
t
,τ)
=
R
(
t
+
T
0
,τ),
(2.1
b
)
where
E
{·}
is the expectation operator,
R
is the autocorrelation
func
tion and
τ
is a
time
delay. Since the mean and autocorrelation are
time
peri
odic, they can be expressed as a
Fourier series
985
A42-3
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
L. Heidt and T. Colonius
E
{
q
(
t
)
}=

k
α
=−∞
ˆ
q
k
α
α
0
e
i2
π
(
k
α
α
0
)
t
,
(2.2
a
)
R
(
t
,τ)
E
{
q
(
t
+
τ/
2
)
q
(
t
τ/
2
)
}=

k
α
=−∞
ˆ
R
k
α
α
0
(τ )
e
i2
π
(
k
α
α
0
)
t
,
(2.2
b
)
where
k
α
Z
and the Fourier series coefficients are given by
ˆ
q
k
α
α
0
1
T
0

T
0
/
2
T
0
/
2
E
{
q
(
t
)
}
e
i2
π
(
k
α
α
0
)
t
d
t
,
(2.3
a
)
ˆ
R
k
α
α
0
(τ )
1
T
0

T
0
/
2
T
0
/
2
R
(
t
,τ)
e
i2
π
(
k
α
α
0
)
t
d
t
,
(2.3
b
)
where
α
0
=
1
/
T
0
is the fundamental cycle frequency and
(
·
)
is the complex conjugate.
The Fourier coefficients
ˆ
R
k
α
α
0
(τ )
are known as the cyclic autocorrelation functions of
q
(
t
)
at cycle frequency
k
α
α
0
. If a process contains non-zero
ˆ
q
k
α
α
0
and/or
ˆ
R
k
α
α
0
(τ )
,itissaid
to exhibit first- and second-order cyclostationarity at cycle frequency
k
α
α
0
, respectively.
Wide-sense stationary processes are the special case where
ˆ
R
k
α
α
0
(τ ) /
=
0for
k
α
=
0 only.
If the process
q
(
t
)
contains a deterministic periodic component at cycle frequency
k
α
α
0
,
it would exhibit both first-order and second-order (and any
higher
order) cyclostationarity
at cycle frequency
k
α
α
0
. Thus, a deterministic component results in a pure first-order
component and an impure (i.e. made up from components of a lower-order) second-order
(or higher) component (Antoni
et al.
2004
). Antoni
et al.
(
2004
) and Antoni (
2009
)showed
that in physical systems, it is crucial to
anal
yse the first- and second-order components
separately, where the second-order component
q

(
t
)
is defined as
q

(
t
)
q
(
t
)
E
{
q
(
t
)
}
,
(2.4)
such that
q
(
t
)
=
E
{
q
(
t
)
}+
q

(
t
)
and the mean
E
{
q
(
t
)
}=
E
{
q
(
t
+
T
0
)
}
is
T
0
periodic.
This approach makes physical sense considering that the first-order component is the
deterministic tonal component that originates from the forcing, while the second-order
component is a stochastic component that represents the underlying turbulence that is
modified by the forcing. The sequential approach is analogous to the triple decomposition
(Hussain & Reynolds
1970
,
1972
) where the underlying flow is separated into the
first-order (phase-averaged) and second-order (turbulent/residual) components. First-order
and second-order cyclostationarity then refer to a modulation of the first-order and
second-order components, respectively.
In this manuscript, we assume that all processes
anal
ysed using second-order analysis
tools are zero-mean processes (or have had their first-order component removed). Thus, by
stating that a process exhibits second-order cyclostationarity at cycle frequency
k
α
α
0
,we
mean that the process exhibits
pure second-order cyclostationarity at
k
α
α
0
.
We must clarify one point of terminology. Considering stationary processes are a subset
of cyclostationary processes, all stationary processes are also cyclostationary. We use the
most restrictive description, i.e. stationary processes are referred to as stationary and not
cyclostationary. By stating that a process exhibits cyclostationarity, we imply that at least
one cycle frequency
k
α
α
0
,
k
α
/
=
0 exists.
985
A42-4
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
Decomposition of harmonically forced flows
2.1.
Second-order cyclostationary analysis tools
In fluid dynamics, we are frequently interested in the correlation between two quantities.
Thus, we will now consider the complex-valued vector-valued process
q
(
x
,
t
)
at time
t
and independent variables (or spatial locations)
x
instead of the scalar process
q
(
t
)
.Two
processes are jointly cyclostationary if their cross-correlation function can be expressed as
a Fourier series, such that
R
(
x
,
x

,
t
,τ)
E
{
q
(
x
,
t
+
τ/
2
)
q
(
x

,
t
τ/
2
)
}
=

k
α
=−∞
ˆ
R
k
α
α
0
(
x
,
x

,τ)
e
i2
π
(
k
α
α
0
)
t
,
(2.5)
where the Fourier series coefficients are given by
ˆ
R
k
α
α
0
(
x
,
x

,τ)
1
T
0

T
0
/
2
T
0
/
2
R
(
x
,
x

,
t
,τ)
e
i2
π
(
k
α
α
0
)
t
d
t
,
(2.6)
and are known as the cyclic cross-correlation functions between
q
(
x
)
and
q
(
x

)
at cycle
frequency
k
α
α
0
. If the only non-zero cycle frequency is
k
α
α
0
=
0, then
q
(
x
)
and
q
(
x

)
are
jointly wide-sense stationary. Similar to the common assumption in stationary analysis,
we assume that all processes are separately and jointly cyclostationary.
A cyclostationary process can be
anal
ysed in the dual-frequency domain via the cyclic
cross-spectral density (CCSD). The CCSD is the generalization of the cross-spectral
density (CSD) for cyclostationary processes and is related to the cyclic cross-correlation
function via the cyclic
Wiener–Khinchin relation (Gardner & Robinson
1989
)
S
k
α
α
0
(
x
,
x

,
f
)
=

−∞
ˆ
R
k
α
α
0
(
x
,
x

,τ)
e
i2
π
f
τ
d
τ.
(2.7)
The CCSD can also be written as
S
k
α
α
0
(
x
,
x

,
f
)
lim

f
0
lim
T
→∞
1
T

T
/
2
T
/
2

fE

ˆ
q
1
/
f

x
,
t
,
f
+
1
2
k
α
α
0

ˆ
q
1
/
f

x

,
t
,
f
1
2
k
α
α
0

d
t
,
(2.8)
where
ˆ
q
W
(
x
,
t
,
f
)

t
+
W
/
2
t
W
/
2
q
(
x
,
t

)
e
i2
π
ft

d
t

is the short-time Fourier transform of
q
(
x
,
t
)
,
f
is the spectral
fre
quency and
k
α
α
0
is the cycle frequency. This shows that the
CCSD represents the time-averaged statistical correlation (with zero lag) of two spectral
components at frequencies
f
+
1
2
k
α
α
0
and
f
1
2
k
α
α
0
as the bandwidth approaches zero
(Napolitano
2019
). Formally, the CCSD is defined as
in (
2.8
) and then the Fourier
transform version (
2.7
) is proved, which is then known as the Gardner relation or
as the cyclic
Wiener–Khinchin relation (Napolitano
2019
). For
k
α
=
0, the CCSD
naturally reduces to the CSD, i.e.
S
0
(
x
,
x

,
f
)
. Correlation between spectral components
in cyclostationary processes is critical in the derivation of CS-SPOD, and for stationary
processes, the lack of correlation between spectral components is why SPOD can
anal
yse
each frequency independently.
The
Wigner–Ville (WV) spectrum (Martin
1982
; Martin & Flandrin
1985
; Antoni
2007
)
shows the spectral information of the process as a function of time (or phase) and, for a
985
A42-5
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
L. Heidt and T. Colonius
cyclostationary process, is given by
WV
(
x
,
x

,
t
,
f
)
=

k
α
=−∞
S
k
α
α
0
(
x
,
x

,
f
)
e
i2
π
(
k
α
α
0
)
t
.
(2.9)
The WV spectrum of the cyclic power-spectral density is determined by setting
x
=
x

, giving
WV
(
x
,
t
,
f
)
=
k
α
=−∞
S
k
α
α
0
(
x
,
f
)
e
i2
π
(
k
α
α
0
)
t
. While
non
-
phys
i
cal, the WV
spectrum may contain negative energy densities due to the negative interaction terms in the
WV spectrum (Flandrin
1998
; Antoni
2007
). However, Antoni (
2007
) showed this could
be arbitrarily reduced with increasing sampling time. The CCSD and WV spectrum can be
integrated with respect to frequency (Gardner
1994
; Randall, Antoni & Chobsaard
2001
),
which results in the instantaneous variance and the cyclic distribution of the instantaneous
variance, respectively
m
(
x
,
t
)
=
E
{
q
(
x
,
t
)
q
(
x
,
t
)
}=

−∞
WV
(
x
,
t
,
f
)
d
f
,
(2.10
a
)
ˆ
m
k
α
α
0
(
x
)
=

−∞
S
k
α
α
0
(
x
,
f
)
d
f
,
(2.10
b
)
where
m
(
x
,
t
)
is the mean-variance of the process and
ˆ
m
k
α
α
0
(
x
)
quantifies the
mean-variance contribution from each cycle frequency
k
α
α
0
.
So far, we have assumed that the cycle frequencies are known, but this may not always
be the case. To determine the cycle frequencies present in the system, all possible cycle
frequencies
α
are explored by rewriting the CCSD as
S
(
x
,
x

,α,
f
)
lim

f
0
lim
T
→∞
1
T

T
/
2
T
/
2

fE
ˆ
q
1
/
f
x
,
t
,
f
+
α
2
ˆ
q
1
/
f
x

,
t
,
f
α
2
d
t
.
(2.11)
A process exhibits cyclostationarity at cycle frequency
α
when
S
(
x
,
x

,α,
f
)/
=
0. For
cyclostationary processes, because the cross-correlation function is periodic, the spectral
correlation becomes discrete in
α
such that
S
(
x
,
x

,α,
f
)
=

k
α
=−∞
S
k
α
α
0
(
x
,
x

,
f
)δ(α
k
α
α
0
),
(2.12)
where
δ
is the Kronecker delta. The cyclic distribution of the instantaneous variance is
rewritten as
ˆ
m
(
x
,α)
=

−∞
S
(
x
,α,
f
)
d
f
,
(2.13)
which similarly becomes discrete for a cyclostationary process.
2.2.
Cycloergodicity
In fluid dynamics, it is laborious to require multiple realizations of a single process,
and we often invoke ergodicity in stationary processes to equate the ensemble average
with a long-time average of a single realization. We can similarly leverage the concept
of cycloergodicity as described
by Boyles & Gardner (
1983
), allowing us to replace the
985
A42-6
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
Decomposition of harmonically forced flows
expectation operator with a suitable time average, specifically, the cycle-averaging operator
(Braun
1975
)
̃
q
(
x
,
t
)
=
E
{
q
(
x
,
t
)
}=
lim
P
→∞
1
P
P

p
=
0
q
(
x
,
t
+
pT
0
),
(2.14)
where
̃
q
(
x
,
t
)
is the mean. The cycle-averaging operator is used when the data
are
phase-locked to the forcing (i.e. sampled at an integer number of samples per cycle) and
is identical to the phase-average used in the triple decomposition (Hussain & Reynolds
1970
; Reynolds & Hussain
1972
). As the cycle-averaging operator is periodic, it can be
expressed as a Poisson sum
̃
q
(
x
,
t
)
=
E
{
q
(
x
,
t
)
}=

k
α
=−∞
e
i2
π
(
k
α
α
0
)
t
lim
s
→∞
1
s

s
/
2
s
/
2
q
(
x
,
t

)
e
i2
π
(
k
α
α
0
)
t

d
t

.
(2.15)
This definition is employed for non-phase-locked data or to filter out first-order
components that are assumed to be statistical noise (Sonnenberger, Graichen & Erk
2000
;
Franceschini
et al.
2022
) and is identical to the harmonic-averaging procedure used by
Mezi
́
c(
2013
) and Arbabi & Mezi
́
c(
2017
) when restricted to a temporally periodic average.
2.3.
Computing the CCSD
There are practical considerations and nuances to computing the CCSD from discrete
data that we discuss in this section. Let the vector
q
k
C
N
represent a flow snapshot,
i.e. the instantaneous state of the process
q
(
x
,
t
)
at time
t
k
on a set of points in a spatial
domain
Ω
. The length of the vector
N
is equal to the number of spatial points multiplied
by the number of state variables. We assume that
these
data
are available for
M
equispaced
snapshots, with
t
k
+
1
=
t
k
+

t
. In addition, we assume that
these
data
are phase-locked,
meaning that there are an integer number of time steps in the fundamental period,
T
0
, and define
N
θ
=
T
0
/
t
. This restriction simplifies and reduces the computational
expense of the calculations but can in principle be relaxed by using the Poisson sum
time-average as in (
2.15
)andthe
non
-
com
pu
ta
tion
ally
effi
cient form of CS-SPOD shown
in
Algorithm 3
. Alternatively, non-phased-locked data can be temporally interpolated to
be phase-locked. Adopting similar notation to Towne
et al.
(
2018
), we estimate the CCSD
tensor
S
(
x
,
x

,α,
f
)
, which represents the spectral correlation between
q
(
x
,
t
)
and
q
(
x

,
t
)
at cycle frequency
α
and spectral frequency
f
. For a cyclostationary process,
S
(
x
,
x

,α,
f
)
is non-zero for
α
=
k
α
α
0
only, and therefore is written as
S
k
α
α
0
(
x
,
x

,
f
)
or equivalently
S
k
α
/
T
0
(
x
,
x

,
f
)
. The
space–time data can now be represented as the data matrix
Q
and
time vector
T
:
Q
=
[
q
1
,
q
2
,...,
q
M
]
C
N
×
M
,
(2.16)
T
=
[
t
1
,
t
2
,...,
t
M
]
R
M
.
(2.17)
Although we have a formula for the CCSD as seen in (
2.8
)and(
2.11
), this does not result
in a consistent estimator of the CCSD, as the variance of the
esti
mate of the CCSD does
not tend to zero as the amount of available data becomes large (Jenkins
1968
; Antoni
2007
;
Napolitano
2019
). Instead, this results in an estimate where the variance in the estimate is
equal to the squared value of the estimate itself. A consistent estimate of the CCSD can be
obtained by employing an appropriate averaging technique. The most common technique is
the time-averaging Welch method (Welch
1967
) due to its high computational efficiency.
985
A42-7
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
L. Heidt and T. Colonius
The Welch method averages a number of CCSDs to obtain a consistent
esti
mate of the
CCSD. From (
2.11
), we see that to compute the CCSD
, the Welch procedure is performed
on two frequency-shifted versions of the data, given by
Q
±
α/
2
=
Q
e
i2
π
(
±
α/
2
)
T
=
[
q
1
,
±
α/
2
,
q
2
,
±
α/
2
,...,
q
M
,
±
α/
2
]
,
(2.18
a
)
=
[
q
1
e
i2
π
(
±
α/
2
)
t
1
,
q
2
e
i2
π
(
±
α/
2
)
t
2
,...,
q
M
e
i2
π
(
±
α/
2
)
t
M
]
,
(2.18
b
)
where
q
k
,
±
α/
2
are the
±
1
2
α
frequency-shifted data matrices corresponding to the
k
th
snapshot, i.e.
q
k
,
±
α/
2
=
q
k
e
i2
π
(
±
α/
2
)
t
k
. Next, we split the two frequency-shifted data
matrices into a number of, possibly overlapping, blocks. Each block is written as
Q
(
n
)
±
α/
2
=
[
q
(
n
)
1
,
±
α/
2
,
q
(
n
)
2
,
±
α/
2
,...,
q
(
n
)
N
f
,
±
α/
2
]
C
N
×
N
f
,
(2.19)
where
N
f
is the number of snapshots in each block and the
k
th entry of the
n
th block
is
q
(
n
)
k
,
±
α/
2
=
q
k
+
(
n
1
)(
N
f
N
0
),
±
α/
2
. The total number of blocks,
N
b
,isgivenby
N
b
=

(
M
N
0
)/(
N
f
N
0
)
, where
represents the floor operator and
N
0
is the number
of snapshots that each block overlaps. The cycloergodicity hypothesis states that each of
these blocks is considered to be a single realization in an ensemble of realizations of this
cyclostationary flow. Subsequently, the discrete Fourier transform (DFT) of each block for
both frequency-shifted matrices is computed using a window
w
, giving
ˆ
Q
(
n
)
±
α/
2
=
[
ˆ
q
(
n
)
1
,
±
α/
2
,
ˆ
q
(
n
)
2
,
±
α/
2
,...,
ˆ
q
(
n
)
N
f
,
±
α/
2
]
,
(2.20)
where
ˆ
q
(
n
)
k
,
±
α/
2
=
1

N
f
N
f

j
=
1
w
j
q
(
n
)
j
,
±
α/
2
e
i2
π
(
k
1
)
[
(
j
1
)/
N
f
]
(2.21)
for
k
=
1
,...,
N
f
and
n
=
1
,...,
N
b
, where
ˆ
q
(
n
)
k
,
±
α/
2
is the
k
th Fourier component of
the
n
th block of the
±
α/
2 frequency-shifted data matrix, i.e.
f
k
,
±
α
0
/
2
. The nodal values
w
j
of a window function are
used to mitigate spectral and cyclic leakage arising from
the non-periodicity of the data within each block. Due to the
±
α/
2 frequency-shifting
applied, the
k
th discrete frequencies of the
±
α/
2 frequency-shifted data matrices represent
a frequency of
f
k
,
±
α/
2
=
f
k
±
α
2
=
k
1
N
f

t
for
k
N
f
/
2
,
k
1
N
f
N
f

t
for
k
>
N
f
/
2
,
±
α
2
.
(2.22)
This shows that the frequency components
f
k
+
α/
2and
f
k
α/
2, as required by (
2.11
),
have the same index
k
in the shifted frequency vectors
f
k
,
±
α/
2
. The CCSD tensor
S
(
x
,
x

,α,
f
)
is then estimated at cycle frequency
α
and spectral frequency
f
k
by
S
f
k
=

t
sN
b
N
b

n
=
1
ˆ
q
(
n
)
k
,α/
2
(
ˆ
q
(
n
)
k
,
α/
2
)
,
(2.23)
where
s
=
N
f
j
=
1
w
2
j
is the normalization constant that accounts for the difference in power
between the windowed and non-windowed signal. This is written compactly by arranging
985
A42-8
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
Decomposition of harmonically forced flows
Algorithm 1
Algorithm to compute the CCSD.
1:
for
Each data block,
n
=
1
,
2
,...,
N
b
do
Compute the frequency-shifted block-data matrices
2:
Q
(
n
)
±
α/
2
[
q
1
+
(
n
1
)(
N
f
N
0
),
±
α/
2
,
q
2
+
(
n
1
)(
N
f
N
0
),
±
α/
2
,...,
q
N
f
+
(
n
1
)(
N
f
N
0
),
±
α/
2
]
Using a (windowed) fast Fourier transform, calculate and store the row-wise
DFT for each frequency-shifted block-data matrix
3:
ˆ
Q
(
n
)
±
α/
2
=
FFT
(
Q
(
n
)
±
α/
2
)
=
[
ˆ
q
(
n
)
1
,
±
α/
2
,
ˆ
q
(
n
)
2
,
±
α/
2
,...,
ˆ
q
(
n
)
N
f
,
±
α/
2
]
The column
ˆ
q
(
n
)
k
,
±
α/
2
contains the
n
th realization of the Fourier mode
at the
k
th discrete frequency
f
k
,
±
α/
2
4:
end for
5:
for
Each frequency
k
=
1
,
2
,...,
N
f
(or some subset of interest)
do
Assemble the matrices of Fourier realizations from the
k
th column of each
ˆ
Q
(
n
)
±
α/
2
6:
ˆ
Q
f
k
,
±
α/
2
κ
[
ˆ
q
(
1
)
k
,
±
α/
2
,
ˆ
q
(
2
)
k
,
±
α/
2
,...,
ˆ
q
(
N
b
1
)
k
,
±
α/
2
,
ˆ
q
(
N
b
)
k
,
±
α/
2
]
Compute the CCSD at spectral frequency
f
k
and cycle frequency
α
7:
S
f
k
=
ˆ
Q
f
k
,α/
2
(
ˆ
Q
f
k
,
α/
2
)
.
8:
end for
the Fourier coefficients at the same index
k
into new frequency-data matrices
ˆ
Q
f
k
,
±
α/
2
=
κ
[
ˆ
q
(
1
)
k
,
±
α/
2
,
ˆ
q
(
2
)
k
,
±
α/
2
,...,
ˆ
q
(
N
b
1
)
k
,
±
α/
2
,
ˆ
q
(
N
b
)
k
,
±
α/
2
]
C
N
×
N
b
,
(2.24)
where
κ
=

t
/
sN
b
.
Then,
S
f
k
is estimated by
S
f
k
=
ˆ
Q
f
k
,α/
2
(
ˆ
Q
f
k
,
α/
2
)
.
(2.25)
This estimate converges, i.e. the bias and variance become zero, as
N
b
and
N
f
are
increased together (Welch
1967
; Antoni
2007
; Bendat & Piersol
2011
). The algorithm
to compute the CCSD from data snapshots is outlined in
Algorithm 1
, from which all
other second-order cyclostationary analysis tools can be computed. For efficient memory
management, variables assigned with ‘
’ can be deleted after each iteration in their
respective loop. Similar to the Welch estimate of the CSD, the estimate of the CCSD
suffers from the standard bias-variance trade-off, and caution should be taken to ensure
sufficiently converged statistics. In the CCSD, a phenomenon similar to spectral leakage
is present and is called cyclic leakage (Gardner
1986
a
) that results in erroneous cycle
frequencies. Using 67 % overlap when using a Hanning or Hamming window results in
excellent cyclic leakage minimization and variance reduction (Antoni
2007
). To reduce the
variance sufficiently,
T

f

1 is required (Antoni
2009
). If one does not know the cycle
frequencies
a
pri
ori
, one must search over all possible cycle frequencies with a resolution
=
1
/
T
(Gardner
1986
a
) to ensure all cycle frequencies are captured.
3. Cyclostationary SPOD
3.1.
Derivation
The objective of CS-SPOD is to find deterministic functions that best approximate,
on average, a zero-mean stochastic process. For clarity, we derive CS-SPOD using an
985
A42-9
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press
L. Heidt and T. Colonius
approach and notation analogous to the SPOD derivation presented
by Towne
et al.
(
2018
) and refer the reader to Brereton & Kodal (
1992
), Towne
et al.
(
2018
)
and Schmidt
& Colonius (
2020
) for detailed discussions on POD and SPOD. Like SPOD, we seek
deterministic modes that depend on both space and time such that we can optimally
decompose the
space–time statistics of the flow. Thus, we assume that each realization
of the stochastic process belongs to a Hilbert space with an inner product

q
1
,
q
2

x
,
t
=

−∞

Ω
q
2
(
x
,
t
)
W
(
x
)
q
1
(
x
,
t
)
d
x
d
t
,
(3.1)
where
q
1
(
x
,
t
),
q
2
(
x
,
t
)
are two realizations of the flow,
W
(
x
)
is a positive-definite
weighting tensor (
while we have chosen a time-independent weighting tensor since this
simplifies the derivations and is appropriate for the example cases shown, a time-periodic
weighting tensor could also be
employed), and
Ω
denotes the spatial domain of interest.
We then seek to maximize
λ
=
E
{|
q
(
x
,
t
),
φ
(
x
,
t
)

x
,
t
|
2
}

φ
(
x
,
t
),
φ
(
x
,
t
)

x
,
t
,
(3.2)
which leads to

−∞

Ω
R
(
x
,
x

,
t
,
t

)
W
(
x

)
φ
(
x

,
t

)
d
x

d
t

=
λ
φ
(
x
,
t
),
(3.3)
where
R
(
x
,
x

,
t
,
t

)
E
{
q
(
x
,
t
)
q
(
x

,
t

)
}
is the two-point
space–time correlation tensor.
Until this stage, no assumptions about the flow
have been made and
it is therefore identical
to the derivation of SPOD (Lumley
1967
,
1970
; Towne
et al.
2018
).
Since cyclostationary flows persist indefinitely, they have infinite energy in the
space–time norm, as shown in (
3.1
). Consequently, the eigenmodes of (
3.3
) do not possess
any of the useful quantities relied upon in POD or SPOD. To solve this, a new eigenvalue
decomposition is obtained in the spectral domain from which modes with the desired
properties are determined. We employ a solution ansatz of
φ
(
x
,
t
)
=

k
f
=−∞
ψ
(
x
+
k
f
α
0
)
e
i2
π
+
k
f
α
0
)
t
.
(3.4)
The set of frequencies present in the solution ansatz
φ
(
x
,
t
)
is called the
γ
set of solution
frequencies
Ω
γ
={
..., γ
2
α
0
α
0
,γ, γ
+
α
0
+
2
α
0
, ...
}
.
In
Appendix A
, we then use theory from §
2
to derive the infinite-dimensional CS-SPOD
eigenvalue problem, written compactly as

Ω
S
(
x
,
x

,γ)
W
(
x

)
Ψ
(
x

,γ)
d
x

=
λ
(γ )
Ψ
(
x
,γ),
(3.5)
985
A42-10
https://doi.org/10.1017/jfm.2024.70
Published online by Cambridge University Press