2D spectral energy density using fft2 - energy in spatial domain unequal to that in wavenumber domain

17 Ansichten (letzte 30 Tage)
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

Akzeptierte Antwort

Dr. Seis
Dr. Seis am 5 Jul. 2012
Bearbeitet: Dr. Seis am 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 Kommentare
Alexander
Alexander am 6 Jul. 2012
Bearbeitet: Alexander am 6 Jul. 2012
It is not unambiguously. One might also take the 2.+3. quadrant instead, so there seems to be no unique domain in 2D. Anyway.
You helped me very much. Thanks Elige.
juntian chen
juntian chen am 16 Dez. 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.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by