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 $\boldsymbol {q}(\boldsymbol {x}, t)$ at time $t$ and independent variables (or spatial locations) $\boldsymbol {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

(2.5)\begin{align} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, \tau) &\equiv E\{\boldsymbol{q}(\boldsymbol{x}, t+\tau/2)\boldsymbol{q}^*(\boldsymbol{x}^\prime, t - \tau/2)\}\nonumber\\ & =\sum_{k_\alpha={-}\infty}^\infty \hat{\boldsymbol{R}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, \tau) \,{\rm e}^{{\rm i}2{\rm \pi}(k_\alpha \alpha_0) t}, \end{align}

where the Fourier series coefficients are given by

(2.6)\begin{equation} \hat{\boldsymbol{R}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, \tau) \equiv \frac{1}{T_0} \int_{{-}T_0/2}^{T_0/2} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, \tau) \,{\rm e}^{-{\rm i}2{\rm \pi} (k_\alpha \alpha_0) t} \,{\rm d}t, \end{equation}

and are known as the cyclic cross-correlation functions between $\boldsymbol {q}(\boldsymbol {x})$ and $\boldsymbol {q}(\boldsymbol {x}^\prime )$ at cycle frequency $k_\alpha \alpha _0$. If the only non-zero cycle frequency is $k_\alpha \alpha _0 = 0$, then $\boldsymbol {q}(\boldsymbol {x})$ and $\boldsymbol {q}(\boldsymbol {x}^\prime )$ 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 analysed 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 Reference Gardner and Robinson1989)

(2.7)\begin{equation} {\boldsymbol{S}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, f) = \int_{-\infty}^\infty \hat{\boldsymbol{R}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, \tau) \,{\rm e}^{-{\rm i}2{\rm \pi} f\tau} \,{\rm d}\tau. \end{equation}

The CCSD can also be written as

(2.8) \begin{align} & {\boldsymbol{S}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, f) \nonumber\\ &\quad \equiv \lim_{\Delta f \rightarrow 0} \lim_{T \rightarrow \infty} \frac{1}{T} \int_{{-}T/2}^{T/2} \Delta f \,E\left\{\hat{\boldsymbol{q}}_{1/\Delta f} \left(\boldsymbol{x}, t, f + \frac{1}{2} k_\alpha \alpha_0\right) \hat{\boldsymbol{q}}^*_{1/\Delta f} \left(\boldsymbol{x}^\prime, t, f - \frac{1}{2} k_\alpha \alpha_0\right) \right\}\,{\rm d}t, \end{align}

where $\hat {\boldsymbol {q}}_{W} (\boldsymbol {x}, t, f) \equiv \int _{t - {W}/{2}}^{t + {W}/{2}} \boldsymbol {q}(\boldsymbol {x}, t^\prime ) \,{\rm e}^{-{\rm i}2{\rm \pi} f t^\prime } \,{\rm d}t^\prime$ is the short-time Fourier transform of $\boldsymbol {q}(\boldsymbol {x}, t)$, $f$ is the spectral frequency and $k_\alpha \alpha _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 + \frac {1}{2} k_\alpha \alpha _0$ and $f - \frac {1}{2} k_\alpha \alpha _0$ as the bandwidth approaches zero (Napolitano Reference Napolitano2019). 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 Reference Napolitano2019). For $k_\alpha = 0$, the CCSD naturally reduces to the CSD, i.e. $\boldsymbol {S}_0(\boldsymbol {x}, \boldsymbol {x}^\prime, 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 analyse each frequency independently.

The Wigner–Ville (WV) spectrum (Martin Reference Martin1982; Martin & Flandrin Reference Martin and Flandrin1985; Antoni Reference Antoni2007) shows the spectral information of the process as a function of time (or phase) and, for a cyclostationary process, is given by

(2.9)\begin{equation} \boldsymbol{WV}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, f) = \sum_{k_\alpha={-}\infty}^\infty \boldsymbol{S}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, f)\,{\rm e}^{{\rm i}2{\rm \pi} (k_\alpha \alpha_0) t}. \end{equation}

The WV spectrum of the cyclic power-spectral density is determined by setting $\boldsymbol {x} = \boldsymbol {x}^\prime$, giving $\boldsymbol {WV}(\boldsymbol {x}, t, f) = \sum _{k_\alpha =-\infty }^\infty \boldsymbol {S}_{k_\alpha \alpha _0}(\boldsymbol {x}, f)\,{\rm e}^{{\rm i}2{\rm \pi} (k_\alpha \alpha _0) t}$. While non-physical, the WV spectrum may contain negative energy densities due to the negative interaction terms in the WV spectrum (Flandrin Reference Flandrin1998; Antoni Reference Antoni2007). However, Antoni (Reference Antoni2007) showed this could be arbitrarily reduced with increasing sampling time. The CCSD and WV spectrum can be integrated with respect to frequency (Gardner Reference Gardner1994; Randall, Antoni & Chobsaard Reference Randall, Antoni and Chobsaard2001), which results in the instantaneous variance and the cyclic distribution of the instantaneous variance, respectively

(2.10a)$$\begin{gather} \boldsymbol{m}(\boldsymbol{x}, t) = E\{\boldsymbol{q}(\boldsymbol{x}, t)\boldsymbol{q}^*(\boldsymbol{x}, t)\} = \int_{-\infty}^\infty \boldsymbol{WV}(\boldsymbol{x}, t, f) \,{\rm d}f, \end{gather}$$

(2.10b)$$\begin{gather}\hat{\boldsymbol{m}}_{k_\alpha \alpha_0}(\boldsymbol{x}) = \int_{-\infty}^\infty \boldsymbol{S}_{k_\alpha \alpha_0}(\boldsymbol{x}, f) \,{\rm d}f, \end{gather}$$

where ${\boldsymbol {m}}(\boldsymbol {x}, t)$ is the mean-variance of the process and $\hat {\boldsymbol {m}}_{k_\alpha \alpha _0}(\boldsymbol {x})$ quantifies the mean-variance contribution from each cycle frequency $k_\alpha \alpha _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 $\alpha$ are explored by rewriting the CCSD as

(2.11)\begin{align} & {\boldsymbol{S}}(\boldsymbol{x}, \boldsymbol{x}^\prime, \alpha, f) \nonumber\\ &\quad \equiv \lim_{\Delta f \rightarrow 0} \lim_{T \rightarrow \infty} \frac{1}{T} \int_{{-}T/2}^{T/2} \Delta f\, E\left\{\hat{\boldsymbol{q}}_{1/\Delta f} \left(\boldsymbol{x}, t, f + \frac{\alpha}{2}\right) \hat{\boldsymbol{q}}^*_{1/\Delta f} \left(\boldsymbol{x}^\prime, t, f - \frac{\alpha}{2}\right) \right\}\,{\rm d}t. \end{align}

A process exhibits cyclostationarity at cycle frequency $\alpha$ when ${\boldsymbol {S}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \alpha, f) \ne 0$. For cyclostationary processes, because the cross-correlation function is periodic, the spectral correlation becomes discrete in $\alpha$ such that

(2.12)\begin{equation} \boldsymbol{S}(\boldsymbol{x}, \boldsymbol{x}^\prime, \alpha, f) = \sum_{k_\alpha={-}\infty}^\infty \boldsymbol{S}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, f) \delta (\alpha - k_\alpha \alpha_0), \end{equation}

where $\delta$ is the Kronecker delta. The cyclic distribution of the instantaneous variance is rewritten as

(2.13)\begin{equation} \hat{\boldsymbol{m}}(\boldsymbol{x}, \alpha) = \int_{-\infty}^\infty \boldsymbol{S}(\boldsymbol{x}, \alpha, f) \,{\rm d}f, \end{equation}

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 (Reference Boyles and Gardner1983), allowing us to replace the expectation operator with a suitable time average, specifically, the cycle-averaging operator (Braun Reference Braun1975)

(2.14)\begin{equation} \tilde{\boldsymbol{q}}(\boldsymbol{x}, t) = E\{\boldsymbol{q}(\boldsymbol{x}, t)\} = \lim_{P \rightarrow \infty} \frac{1}{P} \sum_{p = 0}^P \boldsymbol{q}(\boldsymbol{x}, t + pT_0), \end{equation}

where $\tilde {\boldsymbol {q}}(\boldsymbol {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 Reference Hussain and Reynolds1970; Reynolds & Hussain Reference Reynolds and Hussain1972). As the cycle-averaging operator is periodic, it can be expressed as a Poisson sum

(2.15)\begin{equation} \tilde{\boldsymbol{q}}(\boldsymbol{x}, t) = E\{\boldsymbol{q}(\boldsymbol{x}, t)\} = \sum_{k_\alpha={-}\infty}^\infty {\rm e}^{{\rm i}2{\rm \pi} (k_\alpha \alpha_0) t} \lim_{s\rightarrow \infty} \frac{1}{s}\int_{{-}s/2}^{s/2} \boldsymbol{q}(\boldsymbol{x}, t^\prime)\,{\rm e}^{-{\rm i}2{\rm \pi} (k_\alpha \alpha_0) t^\prime} \,{\rm d}t^\prime. \end{equation}

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 Reference Sonnenberger, Graichen and Erk2000; Franceschini et al. Reference Franceschini, Sipp, Marquet, Moulin and Dandois2022) and is identical to the harmonic-averaging procedure used by Mezić (Reference Mezić2013) and Arbabi & Mezić (Reference Arbabi and Mezić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 \in \mathbb {C}^N$ represent a flow snapshot, i.e. the instantaneous state of the process $\boldsymbol {q}(\boldsymbol {x}, t)$ at time $t_k$ on a set of points in a spatial domain $\varOmega$. 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 + \Delta 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_\theta =T_0/\Delta 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) and the non-computationally efficient form of CS-SPOD shown in . Alternatively, non-phased-locked data can be temporally interpolated to be phase-locked. Adopting similar notation to Towne et al. (Reference Towne, Schmidt and Colonius2018), we estimate the CCSD tensor $\boldsymbol {S}(\boldsymbol {x}, \boldsymbol {x}^\prime, \alpha, f)$, which represents the spectral correlation between $\boldsymbol {q}(\boldsymbol {x}, t)$ and $\boldsymbol {q}(\boldsymbol {x}^\prime, t)$ at cycle frequency $\alpha$ and spectral frequency $f$. For a cyclostationary process, $\boldsymbol {S}(\boldsymbol {x}, \boldsymbol {x}^\prime, \alpha, f)$ is non-zero for $\alpha = k_\alpha \alpha _0$ only, and therefore is written as $\boldsymbol {S}_{k_\alpha \alpha _0}(\boldsymbol {x}, \boldsymbol {x}^\prime, f)$ or equivalently $\boldsymbol {S}_{k_\alpha /T_0}(\boldsymbol {x}, \boldsymbol {x}^\prime, f)$. The space–time data can now be represented as the data matrix Q and time vector T:

