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

13 次查看（过去 30 天）
Alexander
Alexander 2012-7-5
评论： juntian chen 2021-12-16
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

请先登录，再回答此问题。

采纳的回答

Dr. Seis
Dr. Seis 2012-7-5
编辑：Dr. Seis 2012-7-6
[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 个评论
显示 4更早的评论
Alexander
Alexander 2012-7-6
编辑：Alexander 2012-7-6
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 2021-12-16
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by