Fortran Coder

标题: 求助检查傅里叶变换程序 [打印本页]

作者: alohomora100 时间: 2017-12-12 19:56
标题: 求助检查傅里叶变换程序
我有1000个CT数据，每个CT间隔1.2飞秒，现在对其进行傅里叶变换，只考虑实部，W是角频率想从0-1000，T是时间从0到1200飞秒。傅里叶变换后想得到横轴W（频率），纵轴FCT的图。但是结果貌似不是预期的样子，和origin中自行FFT做的图不吻合，还请各位帮忙看一下程序。谢谢诸位！

[Fortran] 纯文本查看 复制代码

PROGRAM FOURIER 

      IMPLICIT NONE

      REAL(4)::DT,T,FCT

      CHARACTER::U

      INTEGER::I,W

      REAL(4),DIMENSION(1001)::CT

      OPEN(12,FILE='VACFHI.txt')

      OPEN(13,FILE='FT.txt')

      DO I=0,1000

        READ(12,*)U,CT(I)

      ENDDO

      DO W=0,1000 

         DT=0.012

         T=0

         FCT=0

        DO I=0,1000

    IF (T.LT.1001) THEN

        FCT=FCT+2*(CT(I)*COS(W*T)*DT)

        T=T+DT

        ENDIF

        ENDDO

        WRITE(13,*)W/(0.012/1000)*33,FCT        

      ENDDO

      CLOSE(12)

      CLOSE(13)

      END

草图.png (70.06 KB, 下载次数: 248)

左上输入文件图，右上傅里叶变换图，左下origin做的FFT

FT.txt

30.3 KB, 下载次数: 3

傅里叶变换后输出

VACFHI.txt

36.17 KB, 下载次数: 1

读取的文件

作者: li913 时间: 2017-12-12 20:27
给出你的结果、期望结果、输入数据。

作者: alohomora100 时间: 2017-12-12 20:52

li913 发表于 2017-12-12 20:27
给出你的结果、期望结果、输入数据。

您好，相关文件已经显示在附件中

作者: fcode 时间: 2017-12-13 10:50
看起来，似乎有2个问题。
1. 数组越界。你的数组是从 1:1001，但是循环从 0 开始到 1000，会越界。你使用的编译器没有检查出这个问题，你可以考虑换一个编译器。
2. 很奇怪的是，我觉得你的结果似乎没什么问题。但是你需要理解实序列的DFT，结果只有一半的频率是有效的。你可以把 X 换做频率，而非数据序数。与 orgin 的结果对比，就会差不多一致了。

另外，给你一个建议。如果你的代码里大量出现 1000 这个数字，对于将来修改，会是很困难的。建议你用一个整形常量表示，比如 Integer , parameter :: N = 1000
这样，将来容易修改，如下代码：

[Fortran] 纯文本查看 复制代码

program fourier

  implicit none

  real(4)::dt,t,fct

  character::u

  integer::i,w

  Integer , parameter :: N = 1000

  real(4),dimension(0:N)::ct

  open(12,file='vacfhi.txt')

  open(13,file='ft.txt')

  do i=0,N

    read(12,*)u,ct(i)

  enddo

  do w=0,N

    dt=0.012

    t=0

    fct=0

    do i=0,N

      !if ( t <(N+1) ) then

        fct=fct+2*(ct(i)*cos(w*t)*dt)

        t=t+dt

      !endif

    enddo

    write(13,*)w/(0.012/N)*33,fct

  enddo

  close(12)

  close(13)

end

作者: alohomora100 时间: 2017-12-13 11:45

fcode 发表于 2017-12-13 10:50
看起来，似乎有2个问题。
1. 数组越界。你的数组是从 1:1001，但是循环从 0 开始到 1000，会越界。你使用的 ...

多谢，多谢

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)