(2.16) $$\begin{gather} \textbf{Q} = [\textbf{q}_1, \textbf{q}_2, \ldots, \textbf{q}_M] \in \mathbb{C}^{N\times M}, \end{gather}$$

(2.17)$$\begin{gather}\textbf{T} = [t_1, t_2, \ldots, t_M] \in \mathbb{R}^{M}. \end{gather}$$

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 estimate of the CCSD does not tend to zero as the amount of available data becomes large (Jenkins Reference Jenkins1968; Antoni Reference Antoni2007; Napolitano Reference Napolitano2019). 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 Reference Welch1967) due to its high computational efficiency. The Welch method averages a number of CCSDs to obtain a consistent estimate 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

(2.18a) \begin{align} \textbf{Q}_{{\pm} \alpha/2} &= \textbf{Q} \,{\rm e}^{-{\rm i} 2{\rm \pi} ({\pm} \alpha/2) \textbf{T}} = [{\textbf{q}}_{1, \pm \alpha/2},{\textbf{q}}_{2, \pm \alpha/2}, \ldots, {\textbf{q}}_{M, \pm \alpha/2} ], \end{align}

(2.18b) \begin{align} &= [{\textbf{q}}_{1} \,{\rm e}^{-{\rm i} 2{\rm \pi} ({\pm} \alpha/2) t_1}, {\textbf{q}}_{2} \,{\rm e}^{-{\rm i} 2{\rm \pi} ({\pm} \alpha/2) t_2}, \ldots, {\textbf{q}}_{M} \,{\rm e}^{-{\rm i} 2{\rm \pi} ({\pm} \alpha/2) t_M}], \end{align}

where ${\textbf{q}}_{k, \pm \alpha /2}$ are the $\pm \frac {1}{2}\alpha$ frequency-shifted data matrices corresponding to the $k$th snapshot, i.e. ${\textbf{q}}_{k, \pm \alpha /2} = {\textbf{q}}_{k} \,{\rm e}^{-{\rm i} 2{\rm \pi} (\pm \alpha /2) t_k}$. Next, we split the two frequency-shifted data matrices into a number of, possibly overlapping, blocks. Each block is written as

(2.19) \begin{equation} \textbf{Q}^{(n)}_{{\pm} \alpha/2} = [\textbf{q}_{1,\pm \alpha/2}^{(n)},\textbf{q}_{2,\pm \alpha/2}^{(n)}, \ldots, \textbf{q}_{N_f,\pm \alpha/2}^{(n)} ] \in \mathbb{C}^{N\times N_f}, \end{equation}

where $N_f$ is the number of snapshots in each block and the $k$th entry of the $n$th block is $\textbf{q}^{(n)}_{k,\pm \alpha /2} = \textbf{q}_{k+(n-1)(N_f - N_0), \pm \alpha /2}$. The total number of blocks, $N_b$, is given by $N_b = \lfloor {(M - N_0)}/{(N_f - N_0)}\rfloor$, where $\lfloor {\cdot } \rfloor$ 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

(2.20) \begin{equation} \hat{\textbf{Q}}^{(n)}_{{\pm} \alpha/2} = [\hat{\textbf{q}}_{ 1, \pm \alpha/2}^{(n)},\hat{\textbf{q}}_{2, \pm \alpha/2}^{(n)}, \ldots, \hat{\textbf{q}}_{N_f, \pm \alpha/2}^{(n)} ], \end{equation}

where

(2.21) \begin{equation} \hat{\textbf{q}}_{k, \pm \alpha/2}^{(n)} = \frac{1}{\sqrt{N_f}} \sum_{j = 1}^{N_f} {w_j}{\textbf{q}}_{j, \pm \alpha/2}^{(n)} \,{\rm e}^{-{\rm i}2{\rm \pi} (k-1)[(j - 1)/N_f]} \end{equation}

for $k = 1, \ldots, N_f$ and $n = 1, \ldots, N_b$, where $\hat {\textbf{q}}^{(n)}_{k, \pm \alpha /2}$ is the $k$th Fourier component of the $n$th block of the $\pm \alpha /2$ frequency-shifted data matrix, i.e. $f_{k, \pm \alpha _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 $\pm \alpha /2$ frequency-shifting applied, the $k$th discrete frequencies of the $\pm \alpha /2$ frequency-shifted data matrices represent a frequency of

(2.22) \begin{equation} f_{k, \pm \alpha/2} = f_{k} \pm \frac{\alpha}{2} = \begin{cases} \dfrac{k-1}{N_f\Delta t} & \text{for } k \leq N_f/2, \\ \dfrac{k-1 - N_f}{N_f \Delta t} & \text{for } k > N_f/2, \end{cases} \pm \frac{\alpha}{2}. \end{equation}

This shows that the frequency components $f_k + \alpha /2$ and $f_k - \alpha /2$, as required by (2.11), have the same index $k$ in the shifted frequency vectors $f_{k, \pm \alpha /2}$. The CCSD tensor $\boldsymbol {S}(\boldsymbol {x}, \boldsymbol {x}^\prime, \alpha, f)$ is then estimated at cycle frequency $\alpha$ and spectral frequency $f_k$ by

(2.23) \begin{equation} \textbf{S}_{f_k, \alpha} = \frac{\Delta t}{s N_b} \sum_{n =1}^{N_b} \hat{\textbf{q}}^{(n)}_{k, \alpha/2} (\hat{\textbf{q}}^{(n)}_{k, -\alpha/2})^*, \end{equation}

where $s = \sum _{j = 1}^{N_f}w_j^2$ is the normalization constant that accounts for the difference in power between the windowed and non-windowed signal. This is written compactly by arranging the Fourier coefficients at the same index $k$ into new frequency-data matrices

(2.24) \begin{equation} \hat{\textbf{Q}}_{f_k, \pm \alpha/2} = \sqrt{\kappa} [\hat{\textbf{q}}_{k, \pm \alpha/2}^{(1)}, \hat{\textbf{q}}_{k, \pm \alpha/2}^{(2)}, \ldots, \hat{\textbf{q}}_{k, \pm \alpha/2}^{(N_b-1)}, \hat{\textbf{q}}_{k, \pm \alpha/2}^{(N_b)}] \in \mathbb{C}^{N\times N_b}, \end{equation}

where $\kappa = {\Delta t}/{s N_b}$. Then, $\boldsymbol {S}_{f_k, \alpha }$ is estimated by

(2.25)\begin{equation} \textbf{S}_{f_k, \alpha} = \hat{\textbf{Q}}_{f_k, \alpha/2} (\hat{\textbf{Q}}_{f_k, -\alpha/2})^*. \end{equation}

This estimate converges, i.e. the bias and variance become zero, as $N_b$ and $N_f$ are increased together (Welch Reference Welch1967; Antoni Reference Antoni2007; Bendat & Piersol Reference Bendat and Piersol2011). The algorithm to compute the CCSD from data snapshots is outlined in , from which all other second-order cyclostationary analysis tools can be computed. For efficient memory management, variables assigned with ‘$\leftarrow$’ 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 Reference Gardner1986a) 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 Reference Antoni2007). To reduce the variance sufficiently, $T\Delta f \gg 1$ is required (Antoni Reference Antoni2009). If one does not know the cycle frequencies a priori, one must search over all possible cycle frequencies with a resolution $\Delta \alpha = 1/T$ (Gardner Reference Gardner1986a) to ensure all cycle frequencies are captured.

Algorithm 1 Algorithm to compute the CCSD.

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 approach and notation analogous to the SPOD derivation presented by Towne et al. (Reference Towne, Schmidt and Colonius2018) and refer the reader to Brereton & Kodal (Reference Brereton and Kodal1992), Towne et al. (Reference Towne, Schmidt and Colonius2018) and Schmidt & Colonius (Reference Schmidt and Colonius2020) 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

(3.1)\begin{equation} \langle \boldsymbol{q}_1, \boldsymbol{q}_2 \rangle_{x,t} = \int_{-\infty}^\infty \int_{\varOmega} \boldsymbol{q}_2^*(\boldsymbol{x}, t) \boldsymbol{W}(\boldsymbol{x}) \boldsymbol{q}_1(\boldsymbol{x}, t) \,\mathrm{d}\kern 0.06em \boldsymbol{x} \,\mathrm{d}t, \end{equation}

where $\boldsymbol {q}_1(\boldsymbol {x}, t),\ \boldsymbol {q}_2(\boldsymbol {x}, t)$ are two realizations of the flow, $\boldsymbol {W}(\boldsymbol {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 $\varOmega$ denotes the spatial domain of interest. We then seek to maximize

(3.2)\begin{equation} \lambda = \frac{E\{|\langle \boldsymbol{q}(\boldsymbol{x}, t), \boldsymbol{\phi}(\boldsymbol{x}, t)\rangle_{x,t}|^2\}}{\langle \boldsymbol{\phi}(\boldsymbol{x}, t), \boldsymbol{\phi}(\boldsymbol{x}, t)\rangle_{x,t}}, \end{equation}

which leads to

(3.3)\begin{equation} \int_{-\infty}^\infty \int_{\varOmega} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, t^\prime) \boldsymbol{W}(\boldsymbol{x}^\prime) \boldsymbol{\phi}(\boldsymbol{x}^\prime, t^\prime) \,\mathrm{d}\kern 0.06em \boldsymbol{x}^\prime \,\mathrm{d}t^\prime = \lambda \boldsymbol{\phi}(\boldsymbol{x}, t), \end{equation}

where $\boldsymbol {R}(\boldsymbol {x}, \boldsymbol {x}^\prime, t, t^\prime ) \equiv E\{\boldsymbol {q}(\boldsymbol {x}, t)\boldsymbol {q}^*(\boldsymbol {x}^\prime, t^\prime )\}$ 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 Reference Lumley1967, Reference Lumley1970; Towne et al. Reference Towne, Schmidt and Colonius2018).

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

(3.4)\begin{equation} \boldsymbol{\phi}(\boldsymbol{x}, t) = \sum_{k_f ={-}\infty}^\infty \boldsymbol{\psi}(\boldsymbol{x}, \gamma + k_f\alpha_0) \,{\rm e}^{{\rm i}2{\rm \pi} (\gamma + k_f\alpha_0)t}. \end{equation}

The set of frequencies present in the solution ansatz $\boldsymbol {\phi }(\boldsymbol {x}, t)$ is called the $\gamma$ set of solution frequencies $\mathcal {\varOmega }_{\gamma } = \{\ldots,\ \gamma - 2\alpha _0,\ \gamma - \alpha _0, \gamma,\ \gamma + \alpha _0,\ \gamma + 2\alpha _0,\ \ldots \}$.

