Fortran Coder

查看: 4822|回复: 3
打印 上一主题 下一主题

[通用算法] 傅里叶逆变换

[复制链接]

23

帖子

9

主题

0

精华

熟手

F 币
134 元
贡献
82 点
跳转到指定楼层
楼主
发表于 2017-12-29 11:39:02 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
最近写了一个只取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
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

796

帖子

2

主题

0

精华

大宗师

F 币
3787 元
贡献
2266 点
沙发
发表于 2017-12-30 09:47:04 | 只看该作者
自己验证一下,逆变换能恢复原信号不。

23

帖子

9

主题

0

精华

熟手

F 币
134 元
贡献
82 点
板凳
 楼主| 发表于 2017-12-31 11:49:48 | 只看该作者
li913 发表于 2017-12-30 09:47
自己验证一下,逆变换能恢复原信号不。

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

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
地板
发表于 2018-1-1 21:15:38 | 只看该作者
你只取实部的话,肯定是不能还原原有信号的。

你可以用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
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-4-20 08:23

Powered by Tencent X3.4

© 2013-2024 Tencent

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