Guide to Spectral Proper Orthogonal Decomposition

A. Theory

An aim of this paper is to make SPOD accessible to practitioners, but a thorough understanding and disambiguation between variants of the approach requires some theoretical background. We refer the reader to several sources [3,6,9,10] and briefly summarize the results here.

POD seeks a deterministic function, $\mathit{\varphi }(z)$, or a set of such functions that, on average, best approximate a zero-mean stochastic process, $\{\mathit{q}(z;\xi )\}$, based on a finite ensemble of samples. Here $\mathit{u}$ represents the random variable, and $z\in \mathrm{\Omega }$ represents the independent variables and their domain, which in general can include both spatial dimensions and time. $\xi$ parameterizes the probability space, which is, in turn, equipped with an expectation operator, $E\{\cdot \}$, so that $\overline{\mathit{q}(z)}=E\{\mathit{q}(z)\}=0$. We drop the explicit reference to $\xi$ in what follows. If the process is not zero-mean, then it must first be centered; that is, we work with fluctuations, ${\mathit{q}}^{\prime }=\mathit{q}-\overline{\mathit{q}}$. In what follows, we assume that this step has been taken and drop the prime.

By assuming that each realization of the process belongs to a Hilbert space with an inner product, typically

The ${\mathit{\varphi }}_k(z)$ are eigenfunctions, with associated eigenvalues ${\lambda }_1\ge {\lambda }_2\ge \cdots \ge 0$ of a Fredholm integral equation whose kernel is the covariance tensor

Aside from the orthogonality of the basis, the expansion gives several desirable properties:

1. (Parseval) The eigenvalues sum to the total variance, ${\sum }_k{\lambda }_k=E\{⟨\mathit{q},\mathit{q}⟩\}$, which in many flow problems represents a kind of energy. For example, if $\mathit{q}$ are velocity fluctuations and $\mathit{W}=1$, then the total variance is (twice) the turbulent kinetic energy integrated over space. Other interpretations are provided below.
2. (Optimality) A truncated expansion with $n$ terms captures more of the total variance than any other orthogonal expansion of the same order.
3. The expansion coefficients are uncorrelated:
   $$E\{a_i{a}_j\}={\lambda }_j{\delta }_{ij}$$(4)

In other words, the decomposition provides a way to express the random field in a series of mutually uncorrelated structures that are optimal for expressing the total variance.

The POD approximates the KL expansion by gathering samples (realizations) for the random process, ${\mathit{q}}^{(k)}$, $k=1,\dots ,N$. The resulting covariance tensor is then estimated:

To apply POD to flowfields, we need to decide on the random variables, $\mathit{q}$, the independent variables, $z$, their domain, $\mathrm{\Omega }$, and of course the means by which we can obtain a sufficiently large ensemble of realizations of the process. The major variants of POD reflect different choices and face different constraints regarding obtaining sufficient data. For example, it is experimentally unrealistic to demand a large number of realizations of a fully three-dimensional, time-dependent flowfield. Likewise, numerical simulation data for the general case may be too voluminous to yield practical computations. The theory discussed in this section is agnostic as to these choices, but they are paramount in determining whether the resulting expansions are able to efficiently represent important flow structures. A takeaway message here is that, regardless of what follows, the most fundamental requirement of POD is that we have a (possibly large) number of independent realizations of a random process so that we hope (and then verify) that $N$ is sufficient so that the resulting eigenvectors converge to the true ones.

B. Symmetries

An important simplification (both theoretically and computationally) is obtained when we impose at the outset certain symmetries, or invariances, that we expect in the random process. In fluid flows, the most relevant spatial invariances are translational (homogeneity), rotational, and reflectional. It is important to understand that these symmetries are statistical in nature—turbulence breaks all of the symmetries that laminar flows enjoy—and yet they reappear in the statistics. Take translational invariance, say, in $x$, as an example: the idea is that a flow structure is as likely to occur at one $x$ as any other. This suggests that the covariance should depend only on the relative position, $\mathit{C}(x,x^{\prime },\dots )\to \mathit{C}(x-x^{\prime },\dots )$. Imposing this restriction on the KL expansion above then yields a set of decoupled problems for each wavenumber in a Fourier-series expansion of the covariance [10,11]. Continuous rotational symmetry behaves similarly and can be cast as independent KL expansions for each azimuthal Fourier mode. Reflectional symmetry yields two separate KL expansions for symmetric and antisymmetric modes. $M$-fold rotational symmetry can be addressed with a Fourier–Floquet ansatz [12], and so on.