In Appendix A, we then use theory from § 2 to derive the infinite-dimensional CS-SPOD eigenvalue problem, written compactly as

(3.5)\begin{equation} \int_{\varOmega} {{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x}, \boldsymbol{x}^\prime, \gamma) \boldsymbol{\mathsf{W}}(\boldsymbol{x}^\prime) \boldsymbol{\varPsi}(\boldsymbol{x}^\prime, \gamma)\, \mathrm{d}\kern 0.06em \boldsymbol{x}^\prime = \lambda(\gamma) \boldsymbol{\varPsi}(\boldsymbol{x}, \gamma), \end{equation}

where

(3.6a) \begin{align} &{{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x}, \boldsymbol{x}^\prime, \gamma) \nonumber\\ &\ =\begin{bmatrix} \ddots & \ddots & \ddots & \ddots & \ddots \\ \ddots & \boldsymbol{S}_0(\boldsymbol{x}, \boldsymbol{x}^\prime, \gamma - \alpha_0) & \boldsymbol{S}_{-\alpha_0}\left(\boldsymbol{x}, \boldsymbol{x}^\prime, \gamma - \dfrac{\alpha_0}{2} \right) & \boldsymbol{S}_{{-}2\alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma) & \ddots\\ \ddots & \boldsymbol{S}_{\alpha_0}\left(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma- \dfrac{\alpha_0}{2}\right) & \boldsymbol{S}_0(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma) & \boldsymbol{S}_{-\alpha_0}\left(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma+ \dfrac{\alpha_0}{2} \right) & \ddots\\ \ddots & \boldsymbol{S}_{2\alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma) & \boldsymbol{S}_{\alpha_0}\left(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma+ \dfrac{\alpha_0}{2} \right) & \boldsymbol{S}_0(\boldsymbol{x}, \boldsymbol{x}^\prime,\gamma+ \alpha_0) & \ddots\\ \ddots & \ddots & \ddots & \ddots & \ddots \end{bmatrix},\nonumber\\ \end{align}

(3.6b,c)\begin{gather} {{\boldsymbol{\mathsf{W}}}}(\boldsymbol{x}) =\begin{bmatrix} \ddots & & & & \\ & \hspace{-2mm}\boldsymbol{W}(\boldsymbol{x}) & & & \\ & & \hspace{-2mm}\boldsymbol{W}(\boldsymbol{x}) & & \\ & & & \hspace{-2mm}\boldsymbol{W}(\boldsymbol{x}) & \\ & & & & \hspace{-2mm}\ddots \end{bmatrix},\quad \boldsymbol{\varPsi}(\boldsymbol{x}, \gamma) =\begin{bmatrix} \vdots\\ \boldsymbol{\psi}(\boldsymbol{x},\gamma- \alpha_{0})\\ \boldsymbol{\psi}(\boldsymbol{x},\gamma)\\ \boldsymbol{\psi}(\boldsymbol{x},\gamma+ \alpha_{0})\\ \vdots \end{bmatrix}. \end{gather}

Here, ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ is the CS-SPOD decomposition tensor, ${{\boldsymbol{\mathsf{W}}}}(\boldsymbol {x})$ is the concatenated weight tensor and $\boldsymbol {\varPsi }(\boldsymbol {x}, \gamma )$ are the CS-SPOD eigenvectors. This frequency-domain version $\boldsymbol {\varPsi }(\boldsymbol {x}, \gamma )$ of the CS-SPOD eigenvectors can be converted to the time-domain version $\boldsymbol {\phi }(\boldsymbol {x}, t)$ using (3.4). In essence, we convert the original problem into the frequency domain and then solve for the Fourier series coefficients $\boldsymbol {\psi }(\boldsymbol {x}, f)$ at each $f \in \mathcal {\varOmega }_{\gamma }$.

This coupling of frequencies in CS-SPOD occurs because frequency components separated by $k_\alpha \alpha _0$ are correlated to each other, as shown in (2.8). In contrast, stationary processes do not exhibit correlation between frequencies, and thus each frequency can be solved independently via SPOD. Due to this coupling, CS-SPOD performed at $\gamma$ and $\gamma + \alpha _0$ solve the same problem, i.e. giving $\mathcal {\varOmega }_{\gamma } = \mathcal {\varOmega }_{\gamma +z \alpha _0}$, where $z \in \mathbb {Z}$. This means that CS-SPOD only contains unique solutions for the frequency sets corresponding to $\gamma \in \varGamma$, where $\varGamma = [-\alpha _0/2,\ \alpha _0/2)$.

In practice, the infinite-dimensional problem can not be solved. Instead, we must restrict the cycle frequencies considered and the frequencies present in the solution ansatz to some limit. We restrict the solution ansatz to

(3.7)\begin{equation} \boldsymbol{\phi}(\boldsymbol{x}, t) = \sum_{k_f ={-}K_f}^{K_f} \boldsymbol{\psi}(\boldsymbol{x}, \gamma + k_f\alpha_0) \,{\rm e}^{{\rm i}2{\rm \pi} (\gamma + k_f\alpha_0)t} \end{equation}

and the cycle frequencies to

(3.8)\begin{equation} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, \tau) = \sum_{k_\alpha ={-}K_\alpha}^{K_\alpha} \widetilde{\boldsymbol{R}}_{k_\alpha \alpha_0}(\boldsymbol{x}, \boldsymbol{x}^\prime, \tau) \,{\rm e}^{{\rm i}2{\rm \pi} (k_\alpha \alpha_0) t}, \end{equation}

where $K_f \in \mathbb {Z}^+$ and $K_\alpha \in \mathbb {Z}^+$. This gives a final solution frequency set of $\mathcal {\varOmega }_{\gamma } = \{-K_f\alpha _0 + \gamma, (-K_f+1)\alpha _0 + \gamma,\ \ldots, \gamma,\ \ldots, (K_f-1)\alpha _0+ \gamma, K_f\alpha _0+ \gamma \}$. In addition, the flow may only exhibit cyclostationarity at $K_\alpha$ harmonics of the fundamental cycle frequency. We employ identical notation to restrict the harmonics used to compute various second-order tools, such as the WV spectrum. These limits result in $2 K_f+1$ coupled equations, resulting in a $2K_f+1 \times 2K_f+1$ block eigensystem that is $2 K_\alpha +1$ banded-block-diagonal. In practice, $K_f$ should be chosen such that $\mathcal {\varOmega }_{\gamma }$ encompasses all frequencies of interest/importance, $K_\alpha$ should be chosen to encompass all the cycle frequencies present in the flow and $K_\alpha \le K_f$. An example for $K_f = 2, K_\alpha = 1$ is (for compactness, we have dropped the explicit dependence on $\boldsymbol {x}$ in this equation)

(3.9) \begin{align} {{\boldsymbol{\mathsf{S}}}}(\gamma) =\begin{bmatrix} \boldsymbol{S}_{0}(\gamma - 2\alpha_0) & \boldsymbol{S}_{-\alpha_0}\left(\gamma - \dfrac{3}{2}\alpha_0\right) & 0 & 0 & 0\\ \boldsymbol{S}_{\alpha_0}\left(\gamma - \dfrac{3}{2} \alpha_0\right) & \boldsymbol{S}_{0}(\gamma - \alpha_0) & \boldsymbol{S}_{-\alpha_0}\left(\gamma - \dfrac{1}{2} \alpha_0\right) & 0 & 0 \\ 0 & \boldsymbol{S}_{\alpha_0}\left(\gamma - \dfrac{1}{2} \alpha_0\right) & \boldsymbol{S}_{0}(\gamma ) & \boldsymbol{S}_{-\alpha_0}\left( \gamma + \dfrac{1}{2} \alpha_0\right) & 0 \\ 0 & 0 & \boldsymbol{S}_{\alpha_0}\left(\gamma + \dfrac{1}{2} \alpha_0\right) & \boldsymbol{S}_{0}(\gamma + \alpha_0) & \boldsymbol{S}_{-\alpha_0}\left(\gamma + \dfrac{3}{2} \alpha_0\right)\\ 0 & 0 & 0 & \boldsymbol{S}_{\alpha_0}\left(\gamma + \dfrac{3}{2}\alpha_0\right) & \boldsymbol{S}_{0}(\gamma + 2\alpha_{0}) \end{bmatrix}. \end{align}

In the limiting case that $K_\alpha = 0$, we obtain a block-diagonal CS-SPOD decomposition matrix where each diagonal block is the standard SPOD eigenvalue problem.

3.2. CS-SPOD properties

Since ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ is compact and finite, Hilbert–Schmidt theory guarantees a number of properties analogous to those for POD and SPOD (Lumley Reference Lumley1967, Reference Lumley1970; Towne et al. Reference Towne, Schmidt and Colonius2018). There is a countably infinite set of eigenfunctions $\boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma )$ at each unique frequency set $\mathcal {\varOmega }_{\gamma }$ that is orthogonal to all other modes at the same frequency set $\mathcal {\varOmega }_{\gamma }$ in the spatial inner norm $\langle \boldsymbol {q}_1, \boldsymbol {q}_2 \rangle _x = \int _{\varOmega } \boldsymbol {q}_2^*(\boldsymbol {x}, t) \boldsymbol {W}(\boldsymbol {x}) \boldsymbol {q}_1(\boldsymbol {x}, t) \mathrm {d}\kern 0.06em \boldsymbol {x}$, i.e. $\langle \boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma ), \boldsymbol {\varPsi }_{\!k}(\boldsymbol {x}, \gamma ) \rangle _x = \delta _{j,k}$. The following concatenated vector of each flow realization at the solution frequencies is optimally expanded as

(3.10a,b)\begin{equation} \hat{\boldsymbol{Q}}(\boldsymbol{x}, \gamma) = \begin{bmatrix} \vdots\\ \hat{\boldsymbol{q}}(\boldsymbol{x}, \gamma - \alpha_0)\\ \hat{\boldsymbol{q}}(\boldsymbol{x}, \gamma)\\ \hat{\boldsymbol{q}}(\boldsymbol{x}, \gamma + \alpha_0)\\ \vdots \end{bmatrix}, \quad \hat{\boldsymbol{Q}}(\boldsymbol{x}, \gamma) = \sum_{j = 1}^\infty a_j(\gamma) \boldsymbol{\varPsi}_{\!j}(\boldsymbol{x}, \gamma), \end{equation}

