fortran 复杂到真正的 fftw 问题

Question

我目前正在做一个需要实现傅立叶变换和逆变换的项目。我正在测试从在线示例修改的程序；打印或写入命令通常用于调试目的：

  program testit
  INCLUDE 'fftw3.f'
  double complex out!, in
  real in
  parameter (N=100)
  dimension in(N), out(N)
  integer*8 p,p2
  integer i,j
  real x
  real fact
  write(*,*)"stuff in data"
  OPEN(UNIT=12, FILE="input.txt", ACTION="write", STATUS="replace")
  OPEN(UNIT=20, FILE="dftoutput.txt", ACTION="write", STATUS="replace")
  x=0
  in = 0
  do i=1,N/2
    in(i)=1
  enddo
  do i=1,N
  write(*,"(f10.2,1x,f10.2)")in(i)
  WRITE(12,*)real(in(i))
  enddo
  write(*,*)"create plans"
  call dfftw_plan_dft_r2c_1d(p ,N,in,out,FFTW_ESTIMATE)
  call dfftw_plan_dft_c2r_1d(p2,N,in,out,FFTW_ESTIMATE)
  write(*,*)"do it"
  call dfftw_execute_dft_r2c(p,in,out)
  do i=1,N
    write(*,"(f12.4,1x,f12.4)")out(i)
    WRITE(20,*)abs(out(i))
  enddo
  write(*,*)"undo it"
  call dfftw_execute_dft_c2r(p2,in,out)
  fact=1.0/N
  do i=1,N
    write(*,)in(i)
    write(*,)out(i)
  enddo
  write(*,*)"clean up"
  call dfftw_destroy_plan(p,in,out)
  call dfftw_destroy_plan(p2,in,out)
  end program

真实到复杂的转换效果很好。逆变换给出了错误的值，并以某种方式修改了输入和输出变量。我不知道问题是什么，我无法在网上找到任何答案。感谢您的帮助。

提前致谢！

乍得·W·弗里尔

编辑：我还想知道在 fftw 包中是否有与 matlab 中的 fftshift() 和 ifftshift() 类似的功能。

Answer 1

存在精度问题：在大多数计算机上，real 指的是单精度 32 位浮点数。 real*8 可用于指定双精度浮点数，与double complex 和dfftw_... 保持一致。

对于 FFTW 实数到复数的变换，如果输入的大小是N real*8，那么size of the output 就是N/2+1 double complex。由于输入是真实的，coefficients of negative frequencies (higher than N/2+1) are conjugate of positive frequencies，因此避免了它们的计算。

根据FFTW关于planner flags的文档，

FFTW_PRESERVE_INPUT 指定异地转换不得更改其输入数组。这通常是默认值，除了 c2r 和 hc2r（即复数到实数）转换，FFTW_DESTROY_INPUT 是默认值...

如果您希望将系数保存在傅立叶空间out 中，请复制数组或尝试使用FFTW_PRESERVE_INPUT。如果您有足够的内存，第一个解决方案似乎是最好的方法。

FFTW 函数的参数顺序始终为origin,destination。由于从out 到in 执行反向变换：

call dfftw_plan_dft_c2r_1d(p2,N,out,in,FFTW_ESTIMATE)
...
call dfftw_execute_dft_c2r(p2,out,in)

这是基于您的代码，它执行向前和向后 fft。由gfortran main.f90 -o main -lfftw3 -lm -Wall编译：

program testit
INCLUDE 'fftw3.f'
double complex out!, in
real*8 in
parameter (N=100)
dimension in(N), out((N/2+1))
integer*8 p,p2
integer i
real x
real fact
write(*,*)"stuff in data"
OPEN(UNIT=12, FILE="input.txt", ACTION="write", STATUS="replace")
OPEN(UNIT=20, FILE="dftoutput.txt", ACTION="write", STATUS="replace")
x=0
in = 0
do i=1,N/2
  in(i)=1
enddo
do i=1,N
  write(*,"(f10.2,1x/)")in(i)
  WRITE(12,*)real(in(i))
enddo
write(*,*)"create plans"
call dfftw_plan_dft_r2c_1d(p ,N,in,out,FFTW_ESTIMATE)
call dfftw_plan_dft_c2r_1d(p2,N,out,in,FFTW_ESTIMATE)
write(*,*)"do it"
call dfftw_execute_dft_r2c(p,in,out)
do i=1,N/2+1
  write(*,"(f12.4,1x,f12.4/)")out(i)
  WRITE(20,*)abs(out(i))
enddo
write(*,*)"undo it"
call dfftw_execute_dft_c2r(p2,out,in)
fact=1.0/N
write(*,*)"input back, except for scaling"
do i=1,N
  write(*,"(f10.2,1x)")in(i)
enddo
write(*,*)"output may be modified by c2r...except if FFTW_PRESERVE_INPUT is set"
do i=1,N/2+1
  write(*,"(f10.2,1x,f10.2/)")out(i)
enddo
write(*,*)"clean up"
call dfftw_destroy_plan(p,in,out)
call dfftw_destroy_plan(p2,out,in)
end program

在前向 fft in->out 和后向 fft out->in 之后，in 与其初始值相同，除了比例因子 N。