2D spectral energy density using fft2 - energy in spatial domain unequal to that in wavenumber domain
15 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Alexander
el 5 de Jul. de 2012
Comentada: juntian chen
el 16 de Dic. de 2021
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 comentarios
Respuesta aceptada
Dr. Seis
el 5 de Jul. de 2012
Editada: Dr. Seis
el 6 de Jul. de 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 comentarios
juntian chen
el 16 de Dic. de 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.
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!