fortran程序 频率域滤波

PROGRAM MAIN
IMPLICIT NONE
COMPLEX(4) C(100000), B(100000)
REAL(8) DATA(100000), AMP(100000), TIME(100000), UP, LOWER, H(100000)
REAL(8) FRE, DT, RREAL(100000), AAIMAG(100000)
INTEGER:: M, K, I, N, NN, MM, KK, GG
OPEN(10,FILE='mb.txt')
! OPEN(20,FILE='K_FRE_H.DAT')
OPEN(30,FILE='T_X.DAT')
N=0
DO
READ(10,*,END=2) TIME(N+1), DATA(N+1)
N=N+1
END DO
2 CONTINUE
CLOSE(10)
! 判断是2的多少次方
DO NN=0,50
IF (2**NN==N) THEN
EXIT
ELSE IF (2**(NN-1).LE.N.AND.2**(NN+1).GE.N) THEN
MM=2**(NN+1)-N
EXIT
END IF
END DO
DT=TIME(2)-TIME(1)
DT=DT
! 自动补零
DO KK=N+1,N+MM
TIME(KK)=TIME(N)+DT*(KK-N)
DATA(N+1)=0.
END DO
N=N+MM
DO M=1,N
C(M)=CMPLX(DATA(M),0.)
END DO
CALL FAST(N,C,N,-1) !傅里叶正变换
I=N
DO K=1,N
C(K)=C(K)/REAL(I)
AMP(K)=CABS(C(K))
END DO
! 生成频率域带通滤波器
WRITE(*,*)"PLEASE INPUT LOWER LIMIT"
READ(*,*) LOWER
WRITE(*,*)"PLEASE INPUT UPPER LIMIT"
READ(*,*) UP
DO GG=0,N/2
FRE=REAL(GG)/(REAL(N)*DT)
IF(FRE.LE.UP.AND.FRE.GE.LOWER) THEN
H(GG)=1.
ELSE
H(GG)=0.
END IF
END DO
DO GG=N/2+1,N-1
H(GG)=H(N-GG)
END DO
DO I=1,N
RREAL(I)=REAL(C(I))*H(I-1)
AAIMAG(I)=AIMAG(C(I))*H(I-1)
END DO
DO M=1,N
B(M)=CMPLX(RREAL(M),AAIMAG(M))
END DO
CALL FAST(N,B,N,+1) !傅里叶正变换

WRITE(30,12) (DT*(K-1),REAL(B(K)),K=1,N)
12 FORMAT(F15.6,TR3,F15.6)
CLOSE(30)
END
! *******************************************
! SUBROUTINE FOR FAST FOURIER TRANSFORM
! *******************************************
! N为复型数据及复型变换值的总数
! X为等间隔复型数据
! ND主程序中X的维数
! IND傅里叶变换时：IND=-1,傅里叶逆变换时：IND=1
! 在进行傅里叶变换时，IND=-1时，变换值被扩大N倍！
SUBROUTINE FAST(N,X,ND,IND)
IMPLICIT NONE
COMPLEX X(ND), TEMP, THETA
INTEGER J, I, N, M, k, KMAX, ISTEP,ND, IND
J=1
DO I=1, N
IF(I.GE.J) GOTO 110
TEMP=X(J)
X(J)=X(I)
X(I)=TEMP
110 M=N/2
120 IF(J.LE.M) GOTO 130
J=J-M
M=M/2
IF(M.GE.2) GOTO 120
130 J=J+M
END DO
KMAX=1
140 IF(KMAX.GE.N) RETURN
ISTEP=KMAX*2
DO K=1,KMAX
THETA=CMPLX(0.0,3.141593*REAL(IND*(K-1))/REAL(KMAX))
DO I=K,N,ISTEP
J=I+KMAX
TEMP=X(J)*CEXP(THETA)
X(J)=X(I)-TEMP
X(I)=X(I)+TEMP
END DO
END DO
KMAX=ISTEP
GOTO 140
END