Bessel的个人博客分享 http://blog.sciencenet.cn/u/Bessel

博文

练习题: FFTW 2d r2c 和c2r数组的结构问题

已有 4467 次阅读 2016-9-7 15:50 |系统分类:科研笔记

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;
    }





https://blog.sciencenet.cn/blog-301704-1001443.html

上一篇:练习题:FFTW计算卷积
下一篇:C编译器中制定linker的路径
收藏 IP: 159.226.49.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 22:21

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部