1.1 Energy optimality or spectral purity?

The fundamental energy based decomposition is the proper orthogonal decomposition (POD). This decomposition originates from the Karhunen–Loéve (KL) theorem (Loève Reference Loève1977), which formulates the optimal (in a least-square sense) representation of a stochastic process as a series of orthogonal functions. This decomposition has set the basis for data processing methods which are known, depending on the field, as principal component analysis (PCA), Hotelling transform, empirical orthogonal function (EOF) or proper orthogonal decomposition (POD). Literature reviews on these applications can be found in the work of Cordier & Bergmann (Reference Cordier, Bergmann, David and Schram2013), and the monographs from Jackson (Reference Jackson1991) and Holmes et al. (Reference Holmes, Lumley, Berkooz and Rowley2012).

The POD for data-driven analysis of turbulent flows was introduced to the fluid mechanics community by Lumley (Reference Lumley, Yaglom and Tatarski1967, Reference Lumley1970), who focused on the spatial structures of the POD modes and their link to coherent structures in turbulent flows. The link between spatial and temporal structures was investigated by Aubry (Reference Aubry1991), who presented a variant of the decomposition, referred to as bi-orthogonal decomposition (BOD, Aubry, Guyonnet & Lima Reference Aubry, Guyonnet and Lima1991), which focused on the temporal structures of the modes. This link was made more explicit by Sirovich (Reference Sirovich1987, Reference Sirovich1989, Reference Sirovich1991), who proposed a simple algorithm, known as ‘snapshot POD’, to compute the entire decomposition.

It is nevertheless only thanks to the link between the POD and the well-known singular value decomposition (SVD) that the calculation of the POD modes of a discrete dataset could be formulated as a simple matrix factorization (see Kunisch & Volkwein Reference Kunisch and Volkwein1999; Chatterjee Reference Chatterjee2000; Fahl Reference Fahl2000), in contrast to formulations based on statistics or dynamical system theory (such as Berkooz Reference Berkooz1992; Berkooz, Holmes & Lumley Reference Berkooz, Holmes and Lumley1993; Holmes et al. Reference Holmes, Lumley, Berkooz, Mattingly and Wittenberg1997; Noack et al. Reference Noack, Afanasiev, Morzyński, Tadmor and Thiele2003). For an introduction to the SVD and a historical account on its development, the reader is referred to Stewart (Reference Stewart1993) and Golub (Reference Golub2013).

This simple matrix factorization framework has made the POD an extremely popular tool for both experimental and numerical fluid mechanics. Typical applications include the identification of coherent structures from experimental data (e.g. Schrijer, Sciacchitano & Scarano Reference Schrijer, Sciacchitano and Scarano2014; Mallor et al. Reference Mallor, Vila, Ianiro and Discetti2018), video analysis for adaptive masking or image pre-processing (Mendez et al. Reference Mendez, Raiola, Masullo, Discetti, Ianiro, Theunissen and Buchlin2017; Mendez, Scelzo & Buchlin Reference Mendez, Scelzo and Buchlin2018b ), flow sensing (Venturi & Karniadakis Reference Venturi and Karniadakis2004; Willcox Reference Willcox2006), flow control (Bergmann & Cordier Reference Bergmann and Cordier2008; Brunton & Noack Reference Brunton and Noack2015), reduced-order modelling (Deane et al. Reference Deane, Kevrekidis, Karniadakis and Orszag1991), experimental validation of computational fluid dynamics data (Blanchard et al. Reference Blanchard, Ng, Lowe and Vandsburger2014), data-driven identification of nonlinear systems (Brunton, Proctor & Kutz Reference Brunton, Proctor and Kutz2016b ; Loiseau, Noack & Brunton Reference Loiseau, Noack and Brunton2018) and more (see Pollard et al. Reference Pollard, Castillo, Danaila and Glauser2017). Moreover, extended versions of the decomposition, in which the dataset to decompose is constructed by assembling different quantities (e.g. velocities and concentrations), have also found application in the correlation analysis of coupled phenomena (Glezer, Kadioglu & Pearlstein Reference Glezer, Kadioglu and Pearlstein1989; Maurel, Borée & Lumley Reference Maurel, Borée and Lumley2001; Borée Reference Borée2003; Duwig & Iudiciani Reference Duwig and Iudiciani2009; Antoranz et al. Reference Antoranz, Ianiro, Flores and Villalba2018).

While the POD optimality guarantees the identification of the most energetic contributions within a few modes, the interpretation of these taken individually, and thus their usage for feature detection, can become difficult. As later illustrated with exemplary test cases, there exist situations in which different phenomena, possibly occurring at largely different frequencies, have very similar energy content. In such cases, an energy-based formulation cannot distinguish the various contributions, which consequently share the same POD structures. Frequency-based methods, conversely, assign to each mode a single frequency. Such methods are data-driven adaptations of the discrete Fourier transform (DFT) and infer the frequency of each mode from the data rather than defining it a priori (like in the DFT) from the time discretization.

The fundamental frequency-based and data-driven method is the dynamic mode decomposition (DMD). This decomposition originates from the Koopman theory (Koopman Reference Koopman1931), which consists in studying a nonlinear system by mapping it onto a linear one of much larger size (generally infinite-dimensional). This theory was then used by Mezić (Reference Mezić2005) to show how finite-dimensional approximations of such a Koopman (linear) dynamical system could be used for MOR applications. The extraction of such a finite approximation from the data was proposed by Rowley et al. (Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009) and Schmid (Reference Schmid2010). An introduction to the Koopman framework is proposed by Brunton et al. (Reference Brunton, Brunton, Proctor and Kutz2016a ), while extensive reviews on its connection to the DMD are proposed by Budišić, Mohr & Mezić (Reference Budišić, Mohr and Mezić2012), Mezić (Reference Mezić2013) and Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014). After the first two formulations, briefly reviewed in this work, many extensions of the algorithm have been developed; examples are the optimized DMD (Chen, Tu & Rowley Reference Chen, Tu and Rowley2012), the sparsity-promoting DMD (Jovanović, Schmid & Nichols Reference Jovanović, Schmid and Nichols2014), the randomized DMD (Erichson, Brunton & Kutz Reference Erichson, Brunton and Kutz2017) or the multi-resolution DMD (Kutz, Fu & Brunton Reference Kutz, Fu and Brunton2016b ). An overview of the DMD and its application is proposed by Kutz et al. (Reference Kutz, Brunton, Brunton and Proctor2016a ).

Besides the data-driven selection of the frequency of each mode, the advantage of the DMD over the DFT is to allow these to grow or decay exponentially, being each frequency complex. This makes the DMD suited for data-driven stability analysis (Schmid & Brandt Reference Schmid and Brandt2014) and feedback control design for linear systems (Rowley & Dawson Reference Rowley and Dawson2017).

The constraint of fully harmonic modes, on the other hand, can pose significant challenges to frequency-based decompositions for nonlinear datasets. A purely harmonic decomposition cannot easily represent frequency modulations, frequency and phase jitter or temporally localized events. A large number of modes required for the harmonic description of these phenomena in one domain (e.g. time) results in a substantial redundancy in the structures in the other domain (e.g. space). Moreover, since harmonics have infinite temporal support (they do not ‘start’ nor ‘finish’) their localization capabilities are particularly weak. Extending the frequency domain to the complex plane, as in the DMD, further amplifies these problems in the presence of nonlinear phenomena such as saturation of growth or decay rates.

The limits of harmonic bases have motivated the development of time–frequency analysis and wavelet theory (Mallat Reference Mallat2009; Kaiser Reference Kaiser2010), which has found many applications in the multi-scale modal analysis for fluid flows (Farge Reference Farge1992; Berkooz et al. Reference Berkooz, Elezgaray, Holmes, Lumley and Poje1994; De Souza et al. Reference De Souza, Delville, Lewalle and Bonnet1999; Rinoshika & Zhou Reference Rinoshika and Zhou2005; Schneider & Vasilyev Reference Schneider and Vasilyev2010). Wavelet decompositions are nevertheless not data-driven decomposition since one defines a priori the set of basis elements which have a given scale (frequency) and extension (duration). As described in § 3 and appendix A, the proposed mPOD borrows several ideas from the wavelet literature and combines them with the optimality constraints of the POD.

1.2 The need for hybrid decompositions

Focusing on data-driven methods, the limits of energy based (POD) or frequency based (DMD) approaches have been recently debated in the literature, and Noack (Reference Noack2016) has recently discussed the quest for hybrid decompositions. The first single-harmonic variant of the POD was proposed by Lumley (Reference Lumley1970) and formulated as a natural extension of the POD for stationary and homogeneous flows (see Glauser, Leib & George Reference Glauser, Leib and George1987; George Reference George and Kline1988). This approach consists in computing the POD structures as eigenvectors of the cross-spectral density matrix (see Citriniti & George Reference Citriniti and George2000; Gordeyev & Thomas Reference Gordeyev and Thomas2000; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011) and is referred to by several authors (e.g. Picard & Delville Reference Picard and Delville2000; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017; Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018) as spectral proper orthogonal decomposition (SPOD). Bourgeois, Noack & Martinuzzi (Reference Bourgeois, Noack and Martinuzzi2013) uses a pre-filtering of the data with a low-pass Morlet filter before computing the POD, to limit the frequency content of the resulting modes. Cammilleri et al. (Reference Cammilleri, Gueniat, Carlier, Pastur, Memin, Lusseyran and Artana2013) proposes a combination of POD and DMD, referred to as Cronos–Koopman analysis, which consists in performing DMD on the temporal basis of the POD, to force purely harmonic temporal modes. Noack et al. (Reference Noack, Stankiewicz, Morzyński and Schmid2016) uses a recursive approach for forcing the orthogonality in the DMD modes while minimizing the loss of spectral purity.

It is therefore evident that both the energy maximization and the spectral purity are too restrictive constraints, and an ideal decomposition should match the two approaches, possibly allowing for switching from one to the other. A decisive step towards this direction was achieved by Sieber, Paschereit & Oberleithner (Reference Sieber, Paschereit and Oberleithner2016), who proposed a novel decomposition blending the POD and the DFT. This decomposition was also named spectral proper orthogonal decomposition (SPOD), although it has nothing to do with the frequency-based POD that is also referred to as SPOD by Picard & Delville (Reference Picard and Delville2000), Taira et al. (Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017) and Towne et al. (Reference Towne, Schmidt and Colonius2018).

Sieber’s SPOD consists of using a low-pass filter along all the diagonals of the correlation matrix, before proceeding the Sirovich’s ‘snapshot POD’. The idea of such diagonal filtering arises from Szegö theorem (Grenander Reference Grenander2001; Gray Reference Gray2005) which states that the eigenbasis of a Toeplitz circulant matrix is a Fourier basis. Therefore, the more the correlation matrix resembles a circulant Toeplitz matrix, the more a POD resembles a DFT. A correlation matrix, however, approaches a Toeplitz circulant form only for a stationary process (Brockwell & Davis Reference Brockwell and Davis1987; Gray Reference Gray2005) and the diagonal filter used in the SPOD artificially forces such a pattern, thus forcing the resulting eigenvectors to approach the DFT basis. Depending on the filter size, one moves from a DFT (very large filter kernel) to a POD (very short filter kernel). Between these two limits, a good compromise between energy optimality and spectral purity can be obtained.

The main limitation of this method, however, is that this filtering operation can significantly alter the correlation matrix of non-stationary processes, compromising its symmetry and thus the orthogonality of its eigenvectors. Such operation can be seen as a smoothing of the nonlinearities (Sieber et al. Reference Sieber, Paschereit and Oberleithner2016; Sieber, Paschereit & Oberleithner Reference Sieber, Paschereit and Oberleithner2017), but the decision on the filter size is not supported by a mathematical link between the diagonal filtering and the spectral content of the resulting eigenvectors.

1.3 Contribution of this work and manuscript organization

Section 2 presents a general matrix factorization that is common to all the decompositions, and sets the POD, DFT and DMD in this generalized form. This framework is then used in § 3 to derive the novel multi-scale proper orthogonal decomposition (mPOD).

The mPOD proposed in this work combines ideas from the windowed SPOD in Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016) and the SPOD in Towne et al. (Reference Towne, Schmidt and Colonius2018). Instead of filtering the correlation matrix, as in Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016), the mPOD decomposes it into the contributions of different scales using the multi-resolution (MRA) architecture from wavelet theory. Instead of computing various eigenbases from the spectra of different portions of the data, as in Towne et al. (Reference Towne, Schmidt and Colonius2018), the mPOD computes eigenbases on the correlation matrix of different scales. The derivation of the mPOD in § 3 is carried out by analysing the impact of the filtering process on the spectra of the correlation matrix and the POD modes. Finally, the mPOD is tested and compared to POD, DFT and DMD on an illustrative synthetic test case (§ 4), a numerical test case (§ 5) and an experimental test case (§ 6). Conclusions and perspectives for future applications are presented in § 7.

2.1 Matrix form of the proper orthogonal decomposition (POD)

The POD temporal structures $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ are the eigenvectors of the temporal correlation matrix $\unicode[STIX]{x1D646}=\unicode[STIX]{x1D63F}^{\dagger }\unicode[STIX]{x1D63F}$ . If $\unicode[STIX]{x1D63F}$ is real, as in all the applications considered in this work, these are orthonormal and real vectors:

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D646}=\unicode[STIX]{x1D733}_{{\mathcal{P}}}\unicode[STIX]{x1D726}_{{\mathcal{P}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}}^{\text{T}}=\mathop{\sum }_{r=1}^{rk(\unicode[STIX]{x1D63F})}\unicode[STIX]{x1D706}_{{\mathcal{P}}_{r}}\unicode[STIX]{x1D74D}_{{\mathcal{P}}_{r}}\unicode[STIX]{x1D74D}_{{\mathcal{P}}_{r}}^{\text{T}}.\end{eqnarray}$$

