Multi-Scale Proper Orthogonal Decomposition of Complex Fluid Flows

Data-driven decompositions are becoming essential tools in fluid dynamics, allowing for tracking the evolution of coherent patterns in large datasets, and for constructing low order models of complex phenomena. In this work, we analyze the main limits of two popular decompositions, namely the Proper Orthogonal Decomposition (POD) and the Dynamic Mode Decomposition (DMD), and we propose a novel decomposition which allows for enhanced feature detection capabilities. This novel decomposition is referred to as Multiscale Proper Orthogonal Decomposition (mPOD) and combines Multiresolution Analysis (MRA) with a standard POD. Using MRA, the mPOD splits the correlation matrix into the contribution of different scales, retaining non-overlapping portions of the correlation spectra; using the standard POD, the mPOD extracts the optimal basis from each scale. After introducing a matrix factorization framework for data-driven decompositions, the MRA is formulated via 1D and 2D filter banks for the dataset and the correlation matrix respectively. The validation of the mPOD, and a comparison with the Discrete Fourier Transform (DFT), DMD and POD are provided in three test cases. These include a synthetic test case, a numerical simulation of a nonlinear advection-diffusion problem, and an experimental dataset obtained by the Time-Resolved Particle Image Velocimetry (TR-PIV) of an impinging gas jet. For each of these examples, the decompositions are compared in terms of convergence, feature detection capabilities, and time-frequency localization.

