# Energy Spectrum !!!

 Register Blogs Members List Search Today's Posts Mark Forums Read

April 20, 2014, 23:58
#41
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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
Code:
```clear all;
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?

April 21, 2014, 04:17
#42
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,290
Rep Power: 67
Quote:
 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

April 21, 2014, 04:33
#43
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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,
Also, I tried another case which only contains 1 line of velocity in z direction.The result is still incorrect
Code:
```clear all;
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))```

April 21, 2014, 05:13
#44
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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

April 21, 2014, 05:18
#45
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,290
Rep Power: 67
Quote:
 Originally Posted by huangxianbei Hi,Filippo Sorry for the second reply.I followed your advice, here is my code in matlab Code: ```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???

April 21, 2014, 05:31
#46
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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?

April 21, 2014, 05:40
#47
Senior Member

Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,290
Rep Power: 67
Quote:
 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

April 21, 2014, 05:56
#48
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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

April 26, 2014, 10:00
#49
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 Originally Posted by ehsankf modes should change from -N/2:N/2.
This should be -N/2:N/2-1 for an odd N.

 April 29, 2014, 16:40 #50 New Member   Amjad Join Date: May 2012 Posts: 21 Rep Power: 12 hello, could i get ur code pls, my Email is amjadju@yahoo.com thanks

April 29, 2014, 22:42
#51
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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

 July 8, 2014, 02:10 #52 Member   SHEIKH NASIRUDDIN Join Date: Aug 2013 Location: DELHI Posts: 52 Rep Power: 11 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.

 October 27, 2014, 20:42 #53 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 10 Can I get the code for energy spectrum?

October 27, 2014, 21:08
#54
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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?

October 27, 2014, 21:08
#55
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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

 October 27, 2014, 21:30 #56 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 10 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

October 28, 2014, 00:53
#57
Senior Member

Huang Xianbei
Join Date: Sep 2013
Location: Yangzhou,China
Posts: 302
Rep Power: 12
Quote:
 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

 October 28, 2014, 05:34 #58 Senior Member   chandra shekhar pant Join Date: Oct 2010 Posts: 216 Rep Power: 15 Please send it to me also .. chandukec@gmail.com

 October 29, 2014, 15:33 #59 New Member   Join Date: Nov 2013 Posts: 3 Rep Power: 11 Would you mind sending me the code as well? bmeyer907@gmail.com Thx in advance!

 October 29, 2014, 18:55 Energy spectra #60 New Member   Purvi Join Date: Oct 2014 Posts: 4 Rep Power: 10 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!!

 Tags energy spectrum, fft