where $\hat {\boldsymbol {q}}(\boldsymbol {x}, f)$ is the temporal Fourier decomposition of each flow realization $\boldsymbol {q}(\boldsymbol {x}, t)$ at frequency $f$ and $a_j(\gamma ) = \langle \hat {\boldsymbol {Q}}(\boldsymbol {x}, \gamma ), \boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma )\rangle _{x}$ are the expansion coefficients, which are uncorrelated, i.e. $E\{a_j(\gamma ) a^*_k(\gamma )\} = \lambda _j(\gamma )\delta _{j, k}$.

Here, ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ is positive semi-definite, meaning that ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ has the following unique diagonal representation:

(3.11)\begin{equation} {{\boldsymbol{\mathsf{S}}}}(\boldsymbol{x}, \boldsymbol{x}^\prime, \gamma) = \sum_{j = 1}^\infty \lambda_j(\gamma) \boldsymbol{\varPsi}_{\!j}(\boldsymbol{x}, \gamma) \boldsymbol{\varPsi}^*_{j}(\boldsymbol{x}^\prime, \gamma), \end{equation}

in which the CS-SPOD modes are its principal components. This shows that CS-SPOD determines the modes that optimally reconstruct the second-order statistics, one frequency set $\mathcal {\varOmega }_{\gamma }$ at a time.

CS-SPOD modes are optimal in terms of their total energy reconstruction of ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ only. Thus, although each of the CCSDs present in ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$ have a diagonal representation, the individual components of $\boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma )$ are, in general, not orthogonal in the space norm, i.e. $\langle \boldsymbol {\psi }_j(\boldsymbol {x}, f), \boldsymbol {\psi }_k(\boldsymbol {x}, f) \rangle _x \ne \delta _{j,k}$. One exception is for stationary processes where the correlation between different frequency components is zero, resulting in a block-diagonal matrix where $\boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma )$ contains just a single non-zero $\boldsymbol {\psi }_j(\boldsymbol {x}, f)$ component, with $\langle \boldsymbol {\psi }_j(\boldsymbol {x}, f), \boldsymbol {\psi }_k(\boldsymbol {x}, f) \rangle _x = \delta _{j,k}$.

Transforming the eigenvectors $\boldsymbol {\varPsi }_{\!j}(\boldsymbol {x}, \gamma )$ back into the time domain, noting the ansatz defined in (3.4), gives $\boldsymbol {\phi }_{\gamma, j}(\boldsymbol {x}, t) = \sum _{k_f = -\infty }^\infty \boldsymbol {\psi }_j(\boldsymbol {x}, \gamma + k_f\alpha _0)\,{\rm e}^{{\rm i}2{\rm \pi} (\gamma + k_f\alpha _0) t}$, which are orthogonal in the space–time inner product integrated over a complete period. Thus, every mode occurring at each frequency set $\mathcal {\varOmega }_{\gamma }$ can be viewed as a unique space–time mode.

The two-point space–time correlation tensor can be written as

(3.12)\begin{equation} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, t^\prime) = \int_{-\alpha_0/2}^{\alpha_0/2} \sum_{j = 1}^\infty \lambda_j(\gamma) \boldsymbol{\phi}_{\gamma, j}(\boldsymbol{x}, t) \boldsymbol{\phi}^*_{\gamma, j}(\boldsymbol{x}^\prime, t^\prime) \,{\rm d}\gamma. \end{equation}

Substituting in the frequency expansion of $\boldsymbol {\phi }_{\gamma, j}(\boldsymbol {x}, t)$ and applying $t^\prime = t - \tau$ gives

(3.13)\begin{align} \boldsymbol{R}(\boldsymbol{x}, \boldsymbol{x}^\prime, t, \tau) &= \int_{-\alpha_0/2}^{\alpha_0/2} \sum_{j = 1}^\infty \lambda_j(\gamma)\nonumber\\ &\quad \times \sum_{k_f={-}\infty}^\infty \sum_{k_f^\prime={-}\infty}^\infty \boldsymbol{\psi}_{j}(\boldsymbol{x}, \gamma + k_f\alpha_0) \boldsymbol{\psi}^*_{j}(\boldsymbol{x}^\prime, \gamma + k_f^\prime\alpha_0) \,{\rm e}^{{\rm i}2{\rm \pi} (k_f - k_f^\prime) \alpha_0 t} \,{\rm e}^{{\rm i}2{\rm \pi} (\gamma + k_f^\prime \alpha_0) \tau} \,{\rm d}\gamma, \end{align}

resulting in a reconstruction that is time-periodic due to ${\rm e}^{{\rm i}2{\rm \pi} (k_f - k_f^\prime ) \alpha _0 t}$, which is why the ansatz defined by (3.4) was chosen.

In summary, for cyclostationary flows, CS-SPOD leads to modes that oscillate at a set of frequencies ($\mathcal {\varOmega }_{\gamma }$) and optimally represent the second-order space–time flow statistics.

3.3. Computing CS-SPOD modes in practice

We now detail how to compute CS-SPOD modes from data along with a technique to reduce the cost and memory requirements to levels similar to those of SPOD. A MATLAB implementation of the presented algorithms is available at https://github.com/CyclostationarySPOD/CSSPOD. Since the dimension of the CCSD is $N\times N$, the overall eigensystem $\boldsymbol{\mathsf{S}}_{\gamma _{k}}$ (which is the discrete approximation of ${{\boldsymbol{\mathsf{S}}}}(\boldsymbol {x}, \boldsymbol {x}^\prime, \gamma )$) becomes $(2K_f+1)N \times (2K_f+1)N$ in size. For common fluid dynamics problems, this can become a dense matrix ${O}(10^6-10^9) \times {O}(10^6-10^9)$ in size, which is computationally intractable to store in memory, let alone compute its eigendecomposition. This is also the dimension of the inversion required in the CSEOF methods by Kim et al. (Reference Kim, North and Huang1996) and Kim & North (Reference Kim and North1997). Thus, we derive a method-of-snapshots approach similar to the technique employed in POD (Sirovich Reference Sirovich1987) and SPOD (Citriniti & George Reference Citriniti and George2000; Towne et al. Reference Towne, Schmidt and Colonius2018) that reduces the size of the eigenvalue problem from $(2K_f+1)N \times (2K_f+1)N$ to $(2K_f+1)N_b \times (2K_f+1)N_b$. Since $N_b << N$, the method-of-snapshots technique makes the eigenvalue problem computationally tractable.

To determine CS-SPOD with a finite amount of discrete data, we substitute in the Welch computational procedure for the CCSD into each term of the frequency-limited version of (3.6a). We numerically evaluate this as

(3.14a,b) \begin{equation} \boldsymbol{\mathsf{S}}_{\gamma_{k}} = \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}} \widetilde{\boldsymbol{\mathsf{Q}}}^*_{\gamma_{k}}, \quad \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}} = \begin{bmatrix} \hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}\\ \vdots \\ \hat{\textbf{Q}}_{\gamma_{k}, 0}\\ \vdots \\ \hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0} \end{bmatrix}, \end{equation}

where

(3.15) \begin{equation} \hat{\textbf{Q}}_{\gamma_{k}, k_f\alpha_0} = \sqrt{\kappa} [\hat{\textbf{q}}_{k, k_f\alpha_0}^{(1)}, \hat{\textbf{q}}_{k, k_f\alpha_0}^{(2)}, \ldots, \hat{\textbf{q}}_{k, k_f\alpha_0}^{(N_b-1)}, \hat{\textbf{q}}_{k, k_f\alpha_0}^{(N_b)}] \in \mathbb{C}^{N\times N_b}. \end{equation}

Here, $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ is called the concatenated frequency-data matrix at the discrete $\mathcal {\varOmega }_{\gamma _{k}}$ set of solution frequencies and $\hat {\textbf{q}}_{k, k_f\alpha _0}^{(n)}$ is the $k$th DFT component of the $n$th block of the $k_f\alpha _0$ frequency-shifted data matrix. As stated previously, the solution frequency sets are only unique for $\gamma \in \varGamma$, thus the corresponding DFT frequencies are

(3.16)\begin{equation} \gamma_{k} = \begin{cases} \dfrac{k-1}{N_f\Delta t} & \text{for } k \leq \left\lfloor\dfrac{\alpha_0 N_f \Delta t}{2}\right\rfloor + 1,\\ \dfrac{k-1 - N_f}{N_f \Delta t} & \text{for } N_f - \left\lceil\dfrac{\alpha_0 N_f \Delta t}{2}\right\rceil + 1 < k \le N_f, \end{cases} \end{equation}

which forms the elements $\gamma _{k} \in \varGamma _{k}$. Expanding (3.14) gives

(3.17) \begin{align} &\boldsymbol{\mathsf{S}}_{\gamma_{k}} =\begin{bmatrix} \hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, 0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}^* \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \hat{\textbf{Q}}_{\gamma_{k}, 0}\hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, 0}\hat{\textbf{Q}}_{\gamma_{k}, 0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, 0}\hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}^* \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, - K_f\alpha_0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, 0}^* & \cdots & \hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}\hat{\textbf{Q}}_{\gamma_{k}, K_f\alpha_0}^* \end{bmatrix}. \end{align}

This expression shows that $\boldsymbol{\mathsf{S}}_{\gamma _{k}}$ contains off-diagonal terms that represent spectral correlations that are not present in the process (i.e. are not cycle frequencies considered in (3.8)). However, as $N_b$ and $N$ are increased together, this system converges and becomes a consistent estimate of (3.6a). Thus, all terms that represent spectral correlations not present in (3.8) converge to zero. Furthermore, the estimate is numerically positive semi-definite resulting in CS-SPOD modes that will inherit the desired properties. We note that for the numerical computation, one can not choose $K_\alpha$; instead, only $K_f$ is chosen and $K_\alpha = K_f$.

Equation (3.14) shows that the final eigenvalue problem can be compactly written as

(3.18a)$$\begin{gather} \boldsymbol{\mathsf{S}}_{\gamma_{k}} \boldsymbol{\mathsf{W}} \boldsymbol{\Psi}_{\gamma_{k}} = \boldsymbol{\Lambda}_{\gamma_{k}} \boldsymbol{\Psi}_{\gamma_{k}}, \end{gather}$$

(3.18b)$$\begin{gather}\widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}} \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}}^* \boldsymbol{\mathsf{W}} \boldsymbol{\Psi}_{\gamma_{k}} = \boldsymbol{\Lambda}_{\gamma_{k}} \boldsymbol{\Psi}_{\gamma_{k}}. \end{gather}$$

The spatial inner weight

(3.19)\begin{equation} \langle \boldsymbol{q}_1, \boldsymbol{q}_2 \rangle_x = \int_{\varOmega} \boldsymbol{q}_2^*(\boldsymbol{x}, t) \boldsymbol{W}(\boldsymbol{x}) \boldsymbol{q}_1(\boldsymbol{x}, t) \,\mathrm{d}\kern 0.06em \boldsymbol{x} \end{equation}