Needless to say, when symmetries exist, they should be imposed on the KL ansatz rather than relying on the data to express them for us. First, we typically obtain a reduction in computational effort (for example, by being able to consider a limited number of two-dimensional problems rather than a single, large three-dimensional one), but, more fundamentally, symmetries will be imperfectly expressed in any finite set of realizations: we can observe a turbulent structure at several or many locations in $x$, but never at an infinite number of them. It is a mathematical property of the KL ansatz that guarantees that the reduced SPOD analysis is still overall optimal for the global problem.

The fact that homogeneity leads to a Fourier basis is seen by some as failing of POD in the sense that one might want to describe localized structures that occur in these flows. Our experience is that the computational simplification provided by one or more symmetries is an advantage rather than a pitfall, but if localized descriptions are desired, alternative multiscale bases such as wavelets [13] can be sought. These may be suboptimal in terms of variance/energy but better performing in reduced-order models. Likewise, modifications to the KL ansatz, particularly the imposition of traveling waves together with reconstruction techniques for the phase speed (e.g., [14]), are also possible, although to our knowledge they have not been applied widely, and they may be difficult to pose for large, multidimensional problems.

As mentioned in the last section, the KL theorem requires a compact, bounded covariance tensor, which is contradicted by translational invariance, because the associated coordinate direction is infinite. Fortunately, though, the covariance computed one wavenumber at a time is compact in the remaining, inhomogeneous directions.

C. Stationarity, Standard POD, and SPOD

Stationarity is the temporal analog of translational invariance (homogeneity). Most turbulent flows are idealized as (weak, or wide-sense) stationary processes, and the vast majority of POD algorithms being used across fluid dynamics are based on a stationary assumption. If a flow is not stationary, we can still perform a space–time POD as outlined above, but we require an ensemble of realizations of the process and this is computationally challenging for three-dimensional flows.

When stationarity is an appropriate assumption, it can be exploited in two distinct ways. The first approach considers both spatial and temporal variables in the KL ansatz, that is, $z=[x,t]$, and the inner product is

Then, following the discussion in the previous section, we recognize that $\mathit{C}(x,x^{\prime },t,t^{\prime })\to \mathit{C}(x,x^{\prime },t-t^{\prime })$ and, by an analogously to homogeneity, we find that we can solve a series of POD problems in Fourier space, that is, one frequency at a time. Taking the (continuous) Fourier transform of the covariance gives the CSD tensor:

This ansatz is the basis of spectral POD; all that remains is to determine an appropriate way to estimate the CSD tensor from flow data—this is discussed beginning in Sec. III.

Before proceeding, we discuss a simpler but more restrictive form of POD that exploits stationarity in a different way: restricting the KL basis to only spatial variables, and computing the (space-only) covariance by sampling spatial fields at specific instances from a time series. Stationarity and ergodicity then imply that (when the statistics converge) this will be equivalent to an ensemble average of spatial fields selected from different realizations. The inner product is likewise restricted to an integration over space, and, as a result, the resulting eigenfunctions are not functions of time. Instead, we may expand the time-dependent field by

Because of the ubiquity of this space-only POD problem in the fluids literature, we will refer to it in what follows as “standard” POD. Standard POD has, as its kernel, the same covariance tensor defined above but evaluated at zero time separation, that is, $\mathit{C}(x,x^{\prime },0)$. As such, it is immediately obvious that any notion of temporal correlation among the resulting structures is lost. This criticism lead Schmid [4] to formulate the DMD, which approximates the eigenmodes of a linear operator that maps the state of the flow from one instance to the next. This is more general than SPOD in that these modes can represent structures that grow and decay in time, but the resulting modes are neither orthogonal nor optimal in representing the statistics. For the special case when the flow is stationary, Towne et al. [6] show that, in fact, SPOD modes are DMD modes (provided that the data are centered) and, in fact, represent a kind of optimal averaging of an ensemble of DMD modes computed for different realizations of the same flow. The optimality is precisely the one discussed above: the SPOD modes optimally represent the second-order space–time statistics (covariance).

D. Discrete Form of POD