The development of data processing algorithms for detecting coherent features is continuously fostered by the availability of datasets with ever-growing resolutions, and by the ongoing big data revolution which permeates any area of applied science. Identifying the relevant features from high-dimensional datasets is the purpose of data-driven decompositions, which lay the foundation of model order reduction (MOR), data compression, filtering, and pattern recognition, and which are nowadays greatly enlarging the toolbox 2 M. A. Mendez et al at the researcher's disposal. Recent reviews on data-driven decomposition with emphasis on fluid dynamics are proposed by Taira et al. (2017) and , while an extensive overview on the impact of the big data revolution in turbulence research is presented by Pollard et al. (2017), while Duriez et al. (2017) give an overview of machine learning methods for flow control.
Data-driven decompositions represent the data as a linear combination of basis elements, referred to as modes, having spatial and temporal structures with a certain energy contribution. These decompositions can be classified into two major classes: energy-based and frequency-based. This work presents a generalized formulation of both approaches, analyzes their respective limits, and proposes a novel hybrid decomposition referred to as Multiscale Proper Orthogonal Decomposition (mPOD).
In the following subsections, we present a literature review on the energy-based and frequency-based data-driven decompositions ( §1.1) and outline the recent developments on hybrid methods ( §1.2). An overview of the main contribution of this study and the manuscript organization is presented in §1.3.

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 1977), 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 (2013), and the monographs from Jackson (1991) and Holmes et al. (2012).
The POD for data-driven analysis of turbulent flows was introduced to the fluid mechanic's community by Lumley (1967Lumley ( , 1970, 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 (1991), who presented a variant of the decomposition, referred to as bi-orthogonal decomposition (BOD, Aubry et al. 1991), which focused on the temporal structures of the modes. This link was made more explicit by Sirovich (1987Sirovich ( , 1989Sirovich ( , 1991, 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 1999;Fahl 2000;Chatterjee 2000), in contrast to formulations based on statistics or dynamical system theory (such as Holmes et al. 1997;Noack et al. 2003;Berkooz et al. 1993;Berkooz 1992). For an introduction to the SVD and a historical account on its development, the reader is referred to Golub (2013) and Stewart (1993).
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. Mallor et al. 2018;Schrijer et al. 2014), video analysis for adaptive masking or image pre-processing (Mendez et al. 2018a(Mendez et al. , 2017, flow sensing (Venturi & Karniadakis 2004;Willcox 2006), flow control (Brunton & Noack 2015;Bergmann & Cordier 2008), reduced order modeling (Deane et al. 1991), experimental validation of CFD data (Blanchard et al. 2014), data-driven identification of non-linear systems Loiseau et al. 2018) and more (see Pollard et al. 2017). Moreover, extended versions of the decomposition, in which the 4 M. A. Mendez et al

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 (2016) has recently discussed the quest for hybrid decompositions. The first single-harmonic variant of the POD was proposed by Lumley (1970) and formulated as a natural extension of the POD for stationary and homogeneous flows (see Glauser et al. 1987;George 1988). This approach consists in computing the POD structures as eigenvectors of the cross-spectral density matrix (see Gordeyev & Thomas 2000;Citriniti & George 2000;Gudmundsson & Colonius 2011) and is referred by several authors (e.g Picard & Delville 2000;Taira et al. 2017;Towne et al. 2018) to as Spectral Proper Orthogonal Decomposition (SPOD). Bourgeois et al. (2013) 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. (2013) 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. (2016) 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 et al. (2016), 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 frequencybased POD that is also referred to as SPOD by Picard & Delville (2000); Taira et al. (2017); Towne et al. (2018).
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 2001;Gray 2005) 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 1987;Gray 2005) and the diagonal filter used in the SPOD artificially forces such 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. 2016(Sieber et al. , 2017, 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.

Contribution of this work and manuscript organization
Section 2 presents a general matrix factorization that is common to all the decompositions, and set the POD, DFT, and DMD in this generalized form. This framework is then used in §3 to derive the novel Multiscale Proper Orthogonal Decomposition (mPOD).
The mPOD proposed in this work combines ideas from the windowed SPOD in Sieber et al. (2016) and the SPOD in Towne et al. (2018). Instead of filtering the correlation matrix, as in Sieber et al. (2016), the mPOD decomposes it into the contributions of different scales using the multiresolution (MRA) architecture from Wavelet theory. Instead of computing various eigenbases from the spectra of different portions of the data, as in Towne et al. (2018), the mPOD computes eigenbases on the correlation matrix of different scales. The derivation of the mPOD in §3 is carried out by analyzing 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.

Data Decompositions as Matrix Factorizations
Let D(x i , t k ) = D[i, k] ∈ R ns×nt be a matrix collecting the set of time realizations (snapshot) of a real variable. In the time domain, the data is sampled with a uniform temporal discretization {t k = (k − 1)∆t} nt k=1 . In the spatial domain, the data is sampled on a Cartesian grid x i ∈ R nx×ny , with i ∈ [1, . . . , n s ] a matrix linear index, n s = n C n x n y the number of points, and n C the number of components of the data in case of vector quantities (e.g. n C = 1 for a pressure field, n C = 3 for a 3D velocity field). After reshaping each snapshot of the data into a column vector d k [i] ∈ R ns×1 , the data matrix reads The criterion followed to reshape each snapshot into a vector is irrelevant, provided that the same is used when reshaping back the results of the decomposition. The scope of any discrete decomposition is to break this matrix into the sum of rank-1 contributions, referred to as modes, written in variable separated form. Each mode has a spatial structure φ r [i], a temporal structure ψ r [k] and an amplitude σ r : (2.2) A data matrix of rank rk(D) = min(n t , n s ) is said to have a full rank, while truncating the summation at r c < rk(D) leads to an approximationD of rank r c . In the MOR community, this approximation is referred to as reduced order model of the data.
In order to let the scalars σ r be fully representative of the energy of each mode, both the spatial and the temporal structures must have unitary energy (norm). The notion of energy is defined by the choice of inner product in the space and in the time domain. Although different definitions have been proposed in the literature (Cordier & Bergmann 2013;Holmes et al. 2012), we herein focus on the classical L 2 inner product in its continuous and discrete forms. We shall moreover assume that any data realization (and any spatial or temporal structure) is the discrete representation (sampling) of a continuous and square integrable function. Therefore, let the vector φ r [i] ∈ C ns×1 collect n s samples of its continuous counterpart φ r (x), over a spatial domain x ∈ Ω discretized with a Cartesian grid x i . An estimation of the average energy content of this function is provided by the Euclidean inner product where the spatial structure is also reshaped as a column vector, the over-line denotes complex conjugation, a, b = a † b is the Euclidean inner product of two vectors a, b, and † represents the Hermitian transpose. Similarly, assuming that the vector ψ r [k] ∈ C nt×1 collects n t samples of the function ψ r (t), defined in a time domain t ∈ [0, T ], the averaged energy content in the time domain reads: Equations (2.3) and (2.4) are numerical approximations of the continuous inner products using the right endpoint rule, and are valuable only if both the spatial and the temporal discretizations are Cartesian. More advanced quadrature methods or nonuniform discretization requires the definition of a weighted inner product (Volkwein 2013;Cordier & Bergmann 2013;Gordeyev & Thomas 2000). To simplify the following derivations, we here focus on Euclidean inner products and uniformly sampled data.
From (2.3)-(2.4), it is clear that the normalization ||φ r || 2 = ||ψ r || 2 = 1, at a discrete level, implies that the amplitudes σ r in (2.2) must be normalized by √ n s n t to produce a grid-independent estimation of the energy of each mode. Before discussing further the link between the energy contributions σ r and the energy of the entire dataset, it is useful to arrange the spatial and the temporal structures into matrices . Any decomposition of the form (2.2) becomes a matrix factorization of the form where Σ = diag[σ 1 , σ 2 , . . . , σ rk(D) ] is the diagonal matrix containing the energy contribution of each mode, and the superscript T denotes matrix transposition. Besides allowing for a compact notation, (2.5) allows for visualizing the operation as a projection of the dataset onto a spatial basis Φ and a temporal basis Ψ ; the first is a base for the columns of D; the second is a base for its rows. It is worth noticing that since Σ is diagonal and Φ and Ψ have normalized columns, the factorization in (2.5) is entirely defined by either Φ or Ψ . In this work, we focus on projections in the time (row) domain of the dataset, i.e., (2.5) is defined by the temporal structures Ψ . From (2.5), a matrix inversion and a column-wise normalization allow for computing the spatial structures and the amplitudes: where σ r = ||C r || 2 and C r is r th column of the matrix C = D Ψ T −1 . Particularly convenient is the case of an orthonormal temporal basis, for which Ψ † Ψ = Ψ T Ψ = I, with I the identity matrix. In this case, the inversion becomes Ψ T −1 = Ψ and the energy contributions, which become σ r = ||D ψ r || 2 = ||σ r φ r || 2 , can be computed independently for each mode.
Among the infinite possible bases Ψ , we focus on the three most popular ones: the POD, the DFT and the DMD. Their matrix form is discussed in the following subsections, using the subscripts P, F and D are used to distinguish the decompositions.
Before proceeding, it is worth commenting on two important aspects of any decompositions: 1) the treatment of the time (column-wise) average, 2) the link between the contributions σ r and the total energy content in D. Concerning the first, one can see from (2.2) and (2.5) that removing the time average is equivalent to taking the normalized constant vector ψ µ = 1/ √ n t ∈ R nt×1 as one of the temporal structures, such that where d µ [i] = D ψ µ is proportional to the average column of the dataset, and the energy content σ µ = ||d µ || is proportional to the averaged correlation level, to be normalized by √ n s n t to provide a grid independent estimation. TheD [i, k] is the zero-mean shifted dataset and the last summation is its modal decomposition (distinguished by an accent •). The zero-mean shifted dataset has rk(D) = n t − 1 and all its temporal structures -if an orthogonal decomposition is considered-must have zero mean. This is naturally the case for the DFT, but not for the POD or the DMD. Concerning the second point, the energy content of the dataset can be defined via different matrix norms. Using the most classical Frobenius norm yields where Tr indicates the trace of the matrices K = D † D ∈ R nt×nt and C = DD † ∈ R ns×ns . These are the finite dimensional estimators of the two-point temporal and spatial correlation tensors (Holmes et al. 2012;Cordier & Bergmann 2013) respectively, using the inner products in (2.3)-(2.4). These symmetric positive definite matrices share the same non zero eigenvalues λ r . Introducing (2.5) in (2.8) gives two ways of measuring the total energy from the contribution of each mode: It is therefore evident that it is not possible to compute the total energy of the data by solely using the energy contributions of the modes, unless both the spatial and the temporal structures are orthonormal (Φ † Φ = Ψ † Ψ = I). This case is described in 2.1.

