離散フーリエ変換のコード

フーリエ変換については下記のURLが分かりやすい。
フーリエ級数展開をベクトルで直観的に理解するblog.physips.com

デジタル信号→離散フーリエ変換(DFT)→フーリエ級数(FS)→元の関数の出力(FA)という感じのコードです．
更新していきますが，とりあえずメモ代わり．
同一ディレクトリにDFT.f90,FS.f90,FA.f90って作って，

gfortran DFT.f90 FS.f90 FA.f90
./a.out

とかでいちおう動きます．

DFT.f90

Subroutine DFT(a,b_re,b_im,N,T)                                                                     
! This program calculates the descrete fourier transform.                                           
! This program transform tha array according to                                                     
! Xk = T / N * \sum^N-1_n=0 xn W^nk_N (k = 0,1,2,\cdots,N-1)                                        
! T and N is defined below.                                                                         
                                                                                                    
! Input                                                                                             
! array a   : the sampling data                                                                     
! double  T : sampling time                                                                         
! integer N : sampling frequency                                                                    
! Output                                                                                            
! array b_re,b_im : the transformed data                                                            
                                                                                                    
! Sample data                                                                                       
! Input                                                                                             
! N = 8, T = 1.d0                                                                                   
! a(i) = i / 8.d0                                                                                   
! Output                                                                                            
!     X(i)     b_re     b_im                                                                        
!     X(0)    0.438    0.000                                                                        
!     X(1)   -0.062    0.151                                                                        
!     X(2)   -0.062    0.063                                                                        
!     X(3)   -0.062    0.026                                                                        
!     X(4)   -0.062    0.000                                                                        
!     X(5)   -0.063   -0.026                                                                        
!     X(6)   -0.063   -0.062                                                                        
!     X(7)   -0.063   -0.151           
                                                                                                    
  implicit none                                                                                     
  real(8), dimension(0:N-1) :: a,b_re,b_im                                                          
  real(8) :: pi,T                                                                                   
  integer :: i,j,k,N                                                                                
                                                                                                    
  pi = 3.14159265358979                                                                             
                                                                                                    
  do i = 0,N-1                                                                                      
   b_re(i) = 0.d0                                                                                   
   b_im(i) = 0.d0                                                                                   
  end do                                                                                            
                                                                                                    
  do i = 0,N-1    
    do j = 0,N-1                                                                                    
      b_re(i) = b_re(i) + (dble(T) / dble(N)) * cos(2.d0*pi*j*i/8.d0) * a(j)                        
      b_im(i) = b_im(i) + (dble(T) / dble(N)) * (-1.d0) * sin(2.d0*pi*j*i/8.d0) * a(j)              
    end do                                                                                          
  end do                                                                                            
                                                                                                    
  write(*,'(3A9)') "X(i)","b_re","b_im"                                                             
  do i = 0,N-1                                                                                      
    write(*,'(A7,I1,A1,2F9.3)') "X(",i,")", b_re(i),b_im(i)                                         
  end do                                                                                            
end subroutine DFT           

FS.f90

subroutine FS(b_re,b_im,x,N,T)                                                                      
! This program reverse the fourier transform data                                                   
  implicit none                                                                                     
  real(8), dimension(0:N-1) :: b_re,b_im                                                            
  integer :: N, it, i, j, k                                                                         
  real(8) :: pi, T, x, dt, rt                                                                       
                                                                                                    
  pi = 3.141592653589793238                                                                         
  dt = 0.01                                                                                         
  x  = 0.d0                                                                                         
                                                                                                    
  open(10,file='fourier.txt',status='unknown')                                                      
  write(*,'(2A10)') "t","x"                                                                         
                                                                                                    
  N = N / 2                                                                                         
  do it = 1,100                                                                                     
   rt = dt * it                                                                                     
   x = b_re(0)                                                                                      
   do i = 1, N                                                                                      
     x = x + (b_re(i)+b_re(8-N)) * cos(2 * pi * i * rt / T) + (-b_im(i)+b_im(8-i)) * sin(2 * pi * i * rt / T)
   end do                                                                                           
   write(*,'(2F10.2)') rt, x                                                                        
   write(10,'(2F10.2)') rt, x                                                                       
  end do                                                                                            
end subroutine FS           

FA.f90

program FA                                                                                          
  implicit none                                                                                     
  real(8), allocatable, dimension(:) :: a, b_re, b_im                                               
  integer :: N,i                                                                                    
  real(8) :: T,x                                                                                    
                                                                                                    
  N = 8                                                                                             
  T = 1.d0                                                                                          
  allocate(a(0:N-1),b_re(0:N-1),b_im(0:N-1))                                                        
                                                                                                    
  do i = 0,7                                                                                        
   a(i) = dble(i) / dble(N)                                                                         
   b_re(i) = 0.d0                                                                                   
   b_im(i) = 0.d0                                                                                   
  end do                                                                                            
                                                                                                    
  call DFT(a,b_re,b_im,N,T)                                                                         
  call FS(b_re,b_im,x,N,T)                                                                          
end program FA