fortran 复杂到真正的 fftw 问题 [英] fortran complex to real fftw issue

问题描述

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

I am currently working on a project where I need to implement a fourier transform and an inverse transform. I was testing a program I modified from an online example; the print or write commands are normally for debugging purposes:

  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

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

The real to complex transformation works just fine. The inverse transformation gives wrong values and somehow it modifies both the in and out variables. I do not know what the problem is and I was not able to figure out any answers online. Help is appreciated.

提前致谢！

乍得·W·弗里尔

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

I was also wondering if there is a similar function to fftshift() and ifftshift() from matlab in the fftw package.

推荐答案

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

There is a precision issue : on most computers, real refers to single precision 32bit floating point numbers. real*8 can be used to specify double precision floating point numbers, to be consistent with double complex and dfftw_....

对于FFTW实数到复数的变换，如果输入的大小是N real*8，那么输出大小是N/2+1双复数.由于输入是实数，负频率系数(高于 N/2+1)是共轭的正频率，并且避免了它们的计算.

For FFTW real to complex transform, if the size of the input is N real*8, then the size of the output is N/2+1 double complex. Since the input is real, coefficients of negative frequencies (higher than N/2+1) are conjugate of positive frequencies, and their computation is avoided.

根据 FFTW 关于 planner flags 的文档，

According to the documentation of FFTW regarding planner flags ,

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

FFTW_PRESERVE_INPUT specifies that an out-of-place transform must not change its input array. This is ordinarily the default, except for c2r and hc2r (i.e. complex-to-real) transforms for which FFTW_DESTROY_INPUT is the default...

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

If you wish to save the coefficients in the Fourier space out, either copy the array, or try to use FFTW_PRESERVE_INPUT. The first solution seems the best way to go if you have enough memory.

FFTW 函数的参数顺序总是 origin,destination.由于从 out 到 in 执行反向变换:

Order of arguments for FFTW functions is always origin,destination. As the backward transform is performed from out to 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 编译:

Here is a code based on yours, which performs forward and backward fft. It is compiled by 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.

After the forward fft in->out and the backward fft out->in , in is identical to its initial value, except for a scale factor N.

这篇关于fortran 复杂到真正的 fftw 问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！