We now consider a stochastic process represented by a discrete set of $M$ observations rather than a continuously varying field, such as flow quantities discretized at points in space and/or time, that is, $q\in {\mathbb{C}}^M$, where $M$ is the total number of spatial points times the total number of variables times (if desired) the total number of time steps. We allow the flow quantities to be complex to reflect the fact that they may have already been Fourier transformed in one or more homogeneous spatial dimensions. The Hilbert space becomes a finite-dimensional vector space with inner product

The covariance matrix $C$ is, as before, estimated from an ensemble of observations. Denoting $q^{(k)}$ as the $k$th member of the ensemble, we form a data matrix

Regardless of which of the several forms of POD (standard or spectral) discussed above is being computed, the maximum number of nonzero eigenvalues of Eq. (12) is $\min (M,N)$. If the number of observations is much smaller than the ensemble size, $M<N$, then it is most efficient to solve this $M\times M$ directly. More typically in fluid dynamics, however, we have $M\gg N$, and it is not practical to build the matrix $C$, let alone find its eigendecomposition. However, a reduction in cost is possible. The range of $C$ is contained within the span of $Q$, and so we can construct the eigenvectors as a linear combination of the data vectors:

In the literature, solving the (typically) smaller eigenvalue problem (17) is often called the method of snapshots [15], owing to its use in the traditional space-only POD where the ensemble of data is taken (under stationarity) as a set of observations at different times, and we then build the POD eigenvectors from these “snapshots.” Although we retain this jargon, we note that the algebra discussed in this section is completely generic: all that is required is an ensemble (however defined) of discrete data vectors and an appropriately defined inner product (through specification of $W$). Then we solve whichever eigenvalue problem, (12) or (17), is smaller. Finally, it should be obvious that ${\varphi }^{(k)}$, ${\psi }^{(k)}$, and ${\lambda }_k$ are the (generalized) right singular vectors, left singular vectors, and singular values (squared) of the data matrix.

A. Welch’s Method

As discussed above, the SPOD modes are eigenvectors of the sample CSD matrix $S$, which must be estimated from the data. Although there are in principle a number of ways in which this can be done, we focus here on estimation via the Welch periodogram method, which constructs an ensemble of realizations of the temporal Fourier transform of the data from a single time series consisting of $N_t$ snapshots by breaking it into $N_{\mathrm{blk}}$ blocks, or segments. Each segment consists of $N_{\mathrm{FFT}}$ snapshots and overlaps with the next segment by $N_{\mathrm{ovlp}}$ snapshots. The process is the same as that used to compute the power spectral density (PSD) of a single time series. An advantage of using this method is that it is well studied [16,17] and we can exploit the existing theory for selecting the several parameters needed to process the data, and we can estimate the uncertainty in the estimate based on the length and noisiness of the data record.

In what follows, we represent discrete flow data in a slightly different notation then used in the previous section. Namely, we separate the time dependence of the data and denote $q^{(k)}(t_j)$ as the $k$th realization of the vector of observations at discrete times $t_j$. We assume that each realization is available over a period of time, $T$, at regular time intervals, $t_j=t_0+jT/N_{\mathrm{FFT}}$, $j=\mathrm{1,2},\dots ,N_{\mathrm{FFT}}$, where $t_0$, the starting time, is arbitrary because the process is assumed stationary. The discrete Fourier transform in time of each realization and its inverse are

If the original data are real, then transformed data at negative frequencies in the above expression are redundant: they are conjugates of the corresponding positive frequencies, and only positive frequencies need be considered. However, the data are often complex, as is the case when the data have already been Fourier transformed in one or more homogeneous spatial directions. In this case there is a subtlety. In sampled data the negative frequencies will not be redundant. However, the true (as opposed to sample) statistics of the positive and negative frequencies should be identical; for every POD mode of positive frequency there should exist an identical one (with identical eigenvalue) for negative frequency. Therefore, it is advisable to average together positive/negative frequencies.

After transforming the ensemble of data, the remaining analysis is performed one frequency at a time. We construct the data matrix