Introducing this definition in (2.9) is particularly revealing: a comparison of (2.10) to the first factorization in (2.9) yields $\unicode[STIX]{x1D726}_{{\mathcal{P}}}=\unicode[STIX]{x1D72E}_{P}\unicode[STIX]{x1D731}_{{\mathcal{P}}}^{\text{T}}\unicode[STIX]{x1D731}_{{\mathcal{P}}}\unicode[STIX]{x1D72E}_{{\mathcal{P}}}$ . Since this eigenvalue matrix is diagonal, one recovers $\unicode[STIX]{x1D731}_{{\mathcal{P}}}^{\text{T}}\unicode[STIX]{x1D731}_{{\mathcal{P}}}=\unicode[STIX]{x1D644}$ and $\unicode[STIX]{x1D726}_{{\mathcal{P}}}=\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{2}$ . Therefore, also the spatial structures are real and orthonormal and thus eigenvectors of the spatial correlation matrix $\unicode[STIX]{x1D63E}=\unicode[STIX]{x1D63F}\unicode[STIX]{x1D63F}^{\dagger }=\unicode[STIX]{x1D731}_{{\mathcal{P}}}\unicode[STIX]{x1D726}_{{\mathcal{P}}}\unicode[STIX]{x1D731}_{{\mathcal{P}}}^{\text{T}}$ . Therefore, the energy contribution of each mode becomes $\unicode[STIX]{x1D70E}_{{\mathcal{P}}_{r}}=\sqrt{\unicode[STIX]{x1D706}_{r}}$ , and the energy of the dataset can be recovered from $\unicode[STIX]{x1D72E}_{{\mathcal{P}}}$ .

For the choice of temporal basis in (2.10), (2.5) yields the well-known singular value decomposition (SVD), here denoted as $\unicode[STIX]{x1D63F}=\unicode[STIX]{x1D731}_{{\mathcal{P}}}\unicode[STIX]{x1D72E}_{{\mathcal{P}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}}^{\text{T}}$ . The optimality of the decomposition is thus guaranteed by the Eckart–Young–Mirsky theorem, which states that the rank $r_{c}<rk(\unicode[STIX]{x1D63F})$ approximation of a matrix $\unicode[STIX]{x1D63F}$ obtained from a truncated SVD minimizes the norm of the error

(2.11) $$\begin{eqnarray}Err(r_{c})=\mathop{\left\Vert \unicode[STIX]{x1D63F}-\mathop{\sum }_{r=1}^{r_{c}}\unicode[STIX]{x1D70E}_{{\mathcal{P}}r}\unicode[STIX]{x1D753}_{{\mathcal{P}}r}\unicode[STIX]{x1D74D}_{{\mathcal{P}}r}^{\text{T}}\right\Vert }\nolimits_{2}=\unicode[STIX]{x1D70E}_{{\mathcal{P}}r_{c}+1}.\end{eqnarray}$$

Therefore, each coefficients $\unicode[STIX]{x1D70E}_{{\mathcal{P}}_{r}}$ represents also the $L^{2}$ error produced by a rank $r-1$ approximation, and the error decay $Err(r_{c})\rightarrow 0$ as $r_{c}\rightarrow rk(\unicode[STIX]{x1D63F})$ is the strongest possible.

The POD from Lumley (Reference Lumley, Yaglom and Tatarski1967, Reference Lumley1970) computes the spatial structures $\unicode[STIX]{x1D731}_{{\mathcal{P}}}$ from the eigendecomposition of $\unicode[STIX]{x1D63E}$ ; the BOD from Aubry (Reference Aubry1991) computes the temporal structures (referred to as BOD modes) $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ from the eigendecomposition of $\unicode[STIX]{x1D646}$ . The POD from Sirovich (Reference Sirovich1987, Reference Sirovich1989, Reference Sirovich1991) computes both $\unicode[STIX]{x1D731}_{{\mathcal{P}}}$ and $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ by first solving (2.10) (less demanding if $n_{t}\ll n_{s}$ ) and then using (2.6) (which reads $\unicode[STIX]{x1D731}_{{\mathcal{P}}}=\unicode[STIX]{x1D63F}\unicode[STIX]{x1D733}_{{\mathcal{P}}}\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{-1}$ ).

Concerning the impact of the mean removal before computing the POD from (2.7), the POD is equipped with a vector of constant in its temporal basis only if such a vector is an eigenvector of $\unicode[STIX]{x1D646}$ . This implies:

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D646}\text{}\underline{1}=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D707}}^{2}\text{}\underline{1}\quad \Longleftrightarrow \quad \mathop{\sum }_{j=1}^{n_{t}}\unicode[STIX]{x1D646}[i,j]=\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D707}}^{2}\quad \forall j\in [1,n_{t}]\end{eqnarray}$$

that is the sum over any row of $\unicode[STIX]{x1D646}$ is equal to the corresponding eigenvalue $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D707}}^{2}$ . This occurs for a stationary process, for which the temporal correlation matrix tends towards a circulant form. For a non-stationary process, the mean removal imposes an additional constraint that alters the connection with the SVD. In such condition, one can still recover the result on optimality by considering the decomposition process as a linear affine transformation (see Miranda, Borgne & Le Bontempi Reference Miranda, Borgne and Le Bontempi2008, for a derivation).

Finally, it is worth highlighting that the energy optimality is linked to the decomposition uniqueness: in the case of repeated singular values, which occurs when different coherent phenomena have a similar energy content (as shown in § 4), there exist infinite possible choices of singular vectors $\unicode[STIX]{x1D731}_{{\mathcal{P}}}$ and $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ .

2.2 Matrix form of the discrete Fourier transform (DFT)

The DFT temporal structures $\unicode[STIX]{x1D733}_{{\mathcal{F}}}$ are defined a priori, regardless of the dataset at hand, as columns of the well-known Fourier matrix

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D733}_{{\mathcal{F}}}=\frac{1}{\sqrt{n_{t}}}\left[\begin{array}{@{}ccccc@{}}1 & 1 & 1 & \cdots \, & 1\\ 1 & w & w^{2} & \cdots \, & w^{n_{t}-1}\\ \vdots & \vdots & \ddots & \vdots \\ 1 & w^{n_{t}-1} & w^{2(n_{t}-1)} & \cdots \, & w^{(n_{t}-1)^{2}}\\ \end{array}\right],\end{eqnarray}$$

where $w=\exp (2\unicode[STIX]{x03C0}\text{j}/n_{t})$ , with $\text{j}=\sqrt{-1}$ . Each temporal structure is thus a complex exponential $\unicode[STIX]{x1D74D}_{r}[k]=\exp (2\unicode[STIX]{x03C0}\text{j}f_{r}(k-1)\unicode[STIX]{x0394}t)=\exp (2\unicode[STIX]{x03C0}\text{j}(r-1)(k-1)/n_{t})$ with real frequency discretization $f_{r}=(r-1)\unicode[STIX]{x0394}f$ , with $\unicode[STIX]{x0394}f=1/(n_{t}\unicode[STIX]{x0394}t)=f_{s}/n_{t}$ the frequency resolution and $f_{s}$ the sampling frequency.

Besides being orthonormal ( $\unicode[STIX]{x1D733}_{{\mathcal{F}}}^{-1}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}^{\dagger }$ ), this Vandermonde matrix is also symmetric ( $\unicode[STIX]{x1D733}_{{\mathcal{F}}}^{\text{T}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}$ ). The resulting spatial structures $\unicode[STIX]{x1D731}_{{\mathcal{F}}}=\unicode[STIX]{x1D63F}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\unicode[STIX]{x1D72E}_{{\mathcal{F}}}^{-1}$ from (2.6) are complex and generally not orthogonal. Finally, since the DFT satisfies (2.7) by construction, the zero-mean shifting has no effect on the decomposition.

2.3 Matrix form of the dynamic mode decomposition (DMD)

The DMD temporal structures $\unicode[STIX]{x1D733}_{{\mathcal{D}}}$ are computed under the assumption that a linear dynamical system can well approximate the dataset. It is thus assumed that a propagator matrix $\unicode[STIX]{x1D64B}\in \mathbb{R}^{n_{s}\times n_{s}}$ maps each column of the data matrix $\unicode[STIX]{x1D63F}$ onto the following one via simple matrix multiplication

(2.14) $$\begin{eqnarray}\boldsymbol{d}_{k}=\unicode[STIX]{x1D64B}\boldsymbol{d}_{k-1}=\unicode[STIX]{x1D64B}^{k-1}\boldsymbol{d}_{1}.\end{eqnarray}$$

As each time step $k$ involves the $k-1$ power of $P$ acting on the initial data $\boldsymbol{d}_{1}$ , the evolution of such linear system depends on the eigendecomposition $\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D726}\unicode[STIX]{x1D64E}^{-1}$ of the propagator, since $\unicode[STIX]{x1D64B}^{k-1}=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D726}^{k-1}\unicode[STIX]{x1D64E}^{-1}$ :

(2.15) $$\begin{eqnarray}\boldsymbol{d}_{k}[\boldsymbol{i}]=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D726}^{k-1}(\unicode[STIX]{x1D64E}^{-1}\boldsymbol{d}_{1}[\boldsymbol{i}])=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D726}^{k-1}a_{0}=\mathop{\sum }_{r=1}^{n_{s}}a_{0}[r]s_{r}[\boldsymbol{i}]\unicode[STIX]{x1D706}_{r}^{k-1}.\end{eqnarray}$$

The columns of the eigenvector matrix $\unicode[STIX]{x1D64E}=[s_{1},s_{2},\ldots ,s_{n_{s}}]$ represent the spatial basis of such evolution. Observe that the possible growth/decay of each mode makes the notion of mode ‘amplitude’ particularly cumbersome in the DMD. Instead, the vector $\boldsymbol{a}_{0}=\unicode[STIX]{x1D64E}^{-1}\boldsymbol{d}_{1}$ is the projection of the initial data $\boldsymbol{d}_{1}$ onto the eigenvectors of the propagator. To arrange (2.15) into the factorization form (2.5), one should first arrange the vector $\boldsymbol{a}_{0}$ into a diagonal matrix $\unicode[STIX]{x1D63C}_{0}=\text{diag}(a_{0})$ and then build the temporal basis from the powers of the eigenvalues. The resulting factorization can be written as

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D63F}=\unicode[STIX]{x1D64E}\unicode[STIX]{x1D63C}_{0}\unicode[STIX]{x1D655}^{\text{T}}\quad \text{with }\unicode[STIX]{x1D655}^{\text{T}}=\left[\begin{array}{@{}ccccc@{}}1 & \unicode[STIX]{x1D706}_{1} & \unicode[STIX]{x1D706}_{1}^{2} & \cdots \, & \unicode[STIX]{x1D706}_{1}^{n_{t}}\\ 1 & \unicode[STIX]{x1D706}_{2} & \unicode[STIX]{x1D706}_{2}^{2} & \cdots \, & \unicode[STIX]{x1D706}_{2}^{n_{t}}\\ \vdots & \vdots & \ddots & \vdots \\ 1 & \unicode[STIX]{x1D706}_{n_{s}} & \unicode[STIX]{x1D706}_{n_{s}}^{2} & \cdots \, & \unicode[STIX]{x1D706}_{n_{s}}^{n_{t}}\\ \end{array}\right].\end{eqnarray}$$

The temporal basis of the DMD is thus obtained by normalizing each column of $\unicode[STIX]{x1D655}$ independently, that is $\unicode[STIX]{x1D74D}_{{\mathcal{D}}r}=Z_{r}/\Vert Z_{r}\Vert$ , and the spatial structures can be computed using the projection in (2.6) with no need for computing the eigenvectors $\unicode[STIX]{x1D64E}$ . In this final step, it is important to observe that the DMD temporal basis is generally not orthogonal and thus the matrix inversion is generally unavoidable.

To compute the DMD, one needs to calculate the eigenvalues of the best (in a $L^{2}$ sense) propagator describing the dataset. This propagator can be easily defined by arranging the data matrix into two shifted portions, $\unicode[STIX]{x1D63F}_{1}=[\boldsymbol{d}_{1},\boldsymbol{d}_{2},\ldots ,\boldsymbol{d}_{n_{t}-1}]$ and $\unicode[STIX]{x1D63F}_{2}=[\boldsymbol{d}_{2},\boldsymbol{d}_{3},\ldots ,\boldsymbol{d}_{n_{t}}]$ , containing $n_{t}-1$ realization:

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D63F}_{2}\approx [\unicode[STIX]{x1D64B}\boldsymbol{d}_{1},\unicode[STIX]{x1D64B}\boldsymbol{d}_{2},\ldots ,\unicode[STIX]{x1D64B}\boldsymbol{d}_{n_{t}-1}]=\unicode[STIX]{x1D64B}\unicode[STIX]{x1D63F}_{1}\Longleftrightarrow \unicode[STIX]{x1D64B}=\unicode[STIX]{x1D63F}_{2}\unicode[STIX]{x1D63F}_{1}^{+},\end{eqnarray}$$

where $\unicode[STIX]{x1D63F}_{1}^{+}$ denotes the pseudo-inverse of the first portion of the dataset. Except for the very special case in which $n_{s}=n_{t}$ and the matrix $\unicode[STIX]{x1D63F}_{1}$ is invertible, the definition of this propagator is not well posed in the classical sense of Hadamard. Moreover, even resorting to the classical Moore–Penrose inverse $\unicode[STIX]{x1D63F}_{1}^{+}=\unicode[STIX]{x1D733}_{{\mathcal{P}}}\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{-1}\unicode[STIX]{x1D731}_{{\mathcal{P}}}^{\text{T}}$ , the size of this matrix is usually prohibitively large ( $\unicode[STIX]{x1D64B}\in \mathbb{R}^{n_{s}\times n_{s}}$ ) and thus the DMD decomposition based on the eigenvalues of $\unicode[STIX]{x1D64B}$ (referred to as exact DMD in Tu et al. Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014) is rarely of practical interest. Many variants of the DMD have been developed in the last few years as reviewed in § 1.2. In this work, we only consider the two original DMD algorithms to compare with the proposed mPOD.

The first DMD algorithm, herein denoted as cDMD, reproduces the action of the propagator $\unicode[STIX]{x1D64B}$ with a much smaller matrix $\boldsymbol{{\mathcal{C}}}\in \mathbb{R}^{(n_{t}-1)\times (n_{t}-1)}$ acting on the right side (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009). This matrix is the companion matrix, defined so that

(2.18) $$\begin{eqnarray}\boldsymbol{{\mathcal{C}}}=\left[\begin{array}{@{}ccccc@{}}0 & 0 & \cdots \, & 0 & c_{1}\\ 1 & 0 & & 0 & c_{2}\\ 0 & 1 & & 0 & \vdots \\ \vdots & & \ddots & & \vdots \\ 0 & 0 & \cdots \, & 1 & c_{n_{t}-1}\\ \end{array}\right]\Longleftrightarrow \unicode[STIX]{x1D63F}_{2}=\unicode[STIX]{x1D64B}\unicode[STIX]{x1D63F}_{1}\approx \unicode[STIX]{x1D63F}_{1}\boldsymbol{{\mathcal{C}}}.\end{eqnarray}$$