is approximated as $\langle \boldsymbol {q}_1, \boldsymbol {q}_2 \rangle _x = \textbf{q}_2^* \textbf{W} \textbf{q}_1$, where $\textbf{W} \in \mathbb {C}^{N\times N}$ is a positive-definite Hermitian matrix that accounts for both the weight and the numerical quadrature of the integral on the discrete grid and $\boldsymbol{\mathsf{W}} \in \mathbb {C}^{(2K_f+1)N \times (2K_f+1)N}$ is the block-diagonal matrix of $\textbf{W}$ (similar to (3.6b)). The CS-SPOD modes are then given by the columns of $\boldsymbol {\Psi }_{\gamma _{k}}$ and are ranked by their corresponding eigenvalues given by the diagonal matrix $\boldsymbol {\Lambda }_{\gamma _{k}}$. These discrete CS-SPOD modes hold analogous properties to all those previously discussed, including that they are discretely orthogonal $\boldsymbol {\Psi }^*_{\gamma _{k}} \textbf{W} \boldsymbol {\Psi }_{\!\gamma _{k}} = {{\boldsymbol{\mathsf{I}}}}$ and optimally decompose the estimated CS-SPOD decomposition matrix $\boldsymbol{\mathsf{S}}_{\gamma _{k}} = \boldsymbol {\Psi }_{\!\gamma _{k}} \boldsymbol {\Lambda }_{\gamma _{k}} \boldsymbol {\Psi }^*_{\gamma _{k}}$ (i.e. the second-order statistics).

At most, $\min (N, N_b )$ number of non-zero eigenvalues can be obtained. Thus, it is possible to show that the following $N_b \times N_b$ eigenvalue problem

(3.20) \begin{equation} \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}}^* \boldsymbol{\mathsf{W}} \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}} \boldsymbol{\Theta}_{\gamma_{k}} = \tilde{\boldsymbol{\Lambda}}_{\gamma_{k}} \boldsymbol{\Theta}_{\gamma_{k}} \end{equation}

contains the same non-zero eigenvalues as (3.18). This approach is known as the method-of-snapshots (Sirovich Reference Sirovich1987). The corresponding eigenvectors are exactly recovered as

(3.21) \begin{equation} \widetilde{\boldsymbol{\Psi}}_{\!\gamma_{k}} = \widetilde{\boldsymbol{\mathsf{Q}}}_{\gamma_{k}} \boldsymbol{\Theta}_{\gamma_{k}} \tilde{\boldsymbol{\Lambda}}_{\gamma_{k}}^{{-}1/2}. \end{equation}

Other than the simple weighting matrix $\boldsymbol{\mathsf{W}}$, only the concatenated data matrix $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ must be determined, which is easily achieved by computing each term ($\hat {\textbf{Q}}_{\gamma _{k}, k_f\alpha _0}$) in $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ using . Once $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ is determined, one computes $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}^* \boldsymbol{\mathsf{W}} \widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ and then performs the eigenvalue decomposition. Typically, only the first few modes are of physical interest, which allows us to employ a truncated decomposition where we determine a limited number of the most energetic CS-SPOD modes using randomized linear algebra methods (Martinsson & Tropp Reference Martinsson and Tropp2020). The total energy can be efficiently evaluated by taking the trace of $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}^* \boldsymbol{\mathsf{W}} \widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$. In Appendix B, we show a practical but computationally inefficient implementation of CS-SPOD. The algorithm requires computing $2K_f+1$ CCSDs, and thus the cost is approximately $2K_f + 1$ times that of the SPOD. The memory requirement scales similarly. This can be prohibitive when analysing large data sets.

However, significant savings are realized since all the terms in $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$ are in the form of $\hat {\textbf{Q}}_{\gamma _{k}, k_f\alpha _0}$, which represent the $k$th frequency component of the temporal Fourier transform of the $k_f\alpha _0$ frequency-shifted data matrix. The temporal Fourier transform of the $n$th realization of the $k_f\alpha _0$ frequency-shifted data is given by

(3.22a)\begin{align} \hat{\textbf{q}}_{k, k_f\alpha_0}^{(n)} &= \frac{1}{\sqrt{N_f}} \sum_{j = 1}^{N_f}w_j{\textbf{q}}_{j, k_f\alpha_0}^{(n)}\,{\rm e}^{-{\rm i}2{\rm \pi} (k-1)[{(j - 1)}/{N_f}]}, \end{align}

(3.22b)\begin{align} &= \frac{1}{\sqrt{N_f}} \sum_{j = 1}^{N_f}w_j{\boldsymbol{q}}_{j}^{(n)} \,{\rm e}^{-{\rm i} 2{\rm \pi} (k_f\alpha_0)[{(j-1) + (n-1)(N_f - N_0)}]\Delta t} \,{\rm e}^{-{\rm i}2{\rm \pi} (k-1)[{(j - 1)}/{N_f}]}, \end{align}

where ${\rm e}^{-{\rm i} 2{\rm \pi} (k_f\alpha _0)[{(j-1) + (n-1)(N_f - N_0)}]\Delta t}$ is the frequency-shifting operation. We separate these components into a phase-shifting component and a zero-phase-shift frequency-shifting component, by

(3.23a)$$\begin{gather} \hat{\textbf{q}}_{k, k_f\alpha_0}^{(n)} = {\rm e}^{-{\rm i} 2{\rm \pi} (k_f\alpha_0) \Delta t[{(n-1)(N_f - N_0)}]} \frac{1}{\sqrt{N_f}} \sum_{j = 1}^{N_f}w_j{\textbf{q}}_{j}^{(n)} \,{\rm e}^{-{\rm i}2{\rm \pi} (k_f\alpha_0\Delta t N_f + k-1)[{(j - 1)}/{N_f}]}, \end{gather}$$

(3.23b)$$\begin{gather}\hat{\textbf{q}}_{k, k_f\alpha_0}^{(n)} = {\rm e}^{-{\rm i} 2{\rm \pi} (k_f\alpha_0) \Delta t[{(n-1)(N_f - N_0)}]} \hat{\textbf{q}}_{\ell(k, k_f)}^{(n)}, \end{gather}$$

where $\ell (k, k_f)$ is the $\ell$th frequency that is a function of $k, k_f$ and will be defined in (3.26). This shows that the $f_k$ discrete frequency of the $k_f\alpha _0$-frequency-shifted data matrix ($\hat {\textbf{Q}}_{f_k,k_f\alpha _0}$) can be exactly computed as a phase-shifted version of the $\,f_{\ell (k, k_f)}$ discrete frequency component of the non-frequency-shifted data matrix ($\hat {\textbf{Q}}_{f_{\ell (k, k_f)}, 0}$). To employ this method, $k_f\alpha _0\Delta t N_f \in \mathbb {Z}$. Since $\alpha _0 \Delta T = 1/N_\theta$, where $N_\theta$ is the number of snapshots per fundamental period, this gives ${k_f N_f}/{N_\theta } \in \mathbb {Z}$, which requires $N_f = c N_\theta$, where $c \in \mathbb {Z}^+$ (i.e. there is a restriction on the length of each realization). This ensures that the change in frequency due to the frequency-shifting operator is equal to an integer change in the index of the frequency vector. With this restriction, the frequency spectrum of the DFT of an $N_f$ length record is

(3.24)\begin{equation} f_{k} = \begin{cases} \dfrac{(k-1)\alpha_0}{c } & \text{for } k \leq \frac{N_f}{2}, \\ \dfrac{(k-1 - N_f)\alpha_0}{c } & \text{for } k > \frac{N_f}{2}, \\ \end{cases} \end{equation}

and the unique frequency sets become

(3.25)\begin{equation} \gamma_{k} = \begin{cases} \dfrac{(k-1)\alpha_0}{c } & \text{for } k \leq \left\lfloor\dfrac{c }{2}\right\rfloor + 1,\\ \dfrac{(k-1 - N_f)\alpha_0}{c } & \text{for } N_f - \left\lceil\dfrac{c }{2} \right\rceil + 1 < k \le N_f. \end{cases} \end{equation}

This demonstrates that a frequency shift of $k_f\alpha _0$ corresponds to an integer change in the frequency index, i.e. the $k$th frequency component of the $k_f\alpha _0$-frequency-shifted data matrix corresponds to the phase-shifted version of the $\ell (k, k_f)$th frequency component ($\,f_{\ell (k, k_f)}$) of the non-frequency-shifted data matrix, i.e. $f_{k, k_f\alpha _0} = f_{\ell (k, k_f)}$, where

(3.26)\begin{equation} \ell(k, k_f) = \begin{cases} \begin{cases} k + k_f c & \text{for } k_f \ge 0 \\ k + k_f c + N_{f} & \text{for } k_f < 0 \end{cases} & \text{for } k \leq \left\lfloor\dfrac{c }{2}\right\rfloor + 1, \\ \begin{cases} k + k_f c - N_{f} & \text{for } k_f > 0 \\ k + k_f c & \text{for } k_f \le 0 \end{cases}& \text{for } N_f - \left\lceil\dfrac{c }{2} \right\rceil + 1 < k \le N_f. \end{cases} \end{equation}

This means that all the data required for CS-SPOD (for all frequency sets $\mathcal {\varOmega }_{\gamma _{k}}$) are contained within the Fourier transform of the original data matrix.

incorporates these savings and requires only a single DFT of the data matrix, making it similar in computational cost and memory requirement to SPOD. The memory usage to compute CS-SPOD for complex input data is $\approx ({1}/{(1 - N_0/N_f)} + 1)\times \text {mem}(\textbf{Q})$, which is the memory required to store the, possibly overlapping, block-data matrix and the original data matrix. Additional memory is required to store the temporary matrix $\widetilde {\boldsymbol{\mathsf{Q}}}_{\gamma _{k}}$, although the size of this matrix is minimal as typically $2K_f + 1 << N_f$. In extreme cases where only a single snapshot can be loaded at a time, a streaming CS-SPOD algorithm could be developed analogous to the streaming SPOD method by Schmidt & Towne (Reference Schmidt and Towne2019).

Algorithm 2 Efficient algorithm to compute CS-SPOD.

5.1. Application to a modified linearized complex Ginzburg–Landau equation

Our first example is the simple and well-understood linearized complex Ginzburg–Landau equation, which has been used as a model for a convectively unstable flow that exhibits non-modal growth (Chomaz, Huerre & Redekopp Reference Chomaz, Huerre and Redekopp1988; Hunt & Crighton Reference Hunt and Crighton1991; Cossu & Chomaz Reference Cossu and Chomaz1997). It can be written in the form of a generic linear forced system

