Spectral proper orthogonal decomposition (SPOD)

Implements the frequency domain form of proper orthogonal decomposition (POD)
2,3K téléchargements
Mise à jour 19 oct. 2022

Spectral Proper Orthogonal Decomposition in Matlab

SPOD() is a Matlab implementation of the frequency domain form of proper orthogonal decomposition (POD, also known as principle component analysis or Karhunen-Loève decomposition) called spectral proper orthogonal decomposition (SPOD). SPOD is derived from a space-time POD problem for stationary flows [1,2] and leads to modes that each oscillate at a single frequency. SPOD modes represent dynamic structures that optimally account for the statistical variability of stationary random processes.

The large-eddy simulation data provided along with this example is a subset of the database of a Mach 0.9 turbulent jet described in [3] and was calculated using the unstructured flow solver Charles developed at Cascade Technologies. If you are using the database in your research or teaching, please include explicit mention of Brès et al. [3]. The test database consists of 5000 snapshots of the symmetric component (m=0) of a round turbulent jet.

spod.m is a stand-alone Matlab function with no toolbox dependencies. All other Matlab files contained in this repository are related to the six examples that demonstrate the functionality of the code (see file descriptions below). A physical interpretation of the results obtained from the examples can be found in [4]. The reference for the frequency-time analsyis is [5].


Using your browser

Repository zip file with examples (81.5 MB): https://github.com/SpectralPOD/spod_matlab/archive/master.zip

Matlab function only (15 KB): https://raw.githubusercontent.com/SpectralPOD/spod_matlab/master/spod.m

Using Git in the terminal

git clone https://github.com/SpectralPOD/spod_matlab.git


File Description
spod.m Spectral proper orthogonal decomposition in Matlab
example_1.m Inspect data and plot SPOD spectrum
example_2.m Plot SPOD spectrum and inspect SPOD modes
example_3.m Specify spectral estimation parameters and use weighted inner product
example_4.m Calculate the SPOD of large data and save results on hard drive
example_5.m Calculate full SPOD spectrum of large data
example_6.m Calculate and plot confidence intervals for SPOD eigenvalues
example_7_FTanalysis.m Frequency-time analysis
example_8_invspod.m Band-pass filtering using (inverse) SPOD
example_9_multitaperWelch SPOD using Multitaper-Welch estimators
tcoeffs.m Time-continuous expansion coefficients via convolution
invspod.m Inversion of SPOD using block-wise expansion coefficients
jet_data/getjet.m Interfaces external data source with SPOD() (examples 4-5)
utils/trapzWeightsPolar.m Integration weight matrix for cylindrical coordinates (examples 3-6)
utils/jetLES.mat Mach 0.9 turbulent jet test database
LICENSE.txt License


[L,P,F] = SPOD(X) returns the spectral proper orthogonal decomposition of the data matrix X whose first dimension is time. X can have any number of additional spatial dimensions or variable indices. The columns of L contain the modal energy spectra. P contains the SPOD modes whose spatial dimensions are identical to those of X. The first index of P is the frequency and the last one the mode number ranked in descending order by modal energy. F is the frequency vector. If DT is not specified, a unit frequency sampling is assumed. For real-valued data, adjusted one-sided eigenvalue spectra are returned. Although SPOD(X) automatically chooses default spectral estimation parameters, the user is encouraged to manually specify problem-dependent parameters on a case-to-case basis.

[L,P,F] = SPOD(X,WINDOW) uses a temporal window. If WINDOW is a vector, X is divided into segments of the same length as WINDOW. Each segment is then weighted (pointwise multiplied) by WINDOW. If WINDOW is a scalar, a Hamming window of length WINDOW is used. If WINDOW is omitted or set as empty, a Hamming window is used. Multitaper-Welch estimates are computed with the syntax SPOD(X,[NFFT BW],...), where WINDOW is an array of two scalars. NFFT is the window length and BW the time-halfbandwidth product. See [4] for details. By default, a number of FLOOR(2*BW)-1 discrete prolate spheroidal sequences (DPSS) is used as tapers. Custom tapers can be specified in a column matrix of tapers as WINDOW.

[L,P,F] = SPOD(X,WINDOW,WEIGHT) uses a spatial inner product weight in which the SPOD modes are optimally ranked and orthogonal at each frequency. WEIGHT must have the same spatial dimensions as X.

[L,P,F] = SPOD(X,WINDOW,WEIGHT,NOVERLAP) increases the number of segments by overlapping consecutive blocks by NOVERLAP snapshots. NOVERLAP defaults to 50% of the length of WINDOW if not specified.

[L,P,F] = SPOD(X,WINDOW,WEIGHT,NOVERLAP,DT) uses the time step DT between consecutive snapshots to determine a physical frequency F.

[L,P,F] = SPOD(XFUN,...,OPTS) accepts a function handle XFUN that provides the i-th snapshot as x(i) = XFUN(i). Like the data matrix X, x(i) can have any dimension. It is recommended to specify the total number of snaphots in OPTS.nt (see below). If not specified, OPTS.nt defaults to 10000. OPTS.isreal should be specified if a two-sided spectrum is desired even though the data is real-valued, or if the data is initially real-valued, but complex-valued for later snaphots.