The ones in the sub-diagonals shift the columns of $\unicode[STIX]{x1D63F}_{1}$ to obtain those of $\unicode[STIX]{x1D63F}_{2}$ ; the last column contains the set of coefficients $\boldsymbol{c}=[c_{1},c_{2},\ldots ,c_{n_{t}-1}]^{\text{T}}$ that approximates the last temporal realization $\boldsymbol{d}_{n_{t}}$ as a linear combination of the previous ones. As an eigenvalue $\unicode[STIX]{x1D706}$ of $\boldsymbol{{\mathcal{C}}}$ , with eigenvector $\boldsymbol{v}$ such that $\boldsymbol{{\mathcal{C}}}\boldsymbol{v}=\unicode[STIX]{x1D706}\boldsymbol{v}$ , is also an eigenvalue of $\unicode[STIX]{x1D64B}$ with eigenvector $\unicode[STIX]{x1D63F}\boldsymbol{v}$ , this formulation of the DMD relies on the eigendecomposition of $\boldsymbol{{\mathcal{C}}}$ and thus the calculation of coefficients $\boldsymbol{c}$ . This is a least square minimization problem, extremely sensitive to noise and the choice of the final dataset. Using a $\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64D}$ factorization $\unicode[STIX]{x1D63F}_{1}=\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64D}$ , the solution can be written as $\boldsymbol{c}=\unicode[STIX]{x1D64D}^{-1}\unicode[STIX]{x1D64C}^{\text{T}}\boldsymbol{d}_{n_{t}}$ . More robust variants of this method are the optimized DMD by Chen et al. (Reference Chen, Tu and Rowley2012) or the sparsity promoting DMD by Jovanović et al. (Reference Jovanović, Schmid and Nichols2014) in which the $L_{2}$ minimization problem is extended to the entire dataset and not just on the last realization.

The second DMD algorithm, herein denoted as sDMD, avoids the calculation of ${\mathcal{C}}$ and projects the problem onto a space of much lower dimension (Schmid Reference Schmid2010), in which hopefully only a few eigenvectors are relevant. This is the space spanned by the first $r_{c}$ dominant POD modes of the dataset.

This propagator is therefore defined as $\boldsymbol{{\mathcal{S}}}=\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}^{\text{T}}\unicode[STIX]{x1D64B}\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}\in \mathbb{R}^{r_{c}\times r_{c}}$ , where the tilde denotes an approximation $\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}=[\unicode[STIX]{x1D753}_{{\mathcal{P}}1},\ldots ,\unicode[STIX]{x1D753}_{{\mathcal{P}}r_{c}}]$ . Writing the full size propagator using the Moore–Penrose inverse, this projected propagator reads

(2.19) $$\begin{eqnarray}\boldsymbol{{\mathcal{S}}}=\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}^{\text{T}}\unicode[STIX]{x1D64B}\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}=\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}^{\text{T}}\unicode[STIX]{x1D63F}_{2}\tilde{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}\tilde{\unicode[STIX]{x1D72E}}_{{\mathcal{P}}}^{-1}.\end{eqnarray}$$

Observe that (2.19) is a similarity transform, and thus $\boldsymbol{{\mathcal{S}}}$ and $\unicode[STIX]{x1D64B}$ share the same eigenvalues, only if $\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}$ is a square $n_{s}\times n_{s}$ matrix (that is the full POD basis plus its $n_{s}-n_{t}$ orthogonal complements are taken into account so that $\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}^{\text{T}}=I$ ). That is never the case, and one heuristically hopes that at least the first $r_{c}$ eigenvalues of $\unicode[STIX]{x1D64B}$ are available in $\boldsymbol{{\mathcal{S}}}$ , which implies that $\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}\tilde{\unicode[STIX]{x1D731}}_{{\mathcal{P}}}^{\text{T}}\approx I$ .

Although numerically more robust, the ill posedness of this POD-based projection is easily illustrated with the test cases in §§ 4 and 5. Specifically, major problems occur when the POD basis captures the energy content of the dataset within a few modes, each of which possibly having several essential frequencies. In this condition, the projected propagator in (2.19) acts on a space which is too small to allow for a sufficient number of frequencies, and the DMD decomposition diverges exponentially fast.

Finally, concerning the mean-shifted form in (2.7), one should notice that the DMD can provide the constant vector $\unicode[STIX]{x1D74D}_{\unicode[STIX]{x1D707}}=\text{}\underline{1}/\sqrt{n_{t}}$ as basis element only if $\unicode[STIX]{x1D706}=1$ is an eigenvalue of the propagator ( $\boldsymbol{{\mathcal{C}}}$ for the cDMD, $\boldsymbol{{\mathcal{S}}}$ for the sDMD). While there is no simple relation to know when this is the case for $\boldsymbol{{\mathcal{S}}}$ , one can see from the characteristic polynomial of $\boldsymbol{{\mathcal{C}}}$ that this condition sets $\langle \text{}\underline{1},\boldsymbol{c}\rangle =0$ , that is the sum of all the coefficients should be null. This condition is not generally imposed in the calculation of the DMD, so the corresponding basis might not necessarily have a mode to represent the temporal average of data. On the other hand, removing the mean from the dataset results in $\boldsymbol{c}=-\text{}\underline{1}$ regardless of the dataset, since the last realization must cancel the summation. In this case, as shown by Chen et al. (Reference Chen, Tu and Rowley2012), the solution of the characteristic polynomials are the $n_{t}-1$ roots of unity, and the DMD reduces to the DFT.

3.1 Mode spectra and correlation spectra

The multi-resolution analysis (MRA) is introduced in § 3.2 in the frequency domain, that is along the 2-D Fourier transform of the correlation matrix. It is thus interesting to first analyse the standard POD in the frequency domain. At the scope, we denote as $\widehat{\unicode[STIX]{x1D63C}}$ the discrete Fourier transform (DFT) of a matrix $\unicode[STIX]{x1D63C}$ , whether this is done along its rows, columns or both. Using (2.13), it is easy to see that the DFT of a column vector is obtained via left multiplication by $\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}$ , while the DFT of a row vector is a right multiplication by $\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}$ .

The three Fourier pairs to be linked in this section are related to the time evolution of the data ( $\unicode[STIX]{x1D63F}$ , row-wise), the temporal structures of the POD modes ( $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ , column-wise) and the temporal correlation matrix ( $\unicode[STIX]{x1D646}$ , over both columns and rows). These are:

(3.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D63F}}=\unicode[STIX]{x1D63F}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\Longleftrightarrow \unicode[STIX]{x1D63F}=\widehat{\unicode[STIX]{x1D63F}}\unicode[STIX]{x1D733}_{{\mathcal{F}}} & \displaystyle\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}}\Longleftrightarrow \unicode[STIX]{x1D733}_{{\mathcal{P}}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}} & \displaystyle\end{eqnarray}$$

(3.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D646}}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\unicode[STIX]{x1D646}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\Longleftrightarrow \unicode[STIX]{x1D646}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\widehat{\unicode[STIX]{x1D646}}\unicode[STIX]{x1D733}_{{\mathcal{F}}}. & \displaystyle\end{eqnarray}$$

A first link between the DFT of the dataset (3.1a ) and that of the POD modes (3.1b ) is identified by the cross-spectral density matrix, defined and linked to the temporal correlation matrix as follows:

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{{\mathcal{F}}}=\widehat{\unicode[STIX]{x1D63F}}^{\dagger }\widehat{\unicode[STIX]{x1D63F}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}[\unicode[STIX]{x1D63F}^{\dagger }\unicode[STIX]{x1D63F}]\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\unicode[STIX]{x1D646}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\Longleftrightarrow \unicode[STIX]{x1D646}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\unicode[STIX]{x1D646}_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{F}}}.\end{eqnarray}$$

This matrix is obtained from the temporal correlation matrix $\unicode[STIX]{x1D646}$ via a similarity transform (since $\unicode[STIX]{x1D733}_{{\mathcal{F}}}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=\unicode[STIX]{x1D644}$ ), and therefore share the same eigenvalues.

To link its eigenvectors to the POD temporal structures, it is important to observe that all the transformed quantities in (3.1a )–(3.1c ) (that is $\unicode[STIX]{x1D63F}$ , $\unicode[STIX]{x1D733}_{{\mathcal{P}}}$ and $\unicode[STIX]{x1D646}$ ) are real, and thus their Fourier transforms are Hermitian symmetric. For a continuous signal $x(t)$ , having Fourier transform ${\mathcal{F}}\{x(t)\}=X(\unicode[STIX]{x1D714})$ , the Hermitian symmetry sets $X(\unicode[STIX]{x1D714})=\overline{X(-\unicode[STIX]{x1D714})}$ . For a discrete signal, this property can be introduced via a matrix permutation $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}$ , which flips the spectra along the zero frequency (first column of $\unicode[STIX]{x1D733}_{{\mathcal{F}}}$ in (2.13)).

This permutation matrix can be obtained by applying the DFT operator twice:

(3.3)

Therefore, for any discrete real vector $\boldsymbol{a}$ with discrete Fourier transform $\hat{\boldsymbol{a}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\boldsymbol{a}$ , the Hermitian symmetry can be written as $\hat{\boldsymbol{a}}=\overline{\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\hat{\boldsymbol{a}}}$ . The fact that the permutation matrix can be obtained by transforming the discrete vector twice reflects the symmetry (or duality) property of the Fourier transform.

Using the Hermitian symmetry on the first two Fourier pairs in (3.1a )–(3.1b ) gives:

(3.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}=\overline{\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}}\Longleftrightarrow \overline{\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}}=\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}}\Longleftrightarrow \unicode[STIX]{x1D733}_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}}=\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}} & \displaystyle\end{eqnarray}$$

(3.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D63F}}=\overline{\widehat{\unicode[STIX]{x1D63F}}\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}}\Longleftrightarrow \overline{\widehat{\unicode[STIX]{x1D63F}}}=\widehat{\unicode[STIX]{x1D63F}}\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\Longleftrightarrow \unicode[STIX]{x1D63F}\unicode[STIX]{x1D733}_{{\mathcal{F}}}=\widehat{\unicode[STIX]{x1D63F}}\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}. & \displaystyle\end{eqnarray}$$

Introducing the eigenvalue decomposition of $\unicode[STIX]{x1D646}$ in (3.2) and using (3.4a ) yields

(3.5) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{{\mathcal{F}}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}[\unicode[STIX]{x1D733}_{{\mathcal{P}}}\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{2}\unicode[STIX]{x1D733}_{{\mathcal{P}}}^{\dagger }]\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=(\unicode[STIX]{x1D733}_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}})\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{2}(\unicode[STIX]{x1D733}_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{P}}})^{\dagger }=(\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}})\unicode[STIX]{x1D72E}_{{\mathcal{P}}}^{2}(\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}})^{\dagger }.\end{eqnarray}$$

This shows that the eigenvectors of $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ are the permuted spectra of the POD modes. While the similarity argument is essential in any spectral formulation of the POD (Glauser et al. Reference Glauser, Leib and George1987; Citriniti & George Reference Citriniti and George2000; Gordeyev & Thomas Reference Gordeyev and Thomas2000; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Sieber et al. Reference Sieber, Paschereit and Oberleithner2017; Towne et al. Reference Towne, Schmidt and Colonius2018), (3.5) shows that care must be taken with the conjugation (or the flipping of) the eigenvectors before the inverse DFT, if one seeks to compute POD modes from the cross-spectral density matrix in (3.2).

The DFT of the correlation matrix $\unicode[STIX]{x1D646}$ in (3.1c ), on the other hand, is not obtained via similarity transform and is not self-adjoint ( $\widehat{\unicode[STIX]{x1D646}}\neq \widehat{\unicode[STIX]{x1D646}}^{\dagger }$ ). Hence, since it does not share the same eigenvalues of $\unicode[STIX]{x1D646}$ and $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ , its diagonalization is of no interest for the purpose of this work. However, because of its relevance on the filtering process, it is important to highlight – by introducing (3.2) into (3.1c ) – its link to $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ :

(3.6) $$\begin{eqnarray}\widehat{\unicode[STIX]{x1D646}}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}[\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}K_{{\mathcal{F}}}\unicode[STIX]{x1D733}_{{\mathcal{F}}}]\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}K_{{\mathcal{F}}}=\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D646}_{{\mathcal{F}}}.\end{eqnarray}$$

3.2 Decomposing data and decomposing correlations

It is now of interest to understand how temporal filtering of the dataset influences the frequency content of the POD modes. In particular, we are interested not in just one filter, but an array of filters, constructed to isolate portions of the frequency spectra.

Figure 1.

Multi-resolution analysis (MRA) of a signal $\boldsymbol{d}_{p}[k]$ along its spectrum $\widehat{\boldsymbol{d}}_{p}[n]$ . Using a filter bank, the MRA splits the spectrum into a set of $M$ scales. This set consists of a large-scale (transfer function $\unicode[STIX]{x1D643}_{{\mathcal{L}}_{1}}$ ), a fine-scale (transfer function $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{M-1}}$ ) and $M-2$ intermediate scales (transfer functions from $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{1}}$ to $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{M-2}}$ ). Only the positive portion of the spectrum is shown, the negative one being mirrored along the ordinates.

To illustrate the idea, let us consider the Fourier transform $\widehat{\boldsymbol{d}}_{p}=\boldsymbol{d}_{p}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}$ of the temporal evolution (row-wise) of the data for a given point spatial location $\boldsymbol{i}_{p}$ , that is $\boldsymbol{d}_{p}[k]=\unicode[STIX]{x1D63F}[\boldsymbol{i}_{p},k]$ . We seek to split the spectra $\widehat{\boldsymbol{d}}_{p}$ into $M$ scales, each retaining a portion with negligible overlapping as pictorially illustrated in figure 1. The frequency bandwidths of these scales $\unicode[STIX]{x0394}f_{m}=f_{m+1}^{c}-f_{m}^{c}$ , with $m\in [1,\ldots ,M-1]$ , are defined by a frequency splitting vector $\boldsymbol{F}_{V}=[f_{1}^{c},f_{2}^{c},\ldots ,f_{M-1}^{c}]$ . This vector is bounded by the minimum ( $f_{0}^{c}=0$ ) and the maximum ( $f_{M}^{c}=f_{s}/2$ ) possible frequencies.

