马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
我是用fortran编的,但是滤波结果与原始数据有较大差别,请高手指点,谢谢了
C NUMERICAL FILTER 滤波
PARAMETER(MN=5000000)
DIMENSION AD(MN),DD(MN),IO(MN),DD2(MN)
CHARACTER FINAME*15,FNA*15,SQUASH*40
REAL*4 REAL(MN),IMAG(MN)
OPEN(2,FILE=SQUASH(FINAME//'.800'))
C OPEN(12,FILE='AD.DAT')
READ(2,*)(AD(I),I=1,NDH*NTOTAL) !读入数据
*****补零使数据长度为2的N次幂**************
MI=1
222 MNT=2**MI
IF(MNT.LE.NTOTAL) THEN
MI=MI+1
GO TO 222
ELSE
DO J=1,NDH
DO I=NTOTAL+1,MNT
K1=(I-1)*NDH+J
K2=(I-2)*NDH+J
AD(K1)=AD(K2)
END DO
END DO
NTOTAL=MNT
END IF
FN=1/((SCAN)*2.0)
NFREQ=NTOTAL/FN
DF=2./NFREQ
*****开始滤波******
IFLT=1 !低通
IF(IFLT.EQ.1) THEN
WRITE(*,'(A\)')' ENTER HIGH CUT OFF FREQ.(Hz) : '
C READ(*,*) FHI
FHI=0.02
NHI=FHI/DF
END IF
IF(IFLT.EQ.2) THEN !高通
WRITE(*,'(A\)')' ENTER LOW CUT OFF FREQ.(Hz) : '
READ(*,*) FLO
NLO=FLO/DF
END IF
IF(IFLT.EQ.3) THEN !带通
WRITE(*,'(A\)')' ENTER HIGH CUT OFF FREQ.(Hz) : '
READ(*,*) FHI
NHI=FHI/DF
WRITE(*,'(A\)')' ENTER LOW CUT OFF FREQ.(Hz) : '
READ(*,*) FLO
FLO=FLO/DF
END IF
****************************************************************
DO J=1,NDH
DO I=1,NTOTAL
K=(I-1)*NDH+J
REAL(I)=AD(K)
IMAG(I)=0.0
END DO
CALL FFT842(0,NTOTAL,REAL,IMAG) !快速fft
IF (IFLT.EQ.1) THEN
DO I=NHI+1,NTOTAL
REAL(I)=0.0
IMAG(I)=0.0
END DO
END IF
IF (IFLT.EQ.2) THEN
DO I=1,NLO
REAL(I)=0.0
IMAG(I)=0.0
END DO
END IF
IF (IFLT.EQ.3) THEN
DO I=NHI+1,NTOTAL
REAL(I)=0.0
IMAG(I)=0.0
END DO
DO I=1,NLO
REAL(I)=0.0
IMAG(I)=0.0
END DO
END IF
CALL FFT842(1,NTOTAL,REAL,IMAG) !fft反变换
DO I=1,NTOTAL
K=(I-1)*NDH+J
DD(K)=REAL(I)
END DO
END DO |