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  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. . 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 .
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
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
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
 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
 Lumley, J. L., Stochastic tools in turbulence, Academic Press, 1970
 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
 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
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
Regarding your first question: increasing the frequency bins by reducing the number of snapshots per block is one option, as you mentioned (you're sacrificing frequency resolution doing that though). A better option is to use more data, if available/computable, to get more realizations of the Fourier transform. Regarding your second question: assuming that the sampling frequency is fixed, your only option here is to increase the number of snapshots per block. It really sounds like you'll need more data to get nicely converged modes...
Hi Oliver. Do you know any form to reduce the noise that could appear on the figures? I'm trying to calculate the first SPOD mode from some frequencies and I can obtain a cleaner figure from low frequencies instead from highest. I reduced the number of snapshots per block, reaching some clean figures but only for certain frequencies.
Also I'm trying to replicate some results that I have previously obtained using a DMD on my data, so I want to calculate the SPOD modes using an array with some frequencies that I founded relevant after calculate the DMD modes. It is possible to introduce manually some frequencies instead of those that the program calculate itself?
Have you tried out the examples? Each of the 6 examples builds on the previous one and introduces new features. The examples also show you how to visualize the results.
Hi Oliver, Do you have any tutorial to work with this program? I appreciate you if you can prepare a tutorial on youtube or anywhere.
Below this message is a code snippet that you can append at the end of example_2.m to animate the modes on screen.
%% Animate the same modes.
% Note how all wavepackets travel at approximately the same phase
% speed c_ph. The reason is that their streamwise wavenumber k_x changes
% with frequency such that c_ph = omega/k_x is approximately constant.
nt = 30;
T = 1/f(10); % period of the 10th frequency
time = linspace(0,T,nt); % animate over one period
count = 1;
for ti = 1:nt
for fi = [10 15 20]
for mi = [1 2]
pcolor(x,r,real(squeeze(P(fi,:,:,mi)*exp(2i*pi*f(fi)*time(ti))))), shading interp, axis equal tight, caxis(max(abs(caxis))*[-1 1])
xlabel('x'), ylabel('r'), title(['$f=' num2str(f(fi),'%.2f$') ', mode ' num2str(mi) ', $\lambda=' num2str(L(fi,mi),'%.2g$')])
xlim([0 10]); ylim([0 2])
count = count + 1;
count = 1;
I wonder if is it possible to see the time evolution or oscilation of one coherent structure or several coherent structures with SPOD.
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!