This splitting requires the definition of a filter bank (Vaidyanathanm Reference Vaidyanathanm1992; Strang & Nguyen Reference Strang and Nguyen1996; Shukla & Tiwari Reference Shukla and Tiwari2013), consisting of a low-pass filter with cutoff $f_{1}^{c}$ , a high-pass filter with cutoff $f_{M-1}^{c}$ and $M-2$ band-pass filters between these. Following the lossless formulation from multi-resolution analysis (MRA), these filters are constructed from the set of $m$ low-pass filters with transfer function $\boldsymbol{{\mathcal{H}}}_{{\mathcal{L}}_{m}}$ . Of these low-pass filters, only the first one ( $\boldsymbol{{\mathcal{H}}}_{{\mathcal{L}}_{1}}$ ) is retained, while the others (shown in figure 1 with dashed line) are only used to build the band-pass filters as complementary differences $\unicode[STIX]{x1D643}_{{\mathcal{H}}m}=\unicode[STIX]{x1D643}_{{\mathcal{L}}_{m+1}}-\unicode[STIX]{x1D643}_{{\mathcal{L}}_{m}}$ , the finest one being $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{M}}=1-\boldsymbol{{\mathcal{H}}}_{{\mathcal{L}}_{M}}$ . The resulting set of transfer functions (with modulus sketched in figure 1 with a continuous line) is by construction such that $\unicode[STIX]{x1D643}_{{\mathcal{L}}_{1}}+\unicode[STIX]{x1D643}_{{\mathcal{H}}_{1}}+\cdots +\unicode[STIX]{x1D643}_{{\mathcal{H}}_{M-1}}=\text{}\underline{1}$ and thus the entire spectrum of the dataset is retained. The design of the low-pass filters in the filter bank used in this work is presented in appendix A; in what follows, we focus on the algebra of the filtering process and the theoretical foundation of the proposed mPOD.

Considering the filtering in the frequency domain, this operation is performed via simple multiplication (entry by entry), which in matrix notation corresponds to the Hadamard product $\odot$ . Moreover, since each of these filters is applied in the entire spatial domain (along each row of $\unicode[STIX]{x1D63F}$ ), one should copy the transfer functions $\unicode[STIX]{x1D643}_{{\mathcal{L}}_{m}}$ or $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{m}}$ (row vectors $\in \mathbb{C}^{1\times n_{t}}$ ) to obtain $n_{s}\times n_{t}$ matrices having all the rows equal to $\unicode[STIX]{x1D643}_{{\mathcal{L}}_{m}}$ or $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{m}}$ . Let $\unicode[STIX]{x1D643}_{m}^{\prime }$ be the general transfer function (whether low pass or band pass), we denote with an apex the bi-dimensional extension of the 1-D transfer function. The contribution of a scale $m$ from figure 1 in each of the spatial points (the filtered portion of the dataset) is

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D63F}_{m}=\overbrace{[\underbrace{(\unicode[STIX]{x1D63F}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}})}_{\widehat{\unicode[STIX]{x1D63F}}}\odot \,\unicode[STIX]{x1D643}_{m}^{\prime }]}^{\widehat{\unicode[STIX]{x1D63F}}_{m}}\unicode[STIX]{x1D733}_{{\mathcal{F}}}.\end{eqnarray}$$

The first multiplication is the DFT along the time domain ( $\widehat{\unicode[STIX]{x1D63F}}$ ), the result of the Hadamard multiplication in the square bracket is the filtered spectrum of the data ( $\widehat{\unicode[STIX]{x1D63F}}_{m}$ , the frequency contribution of the given scale) and the last right multiplication is the inverse DFT from (3.1a ). Using (3.7), the correlation matrix $\unicode[STIX]{x1D646}_{m}=\unicode[STIX]{x1D63F}_{m}^{\dagger }\unicode[STIX]{x1D63F}_{m}$ from each scale contribution reads

(3.8) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{m}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}[(\widehat{\unicode[STIX]{x1D63F}}\odot \unicode[STIX]{x1D643}_{m}^{\prime })^{\dagger }(\widehat{\unicode[STIX]{x1D63F}}\odot \unicode[STIX]{x1D643}_{m}^{\prime })]\unicode[STIX]{x1D733}_{{\mathcal{F}}}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}[(\widehat{\unicode[STIX]{x1D63F}}^{\dagger }\widehat{\unicode[STIX]{x1D63F}})\odot ((\unicode[STIX]{x1D643}_{m}^{\prime })^{\dagger }\odot \unicode[STIX]{x1D643}_{m}^{\prime })]\unicode[STIX]{x1D733}_{{\mathcal{F}}},\end{eqnarray}$$

where the associative and commutative properties of the Hadamard product are used to move the transfer functions. Introducing the cross-spectral density matrix $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ in (3.2) and writing the 2-D transfer function as $\text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}=(\unicode[STIX]{x1D643}_{m}^{\prime })^{\dagger }\odot \unicode[STIX]{x1D643}_{m}^{\prime }$ yields

(3.9) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{m}=\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}\underbrace{[\unicode[STIX]{x1D646}_{{\mathcal{F}}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}]}_{\unicode[STIX]{x1D646}_{{\mathcal{F}}m}}\unicode[STIX]{x1D733}_{{\mathcal{F}}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\underbrace{[\widehat{\unicode[STIX]{x1D646}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}]}_{\widehat{\unicode[STIX]{x1D646}}_{m}}\unicode[STIX]{x1D733}_{{\mathcal{F}}},\end{eqnarray}$$

having introduced (3.6) in the last step.

This important result lays the foundation of the mPOD. First, the expression on the left of (3.9) is a similarity transform: the filtered cross-spectral density $\unicode[STIX]{x1D646}_{{\mathcal{F}}m}=\unicode[STIX]{x1D646}_{{\mathcal{F}}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}$ shares the same eigenvalues of the correlation of filtered data $\unicode[STIX]{x1D646}_{m}=\unicode[STIX]{x1D63F}_{m}^{\dagger }\unicode[STIX]{x1D63F}_{m}$ . Second, introducing the eigenvalue decomposition of $\unicode[STIX]{x1D646}_{m}$ as in (3.5), shows that – up to a permutation – the eigenvectors of $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ are the DFT of the POD modes of $\unicode[STIX]{x1D63F}_{m}$ :

(3.10) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{{\mathcal{F}}m}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}K_{m}\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}[\unicode[STIX]{x1D733}_{{\mathcal{P}}m}\unicode[STIX]{x1D72E}_{{\mathcal{P}}m}^{2}\unicode[STIX]{x1D733}_{{\mathcal{P}}m}^{\dagger }]\overline{\unicode[STIX]{x1D733}}_{{\mathcal{F}}}=(\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}m})\unicode[STIX]{x1D72E}_{{\mathcal{P}}m}^{2}(\unicode[STIX]{x1D64B}_{\unicode[STIX]{x03C0}}\widehat{\unicode[STIX]{x1D733}}_{{\mathcal{P}}m})^{\dagger },\end{eqnarray}$$

having introduced (3.4a ) in the last step. The impact of the filter on the POD modes is revealed by the diagonal entries of $\unicode[STIX]{x1D646}_{{\mathcal{F}}m}$ and $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ , which are summations of non-negative real numbers. From (3.10), using (3.4a ), these are:

(3.11) $$\begin{eqnarray}\unicode[STIX]{x1D646}_{{\mathcal{F}}m}[i,j]=\unicode[STIX]{x1D646}_{{\mathcal{F}}}[i,j]\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}[i,j]=\mathop{\sum }_{r=1}^{n_{t}}\unicode[STIX]{x1D70E}_{{\mathcal{P}}mr}^{2}\overline{\widehat{\unicode[STIX]{x1D74D}}}_{{\mathcal{P}}mr}[i]\overline{\widehat{\unicode[STIX]{x1D74D}}}_{{\mathcal{P}}mr}^{\dagger }[j].\end{eqnarray}$$

If the filter at a given scale $m$ removes one of the entries along the diagonal of $\unicode[STIX]{x1D646}_{{\mathcal{F}}}$ ( $\text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}[i,i]=0$ ), the associated frequency $\widehat{\unicode[STIX]{x1D74D}}_{{\mathcal{P}}mr}[i]$ cannot be present in any of the $r$ POD modes of the filtered data: the POD operating at this scale is constrained within the frequencies allowed by the transfer function.

Moreover, while the expression on the right side of (3.9) is not a similarity transform, it shows that the correlation matrix associated with the filtered data $\unicode[STIX]{x1D646}_{m}$ can be computed by a simple 2-D filtering of the original correlation matrix, saving considerable computational time with respect to $n_{s}$ filters along the rows of $\unicode[STIX]{x1D63F}$ . Equation (3.9) shows that the MRA of the dataset can be carried out via 2-D MRA of the temporal correlation matrix.

Figure 2.

The magnitude of the four possible 2-D transfer functions $\text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}$ that can be obtained by combining the low-pass ( $\unicode[STIX]{x1D643}_{{\mathcal{L}}_{m}}$ ) and the high-pass ( $\unicode[STIX]{x1D643}_{{\mathcal{H}}_{m}}$ ) filters associated with a given $m$ scale of $\unicode[STIX]{x1D646}$ , following the multi-scale architecture in figure 1. The areas where the transfer function is ${\sim}1$ and ${\sim}0$ are indicated with blue dashed lines and full grey respectively.

The four possible structures for the 2-D transfer functions $\text{}\underline{\boldsymbol{{\mathcal{H}}}}_{m}$ that arise from a combination of the 1-D kernels in the time domain are illustrated graphically in figure 2. For a given scale, composed of a low frequency and a high frequency $\unicode[STIX]{x1D63F}_{m}=\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}+\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}$ , these correspond to the contributions of four temporal correlation matrices:

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D646}_{m}=\unicode[STIX]{x1D63F}_{m}^{\dagger }\unicode[STIX]{x1D63F}_{m} & = & \displaystyle \unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}+\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}+\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}+\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D733}_{{\mathcal{F}}}[\widehat{\unicode[STIX]{x1D646}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{{\mathcal{L}}m}+\widehat{\unicode[STIX]{x1D646}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{{\mathcal{H}}m}+\widehat{\unicode[STIX]{x1D646}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{{\mathcal{L}}{\mathcal{H}}m}+\widehat{\unicode[STIX]{x1D646}}\odot \text{}\underline{\boldsymbol{{\mathcal{H}}}}_{{\mathcal{H}}{\mathcal{L}}m}]\unicode[STIX]{x1D733}_{{\mathcal{F}}}.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The ‘pure’ terms $\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}$ and $\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}$ can be obtained by first filtering the data (3.7) and then computing the related correlations; the mixed terms can only be revealed by a MRA of the full correlation matrix $\unicode[STIX]{x1D646}_{m}$ .

It is worth noticing that the filter bank architecture described in this section is a generalized version of the MRA via the discrete (dyadic) wavelet transform (DWT) proposed in previous works (Mendez et al. Reference Mendez, Scelzo and Buchlin2018b ; Mendez, Balabane & Buchlin Reference Mendez, Balabane and Buchlin2018a ). In the 2-D wavelet terminology, the four terms in (3.12), with spectral band-pass region illustrated in figure 2, correspond to the approximation and the diagonal, horizontal and vertical details of each scale.

Particular emphasis should be given to the ‘pure’ terms $\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}$ and $\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}$ . These contributions have no frequency overlapping and their eigenspaces are orthogonal complements by construction: their set of eigenvectors can be used to assemble a complete and orthogonal basis as illustrated in § 3.4. The mixed terms $\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}$ and $\unicode[STIX]{x1D63F}_{{\mathcal{H}}_{m}}^{\dagger }\unicode[STIX]{x1D63F}_{{\mathcal{L}}_{m}}$ , on the other hand, generates eigenspaces that can potentially overlap with those of the ‘pure’ contributions and must be disregarded in the MRA of the correlation matrix.

3.3 The mPOD as bridge between POD and DFT

Neglecting the mixed terms prevents the full reconstruction of the correlation matrix, but not the full reconstruction of the dataset since the final temporal basis is complete. The consequences of this operation can be analysed while studying how the choice of the frequency splitting vector $\boldsymbol{F}_{V}$ allows the mPOD to move from the energy optimality of the POD to the spectral purity of the DFT.

Figure 3 shows the band-pass region from the 2-D spectra of the correlation matrix $\widehat{\unicode[STIX]{x1D646}}$ for four choices of $\boldsymbol{F}_{V}$ . The case (a) corresponds to a single scale mPOD (empty  $\boldsymbol{F}_{V}$ ), that is the standard POD. All the eigenvectors $\unicode[STIX]{x1D74D}_{{\mathcal{P}}_{r}}$ are allowed to span the entire set of discrete frequencies available. This configuration corresponds to perfect reconstruction of the correlation matrix and maximum time localization capabilities: if an impulse $\unicode[STIX]{x1D6FF}(t_{k})$ is produced at a certain time, the POD is free to choose an impulse (which has a frequency spectra spanning the entire frequency axis) as one of its temporal structures. On the other hand, all the POD modes share the full set of frequencies, potentially leading to spectral mixing, and phenomena spanning the entire frequency and energy spectra (such as random noise) are spread over all the modes.

The case (b) corresponds to a two-scale mPOD, that is with $\boldsymbol{F}_{V}=f_{1}^{c}$ . Each scale has its own correlation matrix ( $\unicode[STIX]{x1D646}_{{\mathcal{L}}}$ and $\unicode[STIX]{x1D646}_{{\mathcal{H}}_{1}}$ ) with no common frequencies. Phenomena occurring at these two scales are spectrally separated and necessarily assigned to different modes. Moreover, the amount of random noise in each mode is distributed proportionally to the spectral bandwidth of the corresponding scale and becomes less important. Part of the information on the correlation matrix is lost, but the eigenspaces of the two scales, together, span the entire $\mathbb{R}^{n_{t}}$ space and thus preserve the decomposition convergence. On the other hand, the temporal localization capabilities are reduced: as no eigenvector is allowed to span the full frequency range, no temporal structure can capture impulses nor other phenomena requiring the entire frequency range (e.g. shocks).

As the number of scales is increased, their spectral bandwidth reduces and the decomposition becomes spectrally more localized, at the cost of an increased loss of information in the correlation and a reduction of time localization capabilities. The case (c) has three scales, with $\boldsymbol{F}_{V}=[f_{1}^{c},f_{2}^{c}]$ , while the case (d) corresponds to the limit $\boldsymbol{F}_{V}\rightarrow \{rf_{s}/(\unicode[STIX]{x0394}tn_{t})\}_{r=1}^{n_{t}-1}$ , with $r\in [1,\ldots ,n_{t}-1]$ . In this limiting case, every scale is constrained to a single frequency and the mPOD reduces to the DFT, inheriting its finest spectral resolution and its limits (windowing problems, Gibbs phenomena, no time localization and poor convergence, Harris Reference Harris1978, Daubechies Reference Daubechies1992, Folland Reference Folland2009).

