Issue with FFT in MATLAB and product of signals

3 次查看（过去 30 天）
carlos g
carlos g 2020-5-30
评论： Bjorn Gustavsson 2020-5-30
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);
  1 个评论
Bjorn Gustavsson
Bjorn Gustavsson 2020-5-30
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by