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

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

Answer 1

[Dr. Seis]

[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".

[Alexander]

Thanks Elidge. The same works for me too in case of the uncropped F. But when I compute the energy of the unique* part of F, keeping in mind, that I have to multiply the non-DC, non-Nyquist points by 4, I get this discrepancy between E_x and E_k. Amd I dont know why?

• unique: I think the other 3/4 of F is a complex conjugate of the unique 1/4.

[Dr. Seis]

Oh yeah... you cant do that. You can only take 1/2 of the spectrum, since the other half is the complex conjugate. Example to come.

[Alexander]

Ok, I see my mistake. But the question then is, what is in 2D unique?

In your plot (you have shifted the DC to the centre, between the 2 snake lines) their are two possibilities: 3.+4.quadrant OR 4.+1.quadrant. It's a little bit confusing, because then I have to choose between the intervals [-Nyq_x, Nyq_x] x [0, Nyq_y] OR [0, Nyq_x] x [-Nyq_y, Nyq_y].

My aim is to calculate the 1D spectral energy from the above result F. So it is important, which wavevectors (kx,ky) I should choose for further calculation.

But thanks again.

[Dr. Seis]

See edit above.

[Alexander]

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]

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.