Matrix form of the Proper Orthogonal Decomposition (POD)
The POD temporal structures Ψ P are the eigenvectors of the temporal correlation matrix K = D † D. If D is real, as in all the applications considered in this work, these are orthonormal and real vectors: λ Pr ψ Pr ψ T Pr . (2.10) Introducing this definition in (2.9) is particularly revealing: a comparison of (2.10) to the first factorization in (2.9) yields Λ P = Σ P Φ T P Φ P Σ P . Since this eigenvalue matrix is diagonal, one recovers Φ T P Φ P = I and Λ P = Σ 2 P . Therefore, also the spatial structures are real and orthonormal and thus eigenvectors of the spatial correlation matrix C = D D † = Φ P Λ P Φ T P . Therefore, the energy contribution of each mode becomes σ Pr = √ λ r , and the energy of the data-set can be recovered from the Σ P .
For the choice of temporal basis in (2.10), (2.5) yields the well-known Singular Value 8

M. A. Mendez et al
Decomposition (SVD), here denoted as D = Φ P Σ P Ψ T P . The optimality of the decomposition is thus guaranteed by the Eckart-Young-Mirsky theorem, which states that the rank r c < rk(D) approximation of a matrix D obtained from a truncated SVD minimizes the norm of the error (2.11) Therefore, each coefficients σ Pr represents also the L 2 error produced by a rank r − 1 approximation, and the error decay Err(r c ) → 0 as r c → rk(D) is the strongest possible.
The POD from Lumley (1967Lumley ( , 1970 computes the spatial structures Φ P from the eigendecomposition of C; the BOD from Aubry (1991) computes the temporal structures (referred to as BOD modes) Ψ P from the eigendecomposition of K. The POD from Sirovich (1987Sirovich ( , 1989Sirovich ( , 1991 computes both Φ P and Ψ P by first solving (2.10) (less demanding if n t n s ) and then using (2.6)(which reads Φ P = D Ψ P Σ −1 P ). 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 K. This implies: that is the sum over any row of K is equal to the corresponding eigenvalue σ 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 et al. 2008 for a derivation).
Finally, it is worth highlighting that the energy optimality is linked to the decomposition uniqueness: in 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 Φ P and Ψ P .

Matrix form of the Discrete Fourier Transform (DFT)
The DFT temporal structures Ψ F are defined a priori, regardless of the dataset at hand, as columns of the well known Fourier Matrix where w = exp 2πj/n t , with j = √ −1. Each temporal structure is thus a complex exponential ψ r [k] = exp 2πjf r (k −1)∆t = exp 2πj(r −1)(k −1)/n t with real frequency discretization f r = (r − 1)∆f , with ∆f = 1/(n t ∆t) = f s /n t the frequency resolution and f s the sampling frequency.
Besides being orthonormal (Ψ −1 F = Ψ † F ), this Vandermonde matrix is also symmetric (Ψ T F = Ψ F ). The resulting spatial structures Φ F = DΨ F Σ −1 F from (2.6) are complex and generally not orthogonal. Finally, since the DFT satisfies (2.7) by construction, the zero-mean shifting has no effects on the decomposition.

Matrix form of the Dynamic Mode Decomposition (DMD)
The DMD temporal structures Ψ D are computed under the assumption that a linear dynamical system can well approximate the dataset. It is thus assumed that a propagator matrix P ∈ R ns×ns maps each column of the data matrix D onto the following one via simple matrix multiplication (2.14) As each time step k involves the k − 1 power of P acting on the initial data d 1 , the evolution of such linear system depends on the eigendecomposition P = SΛS −1 of the propagator, since P k−1 = S Λ k−1 S −1 : (2.15) The columns of the eigenvector matrix S = [s 1 , s 2 , . . . , s ns ] 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 a 0 = S −1 d 1 is the projection of the initial data d 1 onto the eigenvectors of the propagator. To arrange (2.15) into the factorization form (2.5), one should first arrange the vector a 0 into a diagonal matrix A 0 = diag(a 0 ) and then build the temporal basis from the powers of the eigenvalues. The resulting factorization can be written as The temporal basis of the DMD is thus obtained by normalizing each column of Z independently, that is ψ Dr = Z r /||Z r ||, and the spatial structures can be computed using the projection in (2.6) with no need for computing the eigenvectors S. 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, containing n t − 1 realization: where D + 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 D 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 D + 1 = Ψ P Σ −1 P Φ T P , the size of this matrix is usually prohibitively large (P ∈ R ns×ns ) and thus the DMD decomposition based on the eigenvalues of P (referred to as exact DMD in Tu et al. 2014) 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 P with a much smaller matrix C ∈ R (nt−1)×(nt−1) acting on the right side (Rowley et al. 2009). This matrix is the Companion matrix, defined so that The ones in the sub-diagonals shift the columns of D 1 to obtain those of D 2 ; the last column contains the set of coefficients c = [c 1 , c 2 . . . , c nt−1 ] T that approximates the last temporal realization d nt as a linear combination of the previous ones. As an eigenvalue λ of C, with eigenvector v such that Cv = λv, is also an eigenvalue of P with eigenvector D v, this formulation of the DMD relies on the eigendecomposition of C and thus the calculation of coefficients c. This is a least square minimization problem, extremely sensitive to noise and the choice of the final dataset. Using a QR factorization D 1 = Q R, the solution can be written as c = R −1 Q T d nt . More robust variants of this method are the optimized DMD by Chen et al. (2012) or the sparsity promoting DMD by Jovanović et al. (2014) 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 C and projects the problem onto a space of much lower dimension (Schmid 2010), 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 S =Φ T P PΦ P ∈ R rc×rc , where the tilde denotes an approximationΦ P = [φ P1 , . . . φ Prc ]. Writing the full size propagator using the Moore-Penrose inverse, this projected propagator reads Observe that (2.19) is a similarity transform, and thus S and P share the same eigenvalues, only ifΦ P is a square n s × n s matrix (that is the full POD basis plus its n s − n t orthogonal complements are taken into account so thatΦ PΦ T P = I). That is never the case, and one heuristically hopes that at least the first r c eigenvalues of P are available in S, which implies thatΦ PΦ T P ≈ 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 ψ µ = 1/ √ n t as basis element only if λ = 1 is an eigenvalue of the propagator (C for the cDMD, S for the sDMD). While there is no simple relation to know when this is the case for S, one can see from the characteristic polynomial of C that this condition sets 1, c = 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 c = −1 regardless of the dataset, since the last realization must cancel the summation. In this case, as shown by Chen et al. (2012), the solution of the characteristic polynomials are the n t − 1 roots of unity, and the DMD reduces to the DFT.

The Multiscale Proper Orthogonal Decomposition (mPOD)
The temporal basis Ψ M for the Multiscale Proper Orthogonal Decomposition (mPOD) proposed in this work is computed by setting frequency constraints to the classical POD.
The algorithm consists of three major steps. First, the temporal correlation matrix K = D † D is split into the contributions of different scales via multi-resolution analysis (MRA). This can be performed via 2D Wavelet Transform or more generally via an appropriate filter bank that preserves the symmetry of the correlation matrix. Second, each of these contributions is diagonalized as in the classical POD, to obtain the orthonormal basis corresponding to each scale. Third, the bases of each scale are merged into the final temporal basis Ψ M , and the decomposition is completed via (2.6).
The theoretical foundation of the algorithm is described in §3.1, §3.2 and §3.3. In particular, §3.1 investigates the link between the frequency content in D, in K and in its eigenvectors Ψ P . Section 3.2 introduces the generalized MRA of K via filter banks, discussing its impact on the eigenvectors of different scales and the necessary conditions to keep these mutually orthogonal. {Section 3.3 discusses the link between POD, mPOD and DFT, while §3.4 presents the proposed mPOD algorithm.

Mode Spectra and Correlation Spectra
The multiresolution analysis (MRA) is introduced in §3.2 in the frequency domain, that is along the 2D Fourier transform of the correlation matrix. It is thus interesting to first analyze the standard POD in the frequency domain. At the scope, we denote as A the Discrete Fourier Transform (DFT) of a matrix A, 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 Ψ F , while the DFT of a row vector is a right multiplication by Ψ F .
The three Fourier pairs to be linked in this section are related to the time evolution of the data (D, row-wise), the temporal structures of the POD modes (Ψ P , column-wise) and the temporal correlation matrix (K, over both columns and rows). These are: 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: This matrix is obtained from the temporal correlation matrix K via a similarity transform (since Ψ F Ψ F = I), 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 D, Ψ P and K) are real, and thus their Fourier transforms are Hermitian symmetric. For a continuous signal x(t), having Fourier transform F{x(t)} = X(ω), the Hermitian symmetry sets X(ω) = X(−ω). For a discrete signal, this property can be introduced via a matrix permutation P π , which flips the spectra along the zero frequency (first column of Ψ F in (2.13)).
This permutation matrix can be obtained by applying the DFT operator twice: Therefore, for any discrete real vector a with Discrete Fourier Transformâ = Ψ F a, the Hermitian symmetry can be written asâ = P πâ . 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: Introducing the eigenvalue decomposition of K in (3.2) and using (3.4a) yields This shows that the eigenvectors of K 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. 1987;Gordeyev & Thomas 2000;Citriniti & George 2000;Gudmundsson & Colonius 2011;Sieber et al. 2017;Towne et al. 2018), (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 K in (3.1c), on the other hand, is not obtained via similarity transform and is not self-adjoint ( K = K † ). Hence, since it does not share the same eigenvalues of K and K 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 K F : (3.6)

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.
To illustrate the idea, let us consider the Fourier transform d p = d p Ψ F of the temporal evolution (row-wise) of the data for a given point spatial location i p , that is d p [k] = D[i p , k]. We seek to split the spectra d p into M scales, each retaining a portion with negligible overlapping as pictorially illustrated in Figure 1. The frequency bandwidths of these scales , are defined by a frequency  Figure 1 with a continuous line) is by construction such that H L1 + H H1 + · · · + H H M −1 = 1 and thus the entire spectra 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 . Moreover, since each of these filters is applied in the entire spatial domain (along each row of D), one should copy the transfer functions H Lm or H Hm (row vectors ∈ C 1×nt ) to obtain n s × n t matrices having all the rows equal to H Lm or H Hm . Let H m be the general transfer function (whether low-pass or band-pass), we denote with an apex the bi-dimensional extension of the 1D 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) The first multiplication is the DFT along the time domain ( D), the result of the Hadamard multiplication in the square bracket is the filtered spectra of the data ( D 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 K m = D † m D m from each scale contribution reads where the associative and commutative properties of the Hadamard product are used to move the transfer functions. Introducing the cross-spectral density matrix K F in (3.2) and writing the 2D transfer function as 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 K F m = K F H m shares the same eigenvalues of the correlation of filtered data K m = D † m D m . Second, introducing the eigenvalue decomposition of K m as in (3.5), shows that -up to a permutation-the eigenvectors of K F are the DFT of the POD modes of D m : having introduced (3.4a) in the last step. The impact of the filter on the POD modes is revealed by the diagonal entries of K F m and K F , which are summations of nonnegative real numbers. From (3.10), using (3.4a), these are: If the filter at a given scale m removes one of the entries along the diagonal of K F (H m [i, i] = 0), the associated frequency ψ Pm r [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 to the filtered data K m can be computed by a simple 2D filtering of the original correlation matrix, saving considerable computational time with respect to n s filters along the rows of D. Equation (3.9) shows that the MRA of the dataset can be carried out via 2D MRA of the temporal correlation matrix.
The four possible structures for the 2D transfer functions H m that arise from a combination of the 1D kernels in the time domain are illustrated graphically in Figure 2. For a given scale, composed of a low frequency and a high frequency D m = D Lm + D Hm , these correspond to the contributions of four temporal correlation matrices: The 'pure' terms D † Lm D Lm and D † Hm D Hm 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 K m .
It is worth noticing that the filter bank architecture described in this section is a generalized version of the MRA via Discrete (Dyadic) Wavelet Transform (DWT) proposed in previous works (Mendez et al. 2018a,b). In the 2D 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 D † Lm D Lm and D † Hm D Hm . 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 D † Lm D Hm and D † Hm D Lm , 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.

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 analyzed while studying how the choice of the frequency splitting vector 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 2D spectra of the correlation matrix K for four choices of F V . The case (a) corresponds to a single scale mPOD (empty Figure 3: Partition of the correlation matrix K for different choices of the splitting vector F V in the MRA step of the mPOD. The same color notation of Fig.2 is used for the band-pass/band-stop regions. Case a) consists of a single scale, thus the mPOD becomes a POD. Case b) and c) consists of two and three scales respectively. Case d) shows the limit of n t /2 scales, in which the mPOD becomes a DFT.
that is the standard POD. All the eigenvectors ψ Pr 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 δ(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 F V = f c 1 . Each scale has its own correlation matrix (K L and K H1 ) 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 R nt space and thus preserves 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 F V = [f c 1 , f c 2 ], while the case (d) corresponds to the limit F V → {r f s /(∆t n t )} nt−1 r=1 , with r ∈ [1, . . . , 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 1978;Daubechies 1992;Folland 2009).

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 K and its Fourier transform K from (3.1c). From the analysis of the frequency content in the correlation spectra K, the third step consists in defining the frequency splitting vector F V and in constructing the set of associated transfer functions. Input: Set of n t snapshots of a dataset D(x i , t k ) on a Cartesian grid x i ∈ R nx×ny 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 centering the band-pass windows of each scale on the dominant peaks in the diagonal of K and splitting the band-pass regions accordingly. In this step, only the 'pure' (approximation and diagonal detail) terms of the 2D transfer functions are constructed, that is H L1 = (H L1 ) † H L1 and H Hm = (H Hm ) † H Hm . Disregarding the mixed terms (horizontal and vertical details) in each scale, the correlation matrix is approximated as: (3.13) Each of these contributions is a symmetric, real and positive definite matrix, equipped with its orthonormal eigenspace, computed in the fourth step (3.14) If the frequency-overlapping 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. (2016b), 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 analyzed in this work.
Because of the limited size of the filter impulse response (see Annex 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 [diag(Σ L1 ), diag(Σ H1 ), . . . , diag(Σ H M −1 )] in decreasing order and storing the required permutation matrix P Σ . This matrix is then POD DFT sDMD cDMD mPOD (M=3) mPOD (M=10) mPOD (M=50) Test Case 1 ( §4) 1.  used to re-arrange the eigenvector matrices Ψ L , . . . , Ψ Hm into one single matrix (temporal basis) regardless of their scale of origin 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. (2016), in which the filtering procedure can eventually result in a nonorthogonal temporal basis, or the mrDMD by Kutz et al. (2016b), 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 QR factorization Ψ 0 M = Ψ M R, so that Ψ M = Ψ 0 M R −1 ∈ R nt×nt . This compensates for the losses of orthogonality and the upper triangular matrix R 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.
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- §5- §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 ∼ 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.

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. (4.1c) 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/||S r T T r || 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 D[i, k] is a matrix of size 65536 × 256 of rank rk(D) = 3.

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 S in (2.19), is severely challenged by the rank deficiency of this dataset. Since rk(D) = 3, one has S ∈ R 3×3 , that is a set of eigenvalues incapable of reproducing the data: each of the dominant frequencies 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).
f 1 , f 2 , f 3 requires two complex conjugate modes and thus at least a 6×6 propagator. Low dimensional propagators 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 where all the modal contributions are first sorted in descending order σ a > σ b , ∀a > b.
For the POD, (4.2) reduces to Err[r c ]/||D|| 2 = σ Prc+1 /σ P1 since the L 2 norm of a matrix is equivalent to its largest singular value (see eq. 2.11). This is not the case for other decompositions, for which the spatial basis is not orthogonal. The convergence results of this test case are shown in figure 5. As expected, the POD produces no reconstruction error when r c 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, F V = [0.4, 10, 20] and leads to a convergence much closer to that of the POD than the harmonic decompositions. To analyze the detection capabilities of these decompositions, several representative modes for each are now considered. First, the energy contributions of the DFT (diag(Σ F )) and the cDMD (diag(Σ 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 |λ r | < 1) and exponentially decay, or outside (if |λ r | > 1) and exponentially grow.
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 phaseshifted 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 ≈ 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 6 with a square marker.
This localization limit, stemming from the Heisenberg uncertainty principle (see Daubechies (1992)) is the main 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., φ F3 ≈ φ F9 , φ F5 ≈ −φ F11 ) associated, in the time domain, to harmonic corrections to represent finite duration events.
The DMD modes, having a complex frequency, give a valid alternative when describing  Figure 9: Results from the POD of the synthetic test case. The first row shows the normalized temporal correlation matrix K = D † D/n t (left) and corresponding DFT K = Ψ F KΨ F /n t (right). The second row shows the three non-zero POD modes, including their spatial structures φ Pr and temporal structures ψ Pr . 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 , 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.
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 K and its DFT K, 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 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 center in both the diagonal and the antidiagonal. The spectral leakage produces cross-like patterns centered 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. (2016). 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 Σ P = 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 redistribution. Yet, the singular values are sufficiently close to the limit σ 1 = σ 2 = σ 3 to approach the condition of wholly undetermined POD.
The eigenvectors of K, computed using the Matlab commands 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, Fig. 10 shows that the temporal evolutions in Fig. 4 are almost eigenvectors of K: the matrix multiplication KT k leads to vectors that are close to simple multiples of the temporal evolution T k . The more K is close to being defective, the larger the number of vectors that are almost its eigenvectors and thus POD structures.
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 F V = [0.4, 10, 20]. The DFT of each scale is shown in the first row of Fig. 11; 24 M. A. Mendez et al the dominant eigenvector for each scale is shown in the second row of Fig. 11. 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.

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 2D Navier-Stokes (NS) equations in the vorticity-stream function formulation. This formulation is common in the study of large-scale geophysical systems (Flór et al. 1995;Trieling et al. 1998;Mariotti et al. 1994), and consists of a set of two equations. The first is the transport equation of the vorticity ω = ∇ × u, where the velocity field u = (u, v) is described by a stream function ζ, defined so that ∂ y ζ = u and ∂ x ζ = −v: where J ω, ζ = ∂ y ζ ∂ x ω − ∂ x ζ ∂ y ω is the Jacobian determinant of the solution vector (ω, ζ); Re is the Reynolds number weighting the importance of advection and diffusion; S ω is the spatially and temporally varying source or sink of vorticity. The second equation is the Poisson problem linking vorticity and stream function: Sources of vorticity S ω 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 1993;Kutz 2013) The source term is treated explicitly, and consists of coherent (S Cω ), pseudo-random (S Rω ) and fully random (S N ω ) contributions: The four spatial S C φ l (x i ) and temporal S C ψ l (t k ) structures of the coherent source S Cω 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 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 ψ 1 . This temporal structure is composed of a set of Gaussian pulses A i exp((t k − t i ) 2 /σ 2 t1 ) of amplitude varying between ±[1, 1.5] and equal widths of σ 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 ψ 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: with frequencies f 2,3 = [2, 5], equal width σ t2,3 = 4 but different phase delay for both the harmonic (α 2,3 = [14, 6]) and the modulation terms (β 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 center of the domain, and a temporal structure consisting of a single Gaussian impulse located a t k = 5.2 and width σ t4 = 0.05. The pseudo-random source term S Rω in (5.3) produces a periodic injection of random vortices. The spatial structure S Rφ consist of a set of 300 Gaussian vortices, randomly  ]. An exemplary set of vortices representing one of these structures is shown in Figure 13 (left). The temporal evolution S Rψ , shown in Figure 13 (right), 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 Nω in (5.3) is introduced in the post-processing step, using a random number generator varying between 0 and 1 to generate the spatiotemporal noise N (x i , t k ).
The coefficients A C l (l ∈ [1, 4]), A R and A 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 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.
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 ≈ 5.15, the injection of the source term S Cψ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 center 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ω , these two dipoles do not travel the same distance and are broken and wholly dissipated by t ≈ 6.5.
At t k > 6, the third coherent source begins its pulsation. At t k ≈ 7 this generates four wakes which shed at a constant frequency. The overall vorticity in the domain increases and the large-scale clock-wise rotation involves the four wakes in a closed circular pattern. Before this large-scale rotation vanishes, at t k ≈ 14, the second coherent source starts to pulse at a lower frequency, producing other four trains of vortices, initially also rotating counterclockwise. 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.

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 a noisy data with nonlinear evolution. As shown in §2.3, the amplitude coefficients σ Dr 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 is these vanishes 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.
The problem is particularly evident for the sDMD with the reduced propagator S expressed in (2.19). Figure 17 shows the eigenvalues distribution in the complex plane for reduced propagators S of different dimensions r c × r c , obtained by projecting the full propagator onto the first r c POD modes. For r c 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.
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 ||d nt − D 1 c|| 2 .
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 (first row) shows the corresponding spatial structures, which are well related to the Figure 17: Eigenvalue distribution for reduced propagators S ∈ C rc×rc in (2.19) obtained by projecting the full propagator P 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 (left), and cDMD (right) 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.
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 the second row of Figure 19. The amplitude of these modes is labeled with a red marker in the complex plane of Figure 18, and is much lower than the dominant ones, associated to 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  08, 1.97, 4.99], that is the ones closer to those dominating the DFT spectra. 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 K and its DFT spectra K, together with the spatial and the temporal structures of six representative POD modes.
The temporal correlation matrix K 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.
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. 2017) 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 (σ P2 = 4272 and σ P2 = 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 K,  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 eigenvector of K L ; the thirteenth is eigenvector of K H3 ), but no significant overlapping between this event and the other harmonics is produced.
6. Example III: PIV of an Impinging Jet Flow

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 2006;Specht 2014;Buchlin 2011) and the POD of PIV data has been successfully used by many authors to identify the resulting coherent structures (Kim et al. 2007;Charmiyan et al. 2017;Pieris et al. 2017;Hammad & Milanovic 2009).
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 et al. 1978).
The experimental set up is shown in Figure 23. The rectangular outlet section of the nozzle has an opening of H = 4mm 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 aluminum plate placed at a distance of Z = 30 ± 1%mm. The experiment is carried out with P N = 30 ± 1% P a leading to an average jet flow velocity of U 0 = 6.5 ± 2%m/s, that is a Reynolds number of Re ≈ U 0 H/ν = 1733 considering a kinematic viscosity of ν = 15cSt.
The equipped with a zoom objective Tamron CZ-735 70 − 300mm to provide an image scaling factor of 24.1 ± 0.1pixel/mm. The investigated Region of Interest (ROI) is a rectangular area of 20mm × 38mm at a distance of 10mm from the nozzle outlet (cf. Figure 23). The video sequence is made of n t = 2728 images pairs of 1000 × 740 pixels, acquired in frame-straddling mode at f s = 2kHz. The separation time between the frames is set to ∆t p = 70µs to produce a maximum displacement, in the jet flow, of about 8 pixels. The PIV evaluation is carried out with the Matlab package PIVlab (Thielicke & Stamhuis 2014), using standard iterative multi-step interrogation (four passes from 128 × 128 to 8 × 8 windows) with spline window deformation, 2D Gaussian sub-pixel interpolation, and vector validation via median test. The resulting vector field consists of n s = 60 × 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 24b 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 1967) 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. In this location, most of the energy is centered around a peak at St ≈ 0.3, in agreement with other studies on the free plane jet in transitional regime (Thomas & Goldschmidt 1986;Suresh et al. 2008). 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 1982).
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. (2017) 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 1985;Hadžiabdic & Hanjalic 2008).   Figures 6 and 18, but for the third test case presented in this section. 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 F V 1 (m) = [0.1, 0.2, 0.5]; the second consist of M = 50 with F V 2 (m) = 0.01 + m∆f with equal bandwidth ∆f = 0.02. As discussed in §3.2, the finer the frequency splitting, the more the mPOD approaches a DFT.

