2D spectral energy density using fft2 - energy in spatial domain unequal to that in wavenumber domain
15 views (last 30 days)
Show older comments
Hi everybody,
I compute the 2D-spectral-(kinetic)-energy density of a 2D field (in my case the zonal wind component u=u(x,y)).
According to Parseval's theorem the energy in the spatial and wavenumber domain are equal.
I checked this and it works fine, when I compute the energy of the full (uncropped) wavenumber domain.
But in fact I just want the unique part of the fft2 - in the case of 2D- one quarter (more or less). I addintionaly multiply the spectrum by 4 before integrating over wavenumber space. Now the resulting energy is no more exactly the same as the energy computed in spatial domain. In my case E_x/E_k is around 1.1 (regardless whether I multiply by 4 or not!)
This is my code:
% u is the 2D zonal wind component matrix
[m,n] = size(u);
% Energy of u
dx= 2.78; % spatial increment in [km]
fs= 1/dx;
E_x= sum(sum(u.^2))*dx^2;
% Fourier transform
dkm= fs/m;
dkn= fs/n;
km= (0:m-1)*dkm;
kn= (0:n-1)*dkn;
FT= fft2(u,m,n);
%number of unique points
nUpm= ceil((m+1) /2);
nUpn= ceil((n+1) /2);
%adapt this
FT= FT(1:nUpm,1:nUpn);
km= km(1:nUpm);
kn= kn(1:nUpn);
% spectrum
sp = (abs(FT) *dx^2) .^2;
%since I dropped 3/4 of the FFT, multiply by 4 to retain the same amount of energy
%but not multiply the DC or Nyquist frequency components
if rem(m, 2) && rem(n, 2) % odd m,n excludes Nyquist
sp(2:end,2:end) = sp(2:end,2:end)*4;
elseif rem(m, 2) && ~rem(n, 2)
sp(2:end,2:end -1) = sp(2:end,2:end -1)*4;
elseif ~rem(m, 2) && rem(n, 2)
sp(2:end-1,2:end) = sp(2:end-1,2:end)*4;
else % m,n even
sp(2:end-1,2:end -1) = sp(2:end-1,2:end -1)*4;
end
%Energy of sp
E_k= sum(sum(sp)) *dkm*dkn;
% end of code
I would really appreciate if someone could have a look on this problem.
Alexander
0 Comments
Accepted Answer
Dr. Seis
on 5 Jul 2012
Edited: Dr. Seis
on 6 Jul 2012
[EDIT 7/06]
M=8;N=16;
N=8;M=16;
dx=0.1;dy=0.2;
f = randn(M,N);
% Energy in time domain
energy_f = sum(sum(f.^2))*dx*dy
dkx=1/(N*dx);dky=1/(M*dy);
F = fftshift(fft2(f))*dx*dy;
% Energy in wavenumber domain
energy_F = sum(sum(abs(F).^2))*dkx*dky
% Energy using "half" the spectrum
energy_F2 = 2*sum(sum(abs(F(2:M/2, :).^2)))*dkx*dky + ...
2*sum(sum(abs(F(M/2+1,2:N/2).^2)))*dkx*dky + ...
sum(sum(abs(F(M/2+1, 1).^2)))*dkx*dky + ...
sum(sum(abs(F(M/2+1,N/2+1).^2)))*dkx*dky + ...
2*sum(sum(abs(F(1 ,2:N/2).^2)))*dkx*dky + ...
sum(sum(abs(F(1 , 1).^2)))*dkx*dky + ...
sum(sum(abs(F(1 ,N/2+1).^2)))*dkx*dky
% Plot spectrum
Nyq_x = 1/2/dx;
Nyq_y = 1/2/dy;
kx = -Nyq_x : dkx : Nyq_x-dkx;
ky = -Nyq_y : dky : Nyq_y-dky;
imagesc(1:N,1:M,abs(F))
set(gca,'YTick',1:M,'YTickLabel',ky);
set(gca,'XTick',1:N,'XTickLabel',kx);
I drew a GREEN "X" on all the pixels that do not have a complex-conjugate pair - every other pixel has a complex-conjugate "twin" somewhere (some indicated with a RED symbol). The total energy in the spectrum was determined adding the individual stuff marked with the GREEN "X" and by doubling the stuff inside the magenta "box".
6 Comments
juntian chen
on 16 Dec 2021
When I do the above calculation on the actual flow rate data velocity(M,N), M is length in time domain, N is length in space domain, why don't I get any law? The picture is disorganized.
More Answers (0)
See Also
Categories
Find more on Fourier Analysis and Filtering in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!