(5.1)\begin{equation} \frac{\partial q(x, t)}{\partial t} - L(x, t) q(x, t) = f(x, t), \end{equation}

where $q(x, t)$ and $f(x, t)$ represent the state and forcing, respectively, with $|q(x\rightarrow \pm \infty, t)| \rightarrow 0$, and $L(x, t)$ is the linear operator

(5.2)\begin{equation} L(x, t) ={-}\nu_1\frac{\partial}{\partial x} + \nu_2 \frac{\partial^2}{\partial x^2} - \mu(x, t). \end{equation}

We use the commonly used form $\mu (x) = \mu _0 - c_{\mu }^2 + \mu _2 ({x^2}/{2})$ (Hunt & Crighton Reference Hunt and Crighton1991; Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009; Chen & Rowley Reference Chen and Rowley2011; Towne et al. Reference Towne, Schmidt and Colonius2018). All constants in (5.1) and (5.2), except for $\mu _0$, use the values of Bagheri et al. (Reference Bagheri, Henningson, Hoepffner and Schmid2009).

Similar to Franceschini et al. (Reference Franceschini, Sipp, Marquet, Moulin and Dandois2022), we construct periodic dynamics by using $\mu _0 = \bar {\mu }_0 + A_{\mu _0} \sin (2{\rm \pi} f_0 t)$, where $\bar {\mu }_0$ is the average value of $\mu _0$, $A_{\mu _0}$ is the amplitude of the periodic modulation of $\mu _0$ and $f_0$ is the frequency of the periodic modulation. For $A_{\mu _0} = 0$, the system has time-invariant dynamics, while for $|A_{\mu _0}| > 0$, the system has time-periodic dynamics, resulting in a stationary and cyclostationary response, respectively. By varying $A_{\mu _0}$, we modify the degree to which the system is cyclostationary. We choose $f_0 = 0.1$, which is substantial compared with the frequencies of interest (${\approx }[-0.5, 0.5]$), meaning that the quasi-steady approach of Franceschini et al. (Reference Franceschini, Sipp, Marquet, Moulin and Dandois2022) cannot be employed. Like Towne et al. (Reference Towne, Schmidt and Colonius2018), we use $\bar {\mu }_0 = 0.23$, which for $A_{\mu _0} = 0$, strongly amplifies external noise due to the non-normality of $L(x, t)$ and results in a degree of low-rankness typically present in turbulent flows. As per Franceschini et al. (Reference Franceschini, Sipp, Marquet, Moulin and Dandois2022), we confirm the stability of the system using Floquet analysis (results not shown). To demonstrate the utility of CS-SPOD and to facilitate its interpretation, we compare CS-SPOD performed at several levels of cyclostationarity $A_{\mu _0} = 0.0$, 0.2 and 0.4.

A pseudo-spectral approach using Hermite polynomials is employed to discretize the equations (Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009; Chen & Rowley Reference Chen and Rowley2011), where the collocation points $[x_1, x_2,\ldots,x_{N_H}]$ correspond to the first $N_H$ Hermite polynomials with scaling factor ${\rm Re}\{(-\mu _2/(2\nu _2))^\frac {1}{4}\}$. Following Bagheri et al. (Reference Bagheri, Henningson, Hoepffner and Schmid2009) and Towne et al. (Reference Towne, Schmidt and Colonius2018), we use $N_H = 221$, leading to a computational domain $x \in [-85.19, 85.19]$, which is large enough to mimic an infinite domain. The boundary conditions are implicitly satisfied by the Hermite polynomials (Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009). For CS-SPOD, the value of the weighting matrix at $x_i$ is determined as the distance between the midpoints of the neighbouring grid points. Temporal integration is performed using the embedded fifth-order Dormand–Prince Runge–Kutta method (Dormand & Prince Reference Dormand and Prince1980; Shampine & Reichelt Reference Shampine and Reichelt1997). After the initial transients have decayed, a total of 40 000 solution snapshots are saved with $\Delta t = 0.5$, giving a Nyquist frequency of $f_{\text {Nyquist}} = 1$.

To mimic a turbulent system, similar to Towne et al. (Reference Towne, Schmidt and Colonius2018), we force our system using spatially correlated band-limited noise, $f(x, t)$. This is performed by constructing spatially correlated noise with the following covariance kernel:

(5.3)\begin{equation} g(x, x^\prime) = \frac{1}{\sqrt{2{\rm \pi}}\sigma_\eta} \exp\left[ -\frac{1}{2} \left(\frac{x - x^\prime}{\sigma_\eta} \right)^2\right] \exp\left[ -{\rm i}2{\rm \pi} \frac{x - x^\prime}{\lambda_\eta}\right], \end{equation}

where $\sigma _\eta$ is the standard deviation of the envelope and $\lambda _\eta$ is the wavelength of the filter, such that the covariance $c(x, x^\prime ) = E\{\,f(x, t)f^*(x^\prime, t)\} = g(x, x^\prime )$. Spatial correlation is introduced by multiplying white noise by the Cholesky decomposition of the covariance kernel. The white noise has a uniformly distributed phase, normally distributed amplitude with unit variance and is generated as by Towne et al. (Reference Towne, Schmidt and Colonius2018). The forcing is spatially restricted to an interior portion of the domain via the window $\exp [-(x/L)^p]$, where $L = 60$, $p = 10$. The spatially correlated noise is low-pass filtered using a tenth-order finite-impulse-response filter with a cutoff frequency equal to $0.6f_{Nyquist}$. This results in a stationary forcing that is approximately constant in amplitude up to the cutoff frequency ($-$6 dB in amplitude at the cutoff frequency) but has non-zero spatial correlation as defined by (5.3). The forcing is then linearly interpolated to the temporal locations required by the temporal integration. To compute the WV spectrum, SPOD and CS-SPOD, we employed a window length $N_f = 10N_\theta$ and an overlap $67\,\%$, resulting in $N_b = 595$ (realizations) of the process and a frequency discretization of $\Delta f= 0.01$.

In analysing the data, we must first determine those frequencies, if any, where the system exhibits cyclostationarity. To do this, we compute the CCSD and search over all possible values of $\alpha$ in the range $\alpha \in [-1, 1]$, noting the $\alpha$ discretization required as discussed in § 2 to ensure no possible cycle frequencies are missed. Figure 4 shows the CCSD and integrated CCSD for the three values of $A_{\mu _0}$ at $x = x^\prime =0$, and confirms that the system is cyclostationary when $A_{\mu _0} > 0$ as high values of the CCSD and the integrated CCSD are seen at $\alpha = 0$, the modulation frequency ($\,f_0$) and an increasing number of harmonics as $A_{\mu _0}$ is further increased. This demonstrates that $\alpha _0 = f_0$.

Figure 4. The CCSD (top) and integrated CCSD (bottom) of the Ginzburg–Landau system at $x = 0$: (a) $A_{\mu _0} = 0$; (b) $A_{\mu _0} = 0.2$ and (c) $A_{\mu _0} = 0.4$.

We show 100 realizations of the process for each $A_{\mu _0}$ along with the WV spectrum at $x = x^\prime = 0$ as a function of the phase of $\alpha _0$ in figure 5. The WV spectrum is computed using $K_\alpha = 5$ to encompass all cycle frequencies present. Figure 5(a) shows that the statistics are almost constant as a function of phase for $A_{\mu _0} = 0$, which is expected given the time-invariant dynamics. The small degree of modulation observed is due to statistical uncertainty. In figure 5(b,c), we observe increasing levels of modulation in the statistics as $A_{\mu _0}$ increases. Furthermore, the peak value of the spectrum also increases due to the increasing non-normality of the system with increasing $\mu _0$. Given that the largest value of $\mu _0$ occurs at $\theta = 0.5{\rm \pi}$ and the peak of the WV spectrum occurs at $\theta \approx 0.95{\rm \pi}$, there is a phase delay of $\approx 0.45{\rm \pi}$ between when the dynamics of the system are the least stable and when the perturbations are, on average, the largest.

Figure 5. Example Ginzburg–Landau sample paths (top) at $x = 0$, where the red line shows a single representative trajectory for clarity and WV spectrum at $x= 0$ (bottom): (a) $A_{\mu _0} = 0$; (b) $A_{\mu _0} = 0.2$ and(c) $A_{\mu _0} = 0.4$.

Based on the preceding analysis and to ensure we encompass all frequencies of interest, we compute CS-SPOD using $K_f = 5$, resulting in a frequency range of $\mathcal {\varOmega }_{\gamma } = [-0.5, 0.5] + \gamma$. We first consider the stationary process with $A_{\mu _0} = 0.0$. Although CS-SPOD modes are theoretically equivalent to SPOD for the stationary case, finite data length leads to differences.

Figure 6 shows the SPOD eigenspectrum for $A_{\mu _0} = 0.0$. Note that the spectrum is not symmetric in $f$ because the Ginzburg–Landau system is complex. We superpose on the SPOD spectra the set of frequencies $f \in \mathcal {\varOmega }_{\gamma }$ for $\gamma = 0.05$, and mark and rank the six intersections with the highest energy. Based on the plot, we should find that the four most dominant CS-SPOD modes correspond to the dominant SPOD mode at a frequency of $\gamma - \alpha _0, \gamma, \gamma + \alpha _0$ and $\gamma + 2\alpha _0$, respectively. Similarly, the fifth and sixth CS-SPOD modes should correspond to the first subdominant SPOD modes at a frequency of $\gamma$ and $\gamma + \alpha _0$, respectively. Figure 7 makes comparisons between SPOD and CS-SPOD (performed assuming a fundamental cycle frequency of $\alpha _0 = f_0$) for the energy and eigenfunctions for each of these six modes. While the results are quite similar in each case, there are differences associated with statistics convergence, and this, as expected, occurs when there is a small energy separation between two distinct modes (e.g. modes 5 and 6).

Figure 6. The SPOD eigenspectrum of the Ginzburg–Landau system with $A_{\mu _0} = 0.0$, showing the three most energetic modes at each discrete frequency $f$. Every other eigenvalue has been omitted to improve readability. The six highest-energy modes occurring at the frequencies present in the CS-SPOD solution frequencies at $\gamma = 0.05$, i.e. $f \in \mathcal {\varOmega }_{\gamma }$, are depicted with the red dots.

Figure 7. Comparison of (a) SPOD and (b) CS-SPOD modes of the Ginzburg–Landau system with ${A_\mu = 0}$. From top to bottom are the six most dominant CS-SPOD modes and the six points identified in figure 6. The contour limits of the CS-SPOD eigenfunctions are set equal to the corresponding SPOD mode $\pm \|{\rm Re}\{\boldsymbol {\phi }_j(\boldsymbol {x}, t)\}\|_\infty$.

