fortran复杂到真正的fftw问题 [英] fortran complex to real fftw issue
问题描述
我目前正在研究一个需要实现傅立叶变换和逆变换的项目。我正在测试我从一个在线示例中修改的程序;打印或写入命令通常用于调试:
程序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有关规划师标志,
FFTW_PRESERVE_INPUT
指定不在场转换不得改变其输入数组。这通常是默认设置,除了c2r和hc2r(即复杂到实际)转换,其中FFTW_DESTROY_INPUT
是默认设置...
如果您希望将系数保存在傅立叶空间 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 $ c
, in $ c $>中的$ 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 ,
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 whichFFTW_DESTROY_INPUT
is the default...
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问题的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!