Figure 3.

Partition of the correlation matrix $\widehat{\unicode[STIX]{x1D646}}$ for different choices of the splitting vector $\boldsymbol{F}_{V}$ in the MRA step of the mPOD. The same colour notation of figure 2 is used for the band-pass/band-stop regions. Case (a) consists of a single scale, thus the mPOD becomes a POD. Cases (b,c) consist of two and three scales, respectively. Case (d) shows the limit of $n_{t}/2$ scales, in which the mPOD becomes a DFT.

3.4 The mPOD algorithm

The steps of the proposed mPOD algorithm are listed in the Algorithm 1. As for any data-driven decomposition, the first step consists in assembling the data matrix as in (2.1). The second step consists in computing the temporal correlation matrix $\unicode[STIX]{x1D646}$ and its Fourier transform $\widehat{\unicode[STIX]{x1D646}}$ from (3.1c ). From the analysis of the frequency content in the correlation spectra $\widehat{\unicode[STIX]{x1D646}}$ , the third step consists in defining the frequency splitting vector $\boldsymbol{F}_{V}$ and in constructing the set of associated transfer functions.

This choice is the fundamental step which will set the spectral constraints to the mPOD modes. The frequency splitting could be user-defined if prior knowledge on the investigated data is available, or automatically computed by centring the band-pass windows of each scale on the dominant peaks in the diagonal of $\widehat{\unicode[STIX]{x1D646}}$ and splitting the band-pass regions accordingly. In this step, only the ‘pure’ (approximation and diagonal detail) terms of the 2-D transfer functions are constructed, that is $\boldsymbol{{\mathcal{H}}}_{{\mathcal{L}}_{1}}=(\unicode[STIX]{x1D643}_{{\mathcal{L}}_{1}}^{\prime })^{\dagger }\odot \unicode[STIX]{x1D643}_{{\mathcal{L}}_{1}}^{\prime }$ and $\boldsymbol{{\mathcal{H}}}_{{\mathcal{H}}_{m}}=(\unicode[STIX]{x1D643}_{{\mathcal{H}}_{m}}^{\prime })^{\dagger }\odot \unicode[STIX]{x1D643}_{{\mathcal{H}}_{m}}^{\prime }$ . Disregarding the mixed terms (horizontal and vertical details) in each scale, the correlation matrix is approximated as:

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D646}\approx \unicode[STIX]{x1D733}_{{\mathcal{F}}}[\widehat{\unicode[STIX]{x1D646}}\odot \boldsymbol{{\mathcal{H}}}_{{\mathcal{L}}_{1}}]\unicode[STIX]{x1D733}_{{\mathcal{F}}}+\mathop{\sum }_{m=1}^{M}\unicode[STIX]{x1D733}_{{\mathcal{F}}}[\widehat{\unicode[STIX]{x1D646}}\odot \boldsymbol{{\mathcal{H}}}_{{\mathcal{H}}_{m}}]\unicode[STIX]{x1D733}_{{\mathcal{F}}}\approx \unicode[STIX]{x1D646}_{{\mathcal{L}}_{1}}+\mathop{\sum }_{m=1}^{M-1}\unicode[STIX]{x1D646}_{{\mathcal{H}}_{m}}.\end{eqnarray}$$

Each of these contributions is a symmetric, real and positive definite matrix, equipped with its orthonormal eigenspace, computed in the fourth step

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D646}\approx \unicode[STIX]{x1D733}_{{\mathcal{L}}_{1}}\unicode[STIX]{x1D72E}_{{\mathcal{L}}_{1}}^{2}\unicode[STIX]{x1D733}_{{\mathcal{L}}_{1}}^{\text{T}}+\mathop{\sum }_{m=1}^{M-1}\unicode[STIX]{x1D733}_{{\mathcal{H}}m}\unicode[STIX]{x1D72E}_{{\mathcal{H}}_{m}}^{2}\unicode[STIX]{x1D733}_{{\mathcal{H}}_{m}}^{\text{T}}.\end{eqnarray}$$

If the frequency overlap is identically null, the eigenspaces of these pure terms are orthogonal complements and therefore there are at most $n_{t}$ non-zero eigenvalues among all the scales. This is a major difference between the mPOD and other multi-scale methods such as continuous wavelet transform (CWT), in which the temporal basis is constructed by shifting and dilating a ‘mother’ function, or the multi-resolution DMD (mrDMD) proposed by Kutz et al. (Reference Kutz, Fu and Brunton2016b ), in which the temporal basis is constructed by performing DMD on different portions of the datasets. These decompositions potentially produce high redundancy and poor convergence since the basis is larger than $n_{t}$ . Moreover, each of the basis elements in the mPOD exists over the entire time domain (although they could be null in an arbitrarily large portion of it), while CWT or mrDMD produce different bases for different portions of the time domain, leading to decompositions more complicated than those analysed in this work.

Because of the limited size of the filter impulse response (see appendix A) of each filter, the number of non-zero eigenvalues among the scales is in practice usually slightly larger than $n_{t}$ . Therefore, only the first $n_{t}$ dominant eigenvectors are selected. This is done in the fifth step, by first sorting the full set of eigenvalues $[\text{diag}(\unicode[STIX]{x1D72E}_{{\mathcal{L}}_{1}}),\text{diag}(\unicode[STIX]{x1D72E}_{{\mathcal{H}}_{1}}),\ldots ,\text{diag}(\unicode[STIX]{x1D72E}_{{\mathcal{H}}_{M-1}})]$ in decreasing order and storing the required permutation matrix $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D72E}}$ . This matrix is then used to re-arrange the eigenvector matrices $\unicode[STIX]{x1D733}_{{\mathcal{L}}},\ldots ,\unicode[STIX]{x1D733}_{{\mathcal{H}}m}$ into one single matrix (temporal basis) regardless of their scale of origin

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D733}_{{\mathcal{M}}}^{0}=[\unicode[STIX]{x1D733}_{{\mathcal{L}}_{1}},\unicode[STIX]{x1D733}_{{\mathcal{H}}_{1}},\unicode[STIX]{x1D733}_{{\mathcal{H}}_{2}},\ldots ,\unicode[STIX]{x1D733}_{{\mathcal{H}}_{M-1}}]\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D72E}}.\end{eqnarray}$$

As discussed in § 3.2, the temporal matrix (3.15) is orthonormal if a perfect spectral separation is achieved. This is a major difference between the mPOD and the SPOD by Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016), in which the filtering procedure can eventually result in a non-orthogonal temporal basis, or the mrDMD by Kutz et al. (Reference Kutz, Fu and Brunton2016b ), in which the loss of orthogonality can be produced by growing/decaying modes.

Due to the finite size of the transition band of the filters, this is generally not the case in practice, and a minor loss of orthogonality usually occurs. The sixth step treats this imperfection, polishing this temporal structure via a reduced $\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64D}$ factorization $\unicode[STIX]{x1D733}_{{\mathcal{M}}}^{0}=\unicode[STIX]{x1D733}_{{\mathcal{M}}}\unicode[STIX]{x1D64D}$ , so that $\unicode[STIX]{x1D733}_{{\mathcal{M}}}=\unicode[STIX]{x1D733}_{{\mathcal{M}}}^{0}\unicode[STIX]{x1D64D}^{-1}\in {\mathcal{R}}^{n_{t}\times n_{t}}$ . This compensates for the losses of orthogonality and the upper triangular matrix $\unicode[STIX]{x1D64D}$ offers an indication of the quality of the filtering process, being close to the identity when the spectral overlapping is negligible. Finally, in the last step, the decomposition is completed via the projection and the normalization in (2.6), as for any orthogonal decomposition.

Table 1.

CPU time (in seconds) for the test cases analysed in this work and the decompositions tested. The analysis is carried out using Matlab R2018b on a Intel(TM) i7-3770 @3.40 GHz processor with a 32 GB RAM.

To conclude the presentation of the mPOD algorithm, table 1 collects the CPU time for all the decompositions used in this work and the test cases in §§ 4–6. The POD is naturally the fastest requiring only an eigenvalue decomposition and a projection, with no need for amplitude normalization and sorting. The DFT is slightly more expensive, mostly due to the normalization and the sorting steps, while the DMD requires the additional cost of a minimization problem or an SVD depending on the chosen algorithm. The computational cost of the mPOD is shown for $M=3,10,50$ and drastically increases with the number of scales considered.

While the cost of the MRA is usually negligible when using an FFT-based formulation in the frequency domain, the major additional cost is in the diagonalizations for each scale. For $M=3$ , the CPU time for the proposed mPOD algorithm is well below the CPU time of the DMD and become comparable to it for $M\sim 10$ . For a much larger number of scales (see table for $M=50$ ), the mPOD requires a more efficient implementation of the diagonalization of the various scales. An improved algorithm which performs this step in the frequency domain is under development and will be presented in future work.

4.1 Dataset description

As a first example, we propose a simple synthetic test case which is already in the decomposed form of (2.5) and yet capable of severely challenging standard decompositions. The spatial domain consists of a square Cartesian grid of $n_{s}=256\times 256$ points over a domain $\boldsymbol{x}_{i}\in [-20,20]\times [-20,20]$ , while the time domain spans $t_{k}\in [0,5.11]$ with a total of $n_{t}=256$ points with sampling frequency $f_{s}=100$ . The dataset is composed of the sum of three modes, having spatial and temporal structures as illustrated in figure 4. These modes consist of identical Gaussians with standard deviations of $\unicode[STIX]{x1D70E}=5$ , located in $\boldsymbol{x}_{1}=[10,-10]$ , $\boldsymbol{x}_{2}=[-10,10]$ and $\boldsymbol{x}_{3}=[0,0]$ , pulsing as

