Issue with FFT in MATLAB and product of signals

I'm having an issue with FFT in MATLAB. This is what I do:

1) 2 signals (Data1 and Data2) are imported into MATLAB. They are signals in the 3D space [NX NY NZ]

2) Directions X and Z are periodic, Y is not. Therefore I would like to FFT in these two directions.

3) I calculate a premultiplied spectra in the z-direction of the product of both signals in Fourier space, for which I average on x direction

4) I plot the spectra.

But I didn't obtain was I supposed to. I have a FORTRAN code which does this using FFTW so points 1, 3 and 4 work therefore I suspect I am not computing the FFTs correctly. Does anyone have a hint on where the mistake is?

clear all

clc


%% Import signals

namefile=['end.h5'];    

nVar1='W';

nVar2='W';

res=double(h5readatt(namefile,'/','Resolution'));

NX=res(1);NY=res(2);NZ=res(3);

retau=2061.465620;rel=1/0.000014;

utau=retau/rel;

deltanu=1/utau/rel;

xl=3;zl=1.5;

ye=129;

ymed=(ye+1)/2;

yref=104;

xgrid=[0:NX-1]*xl/NX;

zgrid=[0:NZ-1]*zl/NZ;


for i=2:NX/2+1

lamx(i)=xl/double(i-1);

kappax(i)=2.0*pi/lamx(i);

end

for k=2:NZ/2+1   

lamz(k)=zl/double(k-1);

kappaz(k)=2.0*pi/lamz(k);

end

kappax(1)=0.0;kappaz(1)=0.0;


gridfile=['grid.h5'];

y=h5read(gridfile,'/grids/y');

yf=0.5*(y(1:end-1)+y(2:end));

ygrid=yf';


Data1=h5read(namefile,['/Timestep/' nVar1],[1 1 1],[NX NY NZ]); 

Data2=h5read(namefile,['/Timestep/' nVar2],[1 1 1],[NX NY NZ]); 


%% Permutation to Z,X,Y and FFT along Z and X directions

Data1=permute(Data1,[3,1,2]);  

Data2=permute(Data2,[3,1,2]);  

for j=1:NY

for i=1:NX
    
for k=1:NZ
       
fData1(k,i,j)=fft2(Data1(k,i,j));     
       
fData2(k,i,j)=fft2(Data2(k,i,j));  
    
end

end

end


%% Recovery of the signals in Fourier space

Fsx=NX/xl;

Fsz=NZ/zl;

fx = Fsx*(0:(NX/2))/NX;

fz = Fsz*(0:(NZ/2))/NZ;

P21 = abs(fData1/NX/NZ); 

P11 = P21(1:NX/2+1,1:NZ/2+1,:); 

P11(2:end-1,2:end-1,:) = 2*P11(2:end-1,2:end-1,:); 

P22 = abs(fData2/NX/NZ);

P12 = P22(1:NX/2+1,1:NZ/2+1,:); 

P12(2:end-1,2:end-1,:) = 2*P12(2:end-1,2:end-1,:);


%%Premultiplied z-spectra of the product of signals

phiz=zeros(NY,NZ/2+1);

for j=1:NY

for k=1:NZ/2+1

for i=1:NX/2+1

phiz(j,k)=phiz(j,k)+(conj(P11(k,i,j))*P12(k,i,j))/(NX/2+1);

end

phikz(j,k)=phiz(j,k)*kappaz(k);

end

end


%%Plot

vecto_prod=zeros((NY+1)/2,NZ/2);

factor=2;

facto=(utau^4)*rel;

for indi=0:(NY+1)/2-1

vecto_prod(indi+1,:)=phikz(indi+1,2:NZ/2+1)/factor+phikz(NY-indi,2:NZ/2+1)/factor;

end

figure1=figure(1)

axes1 = axes('Parent',figure1,'YScale','log','XScale','log');

grid(axes1,'on');

hold(axes1,'on');

[C,h]=contourf(lamz(2:NZ/2+1)*retau,(yf(1:ymed)+1.0)*retau,(2*pi/zl)*retau*kappaz(2:NZ/2+1).*(yf(1:ymed)+1.0).*vecto_prod/facto,[10],'LineColor','none');

clabel(C,h);

[Bjorn Gustavsson]

The first thing you have to fix is your calculation of the Fourier-transforms in the x-z directions. In your first tripple-loop you call fft2 with a single scalar:

fData1(k,i,j)=fft2(Data1(k,i,j));

Fourier-transforming a single scalar is not very productive. I guess you would prefer to do something like this:

for j3 = 1:NY
  fData1(k,i,j)=fft2(Data1(:,:,j3));     
  fData2(k,i,j)=fft2(Data2(:,:,j3));  
end

Then you will at least calculate the 2-D fft of each x-z slice of Data1 and Data2.