The crux of Welch’s method is that we can exploit stationarity to construct the ensemble of realizations from a single, longer time series. In essence, different choices of $t_0$ are equivalent to different realizations of the process provided that $T$ is long enough and the overall length of the signal is long enough that we can obtain a large number, $N$, of independent realizations from it. A schematic that outlines the SPOD algorithm based on Welch’s method is shown in Fig. 1. A difficulty, discussed in more detail below, is determining when one has a sufficient number of realizations, but the algorithm is insensitive to inflating the ensemble with dependent ones—these are simply reflected by zero or very small eigenvalues at the end of the SPOD spectrum.

B. Choice of Norm

The choice of inner product, expressed through the weight matrix $W$, plays a central role in SPOD analysis as it determines the optimality and orthogonality properties of the modes. The weight matrix can also be used to apply SPOD to subdomains or a subset of the variables that comprise the solution vector. This is done by assigning zero weight to the portions of the domain or variables that are to be excluded. Alternatively, these parts can be removed from the data before computing the SPOD. The advantage of retaining them with zero weight is that they remain part of the spatially coherent SPOD modes. This permits their physical interpretation, but does not change the SPOD spectrum. In the following, we discuss the four most common choices of norms for (single-phase) flowfields.

C. Choice of Spectral Estimation Parameters

Obtaining an accurate estimate of the CSD matrix $S$ from data depends on appropriately choosing the data sampling parameters, $\mathrm{\Delta }t$ and $N_t$, and the spectral estimation parameters, $N_{\mathrm{FFT}}$ and $N_{\mathrm{ovlp}}$. Fortunately, there is a substantial literature available on how to do this; we summarize the main results here in order to provide a self-contained reference on using SPOD. Further details on these topics can be found in any signal processing textbook, such as Bendat and Piersol [17] and Manolakis et al. [20].

1. Sampling Biases

The sampling time step $\mathrm{\Delta }t$ determines the sampling frequency,

$$f_s=\frac{1}{\mathrm{\Delta }t}$$(24)

and thereby the maximum resolvable (Nyquist) frequency,

$$f_{\max }=\frac{f_s}{2}$$(25)

The number of time steps used to compute the discrete Fourier transform (18), $N_{\mathrm{FFT}}$, determines the period of each “block” of data,

$$T=N_{\mathrm{FFT}}\mathrm{\Delta }t=\frac{N_{\mathrm{FFT}}}{f_s}$$(26)

and the frequency resolution, or bin size,

$$\mathrm{\Delta }f=\frac{1}{T}=\frac{f_s}{N_{\mathrm{FFT}}}$$(27)

The latter corresponds to the lowest resolvable frequency that is not equal to zero. The zero frequency component itself is expected be close to zero if the block mean is close to the true mean.

Accurate spectral estimation requires that $f_s$ and $\mathrm{\Delta }f$ are sufficiently large and small, respectively, such that sampling reproduces the true underlying spectrum. The error that results from discrepancies between the estimation and the true value is referred to as bias. Frequencies above the Nyquist frequency cannot be resolved and, if present in the data, introduce bias in the form of aliasing error. Their contribution is “folded” back along the Nyquist frequency onto lower frequencies. Common techniques to mitigate aliasing are anti-aliasing filters and oversampling. Anti-aliasing filters are low-pass filters that attenuate frequency content above the Nyquist frequency. Oversampling, on the other hand, refers to the practice of sampling at a higher rate than theoretically necessary and to subsequently disregard frequencies that are affected by aliasing.

On the other hand, if the bin size is not sufficiently small, then the spectrum will be biased as distinct content at nearby frequencies will be smeared over the frequency bin. This bias may in principle be made small by choosing a suitably small bin size, or a large $N_{\mathrm{FFT}}$. However, when Welch’s method is used, there is ultimately a compromise between reducing the bias and increasing the uncertainty in the estimate because all the data come from segmenting a single, suitably long time series of size $N_t$. If we choose $N_{\mathrm{FFT}}$ large, then we will not have very many segments in the ensemble, and the spectrum will not be statistically converged. The uncertainty, or variance of the estimate, can be reduced by segmenting the data into overlapping segments, thus increasing the number of segments. With an overlap of $N_{\mathrm{ovlp}}$ points in each segment, the number of total blocks is

$$N_{\mathrm{blk}}=\lfloor \frac{N_t-N_{\mathrm{ovlp}}}{N_{\mathrm{FFT}}-N_{\mathrm{ovlp}}}\rfloor$$(28)

where $\lfloor \cdot \rfloor$ denotes the floor operator.

