Fortran Coder

标题: 傅里叶逆变换 [打印本页]

作者: alohomora100    时间: 2017-12-29 11:39
标题: 傅里叶逆变换
最近写了一个只取5实部的傅里叶变换,现在想对结果做个逆变换,不知逆变换程序写得对不对,望各位帮助修改,多谢!
[Fortran] 纯文本查看 复制代码
!      傅里叶变换
      PROGRAM FOURIER
      IMPLICIT NONE
      REAL(8)::DT,t,FCT,df,pi,w
      CHARACTER::U
      INTEGER::I,j
      REAL(8),DIMENSION(10001)::CT
      OPEN(12,FILE='VACFHI.txt')
      OPEN(13,FILE='testFT.txt')
      dt=1.2     ! fs
      df=1.d-4   ! 1 fs^-1 = 33000 cm^-1
      pi=4.d0*atan(1.d0)
      DO i=1,10001
        READ(12,*)U,CT(I)
      ENDDO
      DO j=1,12500
         FCT=0
         w=2.d0*pi*dble(j)*df
        DO i=0,10000
          t=dble(i)*DT
          FCT=FCT+2.d0*(CT(I+1)*dCOS(w*t)*DT)
        ENDDO
        WRITE(13,*)w/pi/2.d0*3.3d4,FCT        
      ENDDO
      CLOSE(12)
      CLOSE(13)
      END      
!傅里叶逆变换
      program inverse FT
      implicit none
      integer,parameter::N=12500
      real(8),dimension(N)::FT,u
      integer::i,j
      real(8)::t,dw,CT,w
      open(11,file='testFT.txt')
      open(12,file='inverseFT.txt')
  
      do i=1,N
      read(11,*)u(i),FT(i)
      enddo
      
      do i=1,12500
        CT=0
        t=i
        dw=3.3     
        do j= 1,N
          w=3.3*j
          CT=CT+2.d0*(FT(j)*cos(w*t)*dw)
        enddo
        write(12,*)t,CT
      enddo
      close(11)
      close(12)
      end

作者: li913    时间: 2017-12-30 09:47
自己验证一下,逆变换能恢复原信号不。
作者: alohomora100    时间: 2017-12-31 11:49
li913 发表于 2017-12-30 09:47
自己验证一下,逆变换能恢复原信号不。

不能,所以不知道是不是写的不对,谢谢。

作者: vvt    时间: 2018-1-1 21:15
你只取实部的话,肯定是不能还原原有信号的。

你可以用matlab来验证一下自己的结果。

[Fortran] 纯文本查看 复制代码
PROGRAM FOURIER 
      IMPLICIT NONE
      Real(8) , parameter :: PI = acos(-1.d0)
      REAL(8)::FCT,w
      CHARACTER::U
      INTEGER::I,j
      Integer , parameter :: N = 10001
      Integer , parameter :: M = 12501
      REAL(8),DIMENSION(N)::CT
      OPEN(12,FILE='VACFHI.txt')
      OPEN(13,FILE='testFT.txt')
      DO i=1,N
        READ(12,*)U,CT(I)
      ENDDO
      DO j=0,M-1
         FCT=0
         w=2.d0*pi*j/M
        DO i=0,N-1
          FCT=FCT+(CT(I+1)*COS(w*i))
        ENDDO
        WRITE(13,*)w/pi/2.d0*3.3d4,FCT        
      ENDDO
      CLOSE(12)
      CLOSE(13)
      END





欢迎光临 Fortran Coder (http://bbs.fcode.cn/) Powered by Discuz! X3.2