時域分解FFT算法的FORTRAN實(shí)現(xiàn)

1、!時域分解FFT算法!離散時域信號的序號對應(yīng)的時刻為n，離散頻域信號的序號對應(yīng)的頻率為m/（n）!離散時域信號保存在文件A1.TXT中! 變換后的離散頻域信號保存在文件 FFT_RESULT.TXT!編譯器為cvf6.6PROGRAM FFTPARAMETER(K=4)!COMPLEX用結(jié)構(gòu)體定義復(fù)數(shù)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!大循環(huán)DO J0=1,2*(K-I)!單元個數(shù)M=0DO J1=(J0-1)*2*I,(J0-1)*2*I+2*(I-1)-1 !每個單元的運(yùn)算次數(shù)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 定義復(fù)數(shù)的加法FUNCTION PLUS(Z1,Z2)TYPE(CMPLX0) Z1,Z2,PLUSPLUS.RE0=Z1.RE0+Z2.RE0PLUS.IM0=Z1.IM0+Z2.IM0ENDFUNCTION!SUBTRACTION 定義復(fù)數(shù)的減法FUNCTION SUBT(Z1,Z2)TYPE(CMPLX0) Z1,Z2,SUBTSUBT.RE0=Z1.RE0-Z2.RE0SUBT.IM0=Z1.IM0-Z2.IM0ENDFUNCTION!MULTIPLACATION 定義復(fù)數(shù)的乘法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的二進(jìn)制逆序數(shù)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