scaling issue in tcoeffs resolved; example 7 update

example_7_FTanalysis.m

@@ -6,52 +6,84 @@
 %         Journal of Fluid Mechanics 926, A26, 2021
 %
 % A. Nekkanti ([email protected]), O. T. Schmidt ([email protected])
-% Last revision: 5-Sep-2022 (OTS)
+% 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','x','r','dt');
+load(fullfile('jet_data','jetLES.mat'),'p','p_mean','x','r','dt');
 
 %% SPOD
-%   We will use standard parameters: a Hann window of length 256 and 50%
+%   We will use standard parameters: a Hamming window of length 256 and 50%
 %   overlap.
-nFFT  	= 128;
-nOvlp 	= floor(nFFT/2);
+nDFT  	= 128;
+nOvlp 	= floor(nDFT/2);
 weight 	= trapzWeightsPolar(r(:,1),x(1,:));
 
-[L,P,f] = spod(p,nFFT,weight,nOvlp,dt);
+[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 3 modes
+%   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,P,nFFT,weight,nModes);
+a               = tcoeffs(p(1:nt,:,:),P,nDFT,weight,nModes);
 
 %% Visualize the results.
-nt  = size(p,1);
 t   = (1:nt)*dt;
 
 figure
 
 subplot(2,1,1)
-pcolor(t,f,abs(squeeze(a(:,:,1)))); shading interp
+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)'])
+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,3)))); shading interp
+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:3)))/sum(L(:))*100,'%3.1f') '% of energy)'])
+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)
 
-%% Alternatively, plot for pre-multiplied spectra
-% pcolor(t,f,abs(squeeze(a(:,:,1)).*f')); shading interp
-% set(gca,'YScale','log')
-% daspect([50 1 1])
-% title('premultiplied frequency-time diagram')
+% %% 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

invspod.m

@@ -15,7 +15,7 @@
 %         Journal of Fluid Mechanics 926, A26, 2021
 %
 % O. T. Schmidt ([email protected])
-% Last revision: 16-Jul-2022 
+% Last revision: 7-Oct-2022 
 
 dim     = size(P);
 if ndims(A)==2
@@ -24,6 +24,7 @@
     nBlks   = dim(end);
 end
 nFreqs  = dim(1);
+nModes  = size(A,2);
 if length(window)==1
     nDFT        = window;
     window      = hammwin(window);
@@ -53,7 +54,7 @@
     xHatBlk(:)  = 0;
 
     % reconstruct Fourier realization
-    xHatBlk(1:nFreqs,:)     = squeeze(sum(P.*squeeze(A(:,:,blk_i)),2));
+    xHatBlk(1:nFreqs,:)     = squeeze(sum(P(:,1:nModes,:).*squeeze(A(:,:,blk_i)),2));
     if issymm
         xHatBlk(nDFT/2+2:end,:)  = conj(xHatBlk(nDFT/2:-1:2,:));
     end

tcoeffs.m

@@ -5,7 +5,7 @@
 %   continuously-discrete temporal SPOD mode expansion coefficients of the
 %   leading NMODES modes. P is the data matrix of SPOD modes returned by
 %   SPOD. X, WINDOW and WEIGHT are the same variables as for SPOD. If no
-%   windowing function is specified, a Hann window of length WINDOW will
+%   windowing function is specified, a Hamming window of length WINDOW will
 %   be used. If WEIGHT is empty, a uniform weighting of 1 is used.
 %
 %   Reference:
@@ -14,7 +14,7 @@
 %         Journal of Fluid Mechanics 926, A26, 2021
 %
 % A. Nekkanti ([email protected]), O. T. Schmidt ([email protected])
-% Last revision: 12-Sep-2022 
+% Last revision: 7-Oct-2022 Brandon Yeung <[email protected]>
 
 dims        = size(X);
 nt          = dims(1);
@@ -36,7 +36,7 @@
 end
 
 X           = reshape(X,nt,nGrid);
-X           = X-mean(X,1);              % subtrct mean
+X           = X-mean(X,1);              % subtract mean
 
 ndims       = size(P);
 P           = permute(P,[1 length(ndims) 2:length(ndims)-1]);
@@ -49,28 +49,26 @@
 end
 
 weight     = reshape(weight,1,nGrid);
-% padding the data with zeroes at both ends
+% zero-padding
 X          = [zeros(ceil(nDFT/2),nGrid); X; zeros(ceil(nDFT/2),nGrid);];
-a          = zeros(nFreq,nt,nModes);
-winCorr_fac= nDFT/winWeight;
+a          = zeros(nFreq,nModes,nt);
+winCorr_fac= winWeight/nDFT;
 disp(' ')
 disp('Calculating expansion coefficients')
 disp('------------------------------------')
 for i=1:nt
     X_blk               = fft(X(i:i+nDFT-1,:).*window);
     X_blk               = X_blk(1:nFreq,:);
-    % correction for windowing of zero-padded data
-    if (i<ceil(nDFT/2)-1)
-        corr 	= sqrt(winCorr_fac/sum(window(ceil(nDFT/2)-i+1:nDFT)));
-    elseif (i>=nt-ceil(nDFT/2)+1)
-        corr	= sqrt(winCorr_fac/sum(window(1:nt+ceil(nDFT/2)-i)));
+    % correction for windowing and zero-padding
+    if (i<ceil(nDFT/2)+1)
+        corr 	= 1/(winCorr_fac*sum(window(ceil(nDFT/2)-i+1:nDFT)));
+    elseif (i>nt-ceil(nDFT/2)+1)
+        corr	= 1/(winCorr_fac*sum(window(1:nt+ceil(nDFT/2)-i)));
     else
         corr 	= 1;
     end
-    for j=1:nFreq
-        for l=1:nModes
-            a(j,i,l)    = corr*(squeeze(P(j,l,:)))'*(weight.*X_blk(j,:)).';
-        end
+    for l=1:nModes
+        a(:,l,i) = corr*winCorr_fac*dot(squeeze(P(:,l,:)),weight.*X_blk,2);
     end
     disp(['time ' num2str(i) '/' num2str(nt)])
 end