Fortran Coder

查看: 79|回复: 3

[数学库] intel fortran调用FFTW,傅里叶正变换值N/2+1之后为0

[复制链接]

6

帖子

2

主题

0

精华

入门

F 币
50 元
贡献
22 点
发表于 2019-11-20 11:32:58 | 显示全部楼层 |阅读模式
各位大佬们,请教一个问题, 我用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

运行结果见图。





matlab计算结果

matlab计算结果

fortran计算结果

fortran计算结果
回复

使用道具 举报

415

帖子

1

主题

0

精华

宗师

F 币
2170 元
贡献
1400 点
发表于 2019-11-23 23:33:40 | 显示全部楼层
本帖最后由 li913 于 2019-11-23 23:35 编辑

结果没错,只是这个函数省略了共轭对称的元素,你把它恢复就是。
微信图片_20191123191853.jpg

6

帖子

2

主题

0

精华

入门

F 币
50 元
贡献
22 点
 楼主| 发表于 2019-12-4 11:11:27 | 显示全部楼层
li913 发表于 2019-11-23 23:33
结果没错,只是这个函数省略了共轭对称的元素,你把它恢复就是。

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

13

帖子

0

主题

0

精华

熟手

F 币
158 元
贡献
66 点

规矩勋章

QQ
发表于 2019-12-4 22:40:41 | 显示全部楼层
时间域频率域信息量不变哟~
彼岸,有永恒的守候...
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2019-12-14 05:36

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表