论坛 › 编程工具交流 › intel fortran调用FFTW，傅里叶正变换值N/2+1之后为0

放耗子的绅 发表于 2019-11-20 11:32:58

intel fortran调用FFTW，傅里叶正变换值N/2+1之后为0

各位大佬们，请教一个问题， 我用intel fortran去调用FFTW，但傅里叶正变换值N/2+1之后为0，通过对比matlab算的结果，发现偏差很大；
麻烦问一下，如果我要得到N/2+1以后的数据，应该怎么做呢？代码如下：
MODULE globalvars
IMPLICIT NONE
! Simulation cell parameters:
        integer,parameter:: Nx = 4, Ny = 4
        integer(kind = 8) :: p_up, p_dn
END MODULE globalvars

!!!-------------------傅里叶正变换
subroutine backward(output,input)
USE globalvars
!calculates backward fourier transform (complex to real) using FFTW
        implicit none
        COMPLEX, INTENT(IN) :: output(nx,ny)
        REAL, INTENT(OUT)   :: input(nx,ny)
        COMPLEX :: tempk(nx,ny)
        INTEGER :: i,j
        REAL :: sizescale

        sizescale = 1.0/FLOAT(nx * ny)
        tempk = 0.0
        DO j=1,ny
                DO i=1,nx
                        tempk(i,j) = output(i,j)
                END DO
        END DO

        CALL sfftw_execute_dft_c2r(p_dn,tempk,input)
        DO j=1,ny
                DO i=1,nx
                        input(i,j) = sizescale*input(i,j)
                END DO
        END DO
!        RETURN
END subroutine backward

!!!-------------------傅里叶逆变换
subroutine forward(input,output)
USE globalvars
!calculates forward fourier transform (real to complex) using FFTW
        implicit none
       
        real, intent(IN) :: input(nx,ny)
        complex, intent(OUT) :: output(nx,ny)

        call sfftw_execute_dft_r2c(p_up,input,output)
!        RETURN
end subroutine forward


!---------------------主程序

program main

USE globalvars

        implicit none
        include 'fftw3.f'
        integer :: i,j
        real,allocatable:: eta(:,:), input(:,:)
        complex,allocatable :: etak(:,:), output(:,:)       
        allocate (eta(nx,ny), input(nx,ny))
        allocate (etak(nx,ny), output(nx,ny))       
        eta = reshape([3,0,-1,0,&
                                   -1,-2,1,1,&
                                   -2,2,2,2,&
                                   2,-1,4,3],shape(eta))        
        do i=1,nx
                do j=1,ny               
                        write(*,*) eta(j,i)               
                enddo
        enddo
        CALL sfftw_plan_dft_r2c_2d(p_up, nx, ny, input, output, FFTW_EXHAUSTIVE)
        CALL sfftw_plan_dft_c2r_2d(p_dn, nx, ny, output, input, FFTW_EXHAUSTIVE)       

        call forward(eta,etak)       
        do i=1,nx
                do j=1,ny               
                        write(*,*) etak(j,i)       
                enddo
        enddo
       
        call backward(etak,eta)       
        do i=1,nx
                do j=1,ny               
                        write(*,*) eta(j,i)               
                enddo
        enddo       
        deallocate (eta, input)
        deallocate (etak, output)
end program

运行结果见图。

li913 发表于 2019-11-23 23:33:40

本帖最后由 li913 于 2019-11-23 23:35 编辑

结果没错，只是这个函数省略了共轭对称的元素，你把它恢复就是。

放耗子的绅 发表于 2019-12-4 11:11:27

li913 发表于 2019-11-23 23:33
结果没错，只是这个函数省略了共轭对称的元素，你把它恢复就是。

非常感谢您的解答，我测试了下，结果没问题！

青衣巷 发表于 2019-12-4 22:40:41

时间域频率域信息量不变哟~

查看完整版本: intel fortran调用FFTW，傅里叶正变换值N/2+1之后为0