J. Fluid Mech. (2016), vol. 792, pp. 798–828.
c
Cambridge University Press 2016
doi:10.1017/jfm.2016.103
798
Spectral proper orthogonal decomposition
Moritz Sieber
1,
†, C. Oliver Paschereit
1
and Kilian Oberleithner
1
1
Institut für Strömungsmechanik und Technische Akustik, Hermann-Föttinger-Institut,
Technische Universität Berlin, Müller-Breslau Str. 8, D-10623 Berlin, Germany
(Received 13 August 2015; revised 10 November 2015; accepted 3 February 2016;
first published online 4 March 2016)
The identification of coherent structures from experimental or numerical data is an
essential task when conducting research in fluid dynamics. This typically involves the
construction of an empirical mode base that appropriately captures the dominant flow
structures. The most prominent candidates are the energy-ranked proper orthogonal
decomposition (POD) and the frequency-ranked Fourier decomposition and dynamic
mode decomposition (DMD). However, these methods are not suitable when the
relevant coherent structures occur at low energies or at multiple frequencies, which
is often the case. To overcome the deficit of these ‘rigid’ approaches, we propose
a new method termed spectral proper orthogonal decomposition (SPOD). It is based
on classical POD and it can be applied to spatially and temporally resolved data.
The new method involves an additional temporal constraint that enables a clear
separation of phenomena that occur at multiple frequencies and energies. SPOD
allows for a continuous shifting from the energetically optimal POD to the spectrally
pure Fourier decomposition by changing a single parameter. In this article, SPOD is
motivated from phenomenological considerations of the POD autocorrelation matrix
and justified from dynamical systems theory. The new method is further applied to
three sets of PIV measurements of flows from very different engineering problems.
We consider the flow of a swirl-stabilized combustor, the wake of an airfoil with a
Gurney flap and the flow field of the sweeping jet behind a fluidic oscillator. For
these examples, the commonly used methods fail to assign the relevant coherent
structures to single modes. The SPOD, however, achieves a proper separation of
spatially and temporally coherent structures, which are either hidden in stochastic
turbulent fluctuations or spread over a wide frequency range. The SPOD requires
only one additional parameter, which can be estimated from the basic time scales of
the flow. In spite of all these benefits, the algorithmic complexity and computational
cost of the SPOD are only marginally greater than those of the snapshot POD.
Key words: computational methods, low-dimensional models, turbulent flows
1. Introduction and motivation
1.1. Contemporary methods for data reduction
Today’s high-fidelity computational fluid dynamics (CFD) and high-end experimental
data acquisition systems tend to produce vast amounts of data that are getting harder
† Email address for correspondence: moritz.sieber@fd.tu-berlin.de
Downloaded from https://www.cambridge.org/core. Universitaetsbibliothek, on 26 Oct 2017 at 07:48:17, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2016.103
Spectral proper orthogonal decomposition 799
to interpret and overview. Methods to analyse such data are numerous and are always
developing to stay in line with acquisition and computation systems. The most
challenging data stem from turbulent flows that feature a huge range of temporal and
spatial scales. A key challenge in turbulent flow data mining is the distinction of
deterministic coherent motion from purely stochastic motion. Numerous methods exist
that exploit the periodicity or energetic dominance of these coherent structures. These
methods range from classic Fourier decomposition to dynamic mode decomposition
(DMD) and proper orthogonal decomposition (POD). The most prominent among
them are briefly introduced in the following.
POD has been widely used since its introduction by Lumley (1970) and Sirovich
(1987). It was applied in nearly every fluid dynamic field. Beyond fluid dynamics, this
method is also known as singular value decomposition, principal component analysis
or Karhunen–Loève expansion (Berkooz, Holmes & Lumley 1993). The basic idea
behind this method is to construct an optimal basis that represents most of the data
variance with as few basis functions as possible. In the context of POD the variance
is turbulent kinetic energy. Therefore, the POD searches for the most energetic modes
whereby coherent structures with high energy content are likely to be represented by
POD basis functions (Holmes et al. 2012).
Another classical approach is the linear stochastic estimation introduced by Adrian
& Moin (1988), where the readings of different sensors are related via a linear
mapping. This is closely related to the extended POD (Boree 2003), also described
in a unified framework (observable inferred decomposition) by Schlegel et al. (2012).
In recent extensions of linear stochastic estimation, the use of time delays between
the different sensors and also the use of one sensor at multiple time instances are
pursued to separate periodic coherent structures from turbulent fluctuations (Durgesh
& Naughton 2010; Lasagna, Orazi & Iuso 2013). This approach was also used
to improve the determination of harmonic POD modes from few pressure sensors
(Hosseini, Martinuzzi & Noack 2015). These utilizations of data from various time
instances are also related to the temporal constraint used for the POD extension
proposed in this article.
Targeting the temporal periodicity of the coherent structures, spectral methods such
as discrete Fourier transform (DFT) and the recently introduced DMD (Rowley et al.
2009; Schmid 2010) come into play. These methods commonly span the mode space
according to fixed frequencies, which enables the identification of coherent structures
within small spectral bandwidths. In contrast to DFT, DMD also distinguishes modes
with respect to their linear amplification. The recently introduced extended DMD
(Williams, Kevrekidis & Rowley 2015) tries to overcome the limitations encountered
by the (linear) DMD approach when trying to decompose data from nonlinear
systems. The idea is to use nonlinear functions that create observables of the data,
which are exactly described by a linear system. This approach translocates the
problem towards the identification of these nonlinear functions, which can be solved
using the ‘kernel trick’ that is common in machine learning. This paper presents
an alternative approach, which extends POD to account for temporal dynamics in
addition to energetic optimality.
1.2. Why yet another method?
After this short and incomplete review of data processing methods, one may ask
if there is need for another method. The answer is probably no, so we take the
most used method (POD) and bring it up to date for present research issues. The
Downloaded from https://www.cambridge.org/core. Universitaetsbibliothek, on 26 Oct 2017 at 07:48:17, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2016.103
800 M. Sieber, C. O. Paschereit and K. Oberleithner
approach pursued here includes a simple yet effective extension to the classical POD,
which leads to a more general method comprising POD and also DFT. This approach
unifies existing methods, but also offers possibilities beyond these. From the authors’
experience, the currently available methods often fail when applied to challenging
flow data. These stem from flows with weak coherent structures where the recorded
data have low signal-to-noise ratios, from flows with intermittent dynamics, or from
flows featuring multi-modal interactions leading to frequency modulations, to name
a few. In such cases, much effort is required to optimize the data processing until
satisfactory results are obtained. The usual escape route is to focus on a certain spatial
region or to apply suitable filters to pick out a certain wavelength or frequency range.
This involves trial and error or requires prior knowledge of the investigated flow.
There is also the danger of cutting off a substantial portion of the data, leading to
false interpretations. These procedures can be collected under the heading ‘identifying
symmetries’ as done by Holmes et al. (2012). The drawback of this approach is that
the investigated flow must feature some symmetries and they must be known a priori.
A recent application shows the huge variety of spatial and temporal filtering together
with POD to separate different phenomena into different modes (Bourgeois, Noack &
Martinuzzi 2013), exemplifying the complexity of this approach.
The usage of spectral methods for highly turbulent flows is even more challenging
than POD. The variable frequency of single coherent structures and intermittent
occurrence of different structures with the same frequency hinders a proper
decomposition. In terms of the DFT, averaging of spectra from multiple measurements
or sensors is essential to obtain reliable results. Analogously, for DMD, averaging
over several events is an option to reject noise (Tu et al. 2014). Nonetheless, the
results obtained with DFT and DMD suffer from limiting the temporal dynamics to
single frequencies. Turbulent flows hardly ever feature discrete frequencies and it is
not always valuable to restrict a single mode (flow phenomenon) to a single frequency.
Coherent structures that feature significant phase jitter or frequency modulation are
represented by many modes at similar frequencies. In contrast, the POD puts no
temporal constraint on the modes. This can result in modes that represent flow
phenomena occurring at largely different temporal scales. Thus, it is often hard to
interpret these modes and draw meaningful conclusions from the temporal dynamics.
From our point of view, there is a big gap between the energetically optimal
decomposition of POD and the spectrally clean decomposition of DFT or DMD.
This gap will be bridged with the spectral proper orthogonal decomposition (SPOD)
introduced in this article. This new method not only places itself somewhere in
between these two extrema, but it allows for a continuous shifting from one to the
other. The main idea is to apply a filter operation to the POD correlation matrix,
which will force the POD towards clear temporal dynamics. Depending on the filter
strength we continuously sweep from classic POD to DFT.
The remainder of this article is organized as follows: the proposed method is
described in detail in § 2. The reader is guided from snapshot POD via an in-depth
interpretation of the correlation matrix towards the general description of the SPOD.
In addition, a method is explained to identify coupled mode pairs describing a single
coherent structure, which becomes handy when working with SPOD. In § 3, the
new method is demonstrated on three different experimental data sets. The results
are compared against POD and DFT to point out the benefits of SPOD. In § 4 the
capabilities of SPOD are summarized, based on the findings from the application to
experimental data.
Downloaded from https://www.cambridge.org/core. Universitaetsbibliothek, on 26 Oct 2017 at 07:48:17, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2016.103
Spectral proper orthogonal decomposition 801
2. Description and interpretation of the proposed method
2.1. Classical snapshot POD
To introduce the method and the nomenclature, the snapshot POD approach is
described first. We start off with a decomposition of a data set into spatial modes
and temporal coefficients:
u(x, t) =u(x) +u
0
(x, t) =u(x) +
N
X
i=1
a
i
(t)Φ
i
(x). (2.1)
Note that only the fluctuating part u
0
(x, t) is decomposed. It is split into a sum
of spatial modes Φ
i
and mode coefficients a
i
. A set of M spatial points recorded
simultaneously over N time steps is considered. To calculate the POD, the correlation
matrix of this data set is needed. For data obtained from particle image velocimetry
(PIV) or CFD, the number of spatial points is usually larger than the number of
snapshots. The correlation matrix is then calculated between individual snapshots
(temporal correlation). The alternative approach (spatial correlation) that applies to
M N is detailed in appendix A. The correlation between two snapshots is calculated
from an appropriate inner product h, i, usually defined as the L
2
inner product
hu(x), v(x)i=
Z
V
u(x)v(x) dV, (2.2)
where V specifies the spatial region or volume over which the correlation is integrated.
The elements of the correlation matrix R are given by
R
i,j
=
1
N
hu
0
(x, t
i
), u
0
(x, t
j
)i. (2.3)
Matrix R is of size N ×N.
The temporal coefficients a
i
=[a
i
(t
1
), . . . , a
i
(t
N
)]
T
and mode energies λ
i
are obtained
from the eigenvectors and eigenvalues of the correlation matrix:
Ra
i
=λ
i
a
i
; λ
1
> λ
2
> ··· > λ
N
> 0. (2.4)
The subscript i refers to single eigenvalues, which are sorted in descending order.
Since the a
i
are the eigenvectors of the real symmetric positive-definite matrix R, they
are orthogonal. Moreover, they are scaled with the energy of the single modes such
that
1
N
(a
i
, a
j
) =λ
i
δ
ij
, (2.5)
where ( , ) denotes the scalar product. The spatial modes are obtained from the
projection of the snapshots onto the temporal coefficients:
Φ
i
(x) =
1
Nλ
i
N
X
j=1
a
i
(t
j
)u
0
(x, t
j
). (2.6)
These modes are orthonormal by construction, i.e.
hΦ
i
, Φ
j
i=δ
ij
. (2.7)
The formulation so far is perfectly in line with classical snapshot POD, which can also
be computed by a singular value decomposition. However, since the SPOD requires a
manipulation of the correlation matrix we retain the classical form.
Downloaded from https://www.cambridge.org/core. Universitaetsbibliothek, on 26 Oct 2017 at 07:48:17, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2016.103
802 M. Sieber, C. O. Paschereit and K. Oberleithner
50
50
100
100
150
150
200
200
100 2000
0.5
0
–0.5
1.0
i
j
(a)(b)
FIGURE 1. (Colour online) (a) Pseudo-colour plot of the correlation matrix elements R
i,j
and (b) the corresponding correlation coefficient
b
R. The displayed data are selected from
PIV measurements of a forced turbulent jet.
2.2. Properties of the correlation matrix
The SPOD described in this article is essentially a filter applied to the correlation
matrix R. To offer a better understanding of this approach, the structure of the
correlation matrix R is inspected first.
Figure 1(a) shows the structure of the correlation matrix for the data set of a
forced turbulent jet. The data were acquired with PIV inside a 2-D plane aligned
with the jet axis. The considered flow shows strong vortex shedding at the forcing
frequency (the acquisition frequency is 25 times the forcing frequency). The presence
of these periodic patterns in the flow, and their convection within the observed flow
field, leads to a diagonal, wave-like structure of the matrix. This is closely related to
the periodicity of the autocorrelation coefficient. In fact, if the individual elements of
the correlation matrix R are summed up along the diagonals, we obtain the spatially
averaged autocorrelation coefficient
b
R(τ ) =
Z
T
τ
hu
0
(x, t), u
0
(x, t −τ )idt
Z
T
0
hu
0
(x, t), u
0
(x, t)idt
, (2.8)
where the upper bound T is the length of the measured sequence. This is depicted in
figure 1(b), showing the same periodicity as the correlation matrix. The autocorrelation
coefficient itself represents the spectral content of different time scales and wavelengths
and it is directly related to the power spectral density of the underlying data. However,
it contains no information on the phase of individual frequencies, due to the reference
of the signal to itself. This is why the elements along the diagonals of R look so
similar, as they represent only relative changes with respect to the time step on the
main diagonal. Thus, increased similarity along the diagonals of R is equivalent to an
increased similarity of the dynamics of the underlying signal. This property will be
discussed in more depth in § 2.4. The obvious consequence from these findings is: if
we want to obtain smooth dynamics from the POD, we have to enforce the diagonal
similarity of the correlation matrix. This is where we step into spectral POD.
Downloaded from https://www.cambridge.org/core. Universitaetsbibliothek, on 26 Oct 2017 at 07:48:17, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2016.103
