|
|
再看一下这个：
http://blog.renren.com/blog/247544700/703436782?bfrom=01020100300
!!时域分解FFT算法
!!离散时域信号的序号对应的时刻为nΔ，离散频域信号的序号对应的频率为m/（nΔ）
!!离散时域信号保存在文件A1.TXT中
!! 变换后的离散频域信号保存在文件 FFT_RESULT.TXT
!!编译器为cvf6.6
PROGRAM FFT
PARAMETER(K=4)
!COMPLEX用结构体定义复数
TYPE CMPLX0
REAL RE0,IM0
END TYPE
REAL M
TYPE(CMPLX0)A(0:2**K-1),A0(0:2**K-1)
OPEN(1,FILE='A1.TXT')
DO I=0,2**K-1
READ(1,*) A(I).RE0,A(I).IM0
ENDDO
! PRINT*,A
CLOSE(1)
DO I=0,2**K-1
A0(I)=A(IFIX(CTY_I(I,K)))
ENDDO
! PRINT*,A0
! PAUSE
DO I=1,K !大循环
DO J0=1,2**(K-I) !单元个数
M=0
DO J1=(J0-1)*2**I,(J0-1)*2**I+2**(I-1)-1 !每个单元的运算次数
A(J1)=PLUS(A0(J1),MULT(A0(J1+2**(I-1)),W(M,2.**K)))
A(J1+2**(I-1))=SUBT
&(A0(J1),MULT(A0(J1+2**(I-1)),W(M,2.**K)))
M=M+2**(K-I)
IF(M>=2**(K-1)) M=0
ENDDO
ENDDO
A0=A
ENDDO
OPEN(2,FILE='FFT_RESULT.TXT',ACTION='WRITE')
DO I=0,2**K-1
IF(A(I).RE0==0) THEN
WRITE(2,100) A(I).IM0
ELSEIF(A(I).IM0==0)THEN
WRITE(2,200) A(I).RE0
ELSE
WRITE(2,300)A(I).RE0,A(I).IM0
ENDIF
ENDDO
100 FORMAT(F18.3,'i')
200 FORMAT(F18.3)
300 FORMAT(F18.3,' +',F18.3,'i')
CONTAINS
!PLUS 定义复数的加法
FUNCTION PLUS(Z1,Z2)
TYPE(CMPLX0) Z1,Z2,PLUS
PLUS.RE0=Z1.RE0+Z2.RE0
PLUS.IM0=Z1.IM0+Z2.IM0
ENDFUNCTION
!SUBTRACTION 定义复数的减法
FUNCTION SUBT(Z1,Z2)
TYPE(CMPLX0) Z1,Z2,SUBT
SUBT.RE0=Z1.RE0-Z2.RE0
SUBT.IM0=Z1.IM0-Z2.IM0
ENDFUNCTION
!MULTIPLACATION 定义复数的乘法
FUNCTION MULT(Z1,Z2)
TYPE(CMPLX0) Z1,Z2,MULT
MULT.RE0=Z1.RE0*Z2.RE0-Z1.IM0*Z2.IM0
MULT.IM0=Z1.RE0*Z2.IM0+Z1.IM0*Z2.RE0
ENDFUNCTION
!W(N)**M 定义WN的M次方
FUNCTION W(M,N)
REAL M,N
TYPE(CMPLX0) W
W.RE0=COS(2*3.1415926*M/N)
W.IM0=-SIN(2*3.1415926*M/N)
ENDFUNCTION
!CONTRARY I 求I的二进制逆序数
FUNCTION CTY_I(I0,N)
INTEGER I,A(N),N,IO
REAL CTY_I
A=0
I=I0
CTY_I=0
DO K0=1,N
IF(I-2**(N-K0)<0)THEN
A(K0)=0
ELSE
A(K0)=1
I=I-2**(N-K0)
ENDIF
CTY_I=CTY_I+A(K0)*2**(K0-1)
ENDDO
ENDFUNCTION
END |
|
回复此楼