时域分解FFT算法的FORTRAN实现.docxVIP

!!时域分解FFT算法!!离散时域信号的序号对应的时刻为nΔ，离散频域信号的序号对应的频率为m/（nΔ）!!离散时域信号保存在文件A1.TXT中!!变换后的离散频域信号保存在文件FFT_RESULT.TXT!!编译器为cvf6.6PROGRAM FFTPARAMETER(K=4)!COMPLEX用结构体定义复数TYPE CMPLX0REAL RE0,IM0END TYPEREAL MTYPE(CMPLX0)A(0:2**K-1),A0(0:2**K-1)OPEN(1,FILE=A1.TXT)DO I=0,2**K-1READ(1,*) A(I).RE0,A(I).IM0ENDDO!PRINT*,ACLOSE(1)DO I=0,2**K-1A0(I)=A(IFIX(CTY_I(I,K)))ENDDO!PRINT*,A0!PAUSEDO I=1,K!大循环DO J0=1,2**(K-I)!单元个数M=0DO 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=0ENDDOENDDOA0=AENDDOOPEN(2,FILE=FFT_RESULT.TXT,ACTION=WRITE)DO I=0,2**K-1IF(A(I).RE0==0) THENWRITE(2,100) A(I).IM0ELSEIF(A(I).IM0==0)THENWRITE(2,200) A(I).RE0ELSEWRITE(2,300)A(I).RE0,A(I).IM0ENDIFENDDO100FORMAT(F18.3,i)200FORMAT(F18.3)300FORMAT(F18.3, +,F18.3,i)CONTAINS!PLUS 定义复数的加法FUNCTION PLUS(Z1,Z2)TYPE(CMPLX0) Z1,Z2,PLUSPLUS.RE0=Z1.RE0+Z2.RE0PLUS.IM0=Z1.IM0+Z2.IM0ENDFUNCTION!SUBTRACTION 定义复数的减法FUNCTION SUBT(Z1,Z2)TYPE(CMPLX0) Z1,Z2,SUBTSUBT.RE0=Z1.RE0-Z2.RE0SUBT.IM0=Z1.IM0-Z2.IM0ENDFUNCTION!MULTIPLACATION 定义复数的乘法FUNCTION MULT(Z1,Z2)TYPE(CMPLX0) Z1,Z2,MULTMULT.RE0=Z1.RE0*Z2.RE0-Z1.IM0*Z2.IM0MULT.IM0=Z1.RE0*Z2.IM0+Z1.IM0*Z2.RE0ENDFUNCTION!W(N)**M 定义WN的M次方FUNCTION W(M,N)REAL M,NTYPE(CMPLX0) WW.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,IOREAL CTY_IA=0I=I0CTY_I=0DO K0=1,NIF(I-2**(N-K0)0)THENA(K0)=0ELSEA(K0)=1I=I-2**(N-K0)ENDIFCTY_I=CTY_I+A(K0)*2**(K0-1)ENDDOENDFUNCTIONEND