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
Best Answer