|
| 1 | +%% EXAMPLE 10: Adaptive sine-taper SPOD. |
| 2 | +% The time-resolved particle image velocimetry data provided along with |
| 3 | +% this example are a subset of the database of a Mach 0.6 turbulent open |
| 4 | +% cavity flow by Zhang et al. [1]. If you are using the database in your |
| 5 | +% research or teaching, please include explicit mention of Zhang et al. |
| 6 | +% [1]. The test database consists of 4000 snapshots of the streamwise |
| 7 | +% velocity component of a turbulent open cavity flow. This is a canonical |
| 8 | +% example of a mixed broadband-tonal turbulent flow, which requires sharp |
| 9 | +% peak resolution at tonal frequencies and smooth spectrum estimates of |
| 10 | +% broadband regions. If you are using this code for your research, please |
| 11 | +% cite Yeung & Schmidt [2]. |
| 12 | +% |
| 13 | +% References: |
| 14 | +% [1] Zhang, Y., Cattafesta, L. N., and Ukeiley, L., Spectral analysis |
| 15 | +% modal methods (SAMMs) using non-time-resolved PIV, Exp. Fluids 61, 226, 1–12, 2020, |
| 16 | +% DOI 10.1007/s00348-020-03057-8 |
| 17 | +% [2] Yeung, B. C. Y., Schmidt, O. T. Adaptive spectral proper orthogonal decomposition |
| 18 | +% of broadband-tonal flows, Theor. Comput. Fluid Dyn. 38, 355–374, 2024, |
| 19 | +% DOI 10.1007/s00162-024-00695-0 |
| 20 | +% |
| 21 | +% B. Yeung ([email protected]) and O. T. Schmidt ([email protected]) |
| 22 | +% Last revision: 10-Sep-2024 |
| 23 | + |
| 24 | +clc, clear variables |
| 25 | +addpath('utils') |
| 26 | +load(fullfile('cavity_data','cavityPIV.mat'),'u','x','y','dt'); |
| 27 | + |
| 28 | +%% Standard Welch SPOD of the test database using one sine window. |
| 29 | +% We manually specify a window length of 256 and an overlap |
| 30 | +% of 128 snaphots. A sine function is automatically used as the window. |
| 31 | + |
| 32 | +% SPOD using a sine window of length 256 and 128 snaphots overlap |
| 33 | +nDFT = 256; |
| 34 | +nOvlp = 128; |
| 35 | +[L,P,f] = spod_adapt(u,nDFT,[],nOvlp,dt); |
| 36 | + |
| 37 | +%% Plot the SPOD spectrum and some modes. |
| 38 | +figure |
| 39 | +loglog(f,L) |
| 40 | +xlabel('frequency'), ylabel('SPOD mode energy'), title('Welch SPOD') |
| 41 | + |
| 42 | +%% Multitaper SPOD using six sine windows. |
| 43 | +% We can achieve higher frequency resolution and resolve down to lower |
| 44 | +% minimum frequencies using six sine windows of length equal to the |
| 45 | +% number of snapshots. Specify nWin as a vector of length 4000. |
| 46 | +% Multitaper SPOD neither requires nor benefits from overlap, so we leave |
| 47 | +% nOvlp empty. |
| 48 | + |
| 49 | +% SPOD using six sine windows of length 4000 and 0 overlap |
| 50 | +nWin = 6*ones(size(u,1),1); |
| 51 | +[L,P,f] = spod_adapt(u,nWin,[],[],dt); |
| 52 | + |
| 53 | +%% Plot the SPOD spectrum. |
| 54 | +figure |
| 55 | +loglog(f,L) |
| 56 | +xlabel('frequency'), ylabel('SPOD mode energy'), title('multitaper SPOD') |
| 57 | + |
| 58 | +%% Adaptive multitaper SPOD. |
| 59 | +% To combine both sharp peak resolution and smooth spectrum estimates of |
| 60 | +% broadband regions, we allow the window number to vary with frequency. |
| 61 | +% The adaptive algorithm requires only a tolerance and determines the |
| 62 | +% window number autonomously. No spectral estimation parameters |
| 63 | +% are needed. |
| 64 | + |
| 65 | +% SPOD using the adaptive algorithm with a tolerance of 1e-4 |
| 66 | +% Can take a long time to compute |
| 67 | +clear opts |
| 68 | +opts.adaptive = true; |
| 69 | +opts.tol = 1e-4; |
| 70 | +[L,P,f,~,nWin] = spod_adapt(u,[],[],[],dt,opts); |
| 71 | + |
| 72 | +%% Plot the SPOD spectrum and window number. |
| 73 | +figure |
| 74 | +yyaxis left |
| 75 | +loglog(f,L) |
| 76 | +xlabel('frequency'), ylabel('SPOD mode energy'), title('adaptive multitaper SPOD') |
| 77 | + |
| 78 | +yyaxis right |
| 79 | +semilogx(f,nWin) |
| 80 | +ylabel('number of windows') |
| 81 | + |
| 82 | +%% Plot some SPOD modes (Rossiter modes) |
| 83 | +figure |
| 84 | +count = 1; |
| 85 | +for fi = [470 770 1063] |
| 86 | + for mi = [1 2] |
| 87 | + subplot(3,2,count) |
| 88 | + contourf(x,y,real(squeeze(P(fi,:,:,mi))),11,'edgecolor','none'), axis equal tight, clim(max(abs(clim))*[-1 1]) |
| 89 | + xlabel('x'), ylabel('y'), title(['f=' num2str(f(fi),'%.2f') ', mode ' num2str(mi) ', \lambda=' num2str(L(fi,mi),'%.2g')]) |
| 90 | + count = count + 1; |
| 91 | + end |
| 92 | +end |
| 93 | + |
| 94 | +% %% Multitple SPOD using non-uniform window number. |
| 95 | +% % If the (potentially non-uniform) window number is known a priori, e.g., |
| 96 | +% % based on physical knowledge, it can be provided directly. As an example, |
| 97 | +% % we use the nWin computed previously to recover the same results as before. |
| 98 | +% |
| 99 | +% % SPOD using sine windows of length 4000, 0 overlap, and |
| 100 | +% % frequency-dependent window number |
| 101 | +% if length(nWin)<size(u,1), nWin(size(u,1)) = 0; end % zero-pad nWin to length 4000 |
| 102 | +% [L,P,f] = spod_adapt(u,nWin,[],[],dt); |
| 103 | +% |
| 104 | +% %% Plot the SPOD spectrum. |
| 105 | +% % The spectrum should be identical to the previous figure. |
| 106 | +% figure |
| 107 | +% loglog(f,L) |
| 108 | +% xlabel('frequency'), ylabel('SPOD mode energy'), title('multitaper SPOD with non-uniform window number') |
0 commit comments