In figure 8(a), we now compare the CS-SPOD eigenspectrum for all $\gamma _{k} \in \varGamma _{k}$ for the three different values of $A_{\mu _0}$. As $A_{\mu _0}$ increases, so does the energy, as the disturbances are increasingly amplified by the increasing non-normality of the linear operator at phases corresponding to positive $A_{\mu _0} \sin (2{\rm \pi} f_0t)$, consistent with the trend shown previously in figure 5. A large energy separation between the dominant and sub-dominant CS-SPOD modes is observed, which increases for greater $A_{\mu _0}$, indicating that the process is increasingly low rank. In figure 8(b), for $\gamma = 0.05$, we show the fraction of the total energy ($\lambda _T = \sum _{j} \lambda _j$) that the first $J$ CS-SPOD or SPOD modes recover. As theoretically expected for $A_{\mu _0} =0$, CS-SPOD and SPOD result in an almost identical energy distribution. In contrast, with increasing $A_{\mu _0}$, CS-SPOD captures an increasingly greater amount of energy than SPOD. For example, at $A_{\mu _0} = 0.4$, the first CS-SPOD mode captures $64\,\%$ of the total energy, while the first SPOD mode captures just $45\,\%$. Furthermore, the first three CS-SPOD modes capture $92\,\%$ of the total energy, while seven SPOD modes are required to capture a similar amount of energy. As theoretically expected, the energy captured by SPOD does not exceed the energy captured by CS-SPOD (since SPOD modes are a subset of CS-SPOD modes). Thus, as the statistics become increasingly cyclostationary (i.e. more phase-dependent), CS-SPOD is able to capture an increasingly larger fraction of the phase-dependent statistics present in the process, which SPOD, due to the fundamentally flawed assumption of statistical stationarity, is unable to achieve. Due to the increased complexity of CS-SPOD modes, since they contain several frequency components, it is expected that they capture more energy. However, the critical difference is that SPOD is unable to capture the phase-dependent structure of the statistics (regardless of the number of modes employed).

Figure 8. (a) The CS-SPOD energy spectrum of the Ginzburg–Landau systems. (b) Total fractional energy captured by a truncated set of CS-SPOD and SPOD modes of the Ginzburg–Landau systems at $\gamma = 0.05$.

We now investigate how $A_{\mu _0}$ modifies the dominant CS-SPOD modes, at $\gamma = 0.05$, by showing the real component and the magnitude of the temporal evolution of the modes $\boldsymbol {\phi }_j(\boldsymbol {x}, t)$ in figures 9 and 10, respectively. We note that due to the multiple frequency components ($\mathcal {\varOmega }_{\gamma }$) present in $\boldsymbol {\phi }_j(\boldsymbol {x}, t)$, $\boldsymbol {\phi }_j(\boldsymbol {x}, t)$ can, unlike SPOD, no longer be completely represented by a single snapshot and instead must be displayed as a function of time. Similarly, the amplitude of the mode is periodic in time with period $T_0 = 1/\alpha _0$, unlike SPOD where the amplitude is constant in time. Thus, the amplitude is displayed as a function of phase $\theta$. Similar results are observed for other values of $\gamma$ not shown here. Overall, across all values of $A_{\mu _0}$, the real component of the CS-SPOD modes show a similar structure. However, as $A_{\mu _0}$ is increased, an additional modulation is seen that results in increasingly time/phase-dependent magnitudes.

Figure 9. Real component of the three dominant CS-SPOD modes of the Ginzburg–Landau systems at $\gamma = 0.05$. The contour limits for each CS-SPOD mode are $\pm \|{\rm Re}\{\boldsymbol {\phi }_j(\boldsymbol {x}, t)\}\|_\infty$.

Figure 10. Magnitude of the three dominant CS-SPOD modes of the Ginzburg–Landau systems at $\gamma = 0.05$. The contour limits for each CS-SPOD mode are $[0, \|\boldsymbol {\phi }_j(\boldsymbol {x}, \theta )\|_\infty ]$.

Finally, in figure 11, we investigate which frequency components are the most energetic via the fractional energy of each frequency component $f \in \mathcal {\varOmega }_{\gamma }$ for each CS-SPOD mode, defined as $E_{f, j} \equiv \boldsymbol {\psi }_j(\boldsymbol {x}, f)^* \boldsymbol {W}(\boldsymbol {x}) \boldsymbol {\psi }_j(\boldsymbol {x}, f)$, where $\sum _{f \in \mathcal {\varOmega }_{\gamma }} E_{f, j} = 1$. As $A_{\mu _0}$ increases, the CS-SPOD modes are constructed from an increasing number of non-zero-energy frequency components and at higher energy levels. For example, at $\gamma = 0.05$, the dominant frequency component, $f = 0.05$, contains $\approx$100 %, 83 % and 64 % of the total energy of the corresponding CS-SPOD mode for $A_{\mu _0} = 0$, 0.2 and 0.4, respectively. This occurs because of the increasing amount of correlation present between different frequency components as $A_{\mu _0}$ increases. Alternatively, this phenomenon can be understood as the following; as $A_{\mu _0}$ increases, the statistics become more time dependent, and thus, the amount of interaction between frequency components in $\mathcal {\varOmega }_{\gamma }$ increases such that the summation of these frequency components result in CS-SPOD modes that capture the time-periodic modulation experienced by the flow.

Figure 11. Fractional CS-SPOD modal energy by frequency, $E_{f, j}$, of the Ginzburg–Landau systems at $\gamma = 0.05$.

5.2. Forced turbulent jet

We now consider a forced turbulent, isothermal, subsonic jet for which data are available from a previous study (Heidt et al. Reference Heidt, Colonius, Nekkanti, Schmdit, Maia and Jordan2021). The large eddy simulation (LES) was computed using the Charles solver by Cascade Technologies using a set-up similar to previous, experimentally validated, simulations of turbulent jets (Brès et al. Reference Brès, Ham, Nichols and Lele2017, Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018). The jet has a Mach number of $M_j = U_j/c_j =0.4$ and a Reynolds number of $Re_j = {\rho _j U_j D}/{\mu _j} = 4.5\times 10^5$, where $\rho$ is the density, $\mu$ is the viscosity, $U$ is the velocity, $c$ is the speed of sound, $D$ is the nozzle diameter, and the subscripts j and $\infty$ represent the jet and free stream conditions, respectively. Frequencies are reported with respect to the Strouhal number $St = fD/U_j$, where $f$ is the frequency.

A schematic of the simulation set-up is shown in figure 12. An axisymmetric acoustic forcing is applied at a frequency $St_f = f_fD/U_j=0.3$ and amplitude $a_0/U_j = 0.1$. This forcing was chosen to roughly model the forced jet experiments of Crow & Champagne (Reference Crow and Champagne1971), and we chose $St_f=0.3$ to match what they observed as the frequency that led to the largest amplification by the flow (i.e. the jet preferred mode). We intentionally used a high amplitude of forcing as we wanted to clearly establish cyclostationarity in the resulting turbulence. The forcing is applied in an annular region surrounding the jet up to $r/D = 5$. The acoustic forcing inlet co-flow is defined by

(5.4a)$$\begin{gather} c(r)= 0.5 [1 - \text{erf}(2 (r - 5)) ], \end{gather}$$

(5.4b)$$\begin{gather}u_f(r, t) = c(r) \sin (2 {\rm \pi}f_f t), \end{gather}$$

(5.4c)$$\begin{gather}u_x(r, t) = u_{\infty} + a_0 u_f(r, t), \end{gather}$$

(5.4d)$$\begin{gather}u_r(r, t) = u_\phi(r, t) = 0, \end{gather}$$

(5.4e)$$\begin{gather}\rho(r, t) = \rho_\infty + \rho_\infty (u_x(r, t) - u_\infty)/a_\infty, \end{gather}$$

(5.4f)$$\begin{gather}p(r, t) = p_\infty + a_\infty\rho_\infty (u_x(r, t) - u_\infty). \end{gather}$$

Figure 12. Schematic of the forced Mach 0.4 turbulent jet, adapted from Brès et al. (Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018): (a) Overall domain and (b) nozzle region.

The simulation was run, post-transient, with a time step of $\Delta t D/c_\infty = 0.001$, for $480$ periods of the forcing frequency (or a total time of $t_{sim}D/c_\infty \approx 4000$), during which ${N_\theta = 48}$ snapshots were saved over each cycle of the forcing. For the subsequent analyses, the unstructured LES data were interpolated onto a structured cylindrical grid spanning $x/D \in [0, 30]$, $r/D \in [0, 6]$ and $\phi \in [0, 2{\rm \pi} ]$ with 656, 138 and 128 points in the axial, radial and circumferential dimensions. For the stochastic estimates, we use a window length $N_f = 12N_\theta$ and an overlap of $67\,\%$, resulting in $N_b = 118$ blocks and a non-dimensional frequency discretization of $\Delta St \approx 0.025$.

In figure 13, we plot the instantaneous and phase-averaged (2.14) velocity at four phases of one forcing cycle. Though not shown, we verified that the phase-averaged field is axisymmetric, consistent with the axisymmetric jet forcing. In the phase-averaged field, a large modulation in the axial velocity of the jet is observed with a vortex roll-up occurring around $x/D = 2.0$. The fundamental frequency fluctuation is primarily located in the potential core region and drives the large-scale periodic modulation. In figure 14, we extract the first four frequency components ($St = 0,\ 0.3,\ 0.6,\ 0.9$) of the phase-averaged field. The total fluctuation level, i.e. $2\times {\rm Re} \{\hat {u}_{x, f}/U_j \}$, for each non-zero frequency is $\approx$40 %, 15 % and 8 % thereby indicating that a substantial, nonlinear periodic modulation of the mean occurs. Harmonic generation similarly peaks near $x=2$ where the strong roll-up is occurring.

Figure 13. Top of each pair of images is $\tilde {u}_x(\theta )/U_j$ of the forced Mach 0.4 turbulent jet at $\theta = 0,\ {\rm \pi}/2,\ {\rm \pi},\ 3{\rm \pi} /2$. Bottom of each pair of images is ${u_x^{\prime \prime }(t)}/U_j$ at a time instant corresponding to a forcing phase of ${\theta = 0},\ {{\rm \pi} /2},\ {\rm \pi},\ {3{\rm \pi} /2}$.

Figure 14. ${\rm Re}\{\hat {u}_{x, St}/U_j \}$ of the forced Mach 0.4 turbulent jet at $St = 0$, 0.3, 0.6 and 0.9.

