|
FFTW 2d r2c 将Nx*Ny 实数数组变换后,复数数组的结构为Nx*(Ny/2+1), 花了半天的时间才将FFTW的数组结构整明白。 下面是自己写的一段测试代码:
#include<complex.h>
#include<fftw3.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define Nx 5
#define Ny 5
#define Nz 5
int main()
{
int i,j,u;
int k1,k2,k3;
double *in;
double complex *F;
fftw_complex *out;
fftw_plan p;
in = (double*)malloc(sizeof(double) * Nx * Ny * Nz);
F = (double complex*)malloc(sizeof(double complex) * Nx * Ny * Nz);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * Nx * Ny * (Nz/2+1));
if((in==NULL)||(out==NULL)||(F==NULL))
{
printf("Error:insufficient available memoryn");
}
else
{
for(i=0; i<Nx*Ny*Nz; i++)
{
in[i] = exp(-1.0 * i / Nx / Ny / Nz);
}
}
p = fftw_plan_dft_r2c_3d(Nx,Ny,Nz,in,out,FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
fftw_cleanup();
for(k1=0;k1<Nx;k1++){
for(k2=0;k2<Ny;k2++){
for(k3=0;k3<Nz;k3++){
F[k3+Nz*(k2+Ny*k1)] = 0.0 + 0.0*I;
for(i=0;i<Nx;i++){
for(j=0;j<Ny;j++){
for(u=0;u<Nz;u++){
F[k3+Nz*(k2+Ny*k1)] += in[u+Nz*(j+i*Ny)]*cexp(-2.0*M_PI*i*k1*I/Nx)*cexp(-2.0*M_PI*j*k2*I/Ny)*cexp(-2.0*M_PI*u*k3*I/Nz);
}
}
}
}
}
}
for(i=0;i<Nx;i++)/*OUTPUT*/{
for(j=0;j<Ny;j++){
for(u=0;u<Nz;u++){
if(u<Nz/2+1)
printf("(%f %lfi)n", creal(out[u+(Nz/2+1)*(j+Ny*i)] - F[u+Nz*(j+i*Ny)] ), cimag(out[u+(Nz/2+1)*(j+Ny*i)] - F[u+Nz*(j+i*Ny)]) /*,creal(F[j+Ny*i]), cimag(F[j+Ny*i])*/);
else
printf("(%f %lfi)n", creal(conj(out[(Nz-u)%Nz+(Nz/2+1)*((Ny-j)%Ny + Ny*((Nx-i)%Nx))])-F[u+Nz*(j+i*Ny)]), cimag(conj(out[(Nz-u)%Nz+(Nz/2+1)*((Ny-j)%Ny + Ny*((Nx-i)%Nx))])-F[u+Nz*(j+i*Ny)])/*,creal(F[j+Ny*i]), cimag(F[j+Ny*i])*/);
}
}
}
free(in);
free(F);
fftw_free(out);
return 0;
}
在对复数数组进行逆变换时, 输入的复数数组一定要满足Nx*(Ny/2+1),结构才能得到正确的实空间Nx*Ny
维数组。 实例代码:
#include<complex.h>
#include<fftw3.h>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#define Nx 5
#define Ny 5
#define Nz 5
int main()
{
int i,j,u;
int k1,k2,k3;
double *in, *out;
double complex *F;
fftw_plan p;
in = (double*)malloc(sizeof(double) * Nx * Ny * Nz);
F = (double complex*)malloc(sizeof(double complex) * Nx * Ny * (Nz/2+1));
out = (double*)malloc(sizeof(double) * Nx * Ny * Nz);
if((in==NULL)||(out==NULL)||(F==NULL))
{
printf("Error:insufficient available memoryn");
}
else
{
for(i=0; i<Nx*Ny*Nz; i++)
{
in[i] = exp(-1.0 * i / Nx / Ny / Nz);
}
}
for(k1=0;k1<Nx;k1++){
for(k2=0;k2<Ny;k2++){
for(k3=0;k3<Nz/2+1;k3++){
F[k3+(Nz/2+1)*(k2+Ny*k1)] = 0.0 + 0.0*I;
for(i=0;i<Nx;i++){
for(j=0;j<Ny;j++){
for(u=0;u<Nz;u++){
F[k3+(Nz/2+1)*(k2+Ny*k1)] += in[u+Nz*(j+i*Ny)]*cexp(-2.0*M_PI*i*k1*I/Nx)*cexp(-2.0*M_PI*j*k2*I/Ny)*cexp(-2.0*M_PI*u*k3*I/Nz);
}
}
}
}
}
}
p = fftw_plan_dft_c2r_3d(Nx,Ny,Nz,F,out,FFTW_ESTIMATE);
fftw_execute(p); /* repeat as needed */
fftw_destroy_plan(p);
fftw_cleanup();
for(i=0; i<Nx*Ny*Nz; i++)
{
printf("%lfn", (in[i]-out[i]/Nx/Ny/Nz));
}
free(in);
free(F);
free(out);
return 0;
}
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-27 06:05
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社