numerical recipe 푸리에 트랜스폼 관련 질문이요...-.-;;;

[걸음느린아이]

numerical recipes in fortran 책을 보고 푸리에 트랜스폼을 공부하고 있는데요.

이걸 적용을 어떻게 시켜야 될지 모르겠습니다.

시간에 관한 프로그램과 값은 아래와 같이 구하였고

이것을 subroutine 파일과 연동하려고 하는데..

SUBROUTINE four1(data,nn,isign)이 부분의 data 값에

시간 축에서 나온 값을 넣어야 하는걸로 아는데

어떻게 넣어야 할지 모르겠네요..

프로그램 run 시켜보니 오류가 나구요..-.-;

프로그램은 power station 4.0 쓰고 있습니다.

파일과 관련 자료 첨부하였습니다.

조언 부탁드립니다.

Text5.f

f12-2.pdf

C
C
 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
 REAL*8 K11,K12,K1,K21,K22,K2
 REAL*8 I11,I12,I13,I21,I22,I23,I10,I20,IT
 COMPLEX*8 RIT
      OPEN(UNIT=7,FILE='D:\test four1\four1 test.txt',
     $  STATUS='OLD')
C
      I1=16.6D0
 I2=9.3D0
 TAU1=1.1D0
 TAU2=2.D0
 TAU3=12.D0
 TAU4=37.D0
 N=1.8
C
      DO 818 LT=1,10000,1
      TSCAL=1.D-11
 TIME=DFLOAT(LT)*TSCAL
 TN=TIME/1.D-9 
C
 K11=-(TAU1/TAU2)
      K12=(N*TAU2/TAU1)**(1/N)
 K1=DEXP(K11*K12)
C
 K21=-(TAU3/TAU4)
 K22=(N*TAU4/TAU3)**(1/N)
 K2=DEXP(K11*K12)
C    
      I11=I1/K1
 I12=((TN/TAU1)**N)/(1+((TN/TAU1)**N))
 I13=DEXP(-TN/TAU2)
C
      I21=I2/K2
 I22=((TN/TAU3)**N)/(1+((TN/TAU3)**N))
 I23=DEXP(-TN/TAU4)
C
      I10=I11*I12*I13
 I20=I21*I22*I23
 IT=I10+I20
 RIT=DBLE(IT)
C
C     ****************
 ISIGN=1
 NN=1000
C     ****************
C     
C
      CALL FOUR1(IT,NN,ISIGN)
C
C
      WRITE(6,784) TN,RIT
  784 FORMAT(1X,'TIME='F8.2,3X,'REAL I=',F12.4,3X,'IMAGINE I=',F12.4)
      WRITE(7,884) TN,RIT
  884 FORMAT(1X,'TIME='F8.2,3X,'REAL I=',F12.4,3X,'IMAGINE I=',F12.4)
C
  818 CONTINUE
      CLOSE(UNIT=7,STATUS='KEEP')
 STOP
 END
C
C***********************************************************************
C                                                                      *
C         FAST FOURIER TRANSFORM SUBROUTINE                            *
C                                                                      *
C***********************************************************************

      SUBROUTINE four1(data,nn,isign)
      INTEGER isign,nn
      REAL data(2*nn)

      INTEGER i,istep,j,m,mmax,n
      REAL tempi,tempr
      DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp
      n=2*nn
      j=1
      do 11 i=1,n,2
        if(j.gt.i)then
          tempr=data(j)
          tempi=data(j+1)
          data(j)=data(i)
          data(j+1)=data(i+1)
          data(i)=tempr
          data(i+1)=tempi
        endif
        m=n/2
1       if ((m.ge.2).and.(j.gt.m)) then
          j=j-m
          m=m/2
        goto 1
        endif
        j=j+m
11    continue
      mmax=2
2     if (n.gt.mmax) then
        istep=2*mmax
        theta=6.28318530717959d0/(isign*mmax)
        wpr=-2.d0*sin(0.5d0*theta)**2
        wpi=sin(theta)
        wr=1.d0
        wi=0.d0
        do 13 m=1,mmax,2
          do 12 i=m,n,istep
            j=i+mmax
            tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
            tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
            data(j)=data(i)-tempr
            data(j+1)=data(i+1)-tempi
            data(i)=data(i)+tempr
            data(i+1)=data(i+1)+tempi
12        continue
          wtemp=wr
          wr=wr*wpr-wi*wpi+wr
          wi=wi*wpr+wtemp*wpi+wi
13      continue
        mmax=istep
      goto 2
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 0"0)P.

[happyday19c]

일단 four이라는 서브루틴에서 data는 배열이고, 이를 호출할때는 그냥 REAL이네요. 이 부분이 문제입니다. 데이터 형식 점검을 다시 해보시는 것이 좋겠습니다.
푸리에 변환은 시간에 따른 전체 데이터가 있어야 합니다. 그러므로 이 데이터가 어디에 있는지가 관건입니다.

[Sentimental]

subroutine four1(data,nn,sign)에서 샘플링 complex x(1,4)이면, x(1,1)=(10,0), x(1,2)=(5,0), x(1,3)=(5,0),x(1,4)=(10,0) 이런 형태여야 하고, real data(1:8) 또는 real real IT(1:8) 이렇게 하시면,
N=4
do i=1,N
IT(2*i-1)=real(x(1,i))
IT(2*i)=imag(x(1,i))
end do

call four1(IT,4,+1) ! one-dimensional forward FFT

do i=1,N
x1(1,i)=cmplx(IT(2*i-1),IT(2*i))
end do

이렇게 하셔야 합니다.