[huangxianbei]

Originally Posted by FMDenaro

just consider a 3d Cartesian mesh and (j,k) some fixed positions. This is some lines of my Matlab code for the energy spectrum

% ____________ E_uu (stream-wise)______________
% ff is a dump array for using the FFT of the field u(i,j,k)
for i = 1:np+1, ff(i) = u(i,j,k); end
% Discrete transform: computation of the Fourier coefficient
c=fft(ff,np); for q = 1:np, c(q)=c(q)/np; end
% Computation of the modulus squared. Cycle for all the wavenumbers
for q=1:np, cu(q,k)=c(q)*conj(c(q)); end

Then, you must simply plot your power spectrum.
Generally, the computed spectrum are also computed along the whole homogeneous directions (e.g. for all the k) and then averaged to have a statistically meaning.

Hi,Filippo
Sorry for the second reply.I followed your advice, here is my code in matlab

clear all;
[x y z Ux Uy Uz]=textread('U_constantPlane.raw','%f%f%f%f%f%f','headerlines',2);
 N=32768;
 Lz=6.28;
 
 kz=(2*pi*(1:256)/Lz);
 kzz=round(kz);
 w=reshape(Uz,256,128);
 wk=fft(w);
for i = 1:256,
    for j=1:128,
    wk(i,j)=wk(i,j)/N;
    end
end
for i = 1:256,
    for j=1:128,
    R(i,j,kzz)=wk(i,j)*conj(wk(i,j)); 
    end
end
for j=1:128,
    E(kzz)=0.5.*R(j,j,kzz(j))*4*pi.*kzz(j)^2;
end
plot(kzz,E(kzz))

I have the velocity on a plane and the points are 256*128, so I turned the U into a 256*128 matrix. By doing the whole program, the result is zero! I don't know why, could you please help me?

[FMDenaro]

Originally Posted by huangxianbei

Hi,Filippo
I think the third line should be: cu(q)=c(q)*conj(c(q))
As known, the velocity spectrum tensor is
phi[i,j](k)=cu(k),the energy spectrum is integrated to
E(k)=1/2*cu(k)*dS(k)
where S(k) is the sphere's surface area with radius mag(k)
is this means to calculate E(k)=1/2*cu(k)*4*pi*k^2？

yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis

[huangxianbei]

Originally Posted by FMDenaro

yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis

Hi,Filippo,
Thank you for your reply, could you please have a look at my previous reply for the code?
Also, I tried another case which only contains 1 line of velocity in z direction.The result is still incorrect

clear all;
[z Ux Uy Uz]=textread('Zleft_U.xy','%f%f%f%f');
 N=128;
 Lz=6.28;
 kz=(1:N)*2*pi/Lz;
 k=round(kz);
 c=fft(Uz,128);
  for q = 1:N, 
      c(q)=c(q)/N; 
  end
  for q=1:N,
      cu(q)=c(q)*conj(c(q)); 
  end
  %up to now, the velocity-spectrum is obtained,however, the q~=k,
  %how to match them? 
  for q=1:N,
      E(k)=0.5*cu(q)*4*pi*k.^2;
  end
  plot(k,E(k))

[huangxianbei]

Originally Posted by FMDenaro

yes, the second integer k is just becasue I stored in a matrix for all the points along the z-axis

I managed to match them with this:
E=0.5.*cu.*4*pi.*k.^2;
Again,the result is incorrect

[FMDenaro]

Originally Posted by huangxianbei

Hi,Filippo
Sorry for the second reply.I followed your advice, here is my code in matlab

clear all;
[x y z Ux Uy Uz]=textread('U_constantPlane.raw','%f%f%f%f%f%f','headerlines',2);
 N=32768;
 Lz=6.28;
 
 kz=(2*pi*(1:256)/Lz);
 kzz=round(kz);
 w=reshape(Uz,256,128);
 wk=fft(w);
for i = 1:256,
    for j=1:128,
    wk(i,j)=wk(i,j)/N;
    end
