时域分解FFT算法的FORTRAN实现

1、!时域分解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!PRIN

2、T*,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=

3、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')CONTA

4、INS!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,MULTENDFUNCTION!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