[L,P,F] = SPOD(X,WINDOW,WEIGHT,NOVERLAP,DT,OPTS) specifies options: OPTS.savefft: save FFT blocks to avoid storing all data in memory [{false} | true] OPTS.deletefft: delete FFT blocks after calculation is completed [{true} | false] OPTS.savedir: directory where FFT blocks and results are saved [ string | {'results'}] OPTS.savefreqs: save results for specified frequencies only [ vector | {all} ] OPTS.loadfft: load previously saved FFT blocks instead of recalculating [{false} | true] OPTS.mean: provide a mean that is subtracted from each snapshot [ array of size X | 'blockwise' | {temporal mean of X; 0 if XFUN} ] OPTS.nsave: number of most energtic modes to be saved [ integer | {all} ] OPTS.isreal: complex-valuedity of X or represented by XFUN [{determined from X or first snapshot if XFUN is used} | logical ] OPTS.nt: number of snapshots [ integer | {determined from X; defaults to 10000 if XFUN is used}] OPTS.conflvl: confidence interval level [ scalar between 0 and 1 | {0.95} ] OPTS.normvar: normalize each block by pointwise variance [{false} | true]

[L,PFUN,F] = SPOD(...,OPTS) returns a function PFUN instead of the SPOD data matrix P if OPTS.savefft is true. The function returns the j-th most energetic SPOD mode at the i-th frequency as p = PFUN(i,j) by reading the modes from the saved files. Saving the data on the hard drive avoids memory problems when P is large. The modal energy spectra are returned in L, and saved in a separate file 'spod_energy.mat'.

[L,P,F,Lc] = SPOD(...) returns the confidence interval Lc of L. By default, the lower and upper 95% confidence levels of the j-th most energetic SPOD mode at the i-th frequency are returned in Lc(i,j,1) and Lc(i,j,2), respectively. The OPTS.conflvl*100% confidence interval is returned if OPTS.conflvl is set. For example, by setting OPTS.conflvl = 0.99 we obtain the 99% confidence interval. A chi-squared distribution is used, i.e. we assume a standard normal distribution of the SPOD eigenvalues.

[L,P,F,Lc,A] = SPOD(...) returns the block-wise expansion coefficients in A. INVSPOD(P,A,WINDOW,NOVLP) uses A to reconstruct the original data from the SPOD.


[1] Towne, A., Schmidt, O. T., Colonius, T., Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis, J. of Fluid Mech. 847, 821–867, 2018

[2] Lumley, J. L., Stochastic tools in turbulence, Academic Press, 1970

[3] G. A. Brès, P. Jordan, M. Le Rallic, V. Jaunet, A. V. G. Cavalieri, A. Towne, S. K. Lele, T. Colonius, O. T. Schmidt, Importance of the nozzle-exit boundary-layer state in subsonic turbulent jets, J. of Fluid Mech. 851, 83-124, 2018

[4] Schmidt, O. T., Towne, A., Rigas, G., Colonius, T., Bres, G. A., Spectral analysis of jet turbulence, J. of Fluid Mech. 855, 953–982, 2018

[5] Nekkanti, A. and Schmidt, O. T., Frequency–time analysis, low-rank reconstruction and denoising of turbulent flows using SPOD , J. of Fluid Mech. 926, A26, DOI 10.1017/jfm.2021.681, 2021

[6] Schmidt, O. T., Spectral proper orthogonal decomposition using multitaper estimates, Theor. Comput. Fluid Dyn., 1-14, DOI 10.1007/s00162-022-00626-x, https://rdcu.be/cUtP3, 2022

Citation pour cette source

Schmidt, O. T., Spectral proper orthogonal decomposition using multitaper estimates, Theor. Comput. Fluid Dyn., 1-14, DOI 10.1007/s00162-022-00626-x, https://rdcu.be/cUtP3, 2022

Compatibilité avec les versions de MATLAB
Créé avec R2017b
Compatible avec toutes les versions
Plateformes compatibles
Windows macOS Linux
En savoir plus sur Matrix Decomposition dans Help Center et MATLAB Answers

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Les versions qui utilisent la branche GitHub par défaut ne peuvent pas être téléchargées

Version Publié le Notes de version

See release notes for this release on GitHub: https://github.com/SpectralPOD/spod_matlab/releases/tag/1.1.2


See release notes for this release on GitHub: https://github.com/SpectralPOD/spod_matlab/releases/tag/1.1.1

SPOD(X,[NFFT BW],...) computes Multitaper-Welch estimates; [L,P,F,Lc,A] = SPOD(...) returns block-wise expansion coefficients in A; new examples 8 & 9; minor cosmetic changes

literature updated

Frequency-time analysis example and code added (ref. [5])

-references updated

Pour consulter ou signaler des problèmes liés à ce module complémentaire GitHub, accédez au dépôt GitHub.
Pour consulter ou signaler des problèmes liés à ce module complémentaire GitHub, accédez au dépôt GitHub.