fortran复杂到真正的fftw问题

问题描述

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

 程序testit 
 INCLUDE'fftw3.f'
 $ double $ out !, in 
 real 
 parameter（N = 100）
维数（N），出（N）
整数* 8 p，p2 
整数i，j 
实数x 
实数
写入（*，*）数据中的数据
 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（* （i）
 WRITE（12，*）real（in（i））
 enddo 
 write（*，*）中的（f10.2,1x，f10.2） 创建计划
调用dfftw_plan_dft_r2c_1d（p，N，in，out，FFTW_ESTIMATE）
调用dfftw_plan_dft_c2r_1d（p2，N，in，out，FFTW_ESTIMATE）
 write（*，*）do它
 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 
写（*，*）撤消它
调用dfftw_execute_dft_c2r（p2，输入，输出）
 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）
结束程序
  

真实到复杂的转换工作就好了。逆变换给出了错误的值，并以某种方式修改了输入和输出变量。我不知道问题是什么，我无法在网上找到答案。

预先感谢！

乍得W.弗里尔

编辑：我也想知道在fftw包中是否有与matlab中的fftshift（）和ifftshift（）函数类似的函数。

解决方案

存在精度问题：在大多数计算机上， real 指的是单精度32位浮点数。 real * 8 可用于指定双精度浮点数，以与 double complex 和<$ c一致对于FFTW实数到复数变换，如果输入的大小是。

N real * 8 ，那么为 N / 2 + 1 double complex 。由于输入是真实的，因此负频率系数（高于N / 2 + 1）是共轭的的正频率，并且可以避免它们的计算。

根据FFTW有关，

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

FFTW函数的参数顺序总是 origin，destination 。由于后向变换是在中执行的 out 到：

 调用dfftw_plan_dft_c2r_1d（p2，N，out，in，FFTW_ESTIMATE）
 ... 
调用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 
整数i 
实数x 
实数
写（* ，*）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 
在（i）中写入（*，（f10.2,1x /））
 WRITE（12，* ）
 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_execu te_dft_c2r（p2，out，in）
 fact = 1.0 / N 
 write（*，*）input back，除了缩放
 do i = 1，N 
 write （*，（f10.2,1x））在（i）
 enddo 
 write（*，*）输出可以被c2r修改...除非FFTW_PRESERVE_INPUT被设置为$ b （$，$（$））$（$）$（$）$（$）$（$） *，*）清理
调用dfftw_destroy_plan（p，in，out）
调用dfftw_destroy_plan（p2，out，in）
结束程序
  

在 - > out ， in 中的$ c>和fft out - > c>与其初始值相同，除了比例因子 N 。

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.

Thanks in advance!

Chad W. Freer

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

解决方案

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_....

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.

According to the documentation of FFTW regarding planner flags ,

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.

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)

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

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问题的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！