(4.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{1}(t_{k})=A_{1}\sin (2\unicode[STIX]{x03C0}f_{1}t_{k})\exp [0.5(t_{k}-3)^{20}] & \displaystyle\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{2}(t_{k})=A_{2}\sin (2\unicode[STIX]{x03C0}f_{2}t_{k}-\unicode[STIX]{x03C0}/3) & \displaystyle\end{eqnarray}$$

(4.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{3}(t_{k})=A_{3}\sin (2\unicode[STIX]{x03C0}f_{3}t_{k})(t_{k}-2.55)^{2}. & \displaystyle\end{eqnarray}$$

The first mode has a smoothed square box modulation, the second has a harmonic with period longer than the observation time and the third has a parabolic modulation. The amplitudes $A_{k}$ are chosen so as to equal each contribution, that is $A_{r}=1/\Vert S_{r}T_{r}^{T}\Vert _{2}$ , having reshaped both spatial and temporal structures as column vectors. The three frequencies are $f_{1}=15$ , $f_{2}=0.1$ and $f_{3}=7$ . The resulting dataset $\unicode[STIX]{x1D63F}[\boldsymbol{i},k]$ is a matrix of size $65\,536\times 256$ of rank $rk(\unicode[STIX]{x1D63F})=3$ .

Figure 4.

Spatial (a–c) and corresponding temporal (d–f) structures of the modes considered in the simplified test case. The amplitude of each contribution is computed so as to give the energy content. An animation of this test case is available as supplemental Movie 1 available at https://doi.org/10.1017/jfm.2019.212.

4.2 Modal analysis

In this test case, we compare the performances of the DFT, POD, cDMD and mPOD. The sDMD algorithm, based on the reduced propagator $\boldsymbol{{\mathcal{S}}}$ in (2.19), is severely challenged by the rank deficiency of this dataset. Since $rk(\unicode[STIX]{x1D63F})=3$ , one has $\boldsymbol{{\mathcal{S}}}\in \mathbb{R}^{3\times 3}$ , that is a set of eigenvalues incapable of reproducing the data: each of the dominant frequencies $f_{1},f_{2},f_{3}$ requires two complex conjugate modes and thus at least a $6\times 6$ propagator. Low-dimensional propagators $\boldsymbol{{\mathcal{S}}}$ that are inadequately too small occur for datasets which are close to being rank deficient, in which few coherent patterns occur at different frequencies.

For all the decomposition, the $L^{2}$ convergence is calculated as

(4.2) $$\begin{eqnarray}\frac{Err[r_{c}]}{\Vert \unicode[STIX]{x1D63F}\Vert _{2}}=\frac{\displaystyle \mathop{\left\Vert \unicode[STIX]{x1D63F}-\mathop{\sum }_{r=1}^{r_{c}}\unicode[STIX]{x1D70E}_{r}\unicode[STIX]{x1D753}_{r}\unicode[STIX]{x1D74D}_{r}\right\Vert }\nolimits_{2}}{\Vert \unicode[STIX]{x1D63F}\Vert _{2}},\end{eqnarray}$$

where all the modal contributions are first sorted in descending order $\unicode[STIX]{x1D70E}_{a}>\unicode[STIX]{x1D70E}_{b},\forall a>b$ . For the POD, (4.2) reduces to $\text{Err}[r_{c}]/\Vert \unicode[STIX]{x1D63F}\Vert _{2}=\unicode[STIX]{x1D70E}_{{\mathcal{P}}r_{c+1}}/\unicode[STIX]{x1D70E}_{{\mathcal{P}}_{1}}$ since the $L^{2}$ norm of a matrix is equivalent to its largest singular value (see (2.11)). This is not the case for other decompositions, for which the spatial basis is not orthogonal.

Figure 5.

Convergence of the POD, cDMD, mPOD and DFT for the illustrative synthetic test case in terms of $L^{2}$ relative error, defined as in (4.2).

The convergence results of this test case are shown in figure 5. As expected, the POD produces no reconstruction error when $r_{c}\geqslant 3$ , while the convergence of frequency based methods is considerably weaker. The finite duration of the first mode, the aperiodicity of the second and the nonlinear modulation of the third are not easily derivable as a linear combination of infinite harmonics, and the convergence of the DFT requires up to $150$ modes to bring the relative error below $10\,\%$ . The DMD performs better than the DFT, although yet not comparable to the POD. The mPOD of this synthetic test case is performed with $M=4$ , $\boldsymbol{F}_{V}=[0.4,10,20]$ and leads to a convergence much closer to that of the POD than the harmonic decompositions. To analyse the detection capabilities of these decompositions, several representative modes for each are now considered. First, the energy contributions of the DFT $(\text{diag}(\unicode[STIX]{x1D72E}_{{\mathcal{F}}}))$ and the cDMD ( $\text{diag}(\unicode[STIX]{x1D72E}_{{\mathcal{D}}})$ ) are mapped in the complex plane in figure 6. By definition, all the DFT modes lie in the unit circle and are Hermitian symmetric; DMD eigenvalues can lay inside the unit circle (if $|\unicode[STIX]{x1D706}_{r}|<1$ ) and exponentially decay, or outside (if $|\unicode[STIX]{x1D706}_{r}|>1$ ) and exponentially grow.

Figure 6.

Amplitude spectra of the DFT (a) and cDMD (b) spectra mapped in the complex plane. Both decompositions have temporal evolutions distributed along the unit circle. For the DFT this is true by definition; for the DMD this is the result from the eigenspectra of the companion matrix assembled from the data.

Although the two spectra are similar, this additional degree of freedom gives to the DMD better convergence and reduced windowing problems, especially on the lower frequency. The spatial structures for the first $6$ odd DFT modes (the even one are phase-shifted copies) are shown in figure 7. The amplitude of these modes is indicated with a square marker in figure 6 while their temporal structures, being simple harmonics, are not shown. The frequency resolution of the DFT is of the order of $f_{s}/n_{t}\approx 0.2$ , with the leading modes pulsing at $f_{n}=[0.19,7.03,15.04]$ . The spectral leakage is more evident at the lowest frequency, for which a full period of the oscillation is not captured. These modes have no time localization capabilities: none of these can locate in time the identified structures.

Figure 7.

Spatial structures $\unicode[STIX]{x1D753}_{{\mathcal{F}}}[\boldsymbol{x}_{i}]$ of the first six DFT modes, with amplitude indicated in figure 6 with a square marker.

Figure 8.

Spatial structures $\unicode[STIX]{x1D753}_{{\mathcal{D}}}[\boldsymbol{x}_{i}]$ of the first three cDMD modes, with amplitude indicated in figure 6 with a square marker.

This localization limit, stemming from the Heisenberg uncertainty principle (see Daubechies Reference Daubechies1992) is the main limit responsible for the poor convergence, which in turns results in a high redundancy in the spatial modes: many modes (besides the conjugate pairs) have identical spatial structures (e.g.  $\unicode[STIX]{x1D753}_{{\mathcal{F}}_{3}}\approx \unicode[STIX]{x1D753}_{{\mathcal{F}}_{9}}$ , $\unicode[STIX]{x1D753}_{{\mathcal{F}}_{5}}\approx -\unicode[STIX]{x1D753}_{{\mathcal{F}}_{11}}$ ) associated, in the time domain, with harmonic corrections to represent finite duration events.

The DMD modes, having a complex frequency, give a valid alternative when describing events that vanish from the initial condition, but not for events that occur after the initial step. A natural way of overcoming such limit is to break the data into multiple windows and perform the DMD in each of these as in the mrDMD (Kutz et al. Reference Kutz, Fu and Brunton2016b ), at the cost of significantly increasing the decomposition complexity as discussed in § 3.4.

Even in its simplest formulation, the cDMD improves the convergence and has better frequency resolution than the DFT. The cDMD spectrum highlights the three dominant modes pulsing at $f_{n}=[0.101,7.00,15.09]$ , and the corresponding spatial structures, shown in figure 8, correctly localize in space these pulsations. The time localization, of course, remains impossible as for any frequency-based method.

Figure 9.

Results from the POD of the synthetic test case. (a,b) Shows the normalized temporal correlation matrix $\unicode[STIX]{x1D646}=\unicode[STIX]{x1D63F}^{\dagger }\unicode[STIX]{x1D63F}/n_{t}$ (a) and corresponding DFT $\widehat{\unicode[STIX]{x1D646}}=\unicode[STIX]{x1D733}_{{\mathcal{F}}}\unicode[STIX]{x1D646}\unicode[STIX]{x1D733}_{{\mathcal{F}}}/n_{t}$ (b). (c–e) Shows the three non-zero POD modes, including their spatial structures $\unicode[STIX]{x1D753}_{{\mathcal{P}}_{r}}$ and temporal structures $\unicode[STIX]{x1D74D}_{{\mathcal{P}}_{r}}$ .

By contrast, the POD avoids both the time localization and the convergence problems of DFT and DMD, but produces a severe spectral mixing between different modes. The results of the POD are collected in figure 9, which includes the temporal correlation matrix $\unicode[STIX]{x1D646}$ and its DFT $\widehat{\unicode[STIX]{x1D646}}$ , as well as the three identified modes. The first carries the information on the time localization of each event, as the superimposition of different patterns is visible (the onset of the temporal evolution $\boldsymbol{T}_{1}$ is marked in the figure). The second carries the information on the frequency localization, with the three dominant harmonics producing peaks mirrored along the centre in both the diagonal and the anti-diagonal. The spectral leakage produces cross-like patterns centred in these peaks, with larger extension at the lowest frequencies (affected by the strongest windowing effect).

As this simple test case is not stationary, the correlation matrix is not Toeplitz circulant, as enforced in the SPOD proposed by Sieber et al. (Reference Sieber, Paschereit and Oberleithner2016). On the contrary, its localized patterns are the footprint of the different phenomena in the dataset and give to the POD the time localization capabilities that are not achievable by DMD/DFT.

The temporal correlation matrix of this synthetic test case, however, is almost defective, because it is close to having repeated eigenvalues. The three detected singular values are $\unicode[STIX]{x1D72E}_{{\mathcal{P}}}=\text{diag}(297,290,284)$ : although the energy contribution of the introduced mode is identical, the non-orthogonality of the temporal structures requires a slight energy re-distribution. Yet, the singular values are sufficiently close to the limit $\unicode[STIX]{x1D70E}_{1}=\unicode[STIX]{x1D70E}_{2}=\unicode[STIX]{x1D70E}_{3}$ to approach the condition of wholly undetermined POD.

The eigenvectors of $\unicode[STIX]{x1D646}$ , computed using the Matlab commands $\mathtt{eig}$ are shown in figure 9. These modes are a linear combination of the introduced ones, with a severe mixing of different phenomena over different modes. To further highlight the uniqueness problem, figure 10 shows that the temporal evolutions in figure 4 are almost eigenvectors of $\unicode[STIX]{x1D646}$ : the matrix multiplication $\unicode[STIX]{x1D646}\boldsymbol{T}_{k}$ leads to vectors that are close to simple multiples of the temporal evolution $\boldsymbol{T}_{k}$ . The more $\unicode[STIX]{x1D646}$ is close to being defective, the larger the number of vectors that are almost its eigenvectors and thus POD structures.

Figure 10.

Test for POD uniqueness for the proposed test case: the temporal structures in figure 4 are almost eigenvectors of the temporal correlation matrix, since $\unicode[STIX]{x1D646}\boldsymbol{T}_{k}\sim \unicode[STIX]{x1D706}_{k}\boldsymbol{T}_{k}$ .

Figure 11.

Fourier spectra of the correlation matrices $\unicode[STIX]{x1D646}_{{\mathcal{L}}_{1}}$ , $\unicode[STIX]{x1D646}_{{\mathcal{H}}_{1}}$ and $\unicode[STIX]{x1D646}_{{\mathcal{H}}_{2}}$ for the three scales identified via filter bank in the synthetic test case. (b) Shows the dominant eigenvector for each scale contribution. These are then assembled to construct the temporal structures $\unicode[STIX]{x1D74D}_{{\mathcal{M}}}$ as described in § 3.4.

It is therefore clear that it is possible to obtain an orthonormal basis which is more representative of the coherent patterns in the data by constraining the POD within different frequency intervals (scales). In the proposed mPOD, this is achieved via the MRA of the correlation matrix, splitting the correlation matrix in $M=4$ contributions with $\boldsymbol{F}_{V}=[0.4,10,20]$ . The DFT of each scale is shown in figure 11(a); the dominant eigenvector for each scale is shown in figure 11(b). The spectral separation of the three contributions is achieved and the temporal structures recover the evolution of each mode. As a result, the spatial structures are identical to the introduced ones and are therefore not shown. At the cost of a minor loss in the decomposition convergence, if compared to the POD, the mPOD correctly identifies the coherent structures introduced and assigns a different mode to each of them.

5.1 Dataset description

As a second test case, we propose a complex dataset which features the main physical mechanisms of complex fluid flows (i.e. nonlinear advection and diffusion), yet allowing for sufficient control over its dynamics so as to validate the decompositions. This test case is the numerical simulation of an advection–diffusion problem with random, pseudo-random and coherent sources. The problem is derived from the incompressible 2-D Navier–Stokes (NS) equations in the vorticity–streamfunction formulation. This formulation is common in the study of large-scale geophysical systems (Mariotti, Legras & Dritschel Reference Mariotti, Legras and Dritschel1994; Flór, van Heijst & Delfos Reference Flór, van Heijst and Delfos1995; Trieling, van Wsenbeeck & van Heijst Reference Trieling, van Wsenbeeck and van Heijst1998), and consists of a set of two equations. The first is the transport equation of the vorticity $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ , where the velocity field $\boldsymbol{u}=(u,v)$ is described by a streamfunction $\unicode[STIX]{x1D701}$ , defined so that $\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D701}=u$ and $\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D701}=-v$ :

(5.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D714}=-{\mathcal{J}}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D701})+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714}+S_{\unicode[STIX]{x1D714}},\end{eqnarray}$$

where ${\mathcal{J}}(\unicode[STIX]{x1D714},\unicode[STIX]{x1D701})=\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D701}\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D714}-\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D701}\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D714}$ is the Jacobian determinant of the solution vector $(\unicode[STIX]{x1D714},\unicode[STIX]{x1D701})$ ; $Re$ is the Reynolds number weighting the importance of advection and diffusion; $S_{\unicode[STIX]{x1D714}}$ is the spatially and temporally varying source or sink of vorticity. The second equation is the Poisson problem linking vorticity and stream function:

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D735}\times \boldsymbol{u}=\unicode[STIX]{x2202}_{x}v-\unicode[STIX]{x2202}_{y}u=-\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D701}.\end{eqnarray}$$

Sources of vorticity $S_{\unicode[STIX]{x1D714}}$ in (5.1) are typically the vortex stretching/tilting due to the velocity gradient tensor, baroclinic effects due to the misalignment between pressure and density fields and non-conservative body forces (e.g. the Coriolis effect). For the purpose of this work, we limit our interest to the nonlinear evolution of prescribed sources (regardless of their origin) and to the ability of data decomposition to correctly identify their spatial and temporal structures.

The solution of (5.1)–(5.2) is computed using a standard finite difference approach (Hoffmann & Chiang Reference Hoffmann and Chiang1993; Kutz Reference Kutz2013), combining a fast Poisson solver for (5.2) and a Runge–Kutta time marching scheme for (5.1). The computational domain consists of a square Cartesian grid of $n_{s}=128\times 128$ points over a domain $\boldsymbol{x}_{i}\in [-20,20]\times [-20,20]$ , while the solution vectors is exported in a time domain $t_{k}\in [0,20.48]$ with a total of $n_{t}=2048$ sample points. Periodic boundary conditions are set in all the boundaries.

The source term is treated explicitly, and consists of coherent ( $S_{C\unicode[STIX]{x1D714}}$ ), pseudo-random $(S_{R\unicode[STIX]{x1D714}})$ and fully random $(S_{{\mathcal{N}}\unicode[STIX]{x1D714}})$ contributions:

(5.3) $$\begin{eqnarray}S_{\unicode[STIX]{x1D714}}(\boldsymbol{x}_{\boldsymbol{i}},t_{k})=\overbrace{\mathop{\sum }_{l=1}^{4}A_{C_{l}}S_{C_{\unicode[STIX]{x1D719}_{l}}}(\boldsymbol{x}_{\boldsymbol{i}})S_{C_{\unicode[STIX]{x1D713}_{l}}}(t_{k})}^{S_{C\unicode[STIX]{x1D714}}}+\overbrace{A_{R}S_{R\unicode[STIX]{x1D719}}(\boldsymbol{x}_{i})S_{R\unicode[STIX]{x1D713}}(t_{k})}^{S_{R\unicode[STIX]{x1D714}}}+\overbrace{A_{{\mathcal{N}}}{\mathcal{N}}(\boldsymbol{x}_{i},t_{k})}^{S_{{\mathcal{N}}\unicode[STIX]{x1D714}}}.\end{eqnarray}$$

The four spatial $S_{C_{\unicode[STIX]{x1D719}_{l}}}(\boldsymbol{x}_{\boldsymbol{i}})$ and temporal $S_{C_{\unicode[STIX]{x1D713}_{l}}}(t_{k})$ structures of the coherent source $S_{C\unicode[STIX]{x1D714}}$ are shown in figure 12. The spatial structures are all constructed as sets of Gaussian vortex pairs of different size and position. The general form is

(5.4) $$\begin{eqnarray}S_{C_{\unicode[STIX]{x1D719}_{l}}}(\boldsymbol{x}_{\boldsymbol{i}})=\pm \exp \left(\frac{(\boldsymbol{x}_{\boldsymbol{i}}\pm \boldsymbol{x}_{Cl}\pm \boldsymbol{d}_{l})^{2}}{\unicode[STIX]{x1D70E}_{l}^{2}}\right),\end{eqnarray}$$

where the distance between the vortices in each pair are $\boldsymbol{d}_{l}=(d_{xl},d_{yl})=([0,0.5,0.5,0.5],[5.5,0,0,0])$ , the pair locations are $\boldsymbol{x}_{Cl}=(\boldsymbol{x}_{Cl},\boldsymbol{x}_{Cl})$ with $\boldsymbol{x}_{Cl}=[0,10,5,1]$ and the spreading factors are $\unicode[STIX]{x1D70E}_{l}=[49,0.5,0.5,0.5]$ .

Figure 12.

Spatial ( $S_{C_{\unicode[STIX]{x1D719}_{l}}}$ , a) and temporal ( $S_{C_{\unicode[STIX]{x1D713}_{l}}}$ , b) structures of the coherent source term $S_{C\unicode[STIX]{x1D714}}$ in (5.3).

The first source is composed of a single vortex pair of equal sign, which naturally merges producing a large recirculation, with angular velocity controlled by its inertia and the sign of the associated temporal structure $S_{C_{\unicode[STIX]{x1D713}_{1}}}$ . This temporal structure is composed of a set of Gaussian pulses $A_{i}\exp ((t_{k}-t_{i})^{2}/\unicode[STIX]{x1D70E}_{t1}^{2})$ of amplitude varying between $\pm [1,1.5]$ and equal widths of $\unicode[STIX]{x1D70E}_{t1}=0.2$ . The large size of the vortex gives too much inertia to this structure to follow such sharp variations, but the change of sign in $S_{C_{\unicode[STIX]{x1D713}_{1}}}$ (for $t>12$ ) changes the large-scale rotation from counter-clockwise to clockwise towards the end of the simulation. This process is much slower than all the other co-occurring.