As an example demonstrating the compromise between bias and uncertainty, consider planar, two-component PIV data in a single measurement plane in the wake of the vertical-axis wind turbine (VAWT). Suppose that $N_t=9000$ consecutive (time-resolved) snapshots are available. Choosing $N_{\mathrm{FFT}}=256$‡ and $N_{\mathrm{ovlp}}=N_{\mathrm{FFT}}/2=128$, corresponding to an overlap of 50%, this segments the data into $N_{\mathrm{blk}}=69$ blocks. The resulting SPOD spectra (eigenvalues) are plotted versus (nondimensional) frequency in Fig. 2, along with the 99% confidence interval computed as discussed below. Although the confidence in the estimate for each bin is reasonable, the frequency bins are large and the peaks, which correspond to the blade-passing frequency (BPF) and its harmonics, show substantial leakage into nearby frequency bins. If we increase the block size to $N_{\mathrm{FFT}}=2048$, we obtain much finer bins, and the BPF and its harmonics are well resolved. But the variance in the estimate is now much larger, and the confidence bounds for the spectral energy in any bin are only certain to within about a factor of 2.

The practice of taking a 50% overlap, or $N_{\mathrm{FFT}}/2$, is a commonly accepted best practice that dates back to Welch’s original work [21]. The author showed that an overlap of 50% reduces the variance of the estimate of the PSD, and therefore of the SPOD spectrum, by approximately a factor of $11/18$ for fixed $N_{\mathrm{FFT}}$ and $N_t$. This is intuitive as we expect to get a more accurate estimate when averaging over more realizations. Overlaps larger than 50% do not yield better results as the realizations become increasingly dependent, in practice, but simply result in a larger number of very small eigenvalues without decreasing the variance of the leading eigenvalues any further.

In all of the above, we assumed that the data are given as a single, long time series and we use the ergodicity assumption to segment the data into presumably independent realizations. Obviously, this is not necessary if the data are given in the form of independent realizations to start with. Take the example of a high-speed camera PIV measurement where internal camera memory restricts the maximum number of consecutive movie frames that can be recorded. In this case, each consecutive block of frames constitutes one realization and there is no need to segment the data in the first place.

The tradeoffs discussed above highlight the importance of balancing the different estimation parameters, in particular during the initial data acquisition process. A good starting point to estimate suitable parameters for an SPOD analysis is to study the PSD of one, or a few, highly resolved time series at some representative locations in the flow. Such data are often readily available, for example, in the form of hot wire or pressure probe measurement data in an experiment, or as time series data from a single grid point of a high-fidelity numerical simulation.

A. Numerical Data of a Turbulent Jet

In this example, we restrict our attention to the symmetric component of the pressure from the LES data. The reader is referred to Schmidt et al. [24] for a detailed study of large-scale coherent structures in turbulent jets for Mach numbers representative of the subsonic, transonic, and supersonic regime, and for different azimuthal wave numbers. We calculate SPOD modes that are optimal with respect to the compressible energy norm introduced in Sec. III.B. Because the data are saved on a nonequidistant cylindrical grid, we use trapezoidal quadrature weights. The instantaneous fluctuating pressure and mean pressure fields are visualized in Fig. 3. The wide range of spatial and temporal scales that are active in the jet is apparent from the instantaneous pressure contours. The following SPOD analysis will demonstrate the method’s ability to systematically separate these different scales. In time, this separation is explicitly enforced through the Fourier transform, whereas the separation of spatial scales results from the interdependence with the temporal scales, as manifested, for example, in the form of an underlying dispersion relation. Figure 4 shows the SPOD energy spectrum for the spectral estimation parameters listed in Table 1. We choose this representation of the SPOD eigenvalues for its compactness and ease of interpretation. It is important to note that the representation by a continuous line does not imply that the same physical mechanism is represented by the same mode at neighboring frequencies. A good indicator for the dominance of one mechanism is a large separation between the eigenvalues associated with the first and the second SPOD modes. In this example, we see such a behavior in the region around $fD/U\approx 0.6$. To understand this behaviour, we next inspect the SPOD eigenfunctions, or modes in Fig. 5. First, we focus on the leading SPOD modes for different frequencies (left column). In all three cases, the leading mode takes the form of a wave train, often referred to as a wavepacket. We also observe that these wavepackets are localized further upstream with increasing frequency. This observation is the key to understanding the origin of the large difference between the leading and the second eigenvalue for $fD/U\approx 0.6$ in Fig. 4. For frequencies $fD/U\gtrsim 0.2$, the dominant coherent structure is the well-known Kelvin–Helmholtz-type (KH) spatial instability of the jet shear layer, whereas Orr-type waves in the region downstream of the potential core dominate at lower frequencies [24]. The main difference between the leading and the second, or subdominant, modes (right column) is that the latter reveal a multilobe structure in the axial direction. Even though a physical interpretation of this phenomenon is not obvious, we can understand this patterning as a mathematical way of achieving orthogonality.

