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

[kgoodwin]

Thank you for your code, FMDenaro. However, it's difficult for me to understand your process without your data. Does anyone have the math formulas they used for the autocorrelation? Again, the steps I'm using are:
1. Calculate autocorrelation to get signal into wave number space
2. Take fft
3. Calculate energy
E(k) = E(k) + abs(V(i,j,m))^2 + abs(U(i,j,m))^2 + abs(W(i,j,m))^2;
4. plot energy vs wave number

[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
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？

[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

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