Modal Analysis
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 ||d nt − D 1 c|| 2 /||d nt || 2 ≈ 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.
The spatial structures of two representative DFT and DMD modes are shown in Figure  27. For the DFT, these structures are associated to 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 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 Φ D27 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 Φ D29 captures the K-H structures, vortical structures of about 2mm 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 ψ Pr associated to each spatial structure, rather than their temporal evolution ψ Pr . All the modes in the range r = 4 to r = 11 are related to traveling structures and arise therefore in pairs that have π/2 phase delay in space and time. Only the four representative modes are shown.
The temporal correlation matrix K displays a regular pattern composed of a broad range of frequencies and K displays a band-like structure. Yet, a dominant frequency around 12Hz (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. Being the data statistically stationary, the constant vector ψ µ = 1/ √ 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 K. 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 φ P2 in Figure 28 is mostly related to the large-scale pulsation in the wall jet (labeled with I) and its periodic separation, as expected from the spectra of the correlation matrix K. However, this spatial structure also captures, to a minor extent, coherent structures in the free jet region (labeled 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 spectra ψ P2 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 φ P3 , which is mostly linked to the traveling 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 worse case in the mode φ P9 . In addition to the two mechanisms previously described, this mode also captures the K-H structures produced further upstream, where they originate at St ≈ 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 F V 1 (m) = [0.1, 0.2, 0.5], producing the L 2 convergence labeled 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.

Conclusions
The state of the art on data-driven decomposition has been extensively discussed, with particular emphasis on the need for bridging the gap between the energy optimality of the Proper Orthogonal Decomposition (POD) and the spectral purity of the Dynamic Mode Decomposition (DMD). A generalized algebraic framework is introduced, and all the investigated decompositions are presented as matrix factorizations. This framework is then used to derive a novel decomposition, referred to as Multiscale Proper Orthogonal Decomposition (mPOD).
The derivation of the mPOD starts with an investigation on the link between the frequency content (the DFT) of the correlation matrix and that of its eigenvectors (the POD modes). This link is used to study the impact of a time filter in the dataset on the frequency of the correlation matrix and thus the related POD modes. Specifically, it is shown how a filter bank can be used to perform multi-resolution analysis (MRA) of the correlation matrix, splitting it into the contribution of different scales retaining nonoverlapping portions of the frequency spectra of the data. These scales are equipped with  Figure 28: Same as 9 and 20, but for the third test case discussed in this section and replacing the temporal structures ψ P with their normalized frequency content | ψ P |/max{ ψ P }.
their optimal eigenbasis which can be kept mutually orthogonal and finally assembled into a single mPOD basis. The novel decomposition is tested on three selected test cases and its performances compared to those of a Discrete Fourier Transform (DFT), a standard POD and the DMD in its two standard formulations: the Arnoldi-based approach (cDMD) and the SVD-based approach (sDMD). The first of the proposed test case is a synthetic dataset with three predefined modes. These have equal energy, identical but shifted spatial structure, and temporal evolution with different frequency and time localization. None of the popular decompositions presented proved capable of distinguishing the introduced modes. Harmonic decompositions such as DFT and DMD allows for proper frequency identification, but yields poor convergence, high redundancy in the spatial structures, and a complete lack of temporal localization due to the infinite duration of the time basis. On the contrary, the POD allows for the best convergence, with only three modes required to describe correctly (i.e., with zero L 2 errors) the dataset. However, it produces a 'spectral mixing' between different modes. This mixing results in modes that are linear combinations of the expected ones since their similar energy content results in almost repeated eigenvalues of the correlation matrix and thus almost undefined (non-unique) POD.
The proposed mPOD proved capable of distinguishing the three introduced modes, with no spectral mixing, excellent time localization, and almost optimal convergence.
The second investigated test case is a 2D numerical simulation of the Navier-Stokes equation, with prescribed coherent, random and pseudo-random sources. The complexity of the data, further polluted with white noise, strongly accentuate all the limits high-lighted by the previous test case, with the DMD L 2 relative error of the order of 200% within the first modes and up to three different mechanisms sharing the same POD mode. Neither the DFT/DMD nor the POD are capable of identifying and localizing in space and time the different coherent source terms, while the mPOD succeeds at the scope with a negligible loss of convergence over the POD.
Finally, the third case is the time-resolved flow field of an impinging jet experimentally investigated via TR-PIV. Although the data is stationary and statistically well converged, the spectral mixing of the POD is observed while the harmonic constraints of the DFT and DMD yield much poorer convergence and noisy spatial structures. The mPOD naturally handles both problems.
To conclude, the discussed test cases prove that a minor spectral constraint in the diagonalization of the correlation matrix allows overcoming the frequency limits of the POD, and show that the proposed mPOD provides an excellent compromise between energy optimality and spectral purity. While this work has focused on data-driven decompositions for fluid flows, the possibility of constraining the POD optimality within specific scales of interest gives to the mPOD more flexibility than POD or DMD/DFT for a large variety of applications. In data filtering/compression, for example, the mPOD is more robust against random noise (equally spread over all the POD/DFT spectra) and allows to better focus on particular features (e.g. allowing for different levels of data compression/filtering for different scales). This could find useful applications, for example, in the background/foreground separation problem in video analysis, where the slow motion of large objects is modeled by few dominant POD modes (or slow DMD ones) while the dynamics of smaller objects is captured by the remaining modes (see Oliver et al. 2000;Sobral et al. 2015;Grosek & Kutz 2014;Mendez et al. 2017). The mPOD provides an additional degree of freedom in this separation by also including considerations on the frequencies involved. In the data-driven identification of coherent patterns, the mPOD allows overcoming the POD non-uniqueness problem in case of phenomena having similar energy content (e.g., Mendez et al. 2018a).
More generally, by providing bases that are spectrally cleaner than POD bases but energetically more relevant than DFT/DMD ones, the mPOD can be useful in any 'projection-based' application (from pattern recognition to low order modeling of dynamical systems) in which the constraints of energy optimality or spectral purity are both unnecessarily extreme.
transform of the associated transfer function. The impulse response (or kernel) is the output of the filter when acting on a Dirac impulse δ[k], and its design can be carried out using standard methods presented in classical textbooks on signal analysis (Oppenheim & Schafer 2009;Hayes 2011). Among these, the current mPOD algorithm uses the windowing method, which consists in smoothing an ideal transfer function The linear phase response is provided by taking the window kernel to be symmetric about its midpoint, that is w[n] = w[N − n]. The choice of the window, and most importantly its size N , controls how much the real transfer function differs from the ideal one, in terms of pass-band/stop-band deviation, and transition band (Oppenheim & Schafer 2009;Hayes 2011). Once the window is chosen, the low-pass transfer function corresponding to the impulse response in (A 3) is: with ω = 2πf n /f s the normalized frequency (radiant/samples) and w = Ψ F w the DFT of the window. Observe that (A 4) has an analytical representation and thus the transfer function of appropriate size H Lm ∈ C 1×nt , here constructed as a row vector for later convenience, can be computed from (A 4) on an arbitrary frequency discretization. In this work, the default window chosen is the Hanning window, which gives a stop-band attenuation of −44dB (decibels) and transition width of 3.1/N . Concerning the filter order, one should observe that a sharp frequency cut-off becomes more and more difficult when lower frequencies are considered. This parameter is therefore left as a function of the frequency bandwidth considered, as it is common in the filter bank formulation of the Continuous Wavelet Transform (Mallat 2009;Strang & Nguyen 1996). For the test case in §6, for example, the filter orders N m is A zero-phase implementation of these filters, using for example the Matlab function filtfilt, is based on backward and forward filtering, in order to cancel the phase delay (Gustafsson 1996;Matlab 2017). This operation doubles the filter order, squares the modulus of its transfer function but renders it non-causal.
Finally, it is worth noticing that in previous formulations of the mPOD (see Mendez et al. 2018a,b), the filter bank architecture is presented in terms of Discrete (dyadic) Wavelet Transform (DWT). This formulation replaces the concept of impulse response with that of wavelet (for high-pass) and scale (for low-pass) functions, constructed from dilatation and shifting of reference functions. The link between wavelet theory and filter banks is extensively discussed by Mallat (2009) and Strang & Nguyen (1996).
For the purposes of this work, it suffices to notice that the MRA of K can by carried out via Wavelet toolboxes for image compression/denoising (Gonzalez & Woods 2007;Mallat 2009), based on the pyramid algorithm proposed by Mallat (1989). In this case, however, both the scaling and the shift of the reference wavelet/scale functions are fixed and the frequency splitting vector F V becomes: F V = [2 −M , . . . 2 −3 , 2 −2 , 2 −1 ] f s/2 with f s the sampling frequency and M = log 2 (n t ) the number of scales.