The second and third coherent sources differ, in space, only in the distance between the vortex dipoles and have temporal structures consisting of Gaussian modulated harmonics:

(5.5) $$\begin{eqnarray}S_{C_{\unicode[STIX]{x1D713}_{2,3}}}=\sin [2\unicode[STIX]{x03C0}f_{2,3}(t_{k}-\unicode[STIX]{x1D6FC}_{2,3})]\exp \left[\frac{-(t_{k}-\unicode[STIX]{x1D6FD}_{2,3})^{2}}{\unicode[STIX]{x1D70E}_{t2,3}^{2}}\right],\end{eqnarray}$$

with frequencies $f_{2,3}=[2,5]$ , equal width $\unicode[STIX]{x1D70E}_{t2,3}=4$ but different phase delay for both the harmonic ( $\unicode[STIX]{x1D6FC}_{2,3}=[14,6]$ ) and the modulation terms ( $\unicode[STIX]{x1D6FD}_{2,3}=[17,11]$ ). These differences in the phase delay create sharp variations that produce a regular vortex shedding. Finally, the fourth contribution is spatially similar to the second and third, with dipoles placed much closer to the centre of the domain, and a temporal structure consisting of a single Gaussian impulse located a $t_{k}=5.2$ and width $\unicode[STIX]{x1D70E}_{t4}=0.05$ .

The pseudo-random source term $S_{R\unicode[STIX]{x1D714}}$ in (5.3) produces a periodic injection of random vortices. The spatial structure $S_{R\unicode[STIX]{x1D719}}$ consist of a set of $300$ Gaussian vortices, randomly located with widths in the range $[0,1]$ , and amplitudes in the range $[-1,1]$ . An exemplary set of vortices representing one of these structures is shown in figure 13(a). The temporal evolution $S_{R\unicode[STIX]{x1D713}}$ , shown in figure 13(b), consists of $10$ equal Gaussian impulses with width $0.1$ and unitary amplitude and constant time shift. As the introduced spatial structures change randomly from one pulse to the other, this source contribution creates a chaotic background which has no spatial coherence but exhibits a temporal regularity. Finally, the white noise term $S_{{\mathcal{N}}_{\unicode[STIX]{x1D714}}}$ in (5.3) is introduced in the post-processing step, using a random number generator varying between $0$ and $1$ to generate the spatio-temporal noise ${\mathcal{N}}(\boldsymbol{x}_{\boldsymbol{i}},t_{k})$ .

Figure 13.

Exemplary spatial and temporal structures of the pseudo-random source term in (5.3). At each of the time pulse, a random distribution of 300 Gaussian vortices is introduced.

The coefficients $A_{C_{l}}$ ( $l\in [1,4]$ ), $A_{R}$ and $A_{{\mathcal{N}}}$ in (5.3) set the relative importance of these three contributions and are taken as $A_{l}=[9,610,1100,2500]$ , $A_{R}=300$ and $A_{{\mathcal{N}}}=7$ . These coefficients were estimated, by trial and error, to give approximately the same importance to the four coherent sources, taking into account the strong damping effect produced by diffusion (especially on the smaller and faster scales), obtained by setting the Reynolds number in (5.1) to $Re=1$ . The initial vorticity field is a random set of vortices similar to those introduced by the pseudo-random source term.

Figure 14.

Exemplary snapshots from the advection–diffusion problem in (5.1). These selected snapshots describe some important instants in the evolution of the vorticity field considered in this test case. An animation of this test case is available as supplemental Movie 2.

Figure 15.

Time evolution of the global energy content (norm of each vorticity realization), extracted from the diagonal of the temporal correlation matrix $\unicode[STIX]{x1D646}$ .

To visualize the complexity of the dataset obtained, figure 14 collects $10$ selected snapshots of the vorticity field, each corresponding to the time indicated in the title. The description of the vorticity evolution in these snapshots is well complemented by the evolution of the energy signal in time, that is the norm of each vorticity snapshot in the diagonal of the temporal correlation matrix $K$ . This signal is shown in figure 15, and the occurrence of some key events are indicated with an annotation.

For $t_{k}<2$ , the dataset describes the chaotic interaction of the initial field, which strongly diffuses to almost vanish (cf. figure 15). Before this occurs, the first coherent source gives its first pulse, producing a large-scale counter-clockwise rotation in the background. Until $t_{k}=5$ , the three injections of random vortices produces three noticeable increase in the energy content against diffusion, but no coherent patterns can be identified at this stage.

At $t_{k}\approx 5.15$ , the injection of the source term $S_{C\unicode[STIX]{x1D713}_{4}}$ stands out from the chaotic motion, and the energy in the vorticity field reaches its highest peak. The four vortex pairs immediately merge into two identical dipoles of equal vorticity, which repel from the centre of the domain in opposite direction and speed. However, due to the background rotation produced by the first source and the chaotic interaction with the other random injections from the pseudo-random term $S_{R\unicode[STIX]{x1D714}}$ , these two dipoles do not travel the same distance and are broken and wholly dissipated by $t\approx 6.5$ .

At $t_{k}>6$ , the third coherent source begins its pulsation. At $t_{k}\approx 7$ this generates four wakes which shed at a constant frequency. The overall vorticity in the domain increases and the large-scale clockwise rotation involves the four wakes in a closed circular pattern. Before this large-scale rotation vanishes, at $t_{k}\approx 14$ , the second coherent source starts to pulse at a lower frequency, producing other four trains of vortices, initially also rotating counter-clockwise. Until approximately $t_{k}=17$ the second and third source coexist and produces patterns of similar energy content but different frequencies.

From $t_{k}>17$ , the second source vanishes and the vortex pair introduced by the last three pulses of the first coherent source slowly changes the direction of the rotation. This breaks the circular shedding into two vertical ones before a new circular shedding rotating in the opposite direction is produced. The change of direction and the strong diffusion reduce the energy content of the flow until only a large-scale vortex, rotating clockwise, slowly dissipates.

5.2 Modal analysis

The relative $L^{2}$ error convergence (4.2) for the DFT, cDMD, sDMD, POD and mPOD are shown in figure 16. The cDMD and the sDMD formulations yield very similar (yet not identical) results if the full POD basis is used in the construction of the reduced propagator in (2.19).

This test case strongly highlights the convergence limits of the DMD for noisy data with a nonlinear evolution. As shown in § 2.3, the amplitude coefficients $\unicode[STIX]{x1D70E}_{{\mathcal{D}}r}$ in the DMD are projections of the initial conditions onto the DMD basis: large amplitudes are assigned to modes which are more strongly linked to the initial conditions, even if these vanish immediately afterwards – in this case due to the strong diffusion. As a result, both algorithms produce DMD eigenvalues lying well inside the unit circle and exponentially decaying. Since decaying modes with real frequencies are not orthogonal, the DMD convergence is drastically compromised and yields, for $r_{c}<200$ , an $L^{2}$ approximation error which is larger than the norm of the dataset itself.

Figure 16.

Decomposition convergence for DFT, cDMD, sDMD, POD and mPOD for the vorticity field in the nonlinear advection–diffusion problem considered in this section.

The problem is particularly evident for the sDMD with the reduced propagator $\boldsymbol{{\mathcal{S}}}$ expressed in (2.19). Figure 17 shows the eigenvalues distribution in the complex plane for reduced propagators $\boldsymbol{{\mathcal{S}}}$ of different dimensions $r_{c}\times r_{c}$ , obtained by projecting the full propagator onto the first $r_{c}$ POD modes. For $r_{c}\ll n_{t}$ , most of the eigenvalues are well inside the unit circle, and their contribution to the decomposition is limited to the first few time steps. Since each of these modes is normalized, the shortest is their duration, the largest the associated amplitude. As $r_{c}$ approaches $n_{t}$ and the DMD basis approaches completeness, the number of modes with exponential decays is reduced, and the sDMD approaches the cDMD.

Figure 17.

Eigenvalue distribution for reduced propagators $\boldsymbol{{\mathcal{S}}}\in \mathbb{C}^{r_{c}\times r_{c}}$ in (2.19) obtained by projecting the full propagator $\unicode[STIX]{x1D64B}$ onto the first $r_{c}$ POD modes. The smaller the propagator, the more the decomposition convergence is limited within a smaller number of time steps from the initial data. At $r_{c}=n_{t}-1$ , the sDMD approaches the cDMD.

Figure 18.

Amplitude spectra of the DFT (a), and cDMD (b) mapped in the complex plane. While both frequency based decompositions are symmetric along the imaginary axis, the DMD is sensitive to the high noise level and the strong diffusion of the initial dataset, which results in all the leading modes lying inside the unit circle.

The amplitudes of the modes for the DFT and the cDMD are mapped in the complex plane in figure 18. The DFT spectrum shows that the energy contribution of the mean flow is stronger than in the previous test case, and highlights three dominant frequencies together with a wide range of modes with nearly identical energy content. These are produced by the high random noise level in the dataset and the presence of impulsive events. For the cDMD, the random noise makes the calculation of the companion matrix more difficult, producing an error of $92\,\%$ of the $L^{2}$ norm of the last realization in the minimization $\Vert \boldsymbol{d}_{n_{t}}-\unicode[STIX]{x1D63F}_{1}\boldsymbol{c}\Vert _{2}$ .

Figure 19.

(a) Spatial structures $\unicode[STIX]{x1D753}_{{\mathcal{F}}}[\boldsymbol{x}_{i}]$ of three dominant DFT modes, pulsing, from left to right, at frequencies of $f_{n}=[0.12,2.01,5.02]$ . These corresponds to leading frequencies of the coherent source term. (b) Spatial structures $\unicode[STIX]{x1D753}_{{\mathcal{D}}}[\boldsymbol{x}_{i}]$ of the dominant sDMD modes with frequencies $[0.08,1.97,4.99]$ , that is the ones closer to those dominating the DFT spectra.

Figure 20.

Same as figure 9 but for the second test case discussed in this section.

The three dominant frequencies in the DFT are $f_{n}=[0.12,2.01,5.02]$ ; the first is linked to the large-scale motion produced by the first coherent source, the second and third are related to the frequencies of the second and third sources pulsing at $f_{2,3}=[2,5]$ . Figure 19 (a) shows the corresponding spatial structures, which are well related to the introduced ones. The three cDMD modes with unitary modulus having frequencies closer to these DFT modes – also captured by all the sDMD in figure 17 – are shown in figure 19(b). The amplitude of these modes is labelled with a red marker in the complex plane of figure 18, and is much lower than the dominant ones, associated with the vanishing initial conditions.

As for the synthetic test case in § 4, both the DFT and the DMD have sufficient frequency resolution to correctly localize periodic phenomena in the frequency domain and – thanks to the projection in (2.6) – in the space domain. However, as no temporal localization is possible, the convergence of these decompositions is remarkably poorer than the POD and results in a high redundancy of the spatial basis. This is particularly evident in the description of the impulsive event produced by the fourth coherent source, the footprint of which appears in a large number of DFT and DMD modes.

The POD has complementary advantages and limitations. Figure 20 shows the results of the POD analysis, including the temporal correlation matrix $\unicode[STIX]{x1D646}$ and its DFT spectra $\widehat{\unicode[STIX]{x1D646}}$ , together with the spatial and the temporal structures of six representative POD modes.

The temporal correlation matrix $\unicode[STIX]{x1D646}$ gives information on the time localization of different events, featuring two significant chessboard patterns corresponding to the harmonic contributions of the second and third sources shown in figure 12, and rectangular areas corresponding to the time interval of major activity of the first large-scale coherent source. In the diagonal, a clear peak is located at $t_{k}=5.2$ , corresponding to the introduction of the fourth impulsive source, while the random noise gives a pedestal correlation level along the whole diagonal.

Figure 21.

Results from the MRA of the correlation matrix $\unicode[STIX]{x1D646}$ , split into four contributions by a filter bank with frequency splitting vector $\boldsymbol{F}_{V}=[1,3,8]$ . (a) Shows the temporal correlations; (b) shows their frequency content.

Figure 22.

Spatial and temporal structures for six representative mPOD modes for the nonlinear vortex test case presented in this section.

The first mode captures the large-scale motion produced by the first coherent source, with a temporal structure slightly polluted by white noise. Uncorrelated noise is uniformly distributed over the entire POD spectra (Mendez et al. Reference Mendez, Raiola, Masullo, Discetti, Ianiro, Theunissen and Buchlin2017) and the low noise level present in this mode is due to the large dataset available. The second and third modes show a severe spectral mixing, as their energy contribution is similar ( $\unicode[STIX]{x1D70E}_{{\mathcal{P}}_{2}}=4272$ and $\unicode[STIX]{x1D70E}_{{\mathcal{P}}_{2}}=4433$ ). This problem, also illustrated on the simple synthetic test case in § 4, results in a combination of the features from the second and third coherent sources. Such mixing extends to other POD modes with different phase delays, to cancel with the summation of multiple modes. The fifth mode mixes three phenomena: a large-scale contribution to the rotation produced by the first source, a harmonic contribution due to the second source and a sharp peak due to the impulsive fourth source. The same occurs with the ninth mode, which mostly describes the diffusion of the initial data. Interestingly, the tenth mode well captures the impulsive event, although both the harmonic contribution of the second and third sources are present to a minor extent.

Figure 21 shows the results of the MRA analysis of the temporal correlation matrix $\unicode[STIX]{x1D646}$ , using $\boldsymbol{F}_{V}=[1,3,8]$ and computing the filter orders from (A 5). (a) Displays the four contributions of the identified scales and (b) their DFT. No significant phenomenon appears in the last contributions, which mostly contains a substantial portion of the noise energy, whereas the other three scales isolate the contributions of different sources. Being the overlapping of frequencies negligible, the eigenvectors of these matrices are almost orthogonal complements for $\mathbb{R}^{n_{t}}$ , and the resulting full basis $\unicode[STIX]{x1D733}_{{\mathcal{M}}}^{0}$ (3.15) requires a minor re-orthogonalization via $\unicode[STIX]{x1D64C}\unicode[STIX]{x1D64D}$ factorization.