end
for i = 1:256,
    for j=1:128,
    R(i,j,kzz)=wk(i,j)*conj(wk(i,j)); 
    end
end
for j=1:128,
    E(kzz)=0.5.*R(j,j,kzz(j))*4*pi.*kzz(j)^2;
end
plot(kzz,E(kzz))

I have the velocity on a plane and the points are 256*128, so I turned the U into a 256*128 matrix. By doing the whole program, the result is zero! I don't know why, could you please help me?

w=reshape(Uz,256,128);
wk=fft(w)

this way, w is a matrix???

[huangxianbei]

Originally Posted by FMDenaro

w=reshape(Uz,256,128);
wk=fft(w)

this way, w is a matrix???

yes,because z values are ranged in a periodic way(256 values one loop), so 32768 values corresponding to 128 loops, representing the whole value on a plane. Is this resonable?
I'm now confused about the wavanumber k, is this write as:
(1:N)*2*pi/Lz
each one is connected to a cu value in order?

[FMDenaro]

Originally Posted by huangxianbei

yes,because z values are ranged in a periodic way(256 values one loop), so 32768 values corresponding to 128 loops, representing the whole value on a plane. Is this resonable?
I'm now confused about the wavanumber k, is this write as:
(1:N)*2*pi/Lz
each one is connected to a cu value in order?

therefore, fft(w) returns the Fourier transform of each column of the matrix and store in wc.
Anyway, I suggest checking for a simple test: generate a sin(k*x) function (try for some k value) and store in a vector, after doing the fft plot the spectrum and check if the energy density is located at the right frequency

[huangxianbei]

Originally Posted by FMDenaro

therefore, fft(w) returns the Fourier transform of each column of the matrix and store in wc.
Anyway, I suggest checking for a simple test: generate a sin(k*x) function (try for some k value) and store in a vector, after doing the fft plot the spectrum and check if the energy density is located at the right frequency

All right, I'll try it. Another question is still needed to be solved as I have asked:
what's definition formula of k?
2*pi*n/Lz?

Best reguards

[huangxianbei]

Originally Posted by ehsankf

modes should change from -N/2:N/2.

This should be -N/2:N/2-1 for an odd N.

[Amjad Asad]

hello,
could i get ur code pls, my Email is amjadju@yahoo.com
thanks

[huangxianbei]

Originally Posted by Amjad Asad

hello,
could i get ur code pls, my Email is amjadju@yahoo.com
thanks

you can find it in http://www.cfd-online.com/Forums/mai...tml#post488460

[nasir.iitd]

Dear Harishg,
I 'm also facing the same problem and want to plot energy spectra using matlab. Can you send me the code you have used. My e-mail is:nasir.iitd@gmail.com
Your help is highly appreciated.
nasiruddin.

[p3vdim]

Can I get the code for energy spectrum?

[huangxianbei]

Originally Posted by Amjad Asad

hello,
could i get ur code pls, my Email is amjadju@yahoo.com
thanks

Sorry that I just check my thread. Do u still need it?

[huangxianbei]

Originally Posted by p3vdim

Can I get the code for energy spectrum?

Yes, most of my code is in this thread, if you need, I'll send you a .m file with email

[p3vdim]

Yes, that will be great and if you can send me understanding of it too!!!
Thanks
Purvi
my email address is : p3vdim@yahoo.co.in

[huangxianbei]

Originally Posted by p3vdim

Yes, that will be great and if you can send me understanding of it too!!!
Thanks
Purvi
my email address is : p3vdim@yahoo.co.in

Hi,
I have sent you the code. The explaination is not so much for the limited time, hoping this can help you.

Xianbei

[chandra shekhar pant]

Please send it to me also ..
chandukec@gmail.com

[Leon]

Would you mind sending me the code as well? bmeyer907@gmail.com
Thx in advance!

[p3vdim]

Hi Guys,
I am also having problems with plotting energy spectrum. Could somebody plz send me the code.
My email address is : p3vdim@yahoo.co.in
Thanks in advance!!