Another important concept that can be inferred from SPOD spectra is referred to as low-rank behavior. Low-rank refers to the property that the data can be well represented by one, or a small number, of basis vectors. Mathematically, this corresponds to a large separation of eigenvalues. Physically, we expect low-rank behavior whenever a physically dominating mechanism, such as a hydrodynamic instability, is present in the flow. For the example of the turbulent jet, the KH instability is such a mechanism and explains the large gap between ${\lambda }_1$ and ${\lambda }_2$ in the frequency band around $fD/U\approx 0.6$. If, on the contrary, all eigenvalues at a given frequencies are similar, as for very low and very high frequencies in Fig. 4, one might speak of non-low-rank behavior. In this case, we cannot expect to observe flow structures that are well represented by a single SPOD mode and are better served by taking a mathematical point of view and interpret the SPOD modes as set of basis functions of comparable importance.

B. Experimental Data of a Vertical-Axis Wind Turbine

The instantaneous and mean streamwise velocity fields for the experimental example are shown in Fig. 6. We choose a weighted 2-norm for the state vector $q=[uv{]}^T$ comprising both velocity components to obtain a (partial) measure of the turbulent kinetic energy; see Sec. III.B. As the measurement windows are not synchronized in time, spatial or temporal correlation between different windows is not meaningful, and we analyze each window separately. If the statistics are sufficiently converged, however, we find that it is possible to phase-match SPOD modes from different windows and provide a global view of the dynamics. This technique is demonstrated later in the context of Fig. 7. The SPOD spectra and the PSD of all seven PIV windows are presented in Fig. 7. The most prominent features are a peak at a low frequency of $fD/u_{\infty }=0.238$ that is found in all seven windows, and three distinct peaks at higher frequencies in the spectrum of PIV window 3, which contains the VAWT. The dominant peak at $fD/u_{\infty }=1.59$ corresponds to the BPF, and the following two peaks to its first two harmonics.

The SPOD modes associated with the two dominant peaks in the spectra are shown in Fig. 8. From the inspection of the leading SPOD mode, we can conclude that the low-frequency peak at $fD/u_{\infty }=0.238$ results from a bluff-body-like wake that manifests downstream of the turbine. This finding is in agreement with the typical bluff-body wake frequency of $fD/u_{\infty }\approx 0.2$. The asymmetry of the wake is a result of the rotation of the turbine and can be explained as follows. On the upper part of the turbine, the rotation direction coincides with the direction of the incoming flow, whereas the opposite is true for the lower part. Therefore, a much higher shear is introduced on the lower side, resulting in a significantly stronger wake. Because the higher frequency peak at $fD/u_{\infty }=1.59$ coincides with the rotor passing frequency, we expect to observe periodic blade vortex shedding at this frequency. This conjecture is confirmed by the SPOD analysis, which isolates this phenomenon in the leading SPOD mode. The comparison of the modes with the instantaneous flow visualization depicted in Fig. 6 shows that SPOD also functions as an efficient filter for measurement noise. The reason is that most types of measurement noise are largely uncorrelated in space and/or time. The same reasoning applies in the context of statistical outliers and rare physical events. Such outliers are associated with a broad frequency response. For that reason, they are not well approximated by SPOD modes. A conditional form of POD that is conducted over a finite time horizon may be applied to extract rare events [25].

We stress that at both frequencies shown in the figure, the SPOD modes were computed independently in each PIV window (because the data were collected in each window in successive experiments). To produce the figure, we “stitched” together the modes in an ad hoc way by finding a single (complex) constant that best matched the phase between the successive windows. The close match between the contours at the boundaries between the PIV windows is strong evidence of a global structure, to which the individual PIV windows are mutually converging.