-
Notifications
You must be signed in to change notification settings - Fork 61
/
example_7_FTanalysis.m
89 lines (79 loc) · 3.18 KB
/
example_7_FTanalysis.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
%% EXAMPLE 7: Frequency-time analysis.
%
% Reference:
% [1] A. Nekkanti, O. T. Schmidt, Frequency–time analysis, low-rank
% reconstruction and denoising of turbulent flows using SPOD,
% Journal of Fluid Mechanics 926, A26, 2021
%
% A. Nekkanti (aknekkan@eng.ucsd.edu), O. T. Schmidt (oschmidt@ucsd.edu)
% Last revision: 14-Oct-2022 (OTS)
clc, clear variables
addpath('utils')
disp('Loading the entire test database might take a second...')
load(fullfile('jet_data','jetLES.mat'),'p','p_mean','x','r','dt');
%% SPOD
% We will use standard parameters: a Hamming window of length 256 and 50%
% overlap.
nDFT = 128;
nOvlp = floor(nDFT/2);
weight = trapzWeightsPolar(r(:,1),x(1,:));
[L,P,f,~,A] ...
= spod(p,nDFT,weight,nOvlp,dt);
figure
loglog(f,L)
xlabel('frequency'), ylabel('SPOD mode energy')
%% Frequency-time analysis.
% We will compute the mode expansion coefficients of the first 10 modes
% using a windowing function and weights consistent with the SPOD.
nt = size(p,1);
nModes = 3;
a = tcoeffs(p(1:nt,:,:),P,nDFT,weight,nModes);
%% Visualize the results.
t = (1:nt)*dt;
figure
subplot(2,1,1)
pcolor(t,f,abs(squeeze(a(:,1,:)))); shading interp
daspect([100 1 1])
title(['frequency-time diagram (first mode, ' num2str(sum(L(:,1))/sum(L(:))*100,'%3.1f') '\% of energy)'])
xlabel('time'), ylabel('frequency'), caxis([0 0.75].*caxis)
subplot(2,1,2)
pcolor(t,f,squeeze(abs(sum(a,2)))); shading interp
daspect([100 1 1])
title(['frequency-time diagram (sum of first ' num2str(nModes) ' modes, ' num2str(sum(sum(L(:,1:nModes)))/sum(L(:))*100,'%3.1f') '\% of energy)'])
xlabel('time'), ylabel('frequency'), caxis([0 0.75].*caxis)
% %% Reconstruction from continuously-discrete temporal expansion coefficients
% % The inverse SPOD using the continuously-discrete temporal expansion
% % coefficients yields a rank-reduced data reconstruction. We use
% % invspod() to reconstruct an entire block at a time, but only retain the
% % central time instant. This is for demonstration purposes only and
% % absolutely not recommended. Please refer to example 8 for the 'proper'
% % way of reconstruction/filtering using the block-wise expansion
% % coefficients.
%
% nt = 300;
% p_rec = zeros([nt size(x)]);
% for t_i = 1:nt
% p_block = invspod(P,squeeze(a(:,:,t_i)),nDFT,nOvlp);
% p_rec(t_i,:,:) = p_block(ceil(nDFT/2)+1,:,:);
% end
% %
% figure('name','Reconstruction from continuously-discrete temporal expansion coefficients')
% for t_i=1:nt
% subplot(3,1,1)
% pcolor(x,r,squeeze(p(t_i,:,:))-p_mean); shading interp, axis equal tight
% if t_i==1; pmax = max(abs(caxis)); end
% caxis(0.5*pmax*[-1 1]), colorbar
%
% title('Original data')
% caxis(0.5*pmax*[-1 1]), colorbar
% subplot(3,1,2)
% pcolor(x,r,squeeze(p_rec(t_i,:,:))); shading interp, axis equal tight
% title('Reconstructed data')
%
% caxis(0.5*pmax*[-1 1]), colorbar
% subplot(3,1,3)
% pcolor(x,r,(squeeze(p(t_i,:,:))-p_mean)-squeeze(p_rec(t_i,:,:))); shading interp, axis equal tight
% title('Filtered/removed component')
% caxis(0.5*pmax*[-1 1]), colorbar
% drawnow
% end