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

58 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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!