Next, we analyse the second-order stochastic component to determine the cycle frequencies present to apply CS-SPOD. Similar to the previous example, to determine what cycle frequencies are present in the flow, we interrogate the CCSD and integrated CCSD for $\alpha = [-3, 3]$ (not shown), again noting the $\alpha$ discretization required as discussed in § 2. We confirm that the only cycle frequencies present are harmonics of the forcing frequency (i.e. $\mathbb {Z} f_f$).

Figures 15(a,b) show the CCSD and corresponding WV spectrum, respectively, of the axisymmetric component of the axial velocity at $x/D = 5, r/D = 0.75$. For clarity, the CCSD is only shown for $\alpha /\alpha _0 \in \mathbb {Z}$ since all other values of $\alpha$ are $\approx 0$ (to within statistical convergence). A large modulation occurs for $\alpha /\alpha _0 = 0, \pm 1, \pm 2$. The WV spectrum shows this large modulation of the statistics, where the phase of the high-energy regions corresponds to when the high-velocity regions pass. Overall, it is clear that the forced turbulent jet exhibits cyclostationarity at frequencies equal to the harmonics of the forcing frequency.

Figure 15. (a) Absolute value of the CCSD of $u_x^{\prime \prime }/U_j$ of the forced Mach 0.4 turbulent jet at $x/D = 7$, $r/D = 0.75$. (b) The WV spectrum of $u_x^{\prime \prime }/U_j$ of the forced Mach 0.4 turbulent jet at $x/D = 7$, $r/D = 0.75$.

Finally, we demonstrate the utility of CS-SPOD on a forced turbulent jet. Recalling that both SPOD and CS-SPOD modes are decoupled among the azimuthal modes of the jet (owing to the statistical axisymmetry of the flow), we focus for brevity only on the axisymmetric $m=0$ component of the fluctuations. We seek modes that are orthogonal in the Chu-compressible energy norm (Chu Reference Chu1965) that have been applied in previous SPOD studies (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018):

(5.5) \begin{align} \langle\textbf{q}_{j}, \textbf{q}_{k}\rangle_{E}&=\iiint \textbf{q}_{k}^{*}\,\operatorname{diag}\left(\frac{\hat{T}_0}{\gamma_g \hat{\rho}_0 M^{2}}, \hat{\rho}_0, \hat{\rho}_0, \hat{\rho}_0, \frac{\hat{\rho}_0}{\gamma_g(\gamma_g-1) \hat{T}_0 M^{2}}\right) \textbf{q}_{j} r \,\mathrm{d}\kern 0.06em x \,\mathrm{d} r \,\mathrm{d} \phi\nonumber\\ & = \textbf{q}_k^* \textbf{W} \textbf{q}_j, \end{align}

where $M$ is the Mach number, $\gamma _g$ is the ratio of specific heats, $\hat {\rho }_0$ and $\hat {T}_0$ are the zero-frequency mean density and temperature components, and the matrix $\textbf{W}$ takes into account the energy and domain quadrature weights. To compute CS-SPOD, we choose $K_f = 10$, resulting in a non-dimensional solution frequency set of $\mathcal {\varOmega }_{\gamma _{St}} = [-3, 3] + \gamma _{St}$, which encompasses all frequencies of interest.

We show the CS-SPOD eigenspectrum of the turbulent jet in figure 16(a). A large difference in the energy between the 1st and 2nd, and 2nd and 3rd modes, at each $\gamma _{St}$, across almost all $\gamma _{St}$ is seen. This shows that the jet is low-rank. Since CS-SPOD solves for multiple frequencies at a time, the energy separation will be smaller than with SPOD, in particular, with a flatter spectrum. The spectrum peaks at $\gamma _{St} = 0.025$ and decays as $|\gamma _{St}| \rightarrow 0.15$ which, because the smallest $|St| \in \mathcal {\varOmega }_{\gamma _{St}}$ occurs at $|\gamma _{St}|$, occurs due to the decaying energy spectrum typically present in a turbulent jet. This low-rank behaviour, which is expected based on previous literature on natural turbulent jets (e.g. Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018), is observed in figure 16(b) where we show the fraction of the total energy captured by the first $J$ SPOD and CS-SPOD modes. The first CS-SPOD mode, at $\gamma _{St} = 0.15$, captures $10\,\%$ of the total energy present in the flow at the set of frequencies $\mathcal {\varOmega }_{\gamma _{St}}$, two modes capture $15.2\,\%$, ten modes capture $38\,\%$ and 50 modes capture $84.5\,\%$. Surprisingly, in contrast to the Ginsburg–Landau model, the energy separation between the most energetic CS-SPOD and SPOD modes is not large despite the high level of modulation present. However, despite this small difference, a large variation in the structure and temporal evolution of the most energetic SPOD and CS-SPOD modes is seen, which we explore next.

Figure 16. (a) The CS-SPOD energy spectrum of the forced Mach 0.4 turbulent jet. (b) Total fractional energy captured by a truncated set of CS-SPOD (blue) and SPOD (black) modes of the forced Mach 0.4 turbulent jet at $\gamma _{St} = 0.15$. The difference in the fractional energy captured by CS-SPOD and SPOD modes is also shown (red).

We show the real component and absolute value of the pressure component of the most energetic SPOD and CS-SPOD mode at $\gamma _{St} = 0.15$ in figure 17. The solid and dashed lines in these figures correspond to the contour lines of $\tilde {u}_x(\theta )/U_j = 0.25$, 0.75. SPOD modes are only shown at a single time instance due to their time-invariant evolution, while CS-SPOD modes are shown at several time instances to show their temporal evolution. The most dominant SPOD mode is focused downstream at $x/D \approx [6, 12]$, has a frequency $St = 0.15$ and has a structure typical of the so-termed ‘Orr modes’ previously observed in unforced turbulent jets (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2018; Pickering et al. Reference Pickering, Rigas, Nogueira, Cavalieri, Schmidt and Colonius2020). By construction, the amplitude of the SPOD mode remains constant over time, and the region of maximum amplitude corresponds to $x/D \approx [6, 12]$ and $r/D \approx [0,1]$. The real component of the most energetic CS-SPOD mode has a structure similar to the respective SPOD mode but with an additional modulation localized to the shear layer in regions of high velocity. This is also observed in the amplitude contours, where the amplitude of the mode substantially varies as a function of phase in a region similar to the amplitude profile of SPOD, but the high-amplitude regions always follow the high-velocity regions of the jet. The CS-SPOD modes follow this region since it is where the greatest amount of shear occurs along with the vortex roll-up (as seen in figures 13 and 14).

Figure 17. Comparison of the (a) real component and (b) magnitude of the pressure component of the dominant CS-SPOD mode to the dominant SPOD mode of the forced Mach 0.4 turbulent jet at $\gamma _{St} = 0.15$. All contours are set to $\pm 0.75\|{\rm Re}\{\phi _{p,1}({x}, {r}, t)\}\|_\infty$ and $[0, 0.75\|\phi _{p, 1}({x}, {r}, \theta )\|_\infty ]$ for the real and magnitude contours, respectively. The solid and dashed lines correspond to contour lines of $\tilde {u}_x(\theta )/U_j = 0.25$, 0.75, respectively.

Figure 18 shows the same CS-SPOD mode in a zoomed-in region near the nozzle exit, plotted with lower contour levels since the fluctuation levels are smaller there. At $t = 0$ (i.e. $\theta = 0$), a short wavelength Kelvin–Helmholtz (KH) mode that is located between the $25\,\%$ and $75\,\%$ velocity lines in the $x/D = [0,\ 1]$ region is seen. The KH mode is angled towards the centreline due to the modulation of the mean flow. Next, at $t = T_0/4$, the KH mode has propagated slightly downstream and has become significantly weaker due to the much thinner shear layer at this phase of the motion. From $t = T_0/4$ to $t = 3T_0/4$, the KH mode increases in strength as it continues to propagate downstream due to the increasing thickness of the boundary layer. The KH mode also rotates due to the roll-up induced by the forcing, as seen in figure 13. At $t = 3T_0/4$, the KH mode is substantially stronger than at $t = T_0/4$, is a lower-frequency structure located around the $x/D = [0.6,\ 1]$ region and is angled away from the centreline. A corresponding interrogation of the SPOD mode shows no near-nozzle KH activity at this frequency, highlighting the ability of CS-SPOD to reveal potentially important dynamical effects that are slaved to the forcing frequency.

Figure 18. Real component of the pressure of the dominant CS-SPOD mode of the forced Mach 0.4 turbulent jet at $\gamma _{St} = 0.15$ (zoomed into $x/D = [0,\ 2]$, $r/D = [0,\ 2]$). All contours are set to $\pm 0.25\|{\rm Re}\{\phi _{p, 1}(x= [0,\ 2], r = [0,\ 2], t) \}\|_\infty$. The solid and dashed lines correspond to contour lines of $\tilde {u}_x(\theta )/U_j = 0.25$, 0.75, respectively.

Figure 19(a) shows the normalized energy as a function of phase of the three dominant modes at $\gamma _{St} = 0.15$, defined as the spatial norm of the modes at each $\theta$. The energy, despite the large phase-dependent modulation seen in figure 17, varies by just $\pm$2 % as a function of phase. This demonstrates that, despite the strong phase-dependent structure of the mode and of the statistics present in the jet, on average, over the flow, the total energy contained within these modes is not strongly phase-dependent. Finally, in figure 19(b), we show the fractional energy of each frequency component $St \in \mathcal {\varOmega }_{\gamma _{St}}$ of the CS-SPOD modes. The large amount of frequency interaction previously observed is visible, where for $j = 1$, the eight highest energy frequency components are $\pm$0.15, $\pm$0.45, $\pm$0.75, $\pm$1.05 which contain 45.6 %, 3.6 %, 0.52 %, 0.15 % of the energy, respectively. Thus, a large amount of interaction occurs between the frequency components in $\mathcal {\varOmega }_{\gamma _{St}}$, which results in the large periodicity observed. It is important to note that although a frequency component may only contain a small fraction of the total energy in a CS-SPOD mode, in many cases, it is still a physically important feature, such as the modulated KH mode discussed previously, and thus should be carefully studied.

Figure 19. (a) Normalized energy of the dominant CS-SPOD modes of the forced Mach 0.4 turbulent jet over the phase of the external forcing at $\gamma _{St} = 0.15$. (b) Fractional CS-SPOD modal energy by frequency, $E_{f, j}$, of the forced Mach 0.4 turbulent jet at $\gamma _{St} = 0.15$, shown in $\text {log}_{10}$ scale.

Overall, we see that the forcing clearly results in a large modulation of the KH and Orr modes present, an effect that SPOD is unable to capture. Thus, the utility of CS-SPOD to describe the coherent structures in a forced turbulent jet is demonstrated.