The spatial and temporal structures of six representative mPOD modes are shown in figure 22. These modes are similar to the POD modes in figure 20, but with no frequency overlapping. The first mode is a cleaner version of the first POD mode, with negligible random noise; the third and fourth capture the contributions of the second and third coherent sources with no mixing. The sixth mode describes the diffusion of the initial data while the tenth and the thirteenth are related to the impulsive event. The frequency bounding makes the decomposition stronger than the standard POD versus the random noise since only a minor portion of its contribution remains in each scale. Yet, the capabilities of capturing impulsive phenomena (also distributed over the entire frequency spectra) are reasonably preserved: traces of this impulse event are present in all the scales (the tenth mode is the eigenvector of $\unicode[STIX]{x1D646}_{{\mathcal{L}}}$ ; the thirteenth is eigenvector of $\unicode[STIX]{x1D646}_{{\mathcal{H}}_{3}}$ ), but no significant overlapping between this event and the other harmonics is produced.

6.1 Dataset description

As a third test case, we propose an experimental dataset obtained via time-resolved particle image velocimetry (TR-PIV) on a classical, statistically stationary flow configuration: a plane jet impinging normally onto a flat plate. Impinging jet flows have received considerable attention for their application in drying, cooling or heating processes (Zuckerman & Lior Reference Zuckerman and Lior2006; Buchlin Reference Buchlin2011; Specht Reference Specht, Tsotsas and Mujumdar2014) and the POD of PIV data has been successfully used by many authors to identify the resulting coherent structures (Kim et al. Reference Kim, Min, Oh, An, Seoudi, Chun and Lee2007; Hammad & Milanovic Reference Hammad and Milanovic2009; Charmiyan et al. Reference Charmiyan, Azimian, Shirani, Aloui, Degouet and Michaelis2017; Pieris et al. Reference Pieris, Zhang, Yarusevych and Peterson2017).

These flows are characterized by largely different time scales in three main regions: the free jet, close to the nozzle exit; the stagnation region, close to the flow impact; the wall jets, released parallel to the wall from impact (Gutmark, Wolfshtein & Wygnanski Reference Gutmark, Wolfshtein and Wygnanski1978).

Figure 23.

Schematic of the experimental set-up for TR-PIV analysis of an impinging plane jet, showing the region of interest (ROI). The dimensions are in mm.

The experimental set-up is shown in figure 23. The rectangular outlet section of the nozzle has an opening of $H=4~\text{mm}$ and an aspect ratio of $H/W=62.5$ . The chamber gauge pressure $P_{N}$ is monitored by a piezoresistive pressure transducer AMS5812 connected to a pressure tap and the flat surface consist of a polished aluminium plate placed at a distance of $Z=30\pm 1\,\%~\text{mm}$ . The experiment is carried out with $P_{N}=30\pm 1\,\%~\text{Pa}$ leading to an average jet flow velocity of $U_{0}=6.5\pm 2\,\%~\text{m}~\text{s}^{-1}$ , that is a Reynolds number of $Re\approx U_{0}H/\unicode[STIX]{x1D708}=1733$ considering a kinematic viscosity of $\unicode[STIX]{x1D708}=15~\text{cSt}$ .

The seeding tracers for the PIV measurements are microdroplets of mineral oil (Ondina Shell 917), produced by a Laskin nozzle (PIVTEC45-M). These are introduced both in the stagnation chamber and in a second injection on one side of the nozzle to seed the entrainment flow. The laser source is a dual diode-pumped ND:YLF laser (Quantronix Darwin Duo 527), reaching the test section via a light sheet probe (LSP, from ILA). The videos are acquired by a CMOS Photron FASTCAM SA1, positioned at about 30 cm and equipped with a zoom objective Tamron CZ-735 70–300 mm to provide an image scaling factor of $24.1\pm 0.1~\text{pixel}~\text{mm}^{-1}$ .

The investigated region of interest (ROI) is a rectangular area of $20~\text{mm}\times 38~\text{mm}$ at a distance of 10 mm from the nozzle outlet (cf. figure 23). The video sequence is made of $n_{t}=2728$ images pairs of $1000\times 740$ pixels, acquired in frame-straddling mode at $f_{s}=2~\text{kHz}$ . The separation time between the frames is set to $\unicode[STIX]{x0394}t_{p}=70~\unicode[STIX]{x03BC}\text{s}$ to produce a maximum displacement, in the jet flow, of approximately $8$ pixels. The PIV evaluation is carried out with the Matlab package PIVlab (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014), using standard iterative multi-step interrogation (four passes from $128\times 128$ to $8\times 8$ windows) with spline window deformation, 2-D Gaussian sub-pixel interpolation and vector validation via a median test. The resulting vector field consists of $n_{s}=60\times 114$ spatial realizations sampled on a uniform Cartesian grid.

An exemplary instantaneous flow field is shown in figure 24, together with the location of six probes. The normalized power spectral density computed from the velocity magnitude in these locations is shown in figure 24(b) as a function of the dimensionless frequency $St=f_{n}H/U_{0}$ and normalized by the corresponding maximum value. These spectra are computed using the Welch’s method (Welch Reference Welch1967) with Hanning windows of $n_{t}/4$ width and $25\,\%$ overlapping. The first probe $P_{1}$ is located in the free-jet flow, where the shear layer instabilities develop unaffected by the presence of the impinged wall. At this location, most of the energy is centred around a peak at $St\approx 0.3$ , in agreement with other studies on the free plane jet in the transitional regime (Thomas & Goldschmidt Reference Thomas and Goldschmidt1986; Suresh et al. Reference Suresh, Srinivasan, Sundararajan and Das2008). This corresponds to the formation of large vortices (primary vortices) which originates from a Kelvin–Helmholtz (K–H) instability and evolve into roller-like structures (Mumford Reference Mumford1982).

As the flow approaches the stagnation region (probes $P_{2}$ , $P_{3}$ ), the frequency spectra broaden and moves towards lower frequencies due to the vortex merging and the velocity decay, as also reported by Pieris et al. (Reference Pieris, Zhang, Yarusevych and Peterson2017) for impinging jet flows at higher Reynolds numbers. Far from the stagnation region, the wall jet flow ( $P_{4}$ , $P_{5}$ , $P_{6}$ ) is governed by a much lower frequency and the unsteady separation due to the intermittent passage of the primary vortices (Didden & Ho Reference Didden and Ho1985; Hadžiabdic̀ & Hanjalic̀ Reference Hadžiabdic̀ and Hanjalic̀2008).

Figure 24.

(a) Exemplary snapshot of velocity field from TR-PIV of an impinging jet flow, with contour plot displaying the magnitude of the velocity in $\text{m}~\text{s}^{-1}$ , and location of six representative probes. An animation of this test case is available as supplemental Movie 3. (b) Normalized power spectral densities in the probes in (a), versus Strouhal number.

6.2 Modal analysis

Figure 25 shows the $L^{2}$ convergence in (4.2) for POD, cDMD, DFT and two mPOD with frequency splitting vectors. The first, in terms of Strouhal numbers, consists of $M=4$ with $\boldsymbol{F}_{V1}(m)=[0.1,0.2,0.5]$ ; the second consists of $M=50$ with $\boldsymbol{F}_{V2}(m)=0.01+m\unicode[STIX]{x0394}f$ with equal bandwidth $\unicode[STIX]{x0394}f=0.02$ . As discussed in § 3.2, the finer the frequency splitting, the more the mPOD approaches a DFT.

Figure 25.

$L^{2}$ convergence (4.2) of the POD, cDMD, DFT and two mPOD with different frequency splitting vectors for the TR-PIV field of the impinging jet flow.

The frequency spectra for DFT and cDMD are mapped in figure 26. The investigated jet flow is in statistically stationary conditions and the mean flow accounts for a significant portion of the dataset. For plotting purposes, therefore, the entry corresponding to the zero frequency is removed. The spectra of the DFT shows that the investigated test case is rather broadband, with few frequencies reasonably standing over the others. Despite the significant random noise and turbulence contribution, the presence of a statistically converged time average results in a relative error $\Vert \boldsymbol{d}_{n_{t}}-\unicode[STIX]{x1D63F}_{1}\boldsymbol{c}\Vert _{2}/\Vert \boldsymbol{d}_{n_{t}}\Vert _{2}\approx 6\,\%$ in the minimization defining the companion matrix for the cDMD. The resulting DMD dominant modes identify the relevant frequencies more sharply than the DFT. Nevertheless, several eigenvalues lie inside the unit circle and the convergence of the cDMD, although not dramatically poor as in the previous test case, is below the DFT.

Figure 26.

Same as figures 6 and 18, but for the third test case presented in this section.

Figure 27.

Spatial structures of two representative DFT (a,b) and DMD (c,d) modes, with amplitudes shown in figure 26 with a red marker. The DFT modes are associated with dimensionless frequencies of $St=[0.033,0.299]$ ; for the DMD these are $St=[0.009,0.293]$ .

The spatial structures of two representative DFT and DMD modes are shown in figure 27. For the DFT, these structures are associated with dimensionless frequencies of $St=[0.033,0.299]$ and are thus expected to be linked to the large-scale oscillations of the wall jet and K–H structures produced in the free-jet region respectively (see figure 24). The constraint of a single harmonic is nevertheless too restrictive to isolate these mechanisms, and the corresponding spatial modes are not capable of fully identifying the related coherent patterns. The DMD spatial structures, pulsing at dimensionless frequencies of $St=[0.009,0.293]$ , appear less noisy than the DFT ones, thanks to the better frequency selection. The mode $\unicode[STIX]{x1D731}_{{\mathcal{D}}_{27}}$ describes the large-scale lifting of the wall jet, due to the interaction with the roller structures that in this region lose most of their momentum and mix with the entrainment flow. As expected, the mode $\unicode[STIX]{x1D731}_{{\mathcal{D}}_{29}}$ captures the K–H structures, vortical structures approximately 2 mm in diameter that arise from the primary shear layer instability of the free-jet flow.

The results of a standard POD analysis are collected in figure 28. Since the POD modes are stationary, it is interesting to show the frequency content $\widehat{\unicode[STIX]{x1D74D}}_{{\mathcal{P}}r}$ associated with each spatial structure, rather than their temporal evolution $\unicode[STIX]{x1D74D}_{{\mathcal{P}}r}$ . All the modes in the range $r=4$ to $r=11$ are related to travelling structures and arise therefore in pairs that have $\unicode[STIX]{x03C0}/2$ phase delay in space and time. Only the four representative modes are shown.

The temporal correlation matrix $\unicode[STIX]{x1D646}$ displays a regular pattern composed of a broad range of frequencies and $\widehat{\unicode[STIX]{x1D646}}$ displays a band-like structure. Yet, a dominant frequency around 12 Hz ( $St=0.0074$ ) is visible. This corresponds to the large-scale flow motion in the wall jet flow, which covers a larger portion of the domain and, therefore, gives a more significant contribution to the overall correlation level. As the data are statistically stationary, the constant vector $\unicode[STIX]{x1D74D}_{\unicode[STIX]{x1D707}}=\text{}\underline{1}/\sqrt{n_{t}}$ is almost an eigenvector of the correlation matrix, and thus the first POD mode reproduces closely the mean flow (not shown). For plotting purposes, the mean correlation is removed from $\widehat{\unicode[STIX]{x1D646}}$ .

The spatial structures identified by the POD are considerably less influenced by the uncorrelated portion of the dataset and highlight coherent features in the flow. Besides the mean flow, the dominant contribution $\unicode[STIX]{x1D753}_{{\mathcal{P}}_{2}}$ in figure 28 is mostly related to the large-scale pulsation in the wall jet (labelled with $I$ ) and its periodic separation, as expected from the spectra of the correlation matrix $\widehat{\unicode[STIX]{x1D646}}$ . However, this spatial structure also captures, to a minor extent, coherent structures in the free-jet region (labelled with $II$ ). These are expected to occur at much higher frequency (in the range $St=0.1{-}0.2$ ) as shown by the probes $P_{2}$ – $P_{3}$ in figure 24.

The frequency spectrum $\widehat{\unicode[STIX]{x1D74D}}_{{\mathcal{P}}_{2}}$ gives the final proof of spectral mixing for this POD mode: a minor energetic content is also present in the frequency range belonging to the flow structures of the free-jet flow, and the projection of the data along this temporal evolution propagates this mixing to the spatial domain. The same occurs in the mode $\unicode[STIX]{x1D753}_{{\mathcal{P}}_{3}}$ , which is mostly linked to the travelling vortices in the jet, but – to a minor extent – also to the wall jet separation. All the other modes suffer from the same problem, which reaches its worst case in the mode $\unicode[STIX]{x1D753}_{{\mathcal{P}}_{9}}$ . In addition to the two mechanisms previously described, this mode also captures the K–H structures produced further upstream, where they originate at $St\approx 0.3$ as clearly shown by the probe $P_{1}$ in figure 24. This test case shows that the spectral mixing problems of the POD, highlighted by the two test cases of §§ 4 and 5, can easily occur also in statistically stationary flows if these are governed by largely different scales with similar energy content.

As for the previous two test cases, the spectral separation produced by the MRA is sufficient to isolate the contribution of different scales at the cost of a minor reduction in the decomposition convergence (cf. figure 25). Four representative mPOD modes, including spatial structures and related frequency content, are shown in figure 29. These are constructed by setting $\boldsymbol{F}_{V1}(m)=[0.1,0.2,0.5]$ , producing the $L^{2}$ convergence labelled as mPOD1 in figure 25. The resulting mPOD modes are very similar to the POD counterparts, showing that these spectral constraints have a minor impact on the decomposition. Yet, no spectral mixing is produced. The second mode captures only the dynamics of the separating wall jet, the fourth is related only to the evolution of the primary vortex advected towards the wall, while the ninth captures only the formation of the primary vortices in the free-jet region. The seventh mode describes the large-scale oscillation of the jet which is associated with the stagnation region and the initial stage of the wall jet development.

By increasing the number of frequency bandwidths, one could further narrow the frequency content of each mode, improving the frequency resolution but decreasing the $L^{2}$ convergence of the method as shown in figure 25 by the second mPOD with finer frequency splitting vector. This can alter the ranking of different modes, giving priority to mechanisms that have narrower frequency content, and eventually move the mPOD towards the DFT for unreasonably fine frequency splitting vectors.

Figure 28.

Same as figures 9 and 20, but for the third test case discussed in this section and replacing the temporal structures $\unicode[STIX]{x1D74D}_{{\mathcal{P}}}$ with their normalized frequency content $|\widehat{\unicode[STIX]{x1D74D}}_{{\mathcal{P}}}|/\text{max}\{\widehat{\unicode[STIX]{x1D74D}}_{{\mathcal{P}}}\}$ .

Figure 29.

Selected mPOD modes